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

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

## 多体電子系の厳密対角化：疎行列
ここまで固体物理で登場する様々な問題を解いてきましたが、最後に量子多体問題を厳密に解く方法について紹介します。この最後のトピックスは、Juliaを使って実際の研究の現場で取り扱うような問題をどのように解くか、という例の紹介です。ですので、書かれている内容は若干難しいかもしれませんが、行っていることは線形代数と対角化がメインです。細かな流れを追わなくても構いません。その場合には数式がどのようにコードに落とされていくのか、という観点で読んでみてください。

今回扱う方法は、厳密対角化と呼ばれます。本来であれば電子系の多体問題を扱うには「第二量子化」を用いた生成消滅演算子を用いたハミルトニアンを使うのが便利です。しかし、第二量子化による生成消滅演算子を導入するのにはこの本の趣旨からは外れてしまいます。そこで、生成消滅演算子を「あるルールに従う演算」とみなして、行列として表現することにします。生成消滅演算子を行列として表現できればハミルトニアンも行列として表現できることになり、対角化して固有値を得ることができます。この節では、上で登場した強束縛模型に電子間斥力相互作用が入った模型であるハバード模型と呼ばれる模型の厳密な基底状態の波動関数を数値計算で求めてみます。その際、Juliaの多彩な機能（疎行列、多重ディスパッチ、独自型等）を利用してコードを作成します。多重ディスパッチと独自型を組み合わせることで、他のプログラミング言語でのオブジェクト指向プログラミングと同等の自由度とメンテナンス性のあるプログラミングが可能となります。
### 1サイト問題
電子にはスピンという自由度がありまして、アップスピンとダウンスピンという状態があります。そこで、スピンが上向きの電子（アップスピン電子）の波動関数を$\psi_{\uparrow}(x)$、スピンが下向きの電子（ダウンスピン電子）の波動関数を$\psi_{\downarrow}(x)$としましょう。さらに、連続的な空間だと大変ですので、有限の数の格子点を考えます。これは、上で登場した強束縛模型の考え方で、一つの格子点が一つの原子に対応していると考えることができます。
格子点が一つしかない時、 格子点1つにアップスピン状態とダウンスピン状態がありますから、電子は最大2個入ることができます。
電子が2個ある時、この2電子の波動関数を$\phi(\sigma_1,\sigma_2)$とします。$\sigma_1$と$\sigma_2$は$\uparrow$か$\downarrow$の2通りを取ることができます。しかし、電子はフェルミオンですので、粒子の区別をつけることができず、同じ量子状態を取ることもできませんので、$\phi(\sigma_1,\sigma_2)$は実は$\phi(\uparrow,\downarrow)$という1種類の状態しか取ることができません。この時のシュレーディンガー方程式は
```math
\left( \epsilon_{\uparrow} + \epsilon_{\downarrow} + U \right) \phi(\sigma_1,\sigma_2) = E \phi(\uparrow,\downarrow)
```


