# 計算科学基礎講座 : 熱電研究のための第一原理計算入門
第 3 回 Boltzmann 輸送方程式を用いた熱電特性計算

桂 ゆかり (東京大学)
## 1.はじめに

熱電材料の開発研究はいつも，「あちらを立てればこ ちらが立たず」の連続である。このため，全体像が見え ない状態を膨大な実験で突破することにより，多数のパ ラメータを同時に最適化して特性を改善することを強い られる。

そのような中，第一原理計算による熱電特性予測は，暗い森の中から, ゴールまでの一本道を照らし出す一縷 の光のように感じられる. $Z T$ の値を決定するゼーベッ ク係数 $S$, 電気伝導率 $\sigma$, 電子熱伝導率 $\kappa_{\mathrm{el}}$ は, すべて電子構造とキャリア濃度に依存して決定するパラメータで あり，Boltzmann 輸送方程式によって描写できる。フォ ノン熱伝導率 $\kappa_{\mathrm{ph}}$ の第一原理計算も順調に研究が進んで いる.このため，第一原理計算によって，各熱電材料の $Z T$ の予測や，既存材料より高い $Z T$ を示す物質の発見 ができるのではないかと期待させられる.

ところが，第一原理計算による熱電特性の予測には， その過程で多くの簡略化と不確定性が入り込む．結晶構造を定義する段階で為される簡略化もあれば，現在の 第一原理計算手法に内在する不確定性もある. さらに， Boltzmann 輸送方程式から得られる情報量にも限界があ る.このため現在の計算技術では, 第一原理計算を $Z T$值の予測までつなげることは，困難である。

連載の最終回にしていきなり冒頭から期待を打ち砕く のもどうかと思う．だがこの認識があってこそ，第一原理計算の熱電材料研究への活用が可能になると考えてい る. 第一原理計算から得られる情報は, 結晶構造から得 られる情報である.ゴールまでの一本道は見えなくとも，辺りを明るく照らし出すだけの光量は持っているはずで ある。

本稿では, Boltzmann 輸送方程式を解いて得られる情報と，それらの分析法について説明する。各熱電特性の キャリアドープ量 $n$, 温度 $T$, 電子構造依存性について,実在物質を計算例として解説を行う。また計算結果から 未知変数をできるだけ括り出すことで，その物質の熱電特性の特徴を概観する方法を説明する．さらに，実験値 と計算值の比較から，電子散乱に関する情報をつかむ手段を紹介する．論文ではなく解説記事であるという特徴 を生かして，WIEN2k ${ }^{11}$ などの第一原理計算ソフトから 熱電特性を計算できる，BoltzTraP2) というフリーウェ アについて，具体的に説明する。ただしこの情報は, 現在公開されているバージョン（1.2.5）に基づいており,今後変更される可能性もある。これにより，できるだけ 多くの実験系研究者に第一原理計算を活用していただ き, 熱電研究のさらなる発展に貢献できれば幸いである.

2. 熱電材料の第一原理計算に打ける不確定性

2. 1. 結晶構造の不確定性

多くの熱電材料は，純䊀な物質を測定しても非常に低 い $Z T$ しか得られないことが多い. 何も考えずに合成し ても，自然に形成する格子欠陥などにより，たまたま最適な $n$ の試料が得られる物質も存在する。しかし多く の場合，キャリアドープや $\kappa_{\mathrm{ph}}$ の低減を目的として，さ まざまな元素や格子欠陥が高濃度で導入されて初めて，高い $Z T$ が達成できる。

ところがそれらの高 $Z T$ 試料を呼ふ際に，組成を全て 言うのは煩雑なので，つい関連する母物質の名前で呼ん でしまう。例えば, 我々が $\left[p\right.$ 型の $\left.\mathrm{Bi}_{2} \mathrm{Te}_{3}\right.\rfloor$ と呼ぶ物質には， Bi サイトの $75 \%$ \% Sbで置換した試料も含まれる。 その 電子構造は純桲な $\mathrm{Bi}_{2} \mathrm{Te}_{3}$ の電子構造とは異なると考え られる.

また，キャリアドープを目的に添加された元素は，電気伝導に関与する最も重要なエネルギー領域（フェルミ エネルギー $\epsilon_{\mathrm{F}}$ 近傍の電子状態）に新たなバンドをつく る. 半導体デバイスのような希薄ドープなら無視できて も，熱電材料のような大量ドープでは無視できない. よって,より実像に近い電子構造を得るには, 図 1 (a) の母物質ではなく，（c）のように置換元素のランダム分布まで再現する結晶構造で計算する必要がある. しかし これには多大な計算負荷が必要な上, バンド構造の解釈 も困難になる。このため多くの場合，（b）のような単純 な結晶構造を定義して計算される。しかしこれでは, 置換元素が空間内に規則的に整列していることになってし まう。

(a) parent compound

(b) simple supercell

(c) more complex supercell図 1 固溶体の結晶構造の定義法.

このほか，高 $Z T$ 材料にはさまざまな格子欠陷や析出物，粒界などが含まれる。これらを全て織り込んだ結晶構造を定義することは，非常に困難であると言えよう．

## 2. 2. 電子構造の不確定性

いくら現実に近い結晶構造を定義したとしても，第一原理計算の過程で生じる誤差もある。

その最たるものは，半導体のバンドギャップ $\epsilon_{\mathrm{g}}$ の詋評価である。これは，現在の密度汎関数理論（DFT）に よる第一原理計算の弱点である. LDA P $\mathrm{GGA}^{3}$ ) など,計算速度が速く, よく用いられる交換相関ポテンシャル では， $\epsilon_{\mathrm{g}}$ は過小評価される傾向がある. $\epsilon_{\mathrm{g}}$ が $2 / 3$ 程度 に減少してしまう半導体もあれば， $\epsilon_{\mathrm{g}}=0$ の金属（半金属) と判定されてしまう半導体もある. 熱電材料の多く がナローギャップ半導体であることを考慮すると，この 誤差は深刻である

ただし，より高精度の $\epsilon_{\mathrm{g}}$ が得られる交換相関ポテン シャルを利用すれば， $\epsilon_{\mathrm{g}}$ の計算値は改善する. 計算量 は膨大だが， GW 法4) や，HSE 法5) などがよく知られて おり，WIEN2kでは比較的計算の軽い， TB-mBJ $\left.{ }^{6}\right)$ もよ く用いられている. ただし，いずれも歴史がまだ浅いた め, 有効質量など他の計算値も正確になるのか, 研究動向の見極めが重要である.
LDA+Uなど，電子を各原子に局在させるようなクー ロンポテンシャル $U$ を導入して $\epsilon_{\mathrm{g}}$ を調整する手法では，物質にもよるが有効質量を過大評価する傾向がある。そ れよりも，LDA や GGA で計算した電子構造の伝導帯の $\epsilon$ を単純にシフトし， $\epsilon_{\mathrm{g}}$ を人為的に実験值に合わせた方 が，有効質量も正確になる傾向がある。この方法は $\epsilon_{\mathrm{g}}$ の計算值が 0 では使えないものの, BoltzTraP 内で簡単 に実行できる，最も時間のかからない方法と言える。

このほか，第一原理計算の過程で起こりうる問題とし て，その物質に本来存在すべき相互作用を無視して計算 してしまった場合が考えられる。磁性のある物質をスピ ン偏極なしで計算してしまった場合や，重元素からなる 物質をスピンー軌道相互作用（SO）なしで計算してし まった場合などが該当する。

格子定数や原子座標の最適化が十分でないと，異な る電子構造になる。ただし, 構造最適化の際, van der Waals 結合をもつ物質では，結合距離が間違った値に収束しやすい．格子定数は温度によっても変化する。

十分に細かい $k$ 点メッシュの使用が望ましい. ただし，半導体では, 金属の計算ほど多くの $k$ 点が必要なわけで はない. $k$ メッシュを少しずつ細かくしながら計算して いき, 全エネルギーや目的の計算値に変化やノイズが現

## 3 . 熱電材料の状態密度曲線

状態密度曲線 (DOS) は, 多くの論文に掲載されるデー タであるため，それを読むことでその化合物の熱電材料 としての有望性を推定できれば便利である.

図 2 単位体積当たりの值に規格化した, 熱電材料の状態密度曲線の例. 縦軸の 1 目盛が 1 に相当する. (a) $s p$ 半導体 (b) $3 d$ 遷移金属半導体 (c) 酸化物.
例として，WIEN2k によって計算したさまざまな熱電変換材料の状態密度（DOS）を図 2 に示す. 化合物間で DOS の勾配を比較するため，格子体積で割って規格化 している.

