Qiita Teams that are logged in
You are not logged in to any team

Log in to Qiita Team
Community
OrganizationEventAdvent CalendarQiitadon (β)
Service
Qiita JobsQiita ZineQiita Blog
3
Help us understand the problem. What are the problem?

More than 1 year has passed since last update.

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

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

はじめに

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

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

目的

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

目次

  1. 教育的な例の確認「スピン演算子$\hat{S_z}$とその固有値」(知っている人は飛ばしてOK)
  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()で得られる固有ベクトルは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$の(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 & 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$を求めよう。やることは第2節と一緒である。

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.29~~0.34~~0.63~~0.63]$である。クーロン力が大きい計算条件なので、$\psi_3, \psi_4$が多く含まれている。$\psi_1, \psi_2$の混ざり方は、サイトのエネルギーが違うので1:1ではない。

Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
3
Help us understand the problem. What are the problem?