となります。ここで、$\epsilon_{\sigma}$は相互作用がない時のアップスピンあるいはダウンスピンが持つエネルギーです。そして、$U$はクーロン相互作用です。アップスピンとダウンスピンが同じ格子点にある場合、お互いにクーロン斥力の相互作用が働きます。
電子が1個ある時、1電子の波動関数を$\phi(\sigma_1)$とします。電子は1個しかありませんので、$\sigma_1$は$\uparrow$と$\downarrow$のどちらも取ることができます。もちろん、相互作用する相手がいませんのでクーロン相互作用も効きません。この時のシュレーディンガー方程式は
```math
\left( \epsilon_{\uparrow} \delta_{\sigma_1,\uparrow} + \epsilon_{\downarrow} \delta_{\sigma_1,\downarrow} \right) \phi(\sigma_1) = E \phi(\sigma_1)
```
です。
電子が0個の時、波動関数はありませんが、便宜上$\psi_0$という波動関数を考えてみます。この時のシュレーディンガー方程式は
```math
\left(0 \right) \psi_0 = E \psi_0
```
となります。
電子数ごとにシュレーディンガー方程式を立ててみましたが、$\epsilon_{\uparrow}$など同じものが出てきていますので、一つの方程式にまとめることはできないのでしょうか。ここで${\vec \phi} = (\psi_0, \phi(\uparrow),\phi(\downarrow),\phi(\uparrow,\downarrow))^{T}$というベクトルを考えますと、
```math
\left( \begin{array}{cccc}
 0 & 0 & 0 & 0 \\
 0 & \epsilon_{\uparrow} & 0 & 0 \\
0 & 0 & \epsilon_{\downarrow} & 0 \\
 0 & 0 & 0 & \epsilon_{\uparrow} + \epsilon_{\downarrow} + U 
 \end{array} \right) {\vec \phi} = E {\vec \phi}
```
という形になりまして、$4 \times 4$の行列に関する固有値問題になります。この行列は対角行列ですから対角化をするまでもなく固有値が求まっており、あまり行列で表示した恩恵がありませんね。
ここで登場したベクトル${\vec \phi}$ですが、今後サイト数が増えたときにわかりやすくするために
```math
\psi_0 \rightarrow |00 \rangle
```
```math
\phi(\uparrow) \rightarrow |01 \rangle
```
```math
\phi(\downarrow) \rightarrow |10 \rangle
```
```math
\phi(\uparrow,\downarrow) \rightarrow |11 \rangle
```
というように0と1で表すことにします。このように書くと波動関数を二進法の数字一つで表すことができます。今は電子が最大2個ですので二進法の桁数が2になっています。