$\epsilon_{\mathrm{g}}$ の計算値は選択した交換相関ポテンシャルによっ て大きく異なる．WIEN2k のデフォルトであり，他の計算コードでも広く利用されているGGAで計算した $\epsilon_{g}$ は，実験値を大きく下回った。一方，TB-mBJ を利用した場合は，収束までにより長時間を要するが，実験値と近い 值が得られた。

$\mathrm{Bi}_{2} \mathrm{Te}_{3}, \mathrm{PbTe}$ などの重元素化合物では，スピン軌道相互作用（SO）の有無によって計算結果が大きく変わる. TB-mBJ で SO を入れた場合よりも， GGAでSO を入れ た場合の方が実験値に近い $\epsilon_{\mathrm{g}}$ が得られているが，どの 計算がより正しい電子構造を反映しているのかは不明で ある.

DOS の勾配を化合物間で比較する際は，単位体積当 たりの DOS として規格化する必要がある. $\mathrm{Si}, \mathrm{Mg}_{2} \mathrm{Si}$ な どの $s p$ 半導体は，バンド端の DOS の勾配がゆるやかで ある。一方，遷移金属化合物では，DOS の勾配は非常 に急峻になる．特に，局在性の強い $3 d$ 遷移金属の場合 に顕著である。また，酸化物では，価電子帯（ $\mathrm{p}$ 型側） の DOS の勾配が大きくなる傾向がある．これは， $\mathrm{O}^{2-}$ イオンの電子がほぼ局在しているうえに，同じ電子状態 に置かれた O の数が非常に多く，バンドが縮退してい るためである.

DOS で着目すべき点は 2 点である. 1 つは，バンド ギャップ $\epsilon_{\mathrm{g}}$ の有無とその幅. もう 1 つは，フェルミエ ネルギー $\epsilon_{\mathrm{F}}$ における状態密度（DOS）の勾配. この 2 点が熱電特性を左右する.

## 4. Boltzmann 輸送方程式による熱電特性の計算

ここでは，第一原理計算から得られた結果を使って，熱電特性にかかわる各パラメータを計算する原理につい て紹介する. BoltzTraPではテンソルを用いて各パラメー タが計算されているため，はじめにテンソル計算の基礎 について解説する. 次に，熱電現象を記述する，電流と 熱流の方程式を紹介する。そして，Boltzmann 輸送方程式の線形解として得られる $\sigma, \mathrm{S}, \boldsymbol{\kappa}_{\mathrm{el}}$ の計算式を紹介す る.

4. 1. テンソル

物質の異方性を取り込んで計算を行うには，テンソル を利用する．たとえば，電気伝導率テンソル $\sigma$ は以下の ように 3 行 3 列の行列で表される.

$
\boldsymbol{\sigma}=\left(\begin{array}{lll}
\sigma_{x x} & \sigma_{x y} & \sigma_{x z} \\
\sigma_{y x} & \sigma_{y y} & \sigma_{y z} \\
\sigma_{z x} & \sigma_{z y} & \sigma_{z z}
\end{array}\right)
$

ここで，添字の 1 つ目は出力方向，すなわち発生する 電流の方向を示す. 2 つ目は入力方向，すなわち電界が かけられる方向を示す.

電界ベクトル $\mathbf{E}=\left(E_{x}, E_{y}, E_{z}\right)$ をかけた時に発生する電流密度べクトル $\mathrm{j}_{\mathrm{el}}$ は, 以下のような掛け算で計算され る.

$
\mathbf{j}_{\mathrm{el}}=\sigma \mathbf{E}=\left(\begin{array}{ccc}
\sigma_{x x} & \sigma_{x y} & \sigma_{x z} \\
\sigma_{y x} & \sigma_{y y} & \sigma_{y z} \\
\sigma_{z x} & \sigma_{z y} & \sigma_{z z}
\end{array}\right)\left(\begin{array}{l}
E_{x} \\
E_{y} \\
E_{z}
\end{array}\right)
$

ただし $\sigma$ の場合, 対角項以外はほぼ 0 である. よって

$
\mathbf{j}_{\mathrm{el}}=\left(\begin{array}{c}
j_{x} \\
j_{y} \\
j_{z}
\end{array}\right) \sim\left(\begin{array}{c}
\sigma_{x x} E_{x} \\
\sigma_{y y} E_{y} \\
\sigma_{z z} E_{z}
\end{array}\right)
$

と表してもほぼ同じである.

電場と磁場が共存する場合の $\sigma$ は, 出力である電流方向に加えて，入力として電場と磁場の方向が必要となる ため，各成分が添字 3 つを持った $3 \times 3 \times 3$ のテンソルと なる.この場合, 電流と電場が平行な成分 $\left(\sigma_{x x x}, \sigma_{x x y}\right.$ など)以外に, 添字が全て異なる成分 $\left(\sigma_{z x y}\right.$ など $)$ も有意な値 を持ち，これらは Hall 効果による誘導電流を表現する。

得られたテンソルからスカラーの代表値を得る最も単純な方法は，その物性を表す有意な項の平均を取ること である. BoltzTraPの case.trace ファイルに出力される のはこの平均値である．平均を取った時点で，異方性に 関するデータは失われている。

4.2. 電流と熱流の方程式)

電流 $\mathbf{j}_{\mathrm{el}}$ が流れ，磁場がなく，温度勾配 $\nabla T$ が存在する 試料内に発生する電界べクトル $\mathrm{E}$ は，以下のように表 せる.

$
\mathbf{E}=\sigma^{-1} \mathbf{j}_{\mathrm{el}}+\mathbf{S} \nabla T
$

$
\mathrm{j}_{\mathrm{el}}=\sigma \mathrm{E}-\sigma \mathrm{S} \nabla T
$

このとき，電子の群速度は

$
v=\sigma \mathrm{S}
$

である.

また，Peltier 係数を $\mathrm{S} T$ ，電子の拡散による熱伝導率 を $\boldsymbol{\kappa}_{0}$ とおくと，以下の熱流の方程式が得られる。

$
\mathbf{j}_{Q}=\mathbf{S} T \mathbf{j}_{\mathrm{el}}+\boldsymbol{\kappa}_{0} \nabla T
$

ここに（5）式を代入すると，

$
\mathbf{j}_{Q}=\mathrm{S} \sigma T \mathrm{E}+\left(\boldsymbol{\kappa}_{0}-\mathrm{S} \sigma \mathrm{S} T\right) \nabla T
$

が得られ, ここから, 電子熱伝導率

$
\boldsymbol{\kappa}_{\mathrm{el}}=\boldsymbol{\kappa}_{0}-\mathbf{S} \boldsymbol{\sigma} \mathbf{S} T=\boldsymbol{\kappa}_{0}-\mathbf{S} \boldsymbol{\nu} T
$

が導かれる.これは, 試料内の電場が 0 で, 温度勾配の みが存在するときに電子が運ぶ熱を表す.

第 1 項は，熱拡散した電子が運ぶ熱を表す。しかし，高温側の電子が低温側に拡散する過程は, 低温側への「電流」である．試料内の電場が 0 のときに，試料内に流れ ている電流は，（5）式より $-\sigma \mathrm{S} \nabla T$ である. よって，こ の電流と Peltier 係数 $\mathrm{S} T$ の積が，熱拡散による Peltier 効果であり，その熱伝導率は $-\mathrm{S} \sigma \mathrm{S} T$ と表せる. すなわ ち，拡散した電子が自己 Peltier 効果によって，高温端 の温度を高く，低温端の温度を低くする方向に温度差を 作り出し，熱の移動量が減少しているのである。ただし 熱力学第 2 法則により， $\kappa_{\mathrm{el}}$ が負になることはない.

金属では第 2 項は無視できるほど小さいが，パワー ファクター $P=S^{2} \sigma$ が十分に大きな半導体や熱電材料で は第 2 項も大きいと考えられる. $\kappa_{\mathrm{el}}$ が 0 になるのは $\sigma(\epsilon)$ がデルタ関数であるときであると示されている 8 . この 第 2 項なしでは多くの $Z T>1$ の熱電材料の存在が説明 できないことについて，後程紹介する。

## 5. Boltzmann 輸送方程式の線形解

ある電子の有効質量を $m^{*}$ とおくと，運動量は $p=m^{*} \nu=$ $\hbar k$ ，運動エネルギーは

$
\epsilon=\frac{p^{2}}{2 m^{*}}=\frac{\hbar^{2} k^{2}}{2 m^{*}}
$

