初めての計算物理：量子力学と線形代数をpythonで。

初めての計算物理：物性物理をpythonで

はじめに

量子力学を学んでいくと、線形代数が量子力学をよく表現することを知るだろう。
ハミルトニアンは行列、エネルギーは固有値、固有状態は固有ベクトルで表現される。
そして、シュレディンガー方程式は永年方程式に帰着する。
物性理論の研究室では、その方程式を計算機を用いて解き、物性予測をしている（もちろん他のこともやっている）。

しかし、大学の授業で計算機を用いるのは、研究室に入ってからであろう。
研究室配属前に「実験器具は使うが、計算機は使わない」という現状で、物性理論のことを知ってもらう役に立てばと思う。

目的

今回の記事では、座学で学ぶような「解ける」問題と、計算機を用いて「数値的に解く」問題との橋渡しをすることを目的とする

目次

1. 教育的な例の確認「スピン演算子$\hat{S_z}$とその固有値」（知っている人は飛ばしてＯＫ）
2. 教育的な例を数値的に解く
3. 想定されるシチュエーション

$$
\def\bra#1{\mathinner{\left\langle{#1}\right|}}
\def\ket#1{\mathinner{\left|{#1}\right\rangle}}
\def\braket#1#2{\mathinner{\left\langle{#1}\middle|#2\right\rangle}}
\def\bbraket#1#2#3{\mathinner{\left\langle{#1}\middle|#2\middle|#3\right\rangle}}
$$

1. 教育的な例の確認「スピン演算子Szとその固有値」

スピン演算子$\hat{S_z}$とその固有状態の関係は次のように与えられる。

\hat{S_z} \ket{\uparrow} = \frac{\hbar}{2} \ket{\uparrow} ,~~~~~\hat{S_z} \ket{\downarrow} = -\frac{\hbar}{2} \ket{\downarrow}

これらを行列で表してみると、

\frac{\hbar}{2}
\begin{pmatrix}
1 & 0 \\
0 & -1 
\end{pmatrix} 
\begin{pmatrix}
1 \\
0
\end{pmatrix}
= \frac{\hbar}{2} \begin{pmatrix}
1 \\
0
\end{pmatrix},~~~~~
\frac{\hbar}{2}
\begin{pmatrix}
1 & 0 \\
0 & -1 
\end{pmatrix} 
\begin{pmatrix}
0 \\
1
\end{pmatrix}
= \frac{-\hbar}{2} \begin{pmatrix}
0 \\
1
\end{pmatrix}

のように書き換えられる。この表示から明らかなように、

\begin{pmatrix}
a \\
b
\end{pmatrix} = a \ket{\uparrow} + b \ket{\downarrow}

となっており、$\ket{\uparrow}$が実際にどんな関数になっているかどうかはどうでもよく、その係数だけ計算されている。これは、「正規直交基底を用いれば、任意の新しい基底を表現できる」ことに基づいている。

さて、$\hat{S_z}$の固有状態は$\hat{S_y}$の固有状態だろうか？もちろん違う。では、$\ket{\uparrow}$に$\hat{S_y}$を作用させるとどうなるだろうか？

これは昇降演算子から理解することができる。つまり、次の式

S^{+}\ket{\uparrow} = 0,~~~~~ S^{+}\ket{\downarrow} = \hbar \ket{\uparrow}\\
S^{-}\ket{\uparrow} = \hbar \ket{\downarrow},~~~~~S^{-} \ket{\downarrow} = 0

と、昇降演算子の定義
$$
S^{+} = \hat{S_x}+ i\hat{S_y},~~~~~ S^{-}= \hat{S_x} - i\hat{S_y}
$$
から、$\hat{S_x},\hat{S_y}$について逆に解けば、

\hat{S_x}\ket{\uparrow} = \frac{\hbar}{2} \ket{\downarrow}, ~~~~~\hat{S_x}\ket{\downarrow} = \frac{\hbar}{2} \ket{\uparrow}\\
\hat{S_y} \ket{\uparrow} = i\frac{\hbar}{2}\ket{\downarrow}, ~~~~~\hat{S_y} \ket{\downarrow} = -i\frac{\hbar}{2} \ket{\uparrow}

が得られる。確かに$\ket{\uparrow}$は$\hat{S_y}$の固有状態になっていなかった。

では、$\hat{S_y}$の固有状態は何だろうか？予測は付くだろうが、敢えてシステマティックに求めてみよう。まずは$\hat{S_y}\ket{\uparrow} = i\frac{\hbar}{2} \ket{\downarrow}$の行列表示は、

\frac{\hbar}{2}\begin{pmatrix}
0 & -i\\
i & 0
\end{pmatrix}
\begin{pmatrix}1\\0\end{pmatrix} = i\frac{\hbar}{2}\begin{pmatrix}0\\1\end{pmatrix}

である。ここで左辺の行列はパウリ行列の$\sigma_y$になっている。この$\sigma_y$を対角化することで、$\hat{S_y}$の固有状態を求める。$\sigma_y$を対角化をすると、固有値と固有ベクトルは、

\lambda_1=1, ~~ u_1=\frac{1}{\sqrt{2}}\begin{pmatrix}1 \\ i \end{pmatrix} ・・・(1)\\
\lambda_2=-1, ~~ u_2=\frac{1}{\sqrt{2}}\begin{pmatrix}-1\\i\end{pmatrix}・・・(2)

となる。実際に固有状態になっているか確かめると、

\hat{\sigma}_yu_1=\begin{pmatrix}
0 & -i\\
i & 0
\end{pmatrix}\frac{1}{\sqrt{2}}\begin{pmatrix}1 \\ i \end{pmatrix}=
\frac{1}{\sqrt{2}}\begin{pmatrix}-i^2\\i\end{pmatrix} = 
\frac{1}{\sqrt{2}}\begin{pmatrix}1\\i\end{pmatrix} = \lambda_1\cdot u_1