### 生成消滅演算子の行列表示
1サイトしかない問題は数値計算するまでもないほぼ自明な問題になっていました。本当に考えたい問題はサイト数が多い問題です。そのような問題のシュレーディンガー方程式を行列の形で表現すれば数値計算で解くことができます。しかし、上でやったように毎回シュレーディンガー方程式を書き下してその行列の表示を求めるのは面倒です。何か良い方法はないでしょうか。
ここで、生成消滅演算子というものを導入することで面倒な手続きを回避してみましょう。生成消滅演算子とは生成演算子と消滅演算子という二種類の演算子のことです。生成演算子というのは電子を生成させる演算子で、消滅演算子というのは電子を消滅させる演算子のことです。場の量子論について触れたことがない方は電子が生成する消滅すると言われてもピンとこないかもしれません。しかしここでは場の量子論について知っていなくても問題ありません。ここでは天下り的に生成消滅演算子を定義します。生成消滅演算子はフェルミオンの性質を記述するために導入された演算子でして、「同じ状態には一つしか電子は占有できない」と「粒子を入れ替えた場合には波動関数にマイナスの符号がつく」と性質を満たすように作られています。例えば、上で出てきた2電子状態の波動関数$\phi(\uparrow,\downarrow)$の電子の位置を入れ替えると$\phi(\downarrow,\uparrow)$となりますが、これは$\phi(\uparrow,\downarrow) = - \phi(\downarrow,\uparrow)$となるような波動関数になっています。今は2電子までしか考えていませんから粒子の入れ替えは単純ですが、複数の電子がいる場合には任意のどの電子2つを選んで入れ替えてもマイナスの符号をつくように波動関数を作る必要があります。これは非常に面倒ですが、生成消滅演算子を使うことでかなり見通しよく扱いやすくなります。
あるスピン状態$\sigma$の電子に関する生成演算子$c_{\sigma}^{\dagger}$と消滅演算子$c_{\sigma}$というものを考えます。この二つの演算子は
```math
[c_{\sigma},c_{\sigma'}^{\dagger}]_+ = \delta_{\sigma \sigma'}
```
```math
[c_{\sigma},c_{\sigma'}]_+ =0
```
```math
[c_{\sigma}^{\dagger},c_{\sigma'}^{\dagger}]_+ =0
```
という条件を満たしているとします。ここで、
```math
[A,B]_{+} \equiv AB + BA
```
は反交換関係と呼ばれます。これらの演算子を使って、上で定義した4つの状態を表すことにします。
まず、電子が何もいない状態$|00 \rangle$を考えます。この状態にアップスピンに関する生成演算子$c_{\uparrow}^{\dagger}$を作用させると$|01 \rangle$という状態になるとします($c_{\uparrow}^{\dagger} |00 \rangle = |01 \rangle$)。また、ダウンスピンに関する生成演算子$c_{\downarrow}^{\dagger}$を作用させると$|10 \rangle$になるとします($c_{\downarrow}^{\dagger} |00 \rangle = |10 \rangle$)。また、$|11 \rangle$という状態は、$|00 \rangle$に$c_{\uparrow}^{\dagger}$を作用させてから$c_{\downarrow}^{\dagger}$を作用させたら得られるとします（$c_{\downarrow}^{\dagger} c_{\uparrow}^{\dagger} |00 \rangle = |11 \rangle$）。これで三つの状態を$|00 \rangle$から定義することができました。これらの定義の他に、消滅演算子$c_{\sigma}$を$|00 \rangle$に作用させた場合には0になることも約束しておきます（$c_{\sigma} |00 \rangle = 0$）。
これらの定義と上の反交換関係を用いると、消滅演算子の挙動が決定されます。例えば、$|01 \rangle$という状態に消滅演算子$c_{\uparrow}$を作用させると$c_{\uparrow} c_{\uparrow}^{\dagger} |00 \rangle = (1 - c_{\uparrow}^{\dagger} c_{\uparrow}) |00 \rangle  = |00 \rangle$となり、$|01 \rangle$からアップスピンの電子が消滅して$|00 \rangle$になります。また、同じ演算子が二つついた場合には上の反交換関係から$c_{\sigma}^{\dagger}c_{\sigma}^{\dagger}+c_{\sigma}^{\dagger}c_{\sigma}^{\dagger} = 2 c_{\sigma}^{\dagger}c_{\sigma}^{\dagger} = 0$となりますので、$c_{\uparrow}^{\dagger}|01 \rangle = 0$となります。
これでそれぞれの波動関数を生成消滅演算子を使って表すことができました。2電子状態の波動関数は$|11\rangle = c_{\downarrow}^{\dagger} c_{\uparrow}^{\dagger} |00 \rangle$と定義されていますが、この波動関数の2電子のスピンを入れ替えると$c_{\uparrow}^{\dagger} c_{\downarrow}^{\dagger} |00 \rangle$となります。反交換関係$c_{\downarrow}^{\dagger} c_{\uparrow}^{\dagger} + c_{\uparrow}^{\dagger} c_{\downarrow}^{\dagger} =0$を使うと、$c_{\uparrow}^{\dagger} c_{\downarrow}^{\dagger} |00 \rangle = - c_{\downarrow}^{\dagger} c_{\uparrow}^{\dagger} |00 \rangle = - |11\rangle$となりまして、ちゃんと2つの電子の入れ替えをしたときにマイナスの符号が出ます。今は2電子までしか考えていませんが、3電子、4電子と増やしていっても、生成消滅演算子で波動関数を書いておけば、自動的にフェルミオンの性質を満たす波動関数になるため便利です。
次に、生成消滅演算子を行列で表現してみましょう。行列で表現することで数値計算で簡単に扱うことができる形になります。今、状態は$|00\rangle$、$|01\rangle$、$|10\rangle$、$|11\rangle$の4つあります。波動関数がこれらの線型結合:
```math
|\psi \rangle = a_{00} |00\rangle + a_{01} |01\rangle + a_{10} |10\rangle + a_{11} |11 \rangle
```

で書かれているとしますと、この波動関数に生成演算子$c_{\uparrow}^{\dagger}$をそれぞれかけてみますと、$c_{\uparrow}^{\dagger}|00\rangle = |01 \rangle$、$c_{\uparrow}^{\dagger}|01\rangle = 0$、$c_{\uparrow}^{\dagger}|10\rangle = c_{\uparrow}^{\dagger} c_{\downarrow}^{\dagger}|00 \rangle = -|11 \rangle$、$c_{\uparrow}^{\dagger}|11\rangle = 0$となりますから、
```math
c_{\uparrow}^{\dagger} |\psi \rangle = a_{00} |01\rangle  - a_{10} |11\rangle 
```
となります。よって、
```math
C_{\uparrow}^{\dagger} 
\left( \begin{array}{c}
  a_{00} \\
  a_{01} \\
  a_{10} \\
  a_{11} 
 \end{array} \right)
 =
 \left( \begin{array}{c}
  0 \\
  a_{00} \\
  0 \\
  - a_{10} 
 \end{array} \right)
= 
\left( \begin{array}{cccc}
 0 & 0 & 0 & 0 \\
 1 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\
 0 & 0 & -1 & 0
 \end{array} \right) 
 \left( \begin{array}{c}
  a_{00} \\
  a_{01} \\
  a_{10} \\
  a_{11} 
 \end{array} \right)
```