と表される．よって，これを $k$ で微分することにより，

$
\nu=\frac{1}{\hbar} \frac{d \epsilon}{d k}
$

が得られる. これより， $\nu$ は $\epsilon-k$ 曲線の傾きに対応する ということがわかる.

$v$ はその電子が属するバンド番号 $i$ と，波数 $\mathrm{k}$ によっ て異なる．第一原理計算によって求められた電子（ $i, \mathbf{k})$ のエネルギーを $\varepsilon(i, \mathbf{k})$ とおくと，その速度べクトルは $\boldsymbol{V}_{x}(i, \mathbf{k})=\left(\begin{array}{c}\nu_{x}(i, \mathbf{k}) \\ \nu_{y}(i, \mathbf{k}) \\ \nu_{z}(i, \mathbf{k})\end{array}\right)=\frac{1}{\hbar}\left(\begin{array}{c}\frac{\partial \epsilon(i, \mathbf{k})}{\partial k_{x}} \\ \frac{\partial \epsilon(i, \mathbf{k})}{\partial k_{y}} \\ \frac{\partial \epsilon(i, \mathbf{k})}{\partial k_{z}}\end{array}\right)$

として求められる.

ところで，試料内を一定の電流が流れる状態では，電場による電子の加速が散乱によって完全に打ち消され て，vが時間変化しない．このとき，電子が散乱を受け る平均時間間隔を緩和時間 $\tau_{\mathrm{el}}$ と定義すると， $\sigma$ の各成分は

$
\sigma_{x y}(i, \mathbf{k})=e^{2} \tau_{\mathrm{el}} v_{x}(i, \mathbf{k}) v_{y}(i, \mathbf{k})
$

などと計算できる。ただし, 実際に有意な値をもつのは，

$
\sigma_{x x}(i, \mathbf{k})=e^{2} \tau_{\mathrm{el}}\left(v_{x}(i, \mathbf{k})\right)^{2}=\frac{e^{2} \tau_{\mathrm{el}}}{\hbar^{2}}\left(\frac{\partial \varepsilon(i, \mathbf{k})}{\partial k_{x}}\right)^{2}
$

などの対角項のみである. $\tau_{\mathrm{el}}$ の値は実際には $i, \mathrm{k}$ に依存すると考えられるが, ここでは定数 $\tau_{\mathrm{el}}$ であると仮定 している（緩和時間近似）。

ところで，電子ドープは $\epsilon_{\mathrm{F}}$ の上昇，ホールドープは $\epsilon_{\mathrm{F}}$ の低下として表される. そこで， $\sigma$ の式を $\epsilon_{\mathrm{F}}$ の関数 に書き換えるため, エネルギー $\epsilon_{\mathrm{F}}$ から $\epsilon_{\mathrm{F}}+\Delta \epsilon$ の間に含 まれている状態 $(i, \mathbf{k})$ を全て検出し, それらの $\sigma_{x x}(i, \mathbf{k})$ の和をとる. 上記の $\Delta \epsilon$ の範囲に $\epsilon(i, \mathbf{k})$ が含まれていた 時に 1 を返し，しなかった場合には 0 を返すデルタ関数 $\delta\left(\epsilon_{\mathrm{F}}-\epsilon(i, \mathbf{k})\right)$ を用いると，

$
\sigma_{x x}\left(\epsilon_{\mathrm{F}}\right)=\frac{\tau_{\mathrm{el}}}{N_{\mathrm{k}} \Delta \epsilon} \sum_{i, \mathrm{k}} \sigma_{x x}(i, \mathbf{k}) \delta\left(\epsilon_{\mathrm{F}}-\epsilon(i, \mathbf{k})\right)
$

と表される.この $\boldsymbol{\sigma}\left(\epsilon_{\mathrm{F}}\right)$ は輸送分布関数と呼ばれる. ま

図 3 複数のバンドが存在する系における，（a）DOS: $N\left(\epsilon_{\mathrm{F}}\right)$ (b) スペクトル伝導度 : $\sigma_{x x}\left(\epsilon_{\mathrm{F}}\right)$ の計算法. DOS は，微小エネルギー領域 $d \epsilon$ に含まれる状態 の数そのものを数えたものである. $\sigma\left(\epsilon_{\mathrm{F}}\right)$ はこれ らの状態における $\epsilon-k$ 曲線の勾配の自乗の和から 計算される. このため, $\sigma_{x x}\left(\epsilon_{\mathrm{F}}\right)$ では勾配の急峻 なバンド（band 1）の寄与が大きくなる.

た，エネルギーの関数として表した物理量は「スペクト ル」であるため，スペクトル伝導度とも呼ばれる。

$\sigma\left(\epsilon_{\mathrm{F}}\right)$ の計算精度は $k$ 点数 $N_{\mathrm{k}}$ が多いほど上がる.こ れは， $\Delta \epsilon$ 内に検出される点の数が増えるからである. このため Boltzmann 輸送方程式の解析には通常非常に細 かい $k$ メッシュを用いる．ただし，BoltzTraPではバン ドエネルギーを Fourier 展開によって補間し, $k$ 点数を 5 倍など（case.intransにおける lpfac で指定した值）に 増やしているため, 多少粗い $k$ メッシュでも精度の良い $\sigma\left(\epsilon_{\mathrm{F}}\right)$ が得られる.

ところでここまでの結果は， $T=0$ の場合である．有限温度では，Fermi-Dirac 分布関数

$
f\left(\epsilon_{\mathrm{F}}, T\right)=\left(e^{\frac{\epsilon-\epsilon_{\mathrm{F}}}{k_{\mathrm{B}} T}}+1\right)^{-1}
$

を用いて

$
\sigma\left(\epsilon_{\mathrm{F}}, T\right)=\int_{0}^{\infty} \sigma(\epsilon) f\left(\epsilon_{\mathrm{F}}, T\right) d \epsilon
$

と表される。

この積分値の計算には, $\epsilon=\epsilon-\epsilon_{\mathrm{F}}$ における Taylor 展開

$
\left.\sigma(\epsilon, T)\right|_{\epsilon=\epsilon_{F}}=\sum_{n} \frac{1}{n !} \int_{0}^{\infty}\left(\epsilon-\epsilon_{\mathrm{F}}\right)^{n} \sigma\left(\epsilon_{\mathrm{F}}\right)\left(-\frac{\partial f_{\epsilon_{\mathrm{F}}}(\epsilon, T)}{\partial \epsilon}\right) d \epsilon
$

が用いられる. $\sigma\left(\epsilon_{\mathrm{F}}, T\right)$ はこの 0 次の項から，

$
\left.\sigma(\epsilon, T)\right|_{\epsilon=\epsilon_{F}}=\int_{-\infty}^{\infty} \sigma(\epsilon)\left(-\frac{\partial f_{\epsilon_{\mathrm{F}}}(\epsilon, T)}{\partial \epsilon}\right) d \varepsilon
$

と計算される. 以降の導出過程は複雑なので省略するが，

$
\begin{aligned}
\left.\boldsymbol{v}(\epsilon, T)\right|_{\epsilon=\epsilon_{F}}= & \frac{1}{e T} \int_{-\infty}^{\infty}\left(\epsilon-\epsilon_{\mathrm{F}}\right) \boldsymbol{\sigma}(\epsilon)\left(-\frac{\partial f_{\epsilon_{\mathrm{F}}}(\epsilon, T)}{\partial \epsilon}\right) d \varepsilon, \\
\left.\mathrm{S}(\epsilon, T)\right|_{\epsilon=\epsilon_{F}}= & \left.\int_{-\infty}^{\infty} \sigma^{-1}(\epsilon, T) \boldsymbol{v}(\epsilon, T) d \epsilon\right|_{\epsilon=\epsilon_{F}}, \\
\left.\boldsymbol{\kappa}_{\mathrm{el}}(\epsilon, T)\right|_{\epsilon=\epsilon_{F}}= & \frac{1}{e^{2} T} \int_{-\infty}^{\infty}\left(\epsilon-\epsilon_{\mathrm{F}}\right)^{2} \boldsymbol{\sigma}(\epsilon)\left(-\frac{\partial f_{\epsilon_{\mathrm{F}}}(\epsilon, T)}{\partial \epsilon}\right) d \varepsilon \\
& -\left.T \int_{-\infty}^{\infty} \mathrm{S}(\epsilon, T) \boldsymbol{v}(\epsilon, T) d \epsilon\right|_{\epsilon=\epsilon_{F}}
\end{aligned}
$