となり、固有状態になっていることが確かめられた。

この固有状態$u_1$を用いて、$\sigma_x$の期待値を計算してみる。

\bbraket{u_1}{\sigma_x}{u_1} = \frac{1}{\sqrt{2}}\begin{pmatrix}1 & -i \end{pmatrix}
\begin{pmatrix}
0 & 1 \\
1 & 0 
\end{pmatrix}
\frac{1}{\sqrt{2}}\begin{pmatrix}1\\i\end{pmatrix}
=0

よって、期待値は0である。

2. 教育的な例を数値的に解く

1.の解析計算をpythonで実装する。

1. $\sigma_y$行列を定義
2. 対角化
3. 出力して、前節の解析解と一致するか確認

import numpy as np
sigma_y = np.array([
    [0,-1j],
    [1j,0]
])
eigenvalues, eigenvectors = np.linalg.eigh(sigma_y)
print(eigenvalues[0], end="  ") 
print(eigenvectors[0])
# output : -1.0  [-0.70710678+0.j -0.70710678+0.j]

一致しない。これは固有ベクトルの格納のされ方に依るものである。以下の操作をするとよく分かる。

print(eigenvectors @ sigma_y @ np.linalg.eigh(eigenvectors))
# output : [[-1           0.+1e-16j]
#           [ 0.+1e-16    1        ]]

つまり、np.linalg.eigh()で得られる固有ベクトルは２つの固有ベクトルを縦に並べた形になっていて、

eigenvectors[0] =[u1(x1),u2(x1)]
eigenvectors[1] =[u1(x2),u2(x2)]

となる。というわけで、np.linalg.eigh()で得られた固有ベクトルは転置を取る必要がある。

import numpy as np
sigma_y = np.array([
    [0,-1j],
    [1j,0]
])
eigenvalues, eigenvectors = np.linalg.eigh(sigma_y)
evs = np.transpose(eigenvectors)
print(eigenvalues[0], end="  ") 
print(evs[0])
# output : -1.0   [-0.71+0.j     0.+0.71j]  ( 1/sqrt(2)=0.71 )

これは解析的に計算された$\sigma_y$の(2)の方の固有値と固有ベクトルに一致する。

3. 想定されるシチュエーション

行列の対角化は物理の様々な状況で行われる。ここでは、一例として、「スピンありの１次元２サイトモデルで基底状態を求める」ということをやってみる。

ハミルトニアンとして、

H=H_{kin}+H_{\mu}+H_{U}\\
H_{kin} = -t\sum_{\sigma}(c^{*}_{2,\sigma}c_{1,\sigma}+c^{*}_{1,\sigma}c_{2,\sigma})\\
H_{\mu} =  \sum_{\sigma}\sum_{i=1,2}(-\mu_i)c^*_{i,\sigma}c_{i,\sigma}\\
H_{U} = U\sum_{i=1,2}n_{i,\uparrow}n_{i,\downarrow}

を考え、状態としては「粒子数とスピンが保存し、アップスピン１個、ダウンスピン１個のハーフフィリング」を考える。この時、状態は４つ考えられる。

\psi_1 = |(\uparrow\downarrow)_1,(0)_2>\\
\psi_2 = |(0)_1,(\uparrow\downarrow)_2>\\
\psi_3 = |(\uparrow)_1,(\downarrow)_2>\\
\psi_4 = |(\downarrow)_1,(\uparrow)_2>

それぞれの状態は規格化されているとする。

$\psi_1$は次のような状態である。

これらの状態を用いて表される状態を$\phi$とし、次のように表すことにする。もちろんこの$\phi$は全ての状態を表現することができる。

\phi=a\psi_1+b\psi_2+c\psi_3+d\psi_4=\begin{pmatrix}a \\ b \\ c \\ d \end{pmatrix}

ただし、$|a|^2+|b|^2+|c|^2+|d|^2=1$。このとき、

H\phi=\begin{pmatrix}
-2\mu_1+U & 0        & -t & -t \\
    0     &-2\mu_2+U & -t & -t \\
-t & -t & -\mu_1-\mu_2 & 0 \\
-t & -t & 0 & -\mu_1-\mu_2
\end{pmatrix}\phi

のようにハミルトニアンを行列表示することができる。この行列を対角化することで、固有状態となる$\phi$を求めよう。やることは第２節と一緒である。

import numpy as np
#計算条件
mu_1 = 0.1
mu_2 = 0.4
t = 1.0
U = 3.0

H = np.array([
    [-2*mu_1 + U,       　   0,          -t,            -t],
    [          0,  -2*mu_2 + U,          -t,             -t],
    [         -t,　         -t, -mu_1 - mu_2,            0],
    [         -t,　         -t,            0, -mu_1 - mu_2]
])
eigenvalues, eigenvectors = np.linalg.eigh(H)
evs = np.transpose(eigenvectors)
for i in range(4):
    print(eigenvalues[i], end="  ") 
    print(evs[i])

# output
# -1.51  [ 0.29    0.34   0.63   0.63 ]
# -0.50  [ 1e-17  1e-17  -0.71   0.71 ]
#  2.44  [ 0.54   -0.83   0.10   0.10 ]
#  3.57  [-0.79   -0.44   0.30   0.30 ]

計算結果を見てみると、固有値が一番小さな状態、つまり基底状態は、$[0.290.340.63~~0.63]$である。クーロン力が大きい計算条件なので、$\psi_3, \psi_4$が多く含まれている。$\psi_1, \psi_2$の混ざり方は、サイトのエネルギーが違うので１：１ではない。