塩見淳一郎，「第一原理熱伝導解析の現状と今後の展望」

第一原理熱伝導解析の現状と今後の展望
志賀拓麿，塩見淳一郎（東京大学）
て材料全体の熱伝導を制御することが困難であり，熱
1. 熱電変換における熱伝導解析の重要性
伝導制御においては異なるスケールの構造を導入した
広域に分散する熱を回収し，蓄え，他のエネルギー
階層性構造制御 4) が鍵になる．このことから，できる
形態に変換する熱マネージメント技術体系において，
だけ正確に各フォノンの熱伝導率への寄与を表す熱伝
熱エネルギーから電気エネルギーに直接変換可能な熱
導率スペクトルや，その積分である累積熱伝導率を評
電変換は，熱マネージメント技術の中核を担う重要な
価することが必要になる．
要素技術である．熱電変換効率を決定する無次元性能
微視的かつ正確な熱伝導解析を実現するためには，
指数 1) は材料の熱伝導率の逆数に比例することから，
結晶格子の振動状態を正確に理解しなければならない．
材料構造の制御による低熱伝導率化は有効なアプロー
その際に重要となるのが，結晶中の原子間相互作用を
チの１つである．
記述する原子間力定数(Interatomic Force Constants,
構造制御の中でも，ナノテクノロジーを駆使した構
IFCs)である．これは原子の平衡位置からの変位(u)に対
造制御が近年注目されており，ナノ粒子を一体焼結し
するポテンシャル(V)もしくは力(F)のテイラー展開係
たバルクナノ粒子複合材 2) や，数ナノオーダーのナノ
数で定義される．
粒子を母材中に析出させるナノインクルージョン 3) な
ど様々な構造制御が提案・実施されている(図 1)．
Fi  
多様化する構造制御の背景には主要な熱キャリアで
V
1
   ij u j   ijk u j u k  
u i
2 j ,k
j
(1)
あるフォノンの輸送特性がある．熱伝導率に寄与する
フォノンの周波数，もしくは平均自由行程(mean free
ここで i, j, k は原子のインデックスであり，とはそ
path, MFP)の幅は，フェルミ準位近傍に制限される電
れぞれ調和 IFC と三次の IFC である．
子・ホールのそれに比べ幅広いため，単一構造によっ
熱物性や熱伝導率評価の基礎である調和 IFCs は，第
一原理に基づいた Frozen フォノン法 5) や密度汎関数摂
動論法 6) などの非経験法によって計算可能であり，こ
れらの方法は汎用第一原理計算パッケージに実装され
ていることから，フォノン分散関係や種々の熱物性量
を短時間かつ正確に計算することができる．一方，非
調和 IFCs や熱伝導の計算が盛んになったのは 2010 年
以降である．
第一原理に基づいた調和・非調和 IFCs 計算および熱
伝導率の定量計算のさきがけは Broido ら 7) ，また
Esfarjani ら 8) のシリコンに対する非調和格子動力学
(Anharmonic Lattice Dynamics, ALD)法計算 5) である．そ
1
図 1 ナノ粒子複合材とナノインクルージョン中
のフォノン散乱および累積熱伝導率の模式図．
の後，ALD 法はハーフホイスラー 9) や，鉛カルコゲナ
イト
電変換材料の熱伝導率はあまり高くないことから，
10,11) ，マグネシウムシリサイド 12,13) に適用され熱
SMRT は妥当な近似である．
伝導の再現や理解に役立って来た．また近年では，ビ
スマステルライド
式(3)から明らかなように，緩和時間は 2-フォノン状
14) や，スクッテルダイト 15) などの複
態密度に依存する．式(3)の第一項のデルタ関数は 1 つ
雑系へも展開されている．ALD 法の発展は近年目覚し
のフォノンの放出・吸収に関する二粒子状態密度であ
く，ShengBTE 16) や
ALAMODE 17) などの汎用パッケージ
り，低周波数領域から材料固有のカットオフ周波数ま
も次々に現れていることから，第一原理に基づいた調
で有限の値をもつ．一方，第二項のデルタ関数は 2 つ
和計算と同様に，非調和計算もより身近になりつつあ
のフォノンの生成・消滅に関する二粒子状態密度であ
る．
り，カットオフ周波数近傍からカットオフ周波数の 2
本稿では我々がこれまでに実施してきた第一原理熱
倍の高周波領域で有限の値をもつ．これらはフォノン
伝導解析法の方法論と実際の熱電材料への適用結果に
同士の散乱におけるエネルギー保存によるものである．
ついて述べるとともに，ナノ構造化バルク熱電材料の
散乱においては，結晶運動量の保存も満たされなけ
マルチスケール熱伝導解析を紹介し，第一原理に基づ
ればならないことから，緩和時間の周波数依存性は大
いた熱伝導解析の現状と展望について述べる．
雑把に言えば，結晶運動量保存とエネルギー保存(二粒
子状態密度 ) を満たす散乱位相空間 (Scattering Phase
2. 方法論
Space, SPS) 19) によって定まると言える．散乱の選択則
緩和時間近似におけるフォノン・ボルツマン輸送方程
と三次非調和 IFCs を含む散乱の行列要素 V3 は，SPS
式から，格子熱伝導率は次のように表せる．
上の散乱チャネルを閉じたり，広げたりする役割を担
っており，それによって緩和時間の値が定まると捉え