が得られる. これらの式を構成する各関数の挙動につい

(15) 式からわかるように， $\sigma(\epsilon)$ には未知定数 $\tau_{\mathrm{el}}$ が 含まれている. このため, $\sigma(\epsilon)$ から計算される $\sigma\left(\epsilon_{\mathrm{F}}, T\right)$, $\boldsymbol{v}\left(\epsilon_{\mathrm{F}}, T\right), \boldsymbol{\kappa}_{\mathrm{el}}\left(\epsilon_{\mathrm{F}}, T\right)$ も全て $\tau_{\mathrm{el}}$ を含んでいる. ただし,緩和時間近似が有効であれば，計算過程において $\tau_{\text {el }}$ が 相殺できるため， $\mathrm{S}\left(\epsilon_{\mathrm{F}}, T\right)$ は値として得ることができる.

ただし，これらの計算ではマルチバンド効果を十分考慮できていない，これは，（15）式の時点で，バンド番号 $i$ の情報が失われているからである.
6. BoltzTraP の出力データの意味

6. 1 . $\epsilon_{\mathrm{F}}$

バンド計算では，絶対零度において伝導電子に占有さ れた上端のエネルギーをフェルミエネルギーと呼ぶ.こ れは， $T=0 \mathrm{~K}$ を仮定して計算される。一方， $T>0 \mathrm{~K} に$ おける Fermi-Dirac 分布関数の中心 (状態の占有率が 0.5 になるエネルギー) は化学ポテンシャルと呼ばれる.フェ ルミ準位は化学ポテンシャルに対応するが，人によって 定義が異なることもある.

フェルミエネルギーと化学ポテンシャルは金属ではほ ぼ一致する. しかし半導体では, フェルミエネルギーは 価電子帯の上端, 化学ポテンシャルはバンドギャップ内 のどこかとなり，異なる値をとる. case.trace に打ける 1 コラム目の $\epsilon_{\mathrm{F}}$ は化学ポテンシャ ルに対応する任意変数である. 単位は Rydberg（1 Ry= 1/2 Hartree $=13.605698066 \mathrm{eV}$ ）である.

半導体や不純物ドープを行った金属では， $\epsilon_{\mathrm{F}}$ を第一原理的に求める手段がない. このため他のカラムのデー タ $\left(n, S, R_{\mathrm{H}}\right.$ など）を見ながら，最も実験データを再現 する $\epsilon_{\mathrm{F}}$ を選択しなければならず，ここに計算結果の解釈の任意性が生じる.

## 6. 2 . $T$

仮想的な温度 $[\mathrm{K}]$ であり, Fermi-Dirac 分布関数の変化を再現する．格子定数の温度依存性は考慮されていな い. $\epsilon_{\mathrm{F}}$ の温度依存性を無視できるならば， $\epsilon_{\mathrm{F}}$ が等しく， Tのみが異なる 1 ブロックをそのまま使えばよい。こ のような近似は金属や縮退半導体など， $n$ の大きな系に は有効である．1 温度についての計算を行うには, case. intrans ファイルの 8 行目で $T_{\max }$ と temperature grid を同 じ值に設定する.

## 6. 3 . $N$

単位格子（reduced cell）に入っている過剩な電荷の 数を示す. キャリアドープを目的として不純物ドープを 行った場合に，組成比と対応させやすい值である.

単位体積当たりのキャリア濃度 $n\left[\mathrm{~cm}^{-3}\right]$ に換算す るには, reduced cell の体積 $V_{\text {red }}$ で割ればよい. $V_{\text {red }}$ は case.outputtrans に bohr ${ }^{3}\left(\right.$ a.u. $\left.{ }^{3}\right)$ 単位で書かれており, 1 bohr $=0.52917721092$ Åである. ホールの濃度を $n_{\mathrm{h}}$, 電子の濃度を $n_{\mathrm{e}}$ とすると, $n=n_{\mathrm{h}}-n_{\mathrm{e}}$ である. この $n$ は, Hall 測定から得られるキャリア濃度とは異なり, $n_{\mathrm{h}}=n_{\mathrm{e}}$ であれば 0 となる. $n_{\mathrm{H}}$ がキャリアの熱励起の影響を受 けやすいのに対し， $n$ の温度依存性は小さいと考えられ る.

## 6. 4. DOS

単位格子あたりの状態密度で, 単位は (states /Ry) である. BoltzTraPでは， $\epsilon_{\mathrm{F}}$ 近傍のエネルギー範井に存在するバンドのみを抽出している。

## 6. 5 . $S$

Seebeck 係数の計算值であり，単位は V/Kである．図 4 に， $\epsilon_{\mathrm{F}}$ を変えた $n$ 型 $\mathrm{Si}$ の $S-T$ 曲線を示すが，多くの 熱電材料において， ドーパント添加量を変化させたとき に得られる $S-T$ 曲線とよく似た挙動が見られる。 $\epsilon_{\mathrm{F}}$ が 伝導帯の下端にある場合はほぼ平坦となり， $\epsilon_{\mathrm{F}}$ が伝導帯内にある場合には， $S$ が低温で小さくなるという金属的挙動を示した． $\epsilon_{\mathrm{F}}$ をバンドギャップ内に置いた場合 は, $S$ が低温で発散するという半導体的挙動が見られた。 ただし，この場合は $\epsilon_{\mathrm{F}}$ に大きな温度依存性があるため， これらの $S-T$ 曲線は実験とは一致しないと考えられる.

図 5 は異なる $T$ について計算した $p$ 型 $\mathrm{Si}$ の $S-n$ 曲線 である. $300 \mathrm{~K}$ の $S$ は，低 $n$ 側で非常に大きく，nの上昇とともに減少した. 同じ $n$ では, $T$ の上昇に伴い $S も$大きくなった．これらの傾向は，自由電子モデルに打け る $S$ の式

$
S=\frac{2 k_{\mathrm{B}}{ }^{2}}{3 \mathrm{e} \hbar^{2}}\left(\frac{\pi}{3 n}\right)^{2} m^{*} T
$

から予測される傾向と一致している.

高温では，低 $n$ 領域から $S$ が低下し， $S-n$ 曲線にピー クが現れた.これはキャリアの熱励起により, 電子とホー ルのぺアが大量に生成し，系が金属的になった効果であ る.これは両極性拡散効果 (bipolar conduction) と呼ばれ,高温になるほど高 $n$ 領域の $S$ まで影響が及ぶことがわ

図 $4 n$ 型 $\mathrm{Si}$ について計算した，さまざまな $\epsilon_{\mathrm{F}}$ に対す る $S$ の温度依存性. 図中の数字は, 伝導帯のバ ンド端からのエネルギーを $\mathrm{eV}$ 単位で示したもの である.
かる.

6. 6. $\sigma / \tau_{\text {el }}\left[\Omega^{-1} \mathrm{~m}^{-1} \mathrm{~s}^{-1}=\mathrm{Sm}^{-1} \mathrm{~s}^{-1}\right]$

電気伝導率の計算值であり,未知変数 $\tau_{\mathrm{el}}$ を含む. 図 6 (a) のように，両極性拡散効果が出る領域を除いて，温度依存性はほとんどない。ただし当然ながら， $\sigma$ 自体 には温度依存性がある。これは， $\tau_{\text {el }}$ に温度依存性があ るためである.

$\sigma$ の概算值を得るには $\tau_{\mathrm{el}}=10^{-14} \mathrm{~s}$ を代入すればよいと いう記述をしばしば見かけるが, これは純粋な $\mathrm{Si}$ 結晶

いてそれほど長い值を得ることは難しい。また，同じ $\mathrm{Si}$ 結晶でも $T$ が違えば $\tau_{\mathrm{el}}$ も異なることを考慮しなけれ ばならない.

$\tau_{\mathrm{el}}$ の温度依存性は散乱機構によって異なる. 純粋物質において支配的な「音響的変形ポテンシャル散乱」, すなわち音響フォノンによる電子散乱では， $\tau_{\text {el }}$ は $T^{-1}$ に比例する。 また, 不純物散乱のうち, 「イオン化不純物散乱」では, $\tau_{\mathrm{el}}$ は $T^{3 / 2}$ に比例する. その他, 粒界散乱など多くの散乱機構が関与していると考えられるが，温度依存性がここまで顕著ではないためか，議論される ことは少ない. 複数の散乱が支配的な系の $\tau_{\text {el }}$ は, Matthiesen 則を仮定すると，

$
\frac{1}{\tau_{e l}}=\sum_{s} \frac{1}{\tau_{\mathrm{el}}(s)} \propto \frac{1}{T^{r}}
$

と計算される. 指数 $r$ は通常 -1 から $3 / 2$ の間の值をとる. $\tau_{\mathrm{el}}$ を実験的に測定するには，ARPES スペクトルのバ ンド幅を測定する方法などがある。

6. 7 . $R_{\mathrm{H}}$

Hall 係数の計算值である. 計算過程で $\tau_{\mathrm{el}}$ がキャンセ

キャリアドープ量依存性．交換相関ポテンシャル は TB-mBJ を用い，実験值とほぼ一致する $\epsilon_{\mathrm{g}}=$ $1.169 \mathrm{eV}$ を得ている. $n=10^{20} \mathrm{~cm}^{-3}$ 前後のプロッ トは, Vining らによる $p$ 型 $\mathrm{Si}_{0.8} \mathrm{Ge}_{0.2}$ 試料の $300 \mathrm{~K}$ での実験值 ${ }^{16)}$ である. なお，格子体積の温度依存性は考慮していない。