となり、生成演算子の行列表示$C_{\uparrow}^{\dagger}$が求まります。同様に消滅演算子$c_{\uparrow}$をかけてみると、$c_{\uparrow}|00\rangle = 0$、$c_{\uparrow}|01\rangle = |00\rangle$、$c_{\uparrow}|10\rangle = 0$、$c_{\uparrow}|11\rangle = c_{\uparrow} c_{\downarrow}^{\dagger} c_{\uparrow}^{\dagger} |00 \rangle = - |10 \rangle$となりますので、その行列は
```math
C_{\uparrow}
\left( \begin{array}{c}
  a_{00} \\
  a_{01} \\
  a_{10} \\
  a_{11} 
 \end{array} \right) = 
\left( \begin{array}{cccc}
 0 & 1 & 0 & 0 \\
 0 & 0 & 0 & 0 \\
0 & 0 & 0 & -1 \\
 0 & 0 & 0 & 0
 \end{array} \right)
 \left( \begin{array}{c}
  a_{00} \\
  a_{01} \\
  a_{10} \\
  a_{11} 
 \end{array} \right)
```

となります。Juliaで書くと
```julia 
julia> Cdag = [0 0 0 0
       1 0 0 0 
       0 0 0 0
       0 0 -1 0 ]
4×4 Matrix{Int64}:
 0  0   0  0
 1  0   0  0
 0  0   0  0
 0  0  -1  0
julia> C=[0 1 0 0 
       0 0 0 0 
       0 0 0 -1
       0 0 0 0]
4×4 Matrix{Int64}:
 0  1  0   0
 0  0  0   0
 0  0  0  -1
 0  0  0   0
```
となりますね。行列にすると$[C_{\uparrow}]^{\dagger} = C_{\uparrow}^{\dagger}$となっていることがはっきりわかります。また、$C_{\uparrow}^{\dagger} C_{\uparrow}$は
```julia 
julia> Cdag*C
4×4 Matrix{Int64}:
 0  0  0  0
 0  1  0  0
 0  0  0  0
 0  0  0  1
```
となりますから、
```math
C_{\uparrow}^{\dagger} C_{\uparrow} =\left( \begin{array}{cccc}
 0 & 0 & 0 & 0 \\
 0 & 1 & 0 & 0 \\
0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 1
 \end{array} \right)
```

となりますね。ダウンスピンの生成演算子$c_{\uparrow}^{\dagger}$も同様にやってみると、$c_{\downarrow}^{\dagger}|00\rangle = |10 \rangle$、$c_{\downarrow}^{\dagger}|01\rangle = |11 \rangle$、$c_{\downarrow}^{\dagger}|10\rangle = 0$、$c_{\downarrow}^{\dagger}|11\rangle = 0$となりますから、
```math
C_{\downarrow}^{\dagger} \left( \begin{array}{c}
  a_{00} \\
  a_{01} \\
  a_{10} \\
  a_{11} 
 \end{array} \right)  = 
\left( \begin{array}{cccc}
 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 \\
1 & 0 & 0 & 0 \\
 0 & 1 & 0 & 0
 \end{array} \right) \left( \begin{array}{c}
  a_{00} \\
  a_{01} \\
  a_{10} \\
  a_{11} 
 \end{array} \right) 
```
と書くことができまして、Juliaでは
```julia
julia> Cdagdown=[0 0 0 0
       0 0 0 0
       1 0 0 0
       0 1 0 0]
4×4 Matrix{Int64}:
 0  0  0  0
 0  0  0  0
 1  0  0  0
 0  1  0  0
```

