$$
\def\bm#1{\boldsymbol{#1}}
\def\bra#1{\mathinner{\left\langle{#1}\right|}}
\def\ket#1{\mathinner{\left|{#1}\right\rangle}}
\def\brakett#1#2{\mathinner{\left\langle{#1}\middle|#2\right\rangle}}
\def\braket#1#2#3{\bra{#1}#2\ket{#3}}
\def\vector#1{\bm{#1}}
$$
初投稿です。
今回は、numpyを用いたバンド構造の計算方法について、自分の備忘録も兼ねてまとめます。バンド構造をはじめとする固体物理に関する説明の多くが欠落しておりますが、順次追加する予定です。
バンドを計算するためのモデルとしては、グラフェンやバレーフォトニック結晶の性質を説明するハニカム格子の Tight-binding model を用います。
物理的背景
バンド構造
周期的構造中に存在する、電子や光、音などの"波"が持つ分散関係のことを、バンド構造と呼びます。
例えば、電子の場合であればエネルギー $E$ と波数 $k$ の関係 $E=E(k)$ のことです。
自由電子の場合、波数$k$に対するエネルギー $E$ の式は $E(k)=\hbar k^2/2m$ で与えらえ、放物線になります。一方で、固体結晶のような周期的なポテンシャル $V$ の下での分散関係は一般にこれとは異なり、結晶構造の対称性や形状によって様々な形になります。周期構造の対称性によっては、$E=E(k)$を満たす$k$が存在しないような領域が現れます。これをバンドギャップと呼びます。
結晶構造の形を工夫してバンド構造を変調することで、電子や光に特徴的な性質を持たせることができます。
Tight-binding model
グラフェンなどの物性を説明するモデルです。
下のような六角形状に並んだサイトを考えます。数周期分しか示されていませんが、この構造が無限に続いていると考えてください。
以下、Tight-binding model を考え、Aサイト、Bサイトの波動関数$\ket{\phi_A}$,$\ket{\phi_B}$を基底とした2バンド模型を考えます。各サイトのポテンシャルと、隣接サイト同士の遷移（ホッピング）のみを考慮すると、波数空間における波動関数の振る舞いを記述する Bloch-Hamiltonian は次のようになります。
\begin{align}
H(\vector{k}) &= t_1
\begin{bmatrix}
0 & h(\vector{k}) \\
h^*(\vector{k}) & 0
\end{bmatrix}
+ M
\begin{bmatrix}
1 & 0 \\
0 & -1
\end{bmatrix}
\\
h(\vector{k}) &= \sum_i e^{i\vector{k}\cdot \vector{a_i}}
\end{align}
ただし、
● $t_1$ ...A-B間の結合の強さ
● $M$ ...A,Bのポテンシャル
であり、いずれも実数であるとします。[](gは非エルミート項であり通常の Tight-binding model には入っていませんが、光系などで発生する吸収や利得を表現する項として有用です。) また、$\vector{a_1}=a'[0,-1]^T$,$\vector{a_2}=a'[\sqrt{3}/2,1/2]^T$,$\vector{a_3}=a'[-\sqrt{3}/2,1/2]^T$とします。$a'$は最隣接サイト間の距離で、格子定数を$a$とすれば$a'=a/\sqrt{3}$です。
この固有方程式を解くことで、エネルギー $E$ を固有値として、系の波動関数 $\ket{\Psi}$ を固有ベクトルとして得ることができます。 $I_2$ は2*2の単位行列であるとします。
\begin{align}
0 &= \det (H(\vector{k}) - E I_2)
\nonumber\\
&=
\det
\begin{bmatrix}
M - E & t_1 h(\vector{k}) \\
t_1 h^*(\vector{k}) & -M - E
\end{bmatrix}
\nonumber\\
&= -{M}^2 + E^2 - t_1^2{|h(\vector{k})|}^2
\nonumber\\
E
&= \pm \sqrt{ t_1^2{|h(\vector{k})|}^2 + {M}^2 }
\end{align}
上の式から、一般に $M\neq 0$ のときには $E_+$ と $E_-$ は交わらず、バンドギャップが出現することがわかります。
一方で、$M=0$ 、つまりAサイトとBサイトが対等であるときには、特定の波数でバンドギャップが閉じます。$E_\pm$ の式に $h(\vector{k})$ の具体的な表式を代入することで、次の式を得ます。
\begin{align}
E_\pm
&= \pm {t_1|h(\vector{k})|}
\nonumber\\
&= t_1
\sqrt{
3
+ 2\cos \left(\sqrt{3} k_x a' \right)
+ 4\cos\left( \frac{\sqrt{3}}{2} k_x a' \right) \cos\left( \frac{3}{2} k_y a' \right) }
\end{align}
上式は Brillouin zone 上の K点およびK'点、すなわち $(k_x, k_y)=\left( \pm 4\pi/(3\sqrt{3}a'), 0\right)$ の近傍では、
$$
E_\pm(K+q) \simeq \pm \frac{3}{2}t_1 a'q
$$
と近似され、Dirac cone と呼ばれる線形分散が出現します。
Berry曲率
工事中
コード
前節で述べた行列計算を実装します。
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import seaborn as sns
from matplotlib import rc
# Pi
pi = np.pi
# Definition of Pauli matrix
Pauli_0 = np.array([[1+0j,0],[0,1]])
Pauli_x = np.array([[0j,1],[1,0]])
Pauli_y = np.array([[0,-1j],[1j,0]])
Pauli_z = np.array([[1+0j,0],[0,-1]])
class Haldane:
def __init__(self, M, t1,kx0=0, ky0=0, resolution=120, magnification = 1, save=0, **kwargs):
self.M = M
self.t1 = t1
# target point in k space
self.kx0 = kx0
self.ky0 = ky0
# resolution
self.resolution = resolution
# magnification
self.magnification = magnification
# hopping vector
self.a = [0,0,0]
self.a[0] = np.array([0, -1])/np.sqrt(3)
self.a[1] = np.array([np.sqrt(3)/2, 1/2])/np.sqrt(3)
self.a[2] = np.array([-np.sqrt(3)/2, 1/2])/np.sqrt(3)
# list of wavenumbers
self.kxx, self.kyy = np.meshgrid(np.linspace(self.kx0 - np.pi*2/self.magnification, self.kx0 + np.pi*2/self.magnification, self.resolution), np.linspace(self.ky0 - np.pi*2/self.magnification, self.ky0 + np.pi*2/self.magnification, self.resolution))
self.dkx = self.kxx[0][1] - self.kxx[0][0]
self.dky = self.kyy[0][1] - self.kyy[0][0]
# save
self.save = save
# To calculate Hermitian conjugate
def hermite(self, vec, **kwargs):
return np.conjugate(vec.T)
def build_H(self, kx, ky, **kwargs):
k = np.array([kx,ky])
h0 = sum([ np.exp( 1.j * np.dot(k, self.a[i]) ) for i in range(len(self.a))])
H_0 = self.t1 * np.array([[0,h0],[np.conj(h0),0]])
H_1 = self.M * Pauli_z
H = H_0 + H_1
return H
def gen_data(self, **kwargs):
self.list_kx = np.array([])
self.list_ky = np.array([])
self.E = np.array([]) #eigenvalue
self.psi = np.array([]) #eigenvector
self.num_bands = self.build_H(0,0).shape[0] # Number of bands
#self.dky = kyy[1][0] - kyy[0][0]
for j in range(self.resolution):
for i in range(self.resolution):
kx = self.kxx[i][j]
ky = self.kyy[i][j]
E_n_ij, psi_n_ij = np.linalg.eig( self.build_H(kx, ky) )
self.list_kx = np.append(self.list_kx, kx)
self.list_ky = np.append(self.list_ky, ky)
self.E = np.append(self.E, E_n_ij)
self.psi = np.append(self.psi, psi_n_ij)
self.E = self.E.reshape(self.resolution, self.resolution, self.num_bands)
self.psi = self.psi.reshape(self.resolution, self.resolution, self.num_bands, self.num_bands)
# sort data
for i in range(self.resolution):
for j in range(self.resolution):
sortE = self.E[i][j].real.argsort()
self.E[i][j] = self.E[i][j][sortE]
self.psi[i][j] = self.psi[i][j][:,sortE]
def gen_berry_curvature(self, **kwargs):
# calculation of the berry curvature
self.berry = np.array([])
for bnd in range(self.num_bands):
for kx in range(1,self.resolution-1):
for ky in range(1,self.resolution-1):
A = list([0,0,0,0])
#print(psi[:,:,:,bnd][ky-1][kx-1])
A[0] = np.dot(self.hermite(self.psi[:,:,:,bnd][ky-1][kx-1]), self.psi[:,:,:,bnd][ky-1][kx+1])
A[1] = np.dot(self.hermite(self.psi[:,:,:,bnd][ky-1][kx+1]), self.psi[:,:,:,bnd][ky+1][kx+1])
A[2] = np.dot(self.hermite(self.psi[:,:,:,bnd][ky+1][kx+1]), self.psi[:,:,:,bnd][ky+1][kx-1])
A[3] = np.dot(self.hermite(self.psi[:,:,:,bnd][ky+1][kx-1]), self.psi[:,:,:,bnd][ky-1][kx-1])
berry_ij = np.log(A[0]*A[1]*A[2]*A[3]) / ((2*self.dkx )**2)
self.berry = np.append(self.berry, berry_ij)
self.berry = self.berry.reshape(self.num_bands, self.resolution-2, self.resolution-2)
self.berry = self.berry.imag
def plot_band(self, **kwargs):
print("M :", self.M)
print("t_1 :", self.t1)
fig = plt.figure(figsize=(8,8))
ax0 = fig.add_subplot(111, projection='3d')
E_lim = 3
props={'linewidth':0,'cmap':'viridis','antialiased':False,
'rstride':1, 'cstride':1,'vmin':-E_lim, 'vmax':E_lim, "alpha":0.8}
for i in range(self.num_bands):
#ax0.plot_wireframe(kxx,kyy,z[:,:,n])
ax0.plot_surface(self.kxx, self.kyy, self.E[:,:,i].real,**props)
ax0.set_xlim3d(-pi*2,pi*2)
ax0.set_ylim3d(-pi*2,pi*2)
ax0.set_zlim3d(-E_lim,E_lim)
fs = 20
lprop = {"fontsize":fs}
ax0.set_xlabel(r"$k_x$", **lprop)
ax0.set_ylabel(r"$k_y$", **lprop)
#ax0.set_zlabel(r"E", **lprop)
ax0.set_xticks(np.linspace(-pi,pi,3))
ax0.set_yticks(np.linspace(-pi,pi,3))
#ax0.set_zticks(np.linspace(-3,3,1))
#以下の3つは記号の関係上qiitaの記事がバグるのでコメントアウトしていますが、使用時はONにしてください。
#ax0.set_xticklabels(["-$\pi$", "$0$", "$\pi$"], va="baseline", ha = "center", **lprop)
#ax0.set_yticklabels(["-$\pi$", "$0$", "$\pi$"], va="baseline", ha = "center", **lprop)
#ax0.set_zticklabels(["$-3$","$-2$","$-1$","$0$","$1$","$2$","$3$"],**lprop)
filename = "band3D_M=" + str(self.M) + "_t1=" + str(self.t1)
plt.title(filename, fontsize=fs*0.75)
if self.save == 1:
plt.savefig(filename + ".png")
plt.show()
def plot_berry_curvature(self, band_number,**kwargs):
fs = 20
lprop = {"fontsize":fs}
z = self.berry[band_number]
zlim = max(abs(z.max()), abs(z.min()))
fig = plt.figure()
ax = fig.add_subplot(111)
im = ax.imshow(z,extent=(-2,2,-2,2),interpolation='nearest',cmap='seismic',aspect=1)
ax.set_xticks(np.linspace(-2,2,3))
ax.set_yticks(np.linspace(-2,2,3))
#以下の3つは記号の関係上qiitaの記事がバグるのでコメントアウトしていますが、使用時はONにしてください。
#ax.set_xticklabels(["-$2\pi$", "$0$", "$2\pi$"], **lprop)
#ax.set_yticklabels(["-$2\pi$", "$0$", "$2\pi$"], **lprop)
im2 = fig.colorbar(im)
im2.mappable.set_clim(-zlim,zlim)
ax.set_xlabel("$k_x$", **lprop)
ax.set_ylabel("$k_y$", **lprop)
plt.tick_params(labelsize=fs)
filename = "Berry_Curvature_M=" + str(self.M) + "_t1=" + str(self.t1)
plt.title(filename, fontsize=fs*0.75)
fig.tight_layout()
if self.save == 1:
plt.savefig(filename + ".png")
plt.show()
以下のコードで実行できます。
M = 0.
t1 = 1.
#the center point
kx0 = 0
ky0 = 0
resolution = 200
magnification = 1
save = 0
model = Haldane(M, t1, kx0, ky0, resolution, magnification, save)
model.gen_data()
model.gen_berry_curvature()
model.plot_band()
model.plot_berry_curvature(band_number=0)
各関数の解説
- hermite
- Hermite共役を計算するために設定した関数です。
- build_H
- 入力した変数 $M,t_1$ から Hamiltonian $H(\vector{k})$ を行列として生成します。
- gen_data
- 各波数 $\vector{k}=(k_x,k_y)$ に対して $H(\vector{k})$ の固有値と固有関数を計算し、リスト化します。
- gen_berry_curvature
- 固有関数のリストからBerry曲率を計算します。
- plot_band
- 得られた固有値を、$\vector{k}$ の関数として3Dグラフにプロットします。
- plot_berry_curvature
- 得られたBerry曲率を、$\vector{k}$ の関数としてヒートマップにプロットします。
- 引数の band_number は、上下どちらのバンドのBerry曲率を表示するかを指定します。（下のバンドなら0）
計算結果
まずはじめに、$M=0, t_1=1$のときの結果を示します。
バンド図は次のようになります。
2つのバンドの交点付近に線形分散が現れていることがわかります。
$M=0$において、Berry曲率はK点およびK'点にdelta関数上に現れるため、下に示すように数値解析では正しく計算することができません。一応delta関数っぽくはなっていますが...
続いて、ポテンシャルの差 $M$ を与えて同様の計算を行いました。
$M=0.3, t_1=1$としたときのバンド図とBerry曲率を示します。
線形分散が消え、上下2つのバンドが開いていることがわかります。すなわち、エネルギーの飛びが現れ、バンドギャップが出現しています。
また、Berry曲率もK点およびK'に広がりを持ち、有限の値を取っていることがわかります。
備考
・バンド理論の適用先の多くは電子系ですが、周期的な誘電体や弾性体構造を用いることで、光や音波の分散関係も変化させることができます。
・今回はポテンシャルやホッピング係数が実数である場合を考えました。これらの係数が複素数になる場合、バンドの縮退や Exceptional Point と呼ばれる特異点の出現といった、奇妙な現象が出現することが知られています。複素数係数で表現されるような系は非エルミート系と呼ばれ、近年特に光系で盛んに研究が進められています。
結論
numpy によって Tight-binding model を実装し、バンド図とBerry曲率を正しく計算することができました。