ルするので，実験值との直接比較が可能である．ホール と電子が共存する場合, Hall 測定から求めたキャリア濃度 $n_{\mathrm{H}}=e /\left|R_{\mathrm{H}}\right|$ は $n=n_{\mathrm{h}}-n_{\mathrm{e}}$ よりも， $n_{\mathrm{p}}+n_{\mathrm{e}}$ に近い值をと る傾向がある．正確には，

$
R_{H}=\frac{1}{|\mathrm{e}|}\left(\frac{\mu_{p}^{2} n_{p}-\mu_{n}^{2} n_{e}}{\left(\mu_{p} n_{p}+\mu_{e} n_{e}\right)^{2}}\right)
$

と表される. $\mu_{p}, \mu_{n}$ はホールと電子の移動度である.

よって，ホールと電子が共存する系において，実験と 対応する $\epsilon_{\mathrm{F}}$ を探すには， $n_{\mathrm{H}}$ の実験値と $n$ が一致する $\epsilon_{\mathrm{F}}$ を探すのではなく, $R_{\mathrm{H}}$ の実験值が計算値と一致する $\epsilon_{\mathrm{F}}$ を探した方が正確である。

6. 8. $\kappa_{0} / \tau_{\text {el }}\left[\mathrm{Wm}^{-1} \mathrm{~K}^{-1} \mathrm{~s}^{-1}\right]$

$\kappa_{\mathrm{el}} / \tau_{\mathrm{el}}$ の計算値であり，未知変数 $\tau_{\mathrm{el}}$ を含む. 現在の BoltzTraP から計算される值は， $\kappa_{\mathrm{el}}$ の第 1 項のみであ るため，第 2 項も含めて計算させるようにするには， BoltzTraPのソースファイルから fermiintgrals.f を開き，

thermal $(1: 3,1: 3)=\operatorname{kappa}(1: 3,1: 3)$

と書かれている行の後に，次のループを挿入してコンパ イルすればよい.

$
\mathrm{DO} \mathrm{i}=1,3
$

DO $\mathrm{j}=1,3$

DO ialp $=1,3$

thermal $(\mathrm{i}, \mathrm{j})=\operatorname{thermal}(\mathrm{i}, \mathrm{j})+\operatorname{seebeck}($ ialp, $\mathrm{i}) *$ nu(ialp, $\mathrm{j}) *$ temp

ENDDO

ENDDO

ENDDO

図 6 (b) のように, $n>10^{21} \mathrm{~cm}^{-3}$ の金属的領域では,第 2 項を入れない場合も入れた場合も， $\kappa_{\mathrm{el}} / \tau_{\mathrm{el}}$ の計算値 は一致する。しかし半導体領域では，2つの計算值は大 きく異なる. これより，「 $\kappa_{\mathrm{el}}$ の第 2 項は無視できる」と いう近似は半導体には適用できないことがわかる．熱励起の影響が出てくる高温・低 $n$ 領域においては再び両者は一致する。

7. BoltzTraP の出力データから計算できるパラメータ

7. 1. 有効 Lorenz 数 $L_{\text {eff }}$

Wiedemann-Franz 則は,

$
\kappa_{\mathrm{el}}=L \sigma T
$

で表される関係式で，固体物性理論の根幹をなす方程式 のひとつである。この $L$ としては，自由電子モデルの 理論値

$
L_{0}=\frac{\pi^{2} k_{\mathrm{B}}^{2}}{3 e^{2}} \sim 2.45 \times 10^{-8} \mathrm{~V}^{2} \mathrm{~K}^{-2}
$

が用いられることが多い。半導体でも，值こそ若干異な るものの ( $L_{0}$ の $80 \%$ 程度 $), L$ はほぼ一定となることが 知られている．また， $\epsilon_{\mathrm{g}}$ を超える熱励起があると，熱励起されたキャリアが $\epsilon_{g}$ に相当する追加エネルギーを 高温端から低温端に運ぶことになり， $\kappa_{\mathrm{el}}$ が増大するこ とが知られている。

そこで，BoltzTraP の計算結果を用い，以下のように 計算した有効 Lorenz 数 $L_{\text {eff }}$ を図 6 (c) に示す.

$
L_{\mathrm{eff}}=\frac{\kappa_{\mathrm{el}} / \tau_{e l}}{\left(\sigma / \tau_{e l}\right) T}
$

すると， $\boldsymbol{\kappa}_{\mathrm{el}}$ の第 2 項を入れた場合も入れない場合も， $n>10^{21} \mathrm{~cm}^{-3}$ の金属的極限ならば, $L_{\text {eff }}$ は $L_{0}$ と一致した しかしそれ以下の $n$ では，第 2 項を入れない場合には， $L_{\text {eff }}$ は大きな温度依存性をもって上昇した. だが第 2 項 を入れた場合は, $L_{\text {eff }}$ は $L_{0}$ を若干下回りつつもほぼ一定值を示し, 熱励起が起こる領域になって急激に増大した。 よって，だが第 2 項を含めて計算した $\kappa_{\mathrm{el}}$ は，最もよく 半導体の物性を再現していると解釈できる。 $L_{\text {eff }}$ として $L_{0}$ を使う近似も，熱励起による効果が無視できるなら ば，悪くはない

7. 2 . $Z_{\mathrm{e}} T$

$Z T$ との比較が可能な物理量として,

$
Z_{\mathrm{e}} T \equiv \frac{S^{2} \sigma T}{\kappa_{\mathrm{el}}}=\frac{S^{2}\left(\sigma / \tau_{e l}\right) T}{\left(\kappa_{\mathrm{el}} / \tau_{e l}\right)}=\frac{S^{2}}{L_{\mathrm{eff}}}
$

が挙げられる.ここでは $Z_{\mathrm{e}} T^{10)}$ と記すが， $A$-factor ${ }^{11)}$ と

図 $6 p$ 型 $\mathrm{Si}$ について計算した，（a） $\sigma / \tau_{\mathrm{el}} ， （ \mathrm{~b} ） \kappa_{\mathrm{el}} / \tau_{\mathrm{el}}$, (c) $L_{\text {eff }}=\kappa_{\mathrm{el}} / \sigma T$ (d) $Z_{\mathrm{e}} T$ のキャリアドープ量依存性. $\kappa_{0}$ は $\kappa_{\mathrm{el}}$ の第 1 項, $\kappa_{\mathrm{P}}$ は第 2 項 (Peltier 項) を示す. $L_{0}$ は自由電子モデルから計算される理論的な Lorenz 数である.
同一である. この計算過程で $\tau_{\mathrm{el}}$ が相殺できる. $Z T$ の式 が

$
Z T \equiv \frac{S^{2} \sigma T}{\kappa_{\mathrm{el}}+\kappa_{\mathrm{ph}}}=\frac{S^{2}\left(\sigma / \tau_{e l}\right) T}{\left(\kappa_{\mathrm{el}} / \tau_{e l}\right)+\left(\kappa_{\mathrm{ph}} / \tau_{e l}\right)}
$

と書けることから考えると，これは， $\kappa_{\mathrm{ph}} \ll \kappa_{\mathrm{el}}$ とみなし た極限，もしくは $\left(\kappa_{\mathrm{ph}} / \tau_{\mathrm{el}}\right) \rightarrow 0$ とみなした極限におけ る $Z T$ であることがわかる．この近似は，半導体におい てはかなり無理がある．低 $n$ 側では， $\kappa_{\mathrm{ph}}$ が無視できる という状況にはなりえず，高 $n$ の金属的極限において のみ無視できる.