となります。ですので、$C_{\downarrow}^{\dagger}C_{\downarrow}$は行列で書くと
```math
C_{\downarrow}^{\dagger}C_{\downarrow} = 
\left( \begin{array}{cccc}
 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 \\
0 & 0 & 1 & 0 \\
 0 & 0 & 0 & 1
 \end{array} \right) 
```
となります。
これらの生成消滅演算子の組$C_{\sigma}^{\dagger}C_{\sigma}$の行列表示をみてわかることは、スピン状態$\sigma$がいる状態の場合には1を返し、いない場合には0を返す、という演算になっているということです。したがって、アップスピンの電子がいるときにはエネルギーが$\epsilon_{\uparrow}$増え、アップスピンとダウンスピンの両方の電子がいるときには相互作用エネルギー$U$が加えられる、という表現をこの行列を使って表すことができます。つまり、ハミルトニアンも生成消滅演算子（あるいはその行列表示）で表すことができるということです。
ハミルトニアンを
```math
H = \sum_{\sigma} \epsilon_{\sigma} c_{\sigma}^{\dagger} c_{\sigma} + U c_{\uparrow}^{\dagger} c_{\uparrow} c_{\downarrow}^{\dagger} c_{\downarrow}
```
と書いてみましょう。生成消滅演算子の行列表示を代入して整理すると、上で登場した$4 \times 4$行列のハミルトニアンになることがわかります。生成消滅演算子の行列表示がわかっていれば、その組み合わせで書けるハミルトニアンは必ず行列で表現できます。行列で表現できれば固有値問題をJuliaで解ける行列の問題として定義することができます。ここで扱ったのは1サイトの問題なのでご利益がわかりにくいかもしれませんので、次の節ではサイト数を増やしてみましょう。


### 2サイト問題
次は格子点が2つある問題を考えてみます。1つのサイトにはアップスピンとダウンスピンの状態があり、アップスピンに入る、ダウンスピンに入る、両方に入る、両方とも入らない、の4状態がありました。2つのサイトの場合にはそれぞれに4状態がありますから、$4^2 = 16$個の状態があります。状態を書き下してみると$|0000\rangle$、$|0001\rangle$、$|0010\rangle$、$|0011\rangle$、$|0100\rangle$、$|0101\rangle$、$|0110\rangle$、$|0111\rangle$、$|1000\rangle$、$|1001\rangle$、$|1010\rangle$、$|1011\rangle$、$|1100\rangle$、$|1101\rangle$、$|1110\rangle$、$|1111\rangle$の16状態です。二進法で0から15までを4桁で書くことでこれらの状態を表現できます。ここでは、左の2つをダウンスピン、右の2つをアップスピンの状態を表しているとしましょう。つまり、$|1111\rangle = c_{2,\downarrow}^{\dagger}c_{1,\downarrow}^{\dagger}c_{2,\uparrow}^{\dagger}c_{1,\uparrow}^{\dagger}|0000 \rangle$という形で状態を定義しておきます。例えば、$|1001 \rangle$は$|1001 \rangle = c_{2,\downarrow}^{\dagger}c_{1,\uparrow}^{\dagger}|0000 \rangle$です。
1サイトの時と同じように生成消滅演算子を行列として表現してみましょう。行列表示する時のポイントは、作用させようとする生成演算子の位置です。例えば、1サイトの時には2電子状態は$|11 \rangle = c_{\downarrow}^{\dagger} c_{\uparrow}^{\dagger} |00\rangle$と書けていましたが、これにアップスピンに関する消滅演算子$c_{\uparrow}$を作用させると、$c_{\uparrow} |11 \rangle = c_{\uparrow} c_{\downarrow}^{\dagger} c_{\uparrow}^{\dagger} |00\rangle = - c_{\downarrow}^{\dagger} c_{\uparrow} c_{\uparrow}^{\dagger} |00\rangle$のようにマイナス符号が出ます。このマイナス符号は反交換関係$c_{\uparrow} c_{\downarrow}^{\dagger} + c_{\downarrow}^{\dagger} c_{\uparrow} =0$によって生じています。ある状態に演算子を作用させるときは、必ず一番左にその演算子がありますから、その演算子を反交換関係によって移動させる必要があります。この移動が奇数回であればマイナスが、偶数回であればプラスの符号がつきます。この符号をJuliaで計算できれば、生成消滅演算子を行列で表現できるようになります。
以後は$N$サイト問題を取り扱えるように一般的なコードを書いてみましょう。