られる．したがって，正確な緩和時間や熱伝導率計算
1
 C j (q)v 2j (q) j (q)
3 q, j
(2)
のためには，調和と三次 IFCs の両方を精度よく計算し
なければならない．
ここでは系の体積，Cj(q)と vj(q)は波数ベクトル q お
以上のことから，第一原理熱伝導解析の「第一原理」
よびフォノン分岐 j で定義されるフォノンの熱容量と
たる理由は，IFCs を非経験的に求めることにある．IFCs
群速度ベクトルの大きさである．緩和時間 τ はフォノ
の計算手法は多数あるが，本稿では実空間変位法
ン同士の散乱および他のすべての散乱の影響を含んで
用いた IFCs 計算を紹介する．この方法ではスーパーセ
いるが，本稿では熱伝導に最も支配的な 3-フォノン散
ル法を用いて，対象とする材料中の各原子を平衡位置
乱のみを考慮する．三次非調和ポテンシャルに対する
から様々な非平衡位置に変位させ，原子間力の第一原
一次の摂動論に基づけば，緩和時間は次で表される
5) ．
20) を
理計算を行う．この計算を結晶の対称性を考慮した上
で，幾つか繰り返すことで，変位とそれに対する力の
1
2
 2
 j (q) 
V
2
3
 q q q
1
2 ,G
データベースが得られる．最終的に IFCs はデータベー

スを式(1)にフィッティングすることによって求める．
1, 2
(n1  n 2 ) (  1   2 )



1
　 (n1  n 2  1) (  1  2 )
 2

なお，IFCs は変位・力のデータベースがあれば求まる
(3)
ため，手動で原子を変位させてデータベースを作成す
る以外にも，ランダム変位
算
21) や第一原理分子動力学計
14) を用いることによっても得られる．
V3 は散乱の行列要素であり，三次非調和 IFCs と固有ベ
クトル，フォノンの周波数を含む．また n はボーズ・
3. 熱電変換材料の熱伝導解析
アンシュタイン分布である．なお，この式では単一モ
第一原理的に得られる IFCs を用いた ALD 計算によ
ード緩和時間(Single Mode Relaxation Time, SMRT)近似
ってバルク単結晶の熱伝導率 κ を定量的に求めること
を施してある
5) ．これは散乱における他の二つのフォ
ができるが，計算結果の中で特に重要になるのが，熱
ノンの分布(n)が平衡状態であること仮定しており，高
伝導率の MFP スペクトルもしくは，その積分である累
熱伝導材料(κ>100
Wm-1 K-1)に対しては熱伝導率を過小
積熱伝導率である．ここで累積熱伝導率の式を述べて
18) ．ただ多くの熱
おく．Dames と Chen によって提案された累積熱伝導
評価してしまうことが知られている
2
率
22) は等方的な材料を考えると，
映して，鉛テルルであれば 1-10 nm の短い MFP 領域で
急激に累積熱伝導率が増加し，シリコンであれば緩や
 c ( 0 ) 
1
3
 j (q )  0