$Z_{\mathrm{e}} T$ は「理論的上限値」という性質上, 実際に達成す ることは不可能な $Z T$ 值である. よって，この值がいく ら大きくても，特に意味はないことを初めに理解する必要がある．たとえこの值が「100」など異常に高い值を 示していても，そこに過剩な希望を抱いてはいけないの である.

$Z_{\mathrm{e}} T$ はむしろ，その値が低かった場合に役に立つ情報 である。電子構造と熱電特性計算が妥当であるにもかか わらず， $Z_{\mathrm{e}} T$ が低い（例：1を超えない）場合には，そ の物質ではそれ以上の $Z T$ の達成が困難だということを 意味している.

計算の妥当性の判断は難しいので一概には言えない が，熱電材料候補物質として，このようなフィルタリン グに耐えた物質のみを扱うことは，より効率的な新規熱電材料探索につながると期待できる.

図 6 (d)に $p$ 型 $\mathrm{Si}$ について計算した $Z_{\mathrm{e}} T-n$ 曲線を示す. $\kappa_{\mathrm{el}}$ の第 1 項のみを用いて計算した場合には, $Z_{\mathrm{e}} T$ が 1 を超えないという奇妙な挙動が見られる. 他の多くの半導体に関しても同様となる. $Z_{\mathrm{e}} T$ は $Z T$ の上限値である から，もしこれが正しければ， $Z T>1$ が得られるのはよ ほど特殊な電子構造をもつ物質しかないことになってし まう.

けれども， $\kappa_{\mathrm{el}}$ の第 2 項（Peltier 項 $\kappa_{\mathrm{P}}$ ）も含めて計算 すれば， $Z_{\mathrm{e}} T$ は軽々と 1 を超える. これより, 少なくと も $Z_{\mathrm{e}} T$ の観点からは，普通の半導体でも $Z T>1$ を得る ことは否定されない.

$Z_{\mathrm{e}} T=S^{2} / L_{0}$ と打いて計算しても， $\kappa_{\mathrm{el}}$ の第 2 項を考慮 した場合と近い值が得られる。ただし，熱励起が顕著な 領域においては $Z_{\mathrm{e}} T$ の値に大きな差が出る.

$n \rightarrow 0$ では， $\epsilon_{\mathrm{F}}$ がバンドギャップ中央付近の深い位置にあり，励起されたごく少数のキャリアがバンド端の 狭いエネルギー範囲に集中し，ほぼ等しい $\epsilon$ を持ってい る.このような状況では，（19）～(22) 式の $\left(\epsilon-\epsilon_{\mathrm{F}}\right)$ は $\epsilon_{\mathrm{g}} / 2$ 程度の定数とみなすことができ，積分の外にくく
りだせる. すると，

$
\begin{aligned}
& \kappa_{0}=-\kappa_{\mathrm{P}}=S^{2} \sigma T \\
& =k_{\mathrm{B}}{ }^{2} T\left(\frac{\epsilon_{g}}{2}\right)^{2} \int \sigma(\epsilon)\left(-\frac{\partial f_{\varepsilon \mathrm{F}}(\epsilon, T)}{\partial \epsilon}\right) d \varepsilon
\end{aligned}
$

と, $\kappa_{0},-\kappa_{\mathrm{P}}, S^{2} \sigma T$ が全て等しくなる. $Z_{\mathrm{e}} T=S^{2} \sigma$ $T /\left(\kappa_{0}+\kappa_{\mathrm{P}}\right)$ であることから， $\kappa_{\mathrm{P}}$ を無視した場合は $Z_{\mathrm{e}} T \rightarrow 1$ となり, 無視しなかった場合は分母が 0 となっ て, $Z_{\mathrm{e}} T \rightarrow \infty$ となることがわかる.

8. 熱電特性の $\left(\kappa_{\mathrm{ph}} / \tau_{\mathrm{el}}\right)$, 電子構造依存性

8. 1. $\left(\kappa_{\mathrm{ph}} / \tau_{\mathrm{el}}\right)$ の影響

$\left(\kappa_{\mathrm{ph}} / \tau_{\mathrm{el}}\right)$ の値を導入し， $Z T$ を計算した図を図 7 に示 す. $\left(\kappa_{\mathrm{ph}} / \tau_{\mathrm{el}}\right)$ が大きいほど, $Z T$ は $Z_{\mathrm{e}} T$ の值から大幅に 低下していく，低下幅は低 $n$ 側, 低温側ほど顕著である。例えば $\kappa_{\mathrm{ph}}=1 \mathrm{~W} \mathrm{~m}^{-1} \mathrm{~K}^{-1}, \tau_{\mathrm{el}}=10^{-14} \mathrm{~s}$ ならば $\left(\kappa_{\mathrm{ph}} / \tau_{\mathrm{el}}\right)$ は $1 \times 10^{14} \mathrm{~W} \mathrm{~m}^{-1} \mathrm{~K}^{-1} \mathrm{~s}^{-1}$ と計算される. 実際の $\mathrm{Si}$ では $\kappa_{\mathrm{ph}}$ が大きいので, $\left(\kappa_{\mathrm{ph}} / \tau_{\mathrm{el}}\right)$ はさらに大きな值となると予想される.

よって, 高い $Z T$ の達成には, できるだけ $\left(\kappa_{\mathrm{ph}} / \tau_{\mathrm{el}}\right)$ の小さい材料が良い. Slackによって提唱された, PGEC (Phonon-Glass Electron-Crystal) ${ }^{12)}$ 材料が理想的である. また $\left(\kappa_{\mathrm{ph}} / \tau_{\mathrm{el}}\right)$ の影響を受けにくい，高 $n$ 側で使うこ とも有効である. このためには，高 $n$ まで高い $Z_{\mathrm{e}} T$ を 保つ物質を選択し, 高 $n$ までのキャリアドープを実現 することが有効である. $\mathrm{Na}_{x} \mathrm{CoO}_{2}$ のような高 $S$, 高 $n$ 材

図 $7\left(\kappa_{\mathrm{ph}} / \tau_{\mathrm{el}}\right)$ として任意の値を代入して計算した, 300, $900 \mathrm{~K}$ に打ける $p$ 型 $\mathrm{Si}, n$ 型 $\mathrm{SrTiO}_{3}$ の $Z T$ のキャリアドープ量依存性. 数字は, $\left(\kappa_{\mathrm{ph}} / \tau_{\mathrm{el}}\right)$ $\left[\mathrm{Wm}^{-1} \mathrm{~K}^{-1} \mathrm{~s}^{-1}\right]$ の值を示す.
料はこのような物質であると考えられる。

## 8. 2 . DOS の勾配と $\epsilon_{\mathrm{g}$ 依存性}

電子構造の違いが熱電特性に与える影響の解明は，本講座の目的そのものである。しかし，数値に確証の持て ない状態において，各熱電材料の計算結果を数値として 紹介することは，避けなければならない。

そこで本節では，実在熱電材料の電子構造から架空の 電子構造をつくって説明する. これには $n$ 型の $\mathrm{ZnO}$ と $n$ 型の $\mathrm{FeSi}_{2}$ を利用する. 図 2 からわかるように，バン ド端の DOS の勾配は $\mathrm{ZnO}$ の伝導帯では非常にゆるやか で， $\mathrm{FeSi}_{2}$ の伝導帯では急峻である. これらの $\epsilon_{\mathrm{g}}$ を任意 の値に調整することで， $\epsilon_{\mathrm{g}}$ と電子構造が $S, \sigma / \tau_{\mathrm{el}}, L_{\mathrm{eff}}$, $Z_{\mathrm{e}} T$ の $n$ 依存性へ及ぼす影響を紹介する.

る $n$ 型 $\mathrm{ZnO}$ の $S-n$ 曲線を示す。 $\epsilon_{\mathrm{g}}$ が十分に広ければ， $S$ は低 $n$ ほど大きくなる. しかし $\epsilon_{\mathrm{g}}$ が狭くなると，低 $n$側の $S$ が低下していき, $S-n$ 曲線にピークができる。ピー クの位置は $\epsilon_{\mathrm{g}}$ が狭くなるほど高 $n$ 側に移動していき， $\epsilon_{\mathrm{g}}=0$ になると，かなり低い $S$ しか得られなくなる. 同様の傾向は，図 9 (a) の $\mathrm{FeSi}_{2}$ についても言える。物質 にもよるが, $\epsilon_{\mathrm{g}}<10 k_{\mathrm{B}} T(900 \mathrm{~K}$ で $0.75 \mathrm{eV}$ 程度)になると,熱励起の影響が $n \sim 10^{19} \mathrm{~cm}^{-3}$ 前後まで及ぶ.