### 独自型の利用とNサイト問題
それでは、Juliaでコードを書いていきましょう。まず、状態を独自の型で定義してしまいます:
```julia 
mutable struct Electronbasis{N}
    state::Array{Bool,1}
    index::Int64
    function Electronbasis(N,index)
        state = i2state(index,N)
        return new{N}(state,index)
    end
end

function i2state(i,N)
    @assert 0 <= i <  2^N "index should be in  0 <= i < $(2^N)!"
    state = zeros(Bool,N)
    bit = string(i, base=2,pad=N)
    for i=1:N
        state[i]= parse(Bool,bit[N-i+1])
    end
    return state
end
```
ここで、```Electronbasis{N}```のように```{}```がついている場合は、この型はパラメトリック型と呼ばれます。パラメトリック型を用いますと、型がどんな型かを指定することができます。実はこれまで登場してきた配列```Array```にも```Array{Bool,1}```のような形で```{}```がついており、```Array{Bool,1}```であれば要素の型が```Bool```で```1```次元、```Array{Float64,2}```であれば要素の型が```Float64```で```2```次元、という意味になります。今回独自定義した```Electronbasis{N}```では、```{N}```の```N```がパラメータです。今回の目的では```N```はサイト数が入ることを想定しています。この型の中にある```state```は状態の0と1の並びを保持しておく変数、```index```は状態のラベルを表します。この型の中に入っている```function```の```Electronbasis(N,index)```はコンストラクタと呼ばれるものでして、どうやってこの型を作るか、ということを指定しています。今回の場合は、引数として```N```と```index```を取る関数でこの型を定義します、と言っています。この関数の中には```i2state(index,N)```という関数が使われていますので、その挙動を説明します。```@assert　0 <= i <  2^N```は引数```i```がこの範囲に入っていないとエラーメッセージ```"index should be in  0 <= i < $(2^N)!"```を返すマクロです。想定していない数字が入ってきた時にエラーで止める役割をしています。```string(i, base=2,pad=N)```は```i```という10進数の数字を桁が```N```の2進数の文字列に変換する関数:
```julia 
julia> a = string(12, base=2,pad=5)
"01100"
```
です。あとは```parse(Bool,bit[N-i+1])```を使って、右から```i```番目の文字が0か1かを判定して```state[i]```に格納しています。なお、```Bool```型は```true```か```false```を格納する型なのですが、```true```は1、```false```は0として言い直すことができるようになっています。
次に、作った```Electronbasis```型を表示する関数を作りましょう。表示する関数としてdisplayという関数がJuliaにはあります。REPLで
```julia 
julia> a = [1 2
       3 4]
2×2 Matrix{Int64}:
 1  2
 3  4
```
のように行列を作るとみやすい形で表示してくれます。これは```display(a)```という関数を呼んでいます。ですので、REPLではなく```julia  hoge.jl```の形で呼び出した時には
```julia 
a = [1 2
    3 4]
display(a)
```
のようにすることで、行列をいい感じに表示してくれます。少し気をつけるべき点としては、```display(a)```としたときには最後に改行されませんので、改行をしたければ```println("")```等を次の行に追加してください。
```display(a)```は変数```a```の型に応じていい感じに出力してくれる関数ですが、今自分で定義した```Electronbasis```型にはもちろん対応していません。対応させるには、多重ディスパッチを用いて```display(a)```の```a```が```Electronbasis```型の時の挙動を定義すればよいです。Juliaの多重ディスパッチの利点として、すでにどこかで定義されている関数に新しい型の挙動を入れることが容易であることがあげられます。```display```という関数は何もパッケージを読み込まなくても使える関数ですが、このような関数は```Base```というパッケージに入っていることを意味しています。ですので、
```julia 
function Base.display(electron::Electronbasis{N}) where N
    label = get_label(electron)
    println(label)
    return
end
function get_label(electron::Electronbasis{N}) where N
    label=""
    label*="|"
    for i=N:-1:1
        if electron.state[i]
            label*="1"
        else
            label*="0"
        end
    end
    label*=">"    
end
```
を書きます。これで、```display```という関数の引数に```Electronbasis```型が来た時の挙動が定義されました。なお、ここで```Electronbasis{N}```とすることで```Electronbasis```型がパラメトリック型であることを明示しています。```where N```というのは```N```がどのような型かを指定する文法です。例えば```where N <: Int```とすると```N```が整数型の時だけこの```display```の定義が使われるようにすることができます。```where N```と何も型の情報を指定していない場合には、```N```はなんでもよい（```Any```型）ことを意味します。パラメトリック型の面白いところは、ここで登場した```N```を関数内で普通に使えることです。例えばこの関数では```for i=N:-1:1```とすることで```for```文の範囲を```N```を使って指定しています。
この```display```関数を使いますと、
```julia 
N = 2*2
index = 3
state = Electronbasis(N,index)
display(state)
```
の出力結果は```|1100>```となります。これでどんな状態を取り扱っているかがすぐにわかるようになりました。