C
j
(q)v j (q) j (q)
かに累積熱伝導率が増加する傾向が見られる．
(4)
構造制御の観点から考えたときに重要になのは，ス
q, j
ペクトルの広がりである．顕著な例はシリコンであり，
と表される． 0 はカットオフ平均自由行程である．な
熱伝導に寄与するフォノンの MFP 分布は幅広い．数
お MFP は緩和時間 τj(q)を用いて， j(q)=vj(q)τj(q)と表
kBT のエネルギー幅に制限される電子やホールとは異
される．系が十分に大きければ，式(4)を次のように表
なり，幅広いフォノン伝導を制御するためには，合金
すと便利である
23) ．
化やナノ構造化など，異なるスケールの構造制御を組
み合わせた階層性構造制御が必要であることがスペク
 c ( 0 )   
j
0
0
1
1
 d 
C j v j  j   d
3
 d 
トル分布からよく分かる．
(5)
このように，第一原理的に得られる累積熱伝導率は
材料設計指針を与えるが，ALD 計算は比較的計算コス
式(5)中の被積分関数が熱伝導率スペクトルである．
トが高く，ナノ構造化が有効な材料探索には適さない．
式(4)および(5)は MFP が 0 以下のフォノンの熱伝導
このような背景から我々はさらに累積熱伝導率のモデ
率に対する寄与のみを足し込むことを意味しているた
ル化に取り組んだ．図 3 から分かる通り，累積熱伝導
め， 0=0 の場合には κc=0，一方で 0=では，κc はバ
率は材料に強く依存するため，全体の傾向をモデル化
ルク単結晶の熱伝導率 κ に一致する(図 1)．
するのは難しいが，短 MFP と長 MFP の両極限に着目
図 2 に四種類の材料(シリコン，ハーフホイスラー
すると，増加もしくは収束する傾向には相似性が見ら
ZrCoSb，鉛テルル，マグネシウムシリサイド)の MFP
れる．したがって，両極限におけるフォノン輸送特性
熱伝導率スペクトルを示す．全体の傾向をみると，比
を次のようにモデル化した
較的短い MFP 領域でスペクトルは最大値を取り，その
13) ．
長 MFP 領域では長波長フォノンの輸送特性に着目
後減少すること，また熱伝導率が高くなるに従い，ス
すればよく，Debye モデルと Klemens モデル
24) でよく
近似できる．フォノン分散関係に対する Debye モデル
図 2 室温におけるシリコン，ハーフホイスラー
ZrCoSb，マグネシウムシリサイド，鉛テルルの熱伝
図 3 バルク単結晶の熱伝導率で規格化した累積熱伝
導率スペクトル．
導率．
ペクトル強度とピーク値を取る MFP が大きくなるこ
では，分散は線形であるから，周波数  ，フォノン分
とが分かる．一つ一つのフォノンをみれば，長い MFP
岐 j の状態密度は Dj()は Dj( )=D0,j2 で表される．一
を持つフォノンの寄与が大きいが，MFP 空間における
方，Klemens モデルから温度 T における長波長フォノ
寄与の総量を考えると，状態密度を考慮する必要があ
ンの緩和時間は j( )=A0,j-2T -1 で表される．これらとフ
る．そのため，スペクトル分布は比較的短い MFP 領域
ォノン気体モデルを用いれば，長 MFP 領域における特
で最大値を取る．次に，熱伝導率スペクトルの積分で
性長 Long は次で表される．
ある累積熱伝導率(図 3)を見ると，スペクトル分布を反
3
熱伝導率低減に与えるインパクトや，ナノ構造体の特
 Long  (c /  ) 2
(6)
徴長さに対する指針を得ることができる．
ここで c は材料パラメータであり，分岐 j の音速 v0,j と
4. 第一原理マルチスケール熱伝導解析
ボルツマン定数 kB を用いて，次で与えられる．
前節まで述べた通り，ALD 法はいまや熱電変換材料
の熱伝導解析に欠かせない計算ツールになったが，こ
c
kB
3T 3 / 2
A
j 音響
3/ 2
0, j
の手法は逆格子空間法であるため，合金化・ナノ構造
D0, j v05,/j2
(7)
化材料，また界面におけるフォノン輸送などの計算は
得意ではない．合金化による熱伝導率低減は摂動論的
一方，短い MFP 領域では，短波長・高周波数のフ
に取り扱うことが可能であるが，その適用範囲は限定
ォノンの群速度の平均(速度スケール)vS と，灰色近似を
的であることが示されている
25,26) ．本節では置換原子
施した気体分子運動論に基づいて，短 MFP における特
やナノ構造体のサイズや空間分布など，実空間情報が
性長を次のように定義した．
あらわに必要となる系の熱伝導解析手法を紹介する．
分子動力学法(Molecular Dynamics, MD)は実空間法
 Short  3 / CvS
