1. ovation6868lx

    No comment

    ovation6868lx
Changes in body
Source | HTML | Preview

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

はじめに

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

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

目的

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

目次

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

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

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

\hat{S_z} |\uparrow> = \frac{\hbar}{2} |\uparrow> ,~~~~~\hat{S_z} |\downarrow> = -\frac{\hbar}{2} |\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 |\uparrow> + b |\downarrow> 

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

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

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

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}$について逆に解けば、

\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>$の行列表示は、

\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}\\
\lambda_2=-1, ~~ u_2=\frac{1}{\sqrt{2}}\begin{pmatrix}1\\-i\end{pmatrix}
\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$の期待値を計算してみる。

<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()で得られる固有ベクトルは2つの固有ベクトルを縦に並べた形になっていて、

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

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

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

ハミルトニアンとして、

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つ考えられる。

\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 & -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節と一緒である。

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ではない。