Juliaで量子スピン系のハミルトニアンを構成する

# はじめに
次のようなハミルトニアン

```math
H = \sum_{i,j}\left( 
J_x\sigma^x_i\sigma^x_{j} 
+J_y\sigma^y_i\sigma^y_{j}
+J_z\sigma^z_i\sigma^z_{j}
\right)
```

をJuliaで具体的に構成する方法について考えます。ここで、$ \sigma^x, \sigma^y , \sigma^z$ は[パウリ行列](https://ja.wikipedia.org/wiki/%E3%83%91%E3%82%A6%E3%83%AA%E8%A1%8C%E5%88%97)を表します。また下付き添字は、スピンが位置するサイトの番号を表します。$N$を全サイト数とすると、 $\sigma^z_i$ は

```math
\sigma^z_i = 
\sigma^0 \otimes \cdots \otimes \sigma^0 \otimes \sigma^z
\otimes \sigma^0 \otimes \cdots \otimes \sigma^0
```

で与えられる $2^N \times 2^N$ 行列です。ここで、$\otimes$ は[クロネッカー積](https://ja.wikipedia.org/wiki/%E3%82%AF%E3%83%AD%E3%83%8D%E3%83%83%E3%82%AB%E3%83%BC%E7%A9%8D)、$\sigma^0$ は $2\times 2$ の単位行列です。

# パウリ行列の定義
パウリ行列とは、

```math
\sigma^x = 
\begin{pmatrix}
0 & 1 \\ 1 & 0
\end{pmatrix}, \quad
\sigma^y = 
\begin{pmatrix}
0 & -i \\ i & 0
\end{pmatrix}, \quad
\sigma^z = 
\begin{pmatrix}
1 & 0 \\ 0 & -1
\end{pmatrix}
```

というものでした。具体的に成分が与えられているので、これをプログラムすることは難しくありません。

```julia
using SparseArrays

σ⁰ = sparse([1.0+0.0im 0.0; 0.0 1.0+0.0im])
σˣ = sparse([0.0 1.0+0.0im; 1.0+0.0im 0.0])
σʸ = sparse([0.0 -1.0im; 1.0im 0.0])
σᶻ = sparse([1.0+0.0im 0.0; 0.0 -1.0+0.0im])
```

Juliaでは`im`で虚数単位を表します。ここでは後のことを見越して、疎行列型としてパウリ行列を定義しました。Juliaでは[変数名にUnicode文字を使用できる](https://docs.julialang.org/en/v1/manual/unicode-input/)ので、`σᶻ`のようなギリシャ文字や上付き添字などが使えます。

念の為、代数関係が成り立っていることを確認しておきます。

```shell-session
julia> σᶻ*σᶻ == σ⁰
true

julia> σˣ*σʸ - σʸ*σˣ == 2.0im*σᶻ
true
```

# クロネッカー積
行列$A,\,B$のクロネッカー積とは

```math
A \otimes B =
\begin{pmatrix}
a_{11}B & \dots & a_{1N}B \\
\vdots & \ddots & \vdots \\
a_{N1}B & \dots & a_{NN}B
\end{pmatrix}
```

で定義される演算です。$a_{ij}$ は $A$ の $i,j$ 成分を表します。Juliaでは`kron(A,B)`によって計算することが可能です。あるいは、

```julia
⊗(A, B) = kron(A, B)
```

のように`⊗`を定義することで、`A⊗B`とより直観的に書くこともできます。

クロネッカー積の基本的な性質を確認しておきます。

```shell-session
julia> A, B, C, D = rand(2, 2), rand(2, 2), rand(2, 2), rand(2, 2)
julia> A⊗(B+C) ≈ A⊗B + A⊗C
true
julia> (B+C)⊗A ≈ B⊗A + C⊗A
true
julia> (A⊗B)⊗C ≈ A⊗(B⊗C)
true
julia> (A⊗B)*(C⊗D) ≈ (A*C)⊗(B*D)
true
```

`≈`は右辺と左辺が機械精度の範囲で一致するかどうかを判定する演算子で、このようなテストをするときに重宝します。Juliaに対応したエディタ上で`\approx`と入力すると出ます。

# ハミルトニアンの構成法1（ベタ打ち, 1次元）
以上の準備のもと、1次元Ising模型のハミルトニアン


```math
H = \sum_{i=1}^N 
\sigma^z_i\sigma^z_{i+1}
```

を構成してみます。$N=3,4$ のときの具体形を手計算で求めておくと

```math
\mathrm{diag}(3,-1,-1,-1,-1,-1,-1,3) \\
\mathrm{diag}(4,0,0,0,0,-4,0,0,0,0,-4,0,0,0,0,4)
```

となります。$\mathrm{diag}(\dots)$ は対角成分がその引数で与えられ、それ以外の成分が0であるような行列を表します。

さて、Julia上ではハミルトニアンを

```julia
H3 = σᶻ⊗σᶻ⊗σ⁰ + σ⁰⊗σᶻ⊗σᶻ + σᶻ⊗σ⁰⊗σᶻ
H4 = σᶻ⊗σᶻ⊗σ⁰⊗σ⁰ + σ⁰⊗σᶻ⊗σᶻ⊗σ⁰ + σ⁰⊗σ⁰⊗σᶻ⊗σᶻ + σᶻ⊗σ⁰⊗σ⁰⊗σᶻ
```

と書くことができます。もはや計算メモを書くのと全く変わりません。手計算の結果と一致することは、

```shell-session
julia> using LinearAlgebra
julia> H3 == Diagonal([3,-1,-1,-1,-1,-1,-1,3])
true
julia> H4 == Diagonal([4,0,0,0,0,-4,0,0,0,0,-4,0,0,0,0,4])
true
```

により確認できます。

# ハミルトニアンの構成法2（一般のN, 1次元）
一般の $N$ に対応するため、クロネッカー積を再帰的に計算できるようにしておきます。

```julia
function ⊗(A,B,C...)
    A = kron(A, B)
    for Ci in C
        A = kron(A, Ci)
    end
    return A
end
```

`C...`は可変個の引数を取れることを表します。これは次のように動作します。

```shell-session
julia> ⊗(σᶻ,σᶻ,σ⁰,σ⁰) == σᶻ⊗σᶻ⊗σ⁰⊗σ⁰
true
julia> ⊗([σᶻ,σᶻ,σ⁰,σ⁰]...) == σᶻ⊗σᶻ⊗σ⁰⊗σ⁰
true
```

あとは、`[σᶻ,σᶻ,σ⁰,σ⁰]` に相当する配列を自由に渡せるようにすればOKです。例えば次が考えられます。

```julia
function set_spins(N, sites, σs)
    list_mats = fill(σ⁰, N)
    for (site, σ) in zip(sites, σs)
        list_mats[site] = σ
    end
    return list_mats
end
```
```shell-session
julia> set_spins(4, [1,2], [σᶻ,σᶻ]) == [σᶻ,σᶻ,σ⁰,σ⁰]
true
julia> set_spins(4, [1,4], [σᶻ,σˣ]) == [σᶻ,σ⁰,σ⁰,σˣ]
true
```

これらを用いると、Isingハミルトニアンは

```julia
function HIsing(N)
    H = ⊗(set_spins(N, [N,1], [σᶻ,σᶻ])...)
    for i in 1:N-1
        H += ⊗(set_spins(N, [i,i+1], [σᶻ,σᶻ])...)
    end
    return H
end
```

で構成できます。前節で作ったものと比較すると、

```shell-session
julia> HIsing(3) == H3
true
julia> HIsing(4) == H4
true
```

となります。一般のスピン系の場合も同様に構成することができます。冒頭のHeisenbergハミルトニアンは

```julia
function HHeisenberg(N)
    H = ⊗(set_spins(N, [N,1], [σˣ,σˣ])...)
    H += ⊗(set_spins(N, [N,1], [σʸ,σʸ])...)
    H += ⊗(set_spins(N, [N,1], [σᶻ,σᶻ])...)
    for i in 1:N-1
        H += ⊗(set_spins(N, [i,i+1], [σˣ,σˣ])...)
        H += ⊗(set_spins(N, [i,i+1], [σʸ,σʸ])...)
        H += ⊗(set_spins(N, [i,i+1], [σᶻ,σᶻ])...)
    end
    return H
end
```
で与えられます。

量子系である以上、取りうる $N$ の大きさにはタイトな限界がありますが、1次元ハイゼンベルク模型ならば、パソコン上でも $N=20$ 程度は工夫なしにメモリに乗せることができると思います。

# 2次元正方格子の場合
ここでは、$N \times N$の正方格子上のIsing模型を考えます。格子に以下のように番号を振ったとしましょう。

```julia
SiteIndeces_2d(N) = transpose(reshape(1:N^2, N,N))
```

```shell-session
julia> SiteIndeces_2d(3)
3×3 Transpose{Int64,Base.ReshapedArray{Int64,2,UnitRange{Int64},Tuple{}}}:
 1  2  3
 4  5  6
 7  8  9
```

正方格子上で周期条件のもとにIsing模型を考えると、$N=3$の場合、ボンドは(1-2)(2-3)(3-1)(4-5)(5-6)(6-4)(7-8)(8-9)(9-7)(1-4)(4-7)(7-1)(2-5)(5-8)(8-2)(3-6)(6-9)(9-3)の計18本あります。一般の$N$に対応するためには、これらのボンドを数え上げるコードを書く必要があります。落ち着いて考えれば難しくありませんが、ちょっと面倒だし、どこかでミスをしそうです。そこでJuliaの`CartesianIndex`を利用することを考えます。`CartesianIndex`は以下のように動作します。

```julia
A = [1 2 3; 4 5 6; 7 8 9]
A[CartesianIndex(1,1)] == 1
A[CartesianIndex(1,2)] == 2
A[CartesianIndex(1,3)] == 3
A[CartesianIndex(2,1)] == 4
A[CartesianIndex(2,2)] == 5
A[CartesianIndex(2,3)] == 6
A[CartesianIndex(3,1)] == 7
A[CartesianIndex(3,2)] == 8
A[CartesianIndex(3,3)] == 9
```

これに加えて、周期境界条件に対応するため、

```julia
function pbc(index, N)
    i = index % N
    i > 0 ? i : i+N
end

pbc.(-3:6, 3) == [3,1,2,3,1,2,3,1,2,3] # true
```

というものを用意しておきます。これらを用いて、

```julia
function bonds(N)
    A = SiteIndeces_2d(N)
    v = [(A[CartesianIndex(i,j)], A[CartesianIndex(i,pbc(j+1,N))]) for i in 1:N for j in 1:N]
    h = [(A[CartesianIndex(i,j)], A[CartesianIndex(pbc(i+1,N),j)]) for i in 1:N for j in 1:N]
    return vcat(v, h)
end
```

とすれば、すべての近接2サイトを結ぶボンドのリストが得られます。

```shell_session
julia> bonds(3)
18-element Array{Tuple{Int64,Int64},1}:
 (1, 2)
 (2, 3)
 (3, 1)
 (4, 5)
 (5, 6)
 (6, 4)
 (7, 8)
 (8, 9)
 (9, 7)
 (1, 4)
 (2, 5)
 (3, 6)
 (4, 7)
 (5, 8)
 (6, 9)
 (7, 1)
 (8, 2)
 (9, 3)
```

これで、2次元正方格子上のIsingハミルトニアンを構成する準備ができました。

```julia
function HIsing_2d(N)
    H = spzeros(2^(N*N), 2^(N*N))
    for b in bonds(N)
        H += ⊗(set_spins(N*N, [b[1],b[2]], [σᶻ,σᶻ])...)
    end
    return H
end
```

# おまけ（もっと直感的に書く方法）
上のような方法でも相当直感的で可読性の高いコードが書けると思いますが、せっかく上付き/下付き添字を使ってコードが書けるのだから、`σᶻ₂`と書くだけで自動的に`σ⁰⊗σᶻ⊗σ⁰...`と解釈してくれないものだろうかという欲望が湧いてきます。同様のことが[「Julia で下付き文字をインデックスに変換するマクロを書いてみた。」](https://qiita.com/antimon2/items/5260dd28b6f049cd453f)で考察されていたので、これを参考にマクロを使って所望の機能が実現できないか考えてみます。

`σᶻ₂`という記号が使いたいからといって、いきなり

```shell-session
julia> σᶻ₂
UndefVarError: σᶻ₂ not defined
```

と書いても、当然そんな記号は定義されていないと怒られてしまいます。しかし例えば、

```shell-session
julia> macro to_str(x) string(:($x)) end
julia> @to_str σᶻ₂
"σᶻ₂"
```

のようにマクロを使って評価することはできます。そこで、

```julia
function sub_to_num(sub)
    if sub == '₁'
        return 1
    elseif sub == '₂'
        return 2
        (省略)
    else
        UndefVarError()
    end
end
# e.g. sub_to_num('₁') == 1

function to_siteindex(array)
    n = length(array)
    sum([array[i]*10^(n-i) for i in 1:n])
end
# e.g. to_siteindex(sub_to_num.(['₁','₂'])) == 12

macro _spin(S, N)
    indeces = []
    spins = []
    for str in split(string(:($S)), "σ")[2:end]
        str = [s for s in "σ"*str]
        index = to_siteindex(sub_to_num.(str[3:end]))
        push!(indeces, index)
        @assert index <= N
        spin = eval(Meta.parse(prod(str[1:2])))
        push!(spins, spin)
    end
    ⊗(set_spins(N, indeces, spins)...)
end
```

のようなことをすれば`@_spin σᶻ₂σᶻ₃ 4`が`σ⁰⊗σᶻ⊗σᶻ⊗σ⁰`と等価になります。ただしこのままでは、2項以上の和を書くときに

```julia
(@_spin σᶻ₁σᶻ₂ 3) + (@_spin σᶻ₂σᶻ₃ 3) + (@_spin σᶻ₃σᶻ₁ 3)
```

などと書かねばならず煩わしいので、この式自身を

```julia
macro spin(S, N)
    str = "(@_spin " * replace(string(:($S)), "+"=>string(N)*")+(@_spin") * " " * string(N)* ")"
    Meta.parse(str)
end
```

というマクロで展開するようにすれば、

```julia
# N=3 Ising
H3 = @spin σᶻ₁σᶻ₂+σᶻ₂σᶻ₃+σᶻ₃σᶻ₁ 3

# N=4 Ising
H4 = @spin σᶻ₁σᶻ₂+σᶻ₂σᶻ₃+σᶻ₃σᶻ₄+σᶻ₄σᶻ₁ 4
```

でイジングハミルトニアンが書けます。

ということで、一応マクロにマクロを重ねれば、コードの見た目を整えることはいくらでも可能と言えそうです。ただ実質は`spin("σᶻ₁σᶻ₂+σᶻ₂σᶻ₃+σᶻ₃σᶻ₁", 3)`というような関数を定義して実行するのと変わらないので、

```julia
spin("σᶻ₁σᶻ₂+σᶻ₂σᶻ₃+σᶻ₃σᶻ₁", 3)
@spin σᶻ₁σᶻ₂+σᶻ₂σᶻ₃+σᶻ₃σᶻ₁ 3
```
のどちらを綺麗に感じるかという趣味の問題です。

# まとめ
Juliaで量子スピン系のハミルトニアンを構成する方法を見てきました。サンプルコードは[Gist](https://gist.github.com/Shoichiro-Tsutsui/c32bedb5f81c42b61db357406b7e66a1)で公開しています。またここで定義した関数たちをパッケージとしてまとめたものを[GitHub](https://github.com/Shoichiro-Tsutsui/QuantumSpin.jl)で公開しています。パッケージを自作したい人にとって参考になると思います。パッケージはブラウザ上で[お試し使用](https://mybinder.org/v2/gh/Shoichiro-Tsutsui/QuantumSpin.jl/HEAD?filepath=jupyter%2Ftutorial.ipynb)することができます。（ただし、ビルドにちょっと時間がかかります。）