の代表例であり，合金化によって生じる局所的な力場
(8)
や質量の変化を近似なしで取り込むことができる利点
vS は正確には累積熱伝導率が立ち上がる MFP 領域のフ
をもつ．古典 MD 法を用いる場合は，材料の熱物性・
ォノンの群速度の平均を取るべきであるが，状態密度
熱伝導特性を再現するように力場の精度を高めること
により，すべてのフォノンの群速度の平均に置き換え
が課題であるが，これは式(1)の IFCs を用いることで
ても実質的に変わらない．また式中の「3」は材料の等
解決できる．さらに鉛テルルに見られる局所格子歪み
方性によるものである．
など高次結晶格子非調和性に由来する現象も MD で取
これらの特性長とは別に，図 1 に示すような，累積
り扱えることから
27) ，
IFCs を用いた MD シミュレー
熱伝導率が立ち上がる(オンセット)MFP と，飽和し始
ションは ALD と並んで汎用性の高い熱伝導・熱物性解
める(オフセット)MFP を考える．これらを，累積熱伝
析ツールである
導率 κc がバルク単結晶の熱伝導率 κ の 10%(κc/κ=0.1)と
8,9,28) ．
しかしながら， MD 法で取り扱える計算系の大きさ
90%(κc/κ=0.9)に相当する MFP( 0.1 ， 0.9)と定義すれば，
は通常数ナノからサブマイクロ程度であり，バルクナ
短 MFP および長 MFP の特性長と 0.1/ Short~0.5，
ノ構造化材料(図 1)の熱伝導解析に対して大きな課題
 0.9/ Long =100 の関係があることが示されている
が残っている．このようなメゾスケール系における熱
13) ．
これら特性長に必要な量はバルク単結晶の熱伝導率
伝導解析で威力を発揮するのが，第一原理的に得られ
と調和量のみであるから，多くの第一原理パッケージ
た群速度と緩和時間を入力としたフォノン・ボルツマ
に実装されている調和計算を用いれば，ナノ構造化が
ン輸送方程式をモンテカルロ(Monte Carlo, MC)法を用
いて解くフォノン MC シミュレーションである
N j ( , r)
t
 v j ( )e 
N j (, r)
r

N j ( , r)  n j (, T )
 j (, T )