熱励起は図 8 (b) - (d) に示すように $\sigma / \tau_{\text {el }}$ と $L_{\text {eff }}$ を増大させ， $Z_{\mathrm{e}} T$ を低下させる. この点から，キャリアの 熱励起は，ZT の上昇には有害であることがわかる。こ れより， $\epsilon_{\mathrm{g}} \sim 0$ の擬ギャップ金属よりも，十分に広い $\epsilon_{\mathrm{g}}$

図 $8 n$ 型 $\mathrm{ZnO}$ について, $\epsilon_{\mathrm{g}}[\mathrm{eV}]$ を変えて計算した, $900 \mathrm{~K}$ に打ける (a) $S$, (b) $\sigma / \tau_{\text {el }}$, (c) $L_{\text {eff }}, \quad$ (d) $Z_{\mathrm{e}} T$ のキャリアドープ量依存性. $\mathrm{ZnO}$ の $\epsilon_{\mathrm{g}}$ の実験値 は $3.362 \mathrm{eVである}{ }^{13)}$. $\left(\sim 10 k_{\mathrm{B}} T\right)$ が開いた半導体が熱電材料として有利と考 えられる。

熱電材料の $\epsilon_{\mathrm{g}}$ は $10 k_{\mathrm{B}} T$ が最適であるという仮説は過去に Mahanによって提唱され ${ }^{15)}$ ， “ $10 k_{\mathrm{B}} T$-rule”として 知られている．広すぎる $\epsilon_{\mathrm{g}}$ が良くない理由は， バンド 端近傍に不純物バンドが形成する確率が低くなり，キャ リアドープが難しくなるためとされる. しかし，もしバ ンド端に不純物バンドを形成できる元素が見つかるなら ば， $10 k_{\mathrm{B}} T$ 以上の $\epsilon_{\mathrm{g}}$ を持つ半導体でもいいはずだと考 えられる。

図 8 と図 9 を比較すると，バンド端における DOS 勾配が急峻な $\mathrm{FeSi}_{2}$ の方が，全体的に高い $S$ と $Z_{\mathrm{e}} T$ を出し ていることがわかる. ただし, DOS 勾配と $\sigma(\epsilon)$ の勾配は図 3 で説明した通り, 単純な比例関係にあるので はない.また， $\sigma(\epsilon)$ の勾配に $S$ が直接比例するのは, $\sigma(\epsilon)$ が $\epsilon_{F}$ 数 $k_{\mathrm{B}} T$ の範囲で 1 次関数で表される場合 のみである．高温の半導体などでは，この範囲がバンド 端の $\boldsymbol{\sigma}(\epsilon)$ の立ち上がりや谷間をまたぐため, 1 次関数 では表せない.

## 9. 実験値を用いた未知パラメータの算出

$\tau_{\mathrm{el}}$ の情報が得られると, 実験結果の解釈に大きく役立つ. たとえば，不純物ドープやナノ構造の導入が電子散乱に与える影響を， $n$ が変動する場合についても知る ことができる. また $\tau_{\mathrm{el}}$ の温度依存性から，その材料に おける支配的な電子散乱機構もわかるようになる. さま

図 $9 n$ 型 $\mathrm{FeSi}_{2}$ について， $\epsilon_{\mathrm{g}}[\mathrm{eV}]$ を変えて計算した, $900 \mathrm{~K}$ における（a） $S ，$ (b) $\sigma / \tau_{\text {el }}$ ，（c） $L_{\text {eff }} ，(d)$ $Z_{\mathrm{e}} T$ のキャリアドープ量依存性. $\mathrm{FeSi}_{2}$ の $\epsilon_{\mathrm{g}}$ の実験值は $0.72 \mathrm{eV}^{14)}$ である.

ざまな熱電材料について $\tau_{\text {el }}$ を計算・比較すれば，「長い $\tau_{\text {el }}$ が得られやすい物質」を見つけられる可能性もある.熱電特性の実験値 $\left(S_{\exp }(T), \sigma_{\exp }(T) ， \kappa_{\exp }(T)\right)$ のデー タが存在し, 正確な熱電特性計算が行われていれば,未知パラメータ $\tau_{\mathrm{el}}(T)$ の算出は可能である。これまで Wiedemann-Franz 則から概算されていた $\kappa_{\mathrm{ph}}$ も，より正確に評価できると期待できる.

## 9. 1. 解析方法

具体的な解析手段としては，各 $T$ に関して，以下の 手順を実行すればよい.

(1) $\epsilon_{\mathrm{F}}$ の選択

Hall 測定のデータがある場合は, $R_{\mathrm{H}}(T)$ の計算值と実験值が一致する $\epsilon_{\mathrm{F}}$ を探す. なければ, $S(T)$ の計算値と 実験值が一致する $\epsilon_{\mathrm{F}}$ を探す。もし計算が実際の試料の 電子構造を正しく反映していれば, どちらで判定しても，同じ $\epsilon_{\mathrm{F}}$ に行き着くと考えられる。一致する $\epsilon_{\mathrm{F}}$ が複数見 つかる場合は，単の中から，nの值が最も自然なものを 選択する。

こうして選択された一行のデータを，計算から得られ たデータセット $\left.\{\epsilon_{\mathrm{F}}, n, S_{\text {calc }},\left(\sigma / \tau_{\mathrm{el}}\right)_{\text {calc }},\left(\kappa_{\mathrm{e}} / \tau_{\mathrm{e}}\right)_{\text {calc }}\right.\}$ と解釈 し，以降これを用いて計算を進める。

（2）以下の式に従い， $\tau_{\text {el }}$ を算出する.

図 10 Vining らによるボールミル十ホットプレス $\mathrm{Si}_{0.8} \mathrm{Ge}_{0.2}$ 多結晶体 ( $n, p$ 型: “BM+ $\mathrm{HP}$ ”) と, Dismukes らによる Zone-levelling 法で作製した $\mathrm{Si}_{0.8} \mathrm{Ge}_{0.2}$ 多結晶体 ( $n$ 型: “ZL”) について, $\mathrm{Si}$ の 電子構造と $S_{\text {exp }}, \quad \sigma_{\text {exp }}, \kappa_{\exp }$ を用いて計算した (a) $\Delta \epsilon_{\mathrm{F}}$, (b) $n$, (c) $\tau_{\mathrm{el}}$, (d) $\tau_{\mathrm{el}}{ }^{-1}$, (e) $\kappa_{\mathrm{ph}}$, (f) $Z T$ の温度依存性. $\Delta \epsilon_{\mathrm{F}}$ はバンド端からの $\epsilon_{\mathrm{F}}$ のず れを示す.

$
\tau_{\mathrm{el}}=\frac{\sigma_{\exp }}{\left(\sigma / \tau_{\mathrm{el}}\right)_{\text {calc }}}
$

（3）以下の式に従い， $\kappa_{\mathrm{ph}}$ を算出する.

$
\kappa_{\mathrm{ph}}=\kappa_{\mathrm{exp}}-\tau_{\mathrm{el}}\left(\kappa_{\mathrm{el}} / \tau_{\mathrm{el}}\right)_{\mathrm{calc}}
$

## 9.2. 解析例: $\mathrm{Si_{0.8} \mathrm{Ge}_{0.2}$ 多結晶体}

このようにして計算した, 3 種類の $\mathrm{Si}_{0.8} \mathrm{Ge}_{0.2}$ 多結晶体の $\epsilon_{\mathrm{F}}, n ， \tau_{\mathrm{el}}, \kappa_{\mathrm{ph}}, Z T$ の温度依存性を図 10 に示す. Vining らの試料は, 溶融法で作成後, ボールミルで粉砕 してホットプレスで作製した $p / n$ 型の試料で，微細な粒子から構成される ${ }^{16}$. $n$ 型の試料のみ，粉砕後の酸化を 防ぐ特別な処置を行ってあるなど, 若干の製法の違いが ある. Dismukesらの $n$ 型の試料は, Zone levelling 法に よって作成した試料で, 結晶成長方向に長い結晶から なっている ${ }^{17)}$ 、データは全て, 文献 16 の Fig. 4 から取得 した. $\epsilon_{\mathrm{F}}$ は $S(T)$ との対応から調べた. 電子構造として,純粋な $\mathrm{Si}$ の電子構造を使った. よって, 以下の計算結果の解釈は参考程度であることを予め述べておく。 $\epsilon_{\mathrm{F}}$ の温度依存性を図 10 (a) に示す. $T$ の上昇に伴い バンド端へと近づく振る舞いが見られる. $n$ の温度依存性を図 10 (b) に示す. $n$ 型の試料の $n$ はゆっくりと上昇し， $p$ 型の試料はいったん低下した後再び上昇に転じ ているが， $300 \mathrm{~K}$ の Hall 測定から得た $n$ とは多少のずれ がある.