さて、次は```Electronbasis```型に色々な機能をつけていきましょう。Pythonであればクラスに色々なメソッドをつけるのと同じような感覚で独自型に関する関数を追加していくことができます。まず、あるサイト$i$に電子がいるかどうかの関数を
```julia 
function check_site(electron::Electronbasis,isite)
    return electron.state[isite]
end
```
と定義します。```Electronbasis```のフィールドの```state```は```Array{Bool,1}```という型ですので、要素は```true```か```false```をとります。```isite```番目のサイト(1サイトの系であればアップスピンの電子を1番目、ダウンスピンの電子を2番目と呼ぶことにします)の電子を消滅させる消滅演算子を作用させたときに符号がどうなるかを判定する関数を
```julia 
function check_sign(electron::Electronbasis{N},isite) where N
    sign = 1
    for jsite=N:-1:isite+1
        if check_site(electron,jsite)
            sign *= -1
        end
    end
    return sign
end
```
で定義します。ここで、```for jsite=N:-1:isite+1```のループは、```isite```より左側にあるサイトのうち、電子がいるサイトを数えて符号を決定しています。そして、```isite```番目のサイトの電子を消滅させる関数を
```julia 
function erace!(electron::Electronbasis,isite)
    @assert check_site(electron,isite) "there is no electron at $(isite)-site!"
    electron.state[isite] = false
    electron.index = state2i(electron.state)
end
```
と定義します。```@assert check_site(electron,isite)```というのは、消滅させる電子がちゃんといるかどうかを判定している箇所で、電子がそのサイトにいないのにもかかわらず消滅演算子を作用させようとした場合にはエラーメッセージを出して止まるようにしています。```electron.state[isite] = false```は```isite```番目の電子をいないことにしており、```electron.index = state2i(electron.state)```は新しいラベルを計算しています。この```state2i```は状態が与えられた時にラベルを計算する関数で、
```julia 
function state2i(state)
    N = length(state)
    ii = 0
    for i=1:N
        ii += state[i]*2^(i-1)
    end
    return ii
end
```
で定義しました。これで```isite```番目のサイトの消滅演算子を行列として作成する準備が整いました。行列を作成する関数を
```julia 
function make_C(isite,nsite)
    N = 2*nsite
    ndim = 4^nsite
    mat_C = zeros(Int8,ndim,ndim)
    for index=0:ndim-1 #消した先の状態のラベル
        electron = Electronbasis(N,index)
        if check_site(electron,isite) #isiteに電子がいるかどうかチェック
            sign = check_sign(electron,isite) #消滅演算子の符号の計算
            erace!(electron,isite)
            jindex = electron.index #消えた状態のラベル
            mat_C[jindex+1,index+1] = sign #行列要素を決める
        end
    end
    return mat_C
end
```
と定義すればOKです。ここで、消滅演算子は$c_{\uparrow} |01 \rangle = |00 \rangle$のような形で定義されているので、$|01 \rangle$の係数$a_{01}$が$|00 \rangle$の成分に移る形になっています。したがって、$|00 \rangle$は```jindex```というラベル、$|01 \rangle$が```index```というラベルになっています。生成演算子は消滅演算子の転置で計算できますから、これでハミルトニアンを行列表示する準備が整ったことになります。

part2に続きます。