Phonon¶
1. Introduction¶
1.1 フォノン (phonon)とは¶
フォノンは、一言で表すと周期境界条件をもつ結晶における振動モードです。
詳細の説明に入る前に、Phononを目で見て感覚を掴むためこちらのサイトを見てましょう。
フォノンのアニメーション表示。http://henriquemiranda.github.io/phononwebsite/phonon.html より。
上記Websiteへ行く -> 左上のmaterials という表から Si を選択 -> 右図の“フォノン分散図”から適当に点をクリックしてみましょう。 すると、真ん中のシリコン結晶におけるシリコン原子それぞれが周期的な動きをしていることが確認できます。そして、全体として片側から反対側に波が伝わるような動きとなっていることが見て取れます。
このような結晶における振動を格子振動と呼びます。格子振動は振動数が非常に大きくなるため、通常は量子力学的な効果を考慮して取り扱われます。格子振動に量子効果を含めたものにはフォノンという名前がついています。ただし、以下の説明では複雑さを避けるため、格子振動とフォノンは特別に用語を区別せずに解説しています。より詳しくは末尾のAppendixや固体物理学の教科書を参考にしてください。
波の伝播する方向は波数ベクトル \(\vec q\) で表されます。 フォノンの解析を行う際に重要な問題は、ある波数ベクトル \(\vec q\) に対し、どのような振動モードが許されているのか、その振動モードにおける原子の振動数がどうなるかを求めることです。
フォノンは準粒子であり、エネルギー \(E\)と運動量 \(\mathbf{p}\)をもちます。 固体中のフォノンは熱物性や音波特性に深い関わりがあります。例えば、フォノンの励起は内部エネルギー \(U\)の増加に寄与したり、フォノンの輸送は音や熱の伝導を引き起こします。このように、フォノン解析は固体の性質を探る上で重要なシミュレーションツールとなります。
1.2 フォノン振動数の計算方法¶
前章で紹介した振動解析と同様、フォノンも調和振動子近似を用いて計算されます。 ただし違いとして、各原子はUnit cellをまたいで向こうにいる原子とも相互作用をします。 周期境界条件を考慮して、結晶構造から微小変位をさせた系のポテンシャルエネルギーは以下のように表されます。
ここで、\(r_i^a\) \(a\)番目のスーパーセルに属する、\(i\)番目の原子の座標です。結晶における力定数マトリクス (前章で紹介) は\(3NM \times 3NM\) 行列となります。 \(N\) はユニットセルの原子数、 \(M\) はユニットセルの数です。 このため、力定数マトリクスはとても大きな行列となります。
特定の振動モードにおける、各原子の平衡点からの変位 \(u\)は次の平面波のような形を持つと仮定します。
詳細は省きますが、この仮定のもとで実空間における離散フーリエ変換を適用することで、 大きな力定数マトリクスを\(3N \times 3N\)のdynamical matrixに変換することができます。 Dynamical matrixは波数ベクトル\(\vec q\)に依存しており、その固有値と固有ベクトルがフォノンの振動数と振動モードに相当します。 詳細な導出については、末尾のappendixを参照ください。
次章では、フォノンの計算をMatlantis上で提供しているライブラリであるmatlantis-featuresを用いて行いながら、フォノンの重要なコンセプトを説明します。
2. Siでのフォノン計算¶
2.1 構造準備¶
ここでは、フォノン計算の例としてシリコン結晶を用います。 まずは構造緩和を行います。 フォノン分散を計算する時は、その結晶構造の基本単位胞 (primitive cell)を用いることが多いです。セル内を最小原子で構成できるためです。
[1]:
import numpy as np
np.set_printoptions(precision=3)
from ase.build import bulk
from ase.optimize import BFGS
from ase.constraints import ExpCellFilter
import pfp_api_client
from pfp_api_client.pfp.calculators.ase_calculator import ASECalculator
from pfp_api_client.pfp.estimator import Estimator, EstimatorCalcMode
estimator = Estimator(calc_mode=EstimatorCalcMode.CRYSTAL, model_version="v3.0.0")
calculator = ASECalculator(estimator)
atoms = bulk("Si")
atoms.calc = calculator
exp = ExpCellFilter(atoms)
BFGS(exp).run(fmax=0.01)
Step Time Energy fmax
BFGS: 0 14:31:50 -9.108033 0.496677
BFGS: 1 14:31:50 -9.113258 0.001942
/tmp/ipykernel_17542/260060846.py:17: FutureWarning: Import ExpCellFilter from ase.filters
exp = ExpCellFilter(atoms)
[1]:
True
2.2 力定数マトリクスの計算¶
振動解析同様、フォノンの振動数計算には力定数マトリクスが必要です。 この力定数マトリクスは、ForceConstantFeatureで計算されます。
Siの基本単位胞は2原子が含まれています。 この基本単位胞をたくさん繰り返したスーパーセル構造を作成することが精度の良い結果を得るために必要です。 力定数マトリクスの各成分 \(\frac{\partial^2 V(\mathbf{r_0})}{\partial r_i^a \partial r_j^b}\) は、\(a\)番目のユニットセルの\(i\)番目の原子と \(b\)番目のユニットセルの\(j\)番目の原子との距離で減衰していきますが、スーパーセルが十分大きい場合にのみ力定数マトリクスで無視できないような要素は含まれているとみなせます。
本チュートリアルでは、 Si結晶を 10x10x10 に拡張して計算を行います。
matlantis-featuresのForceConstantFeatureを用いて計算します。
[2]:
from matlantis_features.features.phonon import ForceConstantFeature
fc = ForceConstantFeature(
supercell = (10,10,10),
delta = 0.1,
)
force_constant = fc(atoms)
print(force_constant.force_constant.shape)
/home/jovyan/.py39/lib/python3.9/site-packages/matlantis_features/features/phonon/force_constant.py:95: UserWarning: Model version is not specified, default 'v7.0.0' will be used.Please specify it by providing an estimator_fn or setting the 'MATLANTIS_PFP_MODEL_VERSION' environment variable.
self.check_estimator_fn(estimator_fn)
/home/jovyan/.py39/lib/python3.9/site-packages/matlantis_features/features/phonon/force_constant.py:95: UserWarning: Calculation mode is not specified, default 'EstimatorCalcMode.CRYSTAL_U0' will be used.Please specify it by providing an estimator_fn or setting the 'MATLANTIS_PFP_CALC_MODE' environment variable.
self.check_estimator_fn(estimator_fn)
(1000, 6, 6)
計算結果として得られる力定数マトリクス \(C\) は 1000x6x6 のshapeを持つtensorとなります。
[3]:
print(force_constant.force_constant[1])
[[-0.218 -0.193 -0.092 -0.01 0.015 0.056]
[-0.193 -0.218 -0.092 0.015 -0.01 0.056]
[ 0.092 0.092 0.433 0.056 0.056 -0.096]
[-3.081 -2.114 2.114 -0.218 -0.193 0.092]
[-2.114 -3.081 2.114 -0.193 -0.218 0.092]
[ 2.114 2.114 -3.081 -0.092 -0.092 0.433]]
force_constant.force_constant[1] は 6x6 のtensorです。 これは、1つめと2つめのユニットセルにおける力定数 \(C^{1,2}\)、すなわち $:nbsphinx-math:frac{partial^2 V(mathbf{r_0})}{partial r_i^1 partial r_j^2} $です。 前章の振動解析同様、力定数 \(C^{1,2}\)は有限差分法で近似計算されています。
今回は、10x10x10 のスーパーセルを用いているため、ユニットセルの組み合わせである1000x1000 の数だけ 6x6 matrix が得られるべきです。 しかし、結晶の並進対称性を考慮し、\(C^{a,b} = C^{0, b-a}\) となることを利用することで、力定数マトリクスを 1000x6x6 tensorで表すことができます。
2.3 Dynamical matrixの計算¶
分子の振動数は力定数マトリクス\(C\)の対角化で得られましたが、結晶の振動数は実空間の力定数マトリクスを離散フーリエ変換したdynamical matrix \(D_{\vec q}\) の固有値として得られます。 Dynamical matrixへの変換の際に、波数ベクトル \(\vec q\) というパラメータが追加で導入されます。
ここでは、matlantis-featuresのPhononFreqencyを用いて、シリコン結晶の波数ベクトル \(\vec q = [0.1, 0.0, 0.0]\) に対応するdynamical matrix を計算してみます。
[4]:
from matlantis_features.features.phonon.utils import PhononFrequency
phonon_frequency = PhononFrequency(
force_constant.force_constant, force_constant.supercell, force_constant.unit_cell_atoms
)
dynamic_matrix = phonon_frequency._dynamic_matrix_q(np.array([0.1, 0.0, 0.0])).real
print("Dynamical matrix of q=[0.1, 0.0, 0.0]:")
print(dynamic_matrix)
Dynamical matrix of q=[0.1, 0.0, 0.0]:
[[ 0.443 -0.004 -0.004 -0.419 -0.011 -0.011]
[-0.004 0.443 0.004 -0.011 -0.418 0.011]
[-0.004 0.004 0.443 -0.011 0.011 -0.418]
[-0.419 -0.011 -0.011 0.443 -0.004 -0.004]
[-0.011 -0.418 0.011 -0.004 0.443 0.004]
[-0.011 0.011 -0.418 -0.004 0.004 0.443]]
Dynamical matrixは 6x6 の対称行列となっています。これは、力定数マトリクス \(C\) から以下の式の計算により得られています。
その後、この行列を対角化することで固有値・固有ベクトルを得ます。 分子の振動解析と似たように、固有値が振動数に対応し、固有ベクトルはそのフォノンモードにおける原子の変位と関連付けられます。 式の導出など、詳細はappendixを参照ください。
それでは、シリコン結晶の \(\vec q = [0.1, 0.0, 0.0]\) におけるフォノンの振動数を計算してみましょう。
[5]:
frequency, eigenvector = phonon_frequency._frequency_mode_q(np.array([0.1, 0.0, 0.0]), return_mode=True)
for i in range(6):
print(f"Band {i}")
print(f" Frequency: {frequency[i]:.3f} meV")
print(f" Eigenvector:")
print(eigenvector[i].real)
Band 0
Frequency: 6.116 meV
Eigenvector:
[[-0.109 -0.049 -0.06 ]
[-0.108 -0.049 -0.06 ]]
Band 1
Frequency: 6.116 meV
Eigenvector:
[[ 0.006 0.097 -0.091]
[ 0.006 0.097 -0.091]]
Band 2
Frequency: 11.568 meV
Eigenvector:
[[ 0.073 -0.073 -0.073]
[ 0.077 -0.077 -0.077]]
Band 3
Frequency: 60.273 meV
Eigenvector:
[[-0.072 0.073 0.073]
[ 0.076 -0.078 -0.077]]
Band 4
Frequency: 60.308 meV
Eigenvector:
[[-0.017 -0.093 0.078]
[ 0.019 0.096 -0.078]]
Band 5
Frequency: 60.309 meV
Eigenvector:
[[-0.094 -0.017 -0.076]
[ 0.097 0.019 0.076]]
波数ベクトル :math:`vec q`
フォノンの計算では、波数ベクトル \(\vec q\) が導入されました。 この波数ベクトルの物理的な意味は、その振動モードの伝播する方向と関連します。 また、理論的に波数ベクトル \(\vec q\) は第一ブリルアンゾーン (first Brillouin zone)の中に制限できることが示せます (詳細はappendixを参照)。
ここでは、get_brillouin_zone_3d関数を使い、シリコン結晶における第一ブリルアンゾーンを可視化してみます。 Wikipedia上の以下の画像と比較してみてください
対称性ラベルを付した面心立方格子構造(FCC構造)の第一ブリルアンゾーン。 https://en.wikipedia.org/wiki/Brillouin_zone より。
[6]:
from matlantis_features.utils.visual_utils.brillouin_zone import get_brillouin_zone_3d
fig = get_brillouin_zone_3d(atoms.cell, show_reciprocal_basis=True, show_special_k_points=True)
fig.update_layout(width=600, height=500)
2.4 フォノン分散の計算¶
PostPhononBandFeatureは、ForceConstantFeatureの計算結果を利用します。 対象性の高い点のパスを指定する際はラベルの列を指定します。 ここでは、G -> X -> W -> K -> G -> L -> U -> W -> L -> K と U -> X というパスでのフォノン振動数を計算してフォノン分散を得ています。
シリコン結晶はFCC構造となっており、その第一ブリルアンゾーンは切頂八面体 (truncated octahedron)となります。 緑の点で示されているのが対称性の高い点で、例えば [0.375, 0.375, 0.75]の点は“K” というラベルが付与されています。上記図をマウスオーバーするとラベルが表示されるので確認してみてください。 第一ブリルアンゾーンの中心は特に\(\Gamma\) (or G)点 (ガンマ点)と名付けられています。
許されるすべてのフォノンの振動数を求めるには、第一ブリルアンゾーン内のすべての\(\vec q\)点にて計算を行う必要があります。 しかし実際には、第一ブリルアンゾーン内の代表的な対象性の高い点をつなぐパス(例えば\(\Gamma\)点とK点を結ぶ直線など)上でのみ計算を行って可視化をして確認するということが慣例的に行われています。 対象性の高い点をつないだパス上のフォノンの振動数のプロットは、フォノン分散曲線 (フォノン分散関係・フォノンバンド)と呼ばれます。
ここでは、matlantis-featuresのPostPhononBandFeatureを用いて、シリコン結晶のフォノン分散を計算します。
[7]:
from matlantis_features.features.phonon import PostPhononBandFeature
band = PostPhononBandFeature()
band_results = band(
force_constant,
labels = ['G', 'X', 'W', 'K', 'G', 'L', 'U', 'W', 'L', 'K', '|', 'U', 'X'],
total_n_kpts = 100,
)
fig = band_results.plot()
fig.update_layout(width=700, height=500)
上記の図をみると、シリコンのフォノン分散では６つの線が存在している事がわかります(場所によっては2本の線が縮退して1本に見える箇所もあります)。 つまり、各波数ベクトルに対して６つの振動モードがあるということを表しています。
これらの振動モードは、光学モード(optical)と音響モード(acoustic)に分類できます。 さらに、その振動モードの伝播方向により、縦方向(longitudinal)と横方向(transverse)にも分類されます。
次は実際にこの分類を行っていきます。
2.5 フォノンモード¶
光学モード・音響モード
上記で説明したように、各波数ベクトル \(\vec q\) に対して \(3N\) 個の振動数を得ることができます。 分散曲線には\(3N\)個の線が存在します。 これは、\(3\)つの音響モードと、\(3N-3\)つの光学モードに分けられます。
音響モードは以下のような特徴を持ちます。
音響モードは低い振動数を持つ。
音響モードの振動数 \(\omega\) は、ガンマ点で0となる。\(\frac{\partial \omega} {\partial q}\)は音速を与える。
音響モードの\(q=0\)における格子振動を確認すると、格子全体が動いていることが見て取れます。これは弾性波と類似しています。
光学モードは以下のような特徴を持ちます。
光学モードは高い振動数を持つ。
近接する原子同士は逆方向に動きあう。
複数の元素を含む結晶、例えばNaClなどでは、カチオンとアニオン (イオン参照)の間でおこる振動によりダイポールモーメントが変化します。このダイポールモーメントが外場の電磁波と相互作用します。そのため、光学モードと呼ばれます。
[8]:
from ipywidgets import HBox
from ase.visualize import view
from matlantis_features.features.phonon import PostPhononModeFeature
from pfcc_extras.visualize.povray import traj_to_apng
mode = PostPhononModeFeature()
mode_results = mode(
force_constant,
k_point=(0.1, 0.0, 0.0),
repeat_of_cell=(4, 4, 4),
n_images=30,
amplitude=5.0
)
for i in range(6):
traj_to_apng(mode_results.trajectories[i], f"output/vib.{i}.png", rotation="0x,0y,0z", clean=True, n_jobs=16)
/home/jovyan/.py39/lib/python3.9/site-packages/networkx/utils/backends.py:135: RuntimeWarning:
networkx backend defined more than once: nx-loopback
[Parallel(n_jobs=16)]: Using backend ThreadingBackend with 16 concurrent workers.
[Parallel(n_jobs=16)]: Done 30 out of 30 | elapsed: 27.2s finished
[Parallel(n_jobs=16)]: Using backend ThreadingBackend with 16 concurrent workers.
[Parallel(n_jobs=16)]: Done 30 out of 30 | elapsed: 27.8s finished
[Parallel(n_jobs=16)]: Using backend ThreadingBackend with 16 concurrent workers.
[Parallel(n_jobs=16)]: Done 30 out of 30 | elapsed: 27.7s finished
[Parallel(n_jobs=16)]: Using backend ThreadingBackend with 16 concurrent workers.
[Parallel(n_jobs=16)]: Done 30 out of 30 | elapsed: 26.6s finished
[Parallel(n_jobs=16)]: Using backend ThreadingBackend with 16 concurrent workers.
[Parallel(n_jobs=16)]: Done 30 out of 30 | elapsed: 26.2s finished
[Parallel(n_jobs=16)]: Using backend ThreadingBackend with 16 concurrent workers.
[Parallel(n_jobs=16)]: Done 30 out of 30 | elapsed: 26.5s finished
Acoustic vibration modes
Mode #0
Mode #1
Mode #2
Optical vibration modes
Mode #3
Mode #4
Mode #5
可視化してみると、Mode 0, 1, 2の原子は全体が同じ方向に動いている一方、Mode 3, 4, 5では原子同士が相互に振動していることが見て取れます。 低い周波数の3つのモードが音響モード、残りが光学モードに対応していることがわかります。
横(Transverse)モードと縦(Longitudinal)モード
フォノン分岐は原子の振動の向きと波数ベクトルの向きとの関係により横モードと縦モードに分けられます。 横モードでは原子の変位は波数ベクトル \(\vec q\) に対して垂直なります。 縦モードでは、原子の変位が波数ベクトル \(\vec q\) の方向に対して並行になります。
以下のように計算で横モード・縦モードを確認することもできます。
[9]:
reciprocal_basis = atoms.cell.reciprocal()
q_real = reciprocal_basis.dot(np.array([0.1, 0, 0]))
def get_angles(eigenvector, q_real):
return np.arccos(eigenvector.dot(q_real) / np.linalg.norm(eigenvector, axis=1) / np.linalg.norm(q_real))
for i in range(6):
angles = get_angles(eigenvector[i].real, q_real)
print(f"Band {i}:")
print(f" Angles between q and eigenvectors:", angles)
Band 0:
Angles between q and eigenvectors: [1.571 1.571]
Band 1:
Angles between q and eigenvectors: [1.571 1.571]
Band 2:
Angles between q and eigenvectors: [3.142 3.142]
Band 3:
Angles between q and eigenvectors: [0.008 3.134]
Band 4:
Angles between q and eigenvectors: [1.565 1.576]
Band 5:
Angles between q and eigenvectors: [1.566 1.578]
波数ベクトル\(\vec{q}\)と固有ベクトルが並行となっているBand 2, 3は縦モード 波数ベクトル\(\vec{q}\)と固有ベクトルが垂直となっているBand 0, 1, 4, 5は横モードであることが確認できます。
まとめると、シリコン結晶におけるフォノンの6分岐は以下のように分類できます。
Transverse acoustic branch (TA): branch 0 and 1
Longitudinal acoustic branch (LA): branch 2
Transverse optic branch (TO): branch 4 and 5
Longitudinal optic branch (LO): branch 3
2.6 フォノンの状態密度 (Phonon DOS)¶
フォノンの分散関係は物性値に関する情報を与えてくれますが、3次元の波数空間上の非自明な関数です。
代わりに、それぞれの振動モードにおけるエネルギーごとに集約をして、その物質中にどのようなモードが存在するかを簡易に表す方法があります。 これをフォノンの状態密度(Phonon DOS: density of states)と呼び、\(\omega\) から \(\omega+d\omega\)の間に存在する振動モードの数として以下のように定義されます。
ここでは、PostPhononDOSFeature を用いて、シリコン結晶のフォノンの状態密度を計算します。
[10]:
from matlantis_features.features.phonon import PostPhononDOSFeature, plot_band_dos
dos = PostPhononDOSFeature()
dos_results = dos(
force_constant,
kpts = [20,20,20],
freq_min = -5.0,
freq_max = 70.0,
freq_bin = 1.0,
scheme = "mp",
unit = "meV"
)
[11]:
from plotly.offline import iplot
fig = plot_band_dos(band_results, dos_results)
fig.update_layout(width=1000, height=500)
iplot(fig)
上図では、フォノン分散とフォノンの状態密度の図が同時に表示されています。
PostPhononDOSFeatureでは、フォノンの振動数を第一ブリルアンゾーン内の等間隔メッシュ \(\vec q\) 上で計算しています。 メッシュのとり方は、Monkhorst-Pack法とGamma-centered法がよく使われます。 上図では、\(20 \times 20 \times 20\) Monkhorst-Pack style gridが使用されています。 フォノンの状態密度に関する詳細はappendixをご参照ください。
2.7 熱物性値¶
固体の熱物性値はフォノン分散と密接な関係にあります。 フォノンの状態密度がわかれば、以下の式により、比熱 \(C_V\)・エントロピー \(S\)・内部エネルギー \(U\)・ヘルムホルツ自由エネルギー \(A\)は計算が可能です。
[Note]
量子効果を含めたフォノンは離散的なエネルギーしかとらないため、温度の依存性が古典的な格子振動と変わることが知られています。そのため上記の式に量子効果と関わる物理定数であるディラック定数\(\hbar\)が登場しています。
Matlantis上では PostPhononThermochemistryFeature を用いることで上記の計算が可能です。 PostPhononDOSFeature同様、波数ベクトル\(\vec q\)のメッシュの切り方を定義する必要があります。
PostPhononThermochemistryFeatureは力定数マトリクスからフォノンの状態密度を計算し、上記式を用いて熱物性値を計算します。
[12]:
from matlantis_features.features.phonon import PostPhononThermochemistryFeature
thermo = PostPhononThermochemistryFeature()
thermo_results = thermo(
force_constant,
kpts = [20, 20, 20],
tmin = 0.0,
tmax = 1000.0,
tstep = 10.0,
scheme = "mp",
)
[13]:
import pandas as pd
import plotly.graph_objects as go
thermo_dft = pd.read_csv("../assets/Si_thermo_dft.csv")
thermo_expt = pd.read_csv("../assets/Si_thermo_expt.csv")
fig = thermo_results.plot()
fig.add_trace(
go.Scatter(
x = thermo_dft["T"].loc[:100],
y = thermo_dft["entropy"].loc[:100] / 28.09 / 2.,
name = "S DFT(J/K/g)"
)
)
fig.add_trace(
go.Scatter(
x = thermo_dft["T"].loc[:100],
y = thermo_dft["cv"].loc[:100] / 28.09 / 2.,
name = "Cv DFT(J/K/g)"
)
)
fig.add_trace(
go.Scatter(
x = thermo_dft["T"].loc[:100],
y = thermo_dft["internal_energy"].loc[:100] / 28.09 / 2./1000.,
name = "U DFT (kJ/g)"
)
)
fig.add_trace(
go.Scatter(
x = thermo_dft["T"].loc[:100],
y = thermo_dft["helmholtz_free_energy"].loc[:100] / 28.09 / 2./1000.,
name = "F DFT(kJ/g)"
)
)
fig.add_trace(
go.Scatter(
x = thermo_expt["T"],
y = thermo_expt["cp"],
mode = "markers",
name = "Cp expt(J/k/g)"
)
)
fig.add_trace(
go.Scatter(
x = thermo_expt["T"],
y = thermo_expt["entropy"],
mode = "markers",
name = "S expt(J/k/g)"
)
)
fig.update_layout(width=800, height=500)
上図は計算結果のプロットです。計算結果が、DFTの計算結果や、実験値と比べてもよく一致していることが確認できます。
DFT calculation results from
Experiment results from
“Condensed Phase Heat Capacity Data” by Eugene S. Domalski and Elizabeth D. Hearing in NIST Chemistry WebBook, NIST Standard Reference Database Number 69, Eds. P.J. Linstrom and W.G. Mallard, National Institute of Standards and Technology, Gaithersburg MD, 20899, https://doi.org/10.18434/T4D303, (retrieved May, 2022).
3. 調和振動子近似の限界¶
これまでに見てきたように、結晶のフォノン分散・フォノンの状態密度・熱物性値が調和振動子近似から計算できることがわかりました。 しかし、ポテンシャルエネルギー局面の２次以降の項を無視してしまっていることによりうまく扱うことができない物性値も存在します。 代表的なものとして、熱膨張率と熱伝導率が挙げられます。
3.1 熱膨張率¶
熱膨張率はフォノン解析から直接得ることができません。調和振動子近似を用いるということは、近接原子がバネでつながっているモデルを仮定しており、エネルギー局面を局所的に２次関数として扱っています。 現実のポテンシャルエネルギー局面はもちろん２次関数ではありません(前章の図1参照)。 近接原子がバネでつながっている調和振動子近似では、温度を上げても結合長が増えません。結果として、固体の体積も温度により増加するという現象を扱えず、熱膨張率は0となってしまいます。
しかし、熱膨張率はquasi-harmonic approximationと呼ばれる、調和振動子近似の拡張系で扱うことが可能です。Quasi-harmonic approximationでは、密度の高い・低い系それぞれでフォノン計算を行っていきます。 体積依存するような物性値である、比熱 (Cp)・ギブス自由エネルギー (G)なども同様に計算することができます。
3.2 熱伝導率¶
フォノンはエネルギーや運動量を持ちます。熱の輸送は結晶内でのフォノンの移動として説明することができます。しかし、フォノン解析は調和振動子近似をもとに行われていることから、フォノンの散乱が無視されてしまいます。つまり、フォノンは別のフォノンと衝突せずに無限に離れた場所まで移動することになり、熱伝導率も無限に大きくなってしまいます。 そのため、単純なフォノン解析では熱伝導率を得ることができません。
熱伝導率を扱うためには、非調和性とフォノン散乱を扱うために、テイラー展開の3次の後まで考慮する必要があります。
つまり、3次の力定数である \(\frac{\partial^3{E}}{\partial{u_{k}^{a\alpha}} \partial{u_{k'}^{b\beta}} \partial{u_{k''}^{c\gamma}}}\) を計算する必要があります。 格子の熱伝導率を計算する方法としては、real-time approximationなどいくつかの手法が開発されています。
4. 参考文献¶
[2] http://www-personal.umich.edu/~sunkai/teaching/Winter_2013/chapter5.pdf
[3] http://staff.ustc.edu.cn/~zqj/posts/LinearTetrahedronMethod/
[4] Otero-De-La-Roza, et. al. Comput. Phys. Commun. 182, 2232–2248 (2011).
[5] 「固体物理学入門 上・下」 キッテル
[6] 「初歩から学ぶ 固体物理学」 矢口裕之
5. Appendix¶
5.1 フォノン計算の導出¶
Illustration
式の導出の前に、結晶構造の2次元模式図を示します。各ベクトルの意味は図で示されています。
\(\vec{R^{a}}\): 基準セルの原点からユニットセル\(\mathit{a}\)の原点へのベクトル
\(\vec{R_{k}^{a}}\): 基準セルの原点からユニットセル\(\mathit{a}\)の原子\(\mathit{k}\)へのベクトル
\(\vec{R_{k}}\): 原点から同じセル内の原子\(\mathit{k}\)へのベクトル
調和振動子近似 (Harmonic approximation)
前章の振動解析同様、原子が平衡点から微小変位した際のポテンシャルエネルギーは以下のように表されます。
\(\vec {u^a}\) は原子の変位 \(\vec {R^a} - \vec {R_0^a}\)で、 \(\alpha\) と \(\beta\) は座標軸x, y, zを表します。
[Note] 調和振動子近似は平衡点からの微小変位でのみ成り立つ近似です。そのため、フォノン解析で得られる熱物性は低温領域に限られます。
運動方程式
原子\(\mathit{k}\)に作用する力はポテンシャルエネルギーの微分として以下のように得られます。
原子\(\mathit{k}\)の加速度
よって、運動方程式は以下となります。
各原子に対して、3 つの運動方程式が得られます(x, y, z軸方向)。 結果、全体では \(3N \times M\) 個の方程式となります (\(N\) はユニットセル内の原子数、\(M\) はユニットセル数)。
これらは連立方程式として結合しています。ユニットセル\(a\) の原子\(k\)の変位\(u_{k}^{a\alpha}(t)\)は、他の全ての原子変位\(u_{k'}^{b\beta}\)に依存します。
仮定
上記式の解を得るため２つの仮定を導入します。
\(u_{k}^{a\alpha}(t)\) は以下の形式(平面波)で書けるとします。
\[u_{k}^{a\alpha}(t) = A_{k}^{\alpha} e^{-i\omega t}\]
方程式の解が複数ある可能性を考慮し、添字 \(\mathit{m}\)でそれぞれの解を表します。
周期境界条件を考慮し、解は平面波に似た形式になるとします。
上記の仮定から、方程式の解が新たなパラメータ\(\vec{q}\) に依存することになります。これが波数ベクトルとなります。
方程式を解く
上記仮定から運動方程は以下のように書き換えられます。
ここで、前章の振動解析同様、\(3N \times 3N\)行列となる力定数マトリクスを定義します。 しかし、周期境界条件を考慮して、ユニットセル \(a\) とユニットセル \(b\)の間の力定数マトリクスは以下のようになります。
結晶格子の並進移動対称性により、 \(C^{a,b} = C^{0, b-a}\) となります。
ここまでの式を代入すると、運動方程式の […] 内の部分は、力定数マトリクスの離散フーリエ変換の形になっていることがわかります。
さらに運動方程式を単純化するため、 \(3N \times 3N\) 行列である dynamics matrix を新たに定義します。
すると、波数ベクトル \(q\) に対する運動方程式は、同次線形方程式の形になります。
ここで、 \(\vec{\gamma}_{\vec{q}}\) は固有ベクトルで、各原子の振動に対し、以下のように関連付けられます。
Dynamical matrix \(\boldsymbol{\tilde{D}_{\vec{q}}}\)の対角化を行うことで、波数ベクトル\(\vec{q}\)に対応する \(3N\) 個の振動数 \(\omega\) を得ることができます。
5.2 波数ベクトル \(\vec{q}\) が第一ブリルアンゾーンに制限される理由¶
波数ベクトル \(\vec{q}\) が第一ブリルアンゾーンに制限される理由を簡単な１次元の場合で説明します。 本質的には、結晶上で格子点が離散的に配置されていることから来ています。
1次元格子中の原子\(a\)の変位は以下で表されます。
原子\(a\) の座標が \(x=0.0\)で、格子定数が1.0 としたときは、逆格子空間は\(G \in (... -2\pi, 0, 2\pi, 4\pi, ...)\)であり、第一ブリルアンゾーンは\([-\pi, \pi]\)となります。
次の図は、\(u_{q}^{a}(t)\)が\(q=-\pi/3\) (青線), \(q=-\pi/3+2\pi\) (緑線), \(q=-\pi/3+4\pi\) (黄色線)の場合を表しています。 赤点で示される原子\(a\) の変位は3つの波で全く同じになっていることがわかります。
\(u_{q}^{a}(t) = u_{q+G}^{a}(t)\)
つまり、任意の波数ベクトル\(\vec{q}\)に対して、第一ブリルアンゾーン内に等価な原子変位\(u_{q}^{a}(t)\)を表す波数ベクトル \(\vec{q}'\) が存在します。
5.3 フォノンの状態密度と熱物性値¶
プランク分布 (Planck distribution)
固体中では、たくさんのフォノンが同時に励起しています。 フォノン解析を用いて熱物性を得る際は、フォノンの分布を知る必要があります。 フォノンのエネルギーは\(E=\hbar\omega\)で与えられます。 フォノンのモード (\(\vec q\) と \(\omega\)固定)が与えられたとき、そのフォノンが \(n\) 個存在する確率は\(P(n)=\frac{1}{Z}e^{-\frac{n\hbar\omega}{k_B T}}\)で表されます、ここで \(Z\) は規格化定数です。
よって、温度\(T\)におけるフォノンの個数の平均値は、
と表されます。これはプランク分布として知られています。以下にプランク分布のプロットを示します。
左図から、低振動数となるフォノンは高振動数となるフォノンよりも多く存在している事がわかり、その傾向は低温で特に顕著になることもわかります。
右図から、フォノンの数は温度とともに増加することがわかります。
振動数\(\omega\)のフォノンにより運ばれるエネルギーは、以下のように計算されます。
フォノンの状態密度 (Phonon DOS)
フォノンの状態密度は結晶の物性のなかでも重要な要素で、以下のように定義されます。
\[D(\omega) = \frac{1} {N} \sum_{m, \vec q} \delta (\omega - \omega_{m, \vec q})\]\(m\) はバンド (モード)の添字で、波数ベクトル \(\vec q\) は第一ブリルアンゾーン内で足される必要があります。 DOS \(D(\omega)\) は、次の制約に従います。
\[N =\int_{0}^{\infty} D(\omega) d \omega\]\(N\) はフォノンの総数です。
もしフォノンの状態密度 \(D(\omega)\)がわかっていれば、フォノンが寄与するエネルギーは以下のように計算できます。
\[U' = \int_{0}^{\infty} d {\omega} D(\omega) E_{\omega} = \int_{0}^{\infty} d {\omega} D(\omega) \frac{\hbar \omega}{e^{\frac{\hbar \omega}{k_B T}}-1}\]フォノンの零点エネルギー (zero-point energy)は量子系の調和振動子と同様 \(\frac {\hbar \omega}{2}\)となります。 この零点エネルギーを含めると、結晶の内部エネルギーは、
\[U = \int_{0}^{\infty} d {\omega} D(\omega) \; \left[ \frac{\hbar \omega}{2} + \frac {\hbar \omega}{e^{\frac{\hbar \omega}{k_B T}} -1} \right]\]結晶の多くの物性値は内部エネルギー\(U\)と関連付けられます。例えば、熱容量 \(C_V\) は、
本チュートリアルでは、比熱 \(C_V\)・エントロピー \(S\)・内部エネルギー \(U\)・ヘルムホルツ自由エネルギー \(A\)もフォノンの状態密度から計算できることを紹介しました。
フォノン状態密度の計算
フォノンの状態密度はどのように計算できるのでしょうか？
まず、フォノンの振動数を第一ブリルアンゾーン内の等間隔メッシュ \(\vec q\) 上で計算します。メッシュのとり方は、Monkhorst-Pack法とGamma-centered法がよく使われます。
次に、フォノンの振動数を第一ブリルアンゾーン上で積分します。ここでは、smearing methodかtetrahedral integral methodが代表的に使われています。
Tetrahedron methodでは、メッシュのが四面体で区切られます。それぞれの四面体で、フォノンの振動数が線形補間され、解析的にphonon DOSを得られます。 [3]
Smearing methodでは、メッシュグリッド\(\vec q\) 上で計算された全てのフォノン振動数に対して、gaussian smearing function が適用されます。