1. ovation6868lx

    No comment

    ovation6868lx
Changes in body
Source | HTML | Preview
@@ -1,271 +1,271 @@
# 初めての計算物理:物性物理をpythonで
## はじめに
量子力学を学んでいくと、線形代数が量子力学をよく表現することを知るだろう。
ハミルトニアンは行列、エネルギーは固有値、固有状態は固有ベクトルで表現される。
そして、シュレディンガー方程式は永年方程式に帰着する。
物性理論の研究室では、その方程式を計算機を用いて解き、物性予測をしている(もちろん他のこともやっている)。
しかし、大学の授業で計算機を用いるのは、研究室に入ってからであろう。
研究室配属前に「実験器具は使うが、計算機は使わない」という現状で、物性理論のことを知ってもらう役に立てばと思う。
### 目的
今回の記事では、座学で学ぶような「解ける」問題と、計算機を用いて「数値的に解く」問題との橋渡しをすることを目的とする
### 目次
1. 教育的な例の確認「スピン演算子$\hat{S_z}$とその固有値」(知っている人は飛ばしてOK)
2. 教育的な例を数値的に解く
3. 想定されるシチュエーション
## 1. 教育的な例の確認「スピン演算子Szとその固有値」
スピン演算子$\hat{S_z}$とその固有状態の関係は次のように与えられる。
```math
\hat{S_z} |\uparrow> = \frac{\hbar}{2} |\uparrow> ,~~~~~\hat{S_z} |\downarrow> = -\frac{\hbar}{2} |\downarrow>
```
これらを行列で表してみると、
```math
\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}
```
のように書き換えられる。この表示から明らかなように、
```math
\begin{pmatrix}
a \\
b
\end{pmatrix} = a |\uparrow> + b |\downarrow>
```
となっており、$|\uparrow>$が実際にどんな関数になっているかどうかはどうでもよく、その係数だけ計算されている。これは、「正規直交基底を用いれば、任意の新しい基底を表現できる」ことに基づいている。
さて、$\hat{S_z}$の固有状態は$\hat{S_y}$の固有状態だろうか?もちろん違う。では、$|\uparrow>$に$\hat{S_y}$を作用させるとどうなるだろうか?
これは昇降演算子から理解することができる。つまり、次の式
```math
S^{+}|\uparrow> = 0,~~~~~ S^{+}|\downarrow> = \hbar | \uparrow>\\
S^{-}|\uparrow> = \hbar | \downarrow>,~~~~~S^{-}|\downarrow> = 0
```
と、昇降演算子の定義
$$
S^{+} = \hat{S_x}+ i\hat{S_y},~~~~~ S^{-}= \hat{S_x} - i\hat{S_y}
$$
から、$\hat{S_x},\hat{S_y}$について逆に解けば、
```math
\hat{S_x}|\uparrow> = \frac{\hbar}{2}|\downarrow>, ~~~~~\hat{S_x}|\downarrow> = \frac{\hbar}{2}|\uparrow>\\
\hat{S_y}|\uparrow>=i\frac{\hbar}{2}|\downarrow>, ~~~~~\hat{S_y}|\downarrow> = -i\frac{\hbar}{2}|\uparrow>
```
が得られる。確かに$|\uparrow>$は$\hat{S_y}$の固有状態になっていなかった。
では、$\hat{S_y}$の固有状態は何だろうか?予測は付くだろうが、敢えてシステマティックに求めてみよう。まずは$\hat{S_y}|\uparrow> = i\frac{\hbar}{2}|\downarrow>$の行列表示は、
```math
\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$を対角化をすると、固有値と固有ベクトルは、
```math
\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)
+\lambda_2=-1, ~~ u_2=\frac{1}{\sqrt{2}}\begin{pmatrix}-1\\i\end{pmatrix}・・・(2)
```
となる。実際に固有状態になっているか確かめると、
```math
\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$の期待値を計算してみる。
```math
<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. 出力して、前節の解析解と一致するか確認
```python
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]
```
一致しない。これは固有ベクトルの格納のされ方に依るものである。以下の操作をするとよく分かる。
```python
print(eigenvectors @ sigma_y @ np.linalg.eigh(eigenvectors))
# output : [[-1 0.+1e-16j]
# [ 0.+1e-16 1 ]]
```
つまり、np.linalg.eigh()で得られる固有ベクトルは2つの固有ベクトルを縦に並べた形になっていて、
```
eigenvectors[0] =[u1(x1),u2(x1)]
eigenvectors[1] =[u1(x2),u2(x2)]
```
となる。というわけで、np.linalg.eigh()で得られた固有ベクトルは転置を取る必要がある。
```python
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. 想定されるシチュエーション
行列の対角化は物理の様々な状況で行われる。ここでは、一例として、「スピンありの1次元2サイトモデルで基底状態を求める」ということをやってみる。
ハミルトニアンとして、
```math
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}
```
を考え、状態としては「粒子数とスピンが保存し、アップスピン1個、ダウンスピン1個のハーフフィリング」を考える。この時、状態は4つ考えられる。
```math
\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$は次のような状態である。
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/485722/69701cb7-5b10-58f9-2ffd-85fb9a50aa8d.png" width=30%>
これらの状態を用いて表される状態を$\phi$とし、次のように表すことにする。もちろんこの$\phi$は全ての状態を表現することができる。
```math
\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$。このとき、
```math
H\phi=\begin{pmatrix}
-2\mu_1+U & -2t & 0 & 0 \\
-2t &-2\mu_2+U & 0 & 0 \\
-t & -t & -\mu_1-\mu_2 & 0 \\
-t & -t & 0 & -\mu_1-\mu_2
\end{pmatrix}\phi
```
のようにハミルトニアンを行列表示することができる。この行列を対角化することで、固有状態となる$\phi$を求めよう。やることは第2節と一緒である。
```python
import numpy as np
#計算条件
mu_1 = 0.1
mu_2 = 0.5
t = 1.0
U = 5.0
H = np.array([
[-2*mu_1 + U,  -2*t, 0, 0],
[ -2*t, -2*mu_2 + U, 0, 0],
[ -t,  -t, -mu_1 - mu_2, 0],
[ -t,  -t, 0, -mu_1 - mu_2]
])
eigenvalues, eigenvectors = np.linalg.eigh(H)
evs = np.transpose(eigenvectors)
print(eigenvalues[0], end=" ")
print(evs[0])
# output
# -1.604 [-0.301 -0.333 -0.632 -0.632 ]
# -0.6 [0.0 0.0 -0.707 0.707 ]
# 3.358 [0.544 0.709 -0.317 -0.317 ]
# 6.446 [-0.783 0.621 0.023 0.023 ]
```
計算結果を見てみると、固有値が一番小さな状態、つまり基底状態は、$[-0.301 -0.333 -0.632 -0.632]$である。クーロン力が大きい計算条件なので、$\psi_3, \psi_4$が多く含まれている。$\psi_1, \psi_2$の混ざり方は、サイトのエネルギーが違うので1:1ではない。