29) ．
(9)
ここで Nj(,r)は位置 r，角周波数 を有する分岐 j のフ
ォノンの非平衡分布，e はフォノンの伝搬方向の単位ベ
クトルである．例えば，我々はフォノン MC を数ナノ
オーダーの粒子が鉛テルル母材中に析出した系に適用
し，ナノインクルージョンのサイズや体積比に対する
熱伝導率の低減や，熱伝導率スペクトルの変化を調べ
おり，フォノンを効果的かつ選択的に散乱することが
図 4 様々な材料のオンセットとオフセット MFP．
できるナノ構造体の設計指針を得ている(図 5) (なお，
4
本計算ではワーストケース・シナリオとしてナノ粒子
なる(図 6)．
はボイドに置き換えている)29) ．
フォノン MC 法をボイドではなく実際のナノ構造
5. 結言
体に対して適用する際には，バルクのフォノン輸送情
本稿では，近年急速に発展している第一原理に基づ
報に加えて，界面におけるフォノン透過関数が必要と
いた非調和格子動力学法の概要と，実際の熱電変換材
なる．フォノン透過関数はこれまで音響整合や拡散整
料への適用例を挙げ，熱電変換における第一原理熱伝
30) が幅広く用いられているが，実空間におけ
導解析の有用性について述べた．また，非調和格子動
合モデル
る界面フォノン輸送解析法
31) を，IFCs
を用いた MD 法
力学法だけではなく，分子動力学法やモンテカルロシ
に応用することで界面フォノン透過関数も非経験的に
ミュレーションなど，実際のナノ構造化バルク熱電材
求めることが可能となっている
料の熱伝導率を解析するにあたって必要となる他の各
32) ．
以上のように，実際の熱電変換材料の熱伝導解析に
手法の有効性と適用例について述べた．
おいては，単一手法にとらわれるのではなく，異なる
本稿で紹介したマルチスケール熱伝導解析法が，多
スケールで効果的な手法(ALD や，MD，MC)を適材適
くの実験的手法や，電気伝導解析と有機的につながり，
所に組み合わせたマルチスケール熱伝導解析が重要と
高性能熱電材料や熱電デバイスの創成に貢献していく
ことを期待する．
6. 参考文献
1) Goldsmid H. J., Introduction to Thermoelectricity
(Springer, 2009).
2) Poudel B. et al.: Science 320, 634 (2008).
3) Biswas K. et al.: Nat Chem 3, 160 (2011).
4) Biswas K. et al.: Nature 489, 414 (2012).
5) Srivastava G. P., The Physics of Phonons (Taylor &
Francis, 1990).
図 5 ナノインクルージョンの体積比に対するナノ
6) Baroni S. et al.: Rev. Mod. Phys. 73, 515 (2001).
構造化鉛テルルの熱伝導率依存性．a はナノインク
7) Broido D. A. et al.: Appl. Phys. Lett. 91, 231922 (2007).
ルージョン(ナノボイド)の一辺の長さである．挿図
8) Esfarjani K. et al.: Phys. Rev. B 84, 085204 (2011).
は体積比 3.7%，a=1 nm における熱伝導率の周波数
9) Shiomi J. et al.: Phys. Rev. B 84, 104302 (2011).
スペクトル．
10) Shiga T. et al.: Phys. Rev. B 85, 155203 (2012).
11) Tian Z. et al.: Phys. Rev. B 85, 184303 (2012).
12) Li W. et al.: Phys. Rev. B 86, 174307 (2012).
13) Aketo D. et al.: Appl. Phys. Lett. 105 (2014).
14) Hellman O. et al.: Phys. Rev. B 90, 134309 (2014).
15) Li W. et al.: Phys. Rev. B 89, 184304 (2014).
16) Li W. et al.: Comput. Phys. Commun. 185, 1747 (2014).
17) Tadano T. et al.: J. Phys.: Condens. Matter 26, 225402
(2014).
18) Ward A. et al.: Phys. Rev. B 80, 125203 (2009).
19) Lindsay L. et al.: J. Phys.: Condens. Matter 20, 165209
(2008).
20) Esfarjani K. et al.: Phys. Rev. B 77, 144112 (2008).
21) Murakami T. et al.: EPL (Europhysics Letters) 102,
46002 (2013).
図 6 マルチスケール熱伝導解析のフローチャート．
22) Dames C. et al., in Thermoelectrics Handbook: Macro to
5
Nano, edited by D. Rowe (CRC Press, 2005).
(Oxford University Press, USA, 2005).
23) Yang F. et al.: Phys. Rev. B 87, 035437 (2013).
31) Chalopin Y. et al.: Appl. Phys. Lett. 101, 221903 (2012).
24) Klemens P. G., Theory of the Thermal Conductivity of
32) Murakami T. et al.: Appl. Phys. Express 7, 121801
Solids 1958), Vol. 7, in Solid State Physics,,
(2014).
p.^pp. 1.
――――――――――――――――――――
25) Larkin J. M. et al.: J. Appl. Phys. 114, 023507 (2013).
26) Shiga, T. et al.: J. J. Appl. Phys. 53, 021802 (2014).
著者連絡先
27) Shiga, T. et al.: Appl. Phys. Express 7, 041801 (2014).
塩見淳一郎：E-mail
[email protected]
28) Murakami T. et al.: Eur. Phys. Lett. 102, 46002 (2013).
TEL
03-5841-6283
29) Hori T. et al.: Appl. Phys. Lett. 104, 021915 (2014).
FAX
03-5841-8091
30) Chen G., Nanoscale Energy Transport and Conversion
URL
6
http://www.phonon.t.u-tokyo.ac.jp/

Download Report