$\tau_{\text {el }}$ の温度依存性を図 10 (c) (d) に示す. いずれの試料でも $\tau_{\mathrm{el}}$ は $10^{15} \mathrm{~s}$ 台となったが, $\mathrm{ZL}$ 試料は粒界が少 ないためか, 長い $\tau_{\mathrm{el}}$ が得られた. $n$ 型の 2 試料では, $\tau_{\mathrm{el}}^{-1}$ が原点を通る直線でフィッティングできたことか ら，音響フォノンによる電子散乱が支配的であると示唆 される。 $p$ 型の試料ではこのようにはならず，計算誤差 または別の散乱機構の存在が示唆される.

$\kappa_{\mathrm{ph}}$ の温度依存性を図 10 (e) に示す. Vining の試料 の $\kappa_{\mathrm{ph}}$ がともに低いのは，結晶粒の微細化の効果である と考えられる. 高温でキャリアの熱励起が顕著になっ た状態でも， $\kappa_{\mathrm{ph}}$ の増加が低く抑えられているのは， Wiedemann-Franz 則に頼らずに $\kappa_{\mathrm{el}}=\kappa_{0}+\kappa_{\mathrm{P}}$ として計算 した効果である.

得られた $\left(\kappa_{\mathrm{ph}} / \tau_{\mathrm{el}}\right)$ を用いて計算した $Z T$ の温度依存性を図 10 (f) に示す. $Z_{\mathrm{e}} T$ は, 実験值の 5 倍程度の值となつ ている. $Z T$ の計算値は，実験值に比べて若干過大評価 されている.もし $S$ の計算值が実験値を上回っていた 場合， $Z_{\mathrm{e}} T$ は過大評価され， $Z T$ の過大評価につながる. だが同時に， $n$ の過大評価によって $\tau_{\text {el }}$ が過小評価され， $Z T$ の過小評価につながる. よって, この $Z T$ の過大評
価の原因については一概には言えないものの，正確な電子構造計算が行われていれば，このような誤差は出ない ものと考えられる。

10.まとめ

全 3 回の講座の最終回となる本稿では, Boltzmann 輸送方程式を利用して，電子構造か $\zeta S, \sigma / \tau_{\mathrm{el}}, \kappa_{\mathrm{el}} / \tau_{\mathrm{el}}$, $L_{\text {eff }}$ と, $Z T$ の理論的な最大値となる $Z_{\mathrm{e}} T$ を計算する方法について解説した。実在する物質の電子構造をべース に, $\log n$ に対する各熱電特性の変化を紹介した. $\kappa_{\mathrm{e}} / \tau_{\mathrm{el}}$ の計算における第 2 項は， $Z_{\mathrm{e}} T$ の計算に不可欠であるこ とがわかった。また，バンド端のDOS の勾配が大きい ほど高 $n$ まで $Z_{\mathrm{e}} T$ が大きくなること, 狭い $\epsilon_{\mathrm{g}}<10 k_{\mathrm{B}} T$ や大きな $\left(\kappa_{\mathrm{ph}} / \tau_{\mathrm{el}}\right)$ が， $Z_{\mathrm{e}} T(Z T)$ を低 $n$ 側から低下させ ていくことを示した。さらに，実験值 $(S, \sigma, \kappa)$ を 合わせた解析から， $\epsilon_{\mathrm{F}}, n, \tau_{\mathrm{el}}, \kappa_{\mathrm{ph}}, Z T$ が概算できるこ とを示した。

ただしこれらの計算結果には，多数の不確定性が含ま れている. 不純物元素など結晶構造の定義における不確定性， $\epsilon_{\mathrm{g}}$ の値など電子構造計算における不確定性，異方性やマルチバンド効果など熱電特性計算における不確定性， $\tau_{\mathrm{el}}$ などの未知定数の存在が挙げられる.

このデータ選択の任意性を利用して，「実験とよく合致する」計算結果を得ることは容易である。その一方，第一原理計算によって正確に熱電特性を予測することは 非常に難しい，それでも，第一原理計算によって我々の 情報量は明らかに増える。よって, 計算を正しく使えば， より合理的で有望な新規熱電材料設計や，既存材料の特性改善指針を得ることにつながると期待できる。

## 11. あとがき

熱電研究の複雑性に立ち向かうべく，第一原理計算と いうツールを身に着けて挑んだものの，結局最後は歯が 立たず，余計に複雑にしてしまったような気がするのは 熱電研究の宿命であらうか. しかし，そんな複雑性には 慣れているのが，熱電研究に携わる皆様だと思う。第一原理計算にどれだけの情報とどれだけの不確定性がある のか，きっとご理解いただけると信じている.

計算をした人が決して取ってはならないのは，「実験 なんてしなくても，計算すればすべてわかるのです」と いう態度である．確かにそれは計算を頑張るモチべー ションであり，計算の目的を簡潔に人に伝える手段でも ある。しかし，それは現実的に未達成である。

実は計算を行った人は多くの場合，計算が完璧ではな
いことを理解している。しかし，成果を強調しなければ ならない論文という場ではそれが伝わりにくい．誤差の 範井が単純なエラーバーで示せればいいのだが，そうも いかない.

実験家の側も，誤差要因の指摘によって，計算全体を 否定することは避けなければならない。第一原理計算に 完璧を求め，完璧でないことを理由にすべてを否定して は，得られるべき情報も得られなくなるからである。

計算と実験は，互いに競争するのではなく，得意分野 を生かして協力するのがべストだと考えられる. 計算は，物事の構成要素を减らして，本質を明らかにすることに 向いている. このため，傾向を調べるのは得意だが，実際の数值の予測には向いていない. 逆に実験は，多数の 相互作用が絡み合って得られる結果や，考慮されていな かった要素のあぶり出しに向いている.

理論家の使うシンプルなモデルと, 実在物質の結晶構造，そして現実の熱電特性をつなぐためのリンクとなる ように，本講座を執筆させていただいた。ここで紹介し た解析方法を，ぜひ応用・改良していただき，皆様の熱電材料の理解に役立てていただければ幸いである.

## 12. 謝辞

本記事の内容について有意義なご意見をくださりまし た，豊田工業大学の竹内恒博先生に深く感謝申し上げま す.

13. 参考文献

1) Blaha P. et al.: WIEN2k. An augmented plane wave plus local orbitals program for calculating crystal properties, Vienna University of Technology, Austria (2001).

2) Madsen G.K.H. et al.: Comp. Phys. Comm. 175, 67 (2006).

3) Perdew J.P. et al.: Phys. Rev. Lett. 77, 3865-3868

4) Aryasetiawan F. et al.: Rep. Prog. Phys. 61, 237 (1998). BIC 5) Heyd J. et al.: J. Chem. Phys. 118, 8207 (2003).

6) Tran F. et al.: Phys. Rev. Lett. 102, 226401 (2009).

7) Ashcroft N.W. et al.: Introduction to Solid State Physics, Saunders, Philadelphia (1976).

8) Mahan G.D. et al.: Proc. Natl. Acad. Sci. 93, 7436 (1996).

9) 竹内恒博，日本熱電学会誌 $8,1,17$ (2011).

10) Stiewe C. et al.: J. Appl. Phys. 97, 4, 044317 (2005).

11) Takeuchi T., Mater. Trans. 50, 2359 (2009).

12) Slack G. A., CRC Handbook of Thermoelectrics, 407

(1995).

13) Reynolds D.C. et al.: Solid State Comm. 99873 (1996).

14) Olk C.H. et al.: Phys. Rev. B 52, 1692 (1995).

15) Mahan J. E. et al.: Phys. Rev. B 54, 16965 (1996).

16) Vining C.B. et al.: J. Appl. Phys. 69, 4333 (1991).

17) Dismukes J. P. et al.: J. Appl. Phys. 10, 2899 (1976).著者連絡先

東京大学 大学院工学系研究科 物理工学専攻

押山研究室 特任研究員

桂 ゆかり：E-mailykatsura@comas.t.u-tokyo.ac.jp
the Thermoelectrics Society of Japan