1 本研究の動機および目的

熱伝導率はマイクロ・ナノスケール材料の性能を左右する重要な物理量であることから近年注目を浴びている。例えば、微細化が進む集積回路では発熱がデバイスの信頼性や効率を低下させることから、より効果的な冷却システムが求められている。そのための材料として、例えば熱伝導率が高く熱力学的にも安定なカーボンナノチューブに期待が集まっている。また、次世代のエネルギー問題に対する取り組みの一つとして、高効率の熱電変換材料を求める研究が盛んに行われているが、熱電変換材料の性能は電子の輸送係数のみならず、格子熱伝導率が極めて重要な役割を担っていることが知られている。格子熱伝導率が低いほど性能指数ZT が向上する傾向があるため、表面や界面構造の導入や非調和性の高い構造を用いた格子熱伝導率の低減とそれに伴うZT の向上が試みられている。

上述のような新奇材料をより効率的に探索するためには、実験だけではなく理論的なアプローチも欠かせない。特に密度汎関数理論(DFT) と近年の計算機の目覚ましい発展により、さまざまな輸送係数を非経験的に見積もることが可能になりつつある。ところで、電子系の輸送係数は比較的複雑な系においても見積もりが可能である一方、フォノン(格子振動) 伝導率の見積もりは、その高い計算コストから、単純な系に限られていた。そのため、例えば理論的にZT を見積もる場合においても、格子熱伝導率はコントロールパラメータとして扱われる事が多く、その結果理論と実験結果の関連づけを困難にしている。

以上のような課題を乗り越え、理論的な性能予測の信頼性を高めるため、我々は固体のフォノン伝導率を第一原理的に計算するための枠組みを開発した。我々の手法では、まず比較的小規模の第一原理計算の結果を参照として汎用的なモデル関数を構築する。次に、そのモデルと多体摂動論あるいは分子動力学法を組み合わせて格子熱伝導率の見積もりを行う。本手法の妥当性はバルクシリコンで確認し、その後、非常に複雑な構造をもつ熱電材料I 型クラスレートにおいてもその有用性を示した。

2 計算手法

固体中で各原子はおのおの平衡点を中心に振動していると仮定する。原子の変位{ui} が十分小さければ、系のポテンシャルエネルギーを以下ようにTaylor 展開で近似することが可能である。

ここで添字i には原子の番号と変位方向の情報が含まれている。各次数の結合定数{Φ} は調和および非調和の原子間力定数(IFC) であり、これらは原子種や結晶構造に依存するパラメータである。上記のモデルを我々は非調和格子モデル(ALM) と呼び、そのパラメータを第一原理的に決定する。具体的には、DFT の枠組みで様々な変位において原子に働く力をサンプルし、ALM による力がDFT の結果を最も上手く再現するようにIFC を決定する。ALM はパラメータに関して線形関数なので、パラメータの推定は線形最小自乗問題に帰着できる。我々の手法では、任意次数のパラメータを推定することが可能であり、DFT によるサンプリングをより効率よく行うために第一原理分子動力学(FPMD) を用いている。

得られた調和および非調和IFC をもとに、我々はBoltzmann 輸送方程式(BTE) もしくは非平衡分子動力学法(NEMD) を用いて格子熱伝導率の見積もりを行う。BTE では緩和時間近似を採用し、フォノン緩和時間として3 次の非調和項に起因する3 フォノン相互作用のみを考慮する。BTE では低温で重要になるフォノンの量子統計性を考慮することが可能である一方、高温で重要な4 次以上の非調和項の効果を考慮することは計算コストの観点から困難である。そのため高次非調和項の効果が重要になる場合には、NEMD に代表される分子動力学法が有用となる。NEMD 法は表面や界面系への適用も容易である一方、BTE にくらべてサイズ依存性が顕著で特に熱伝導率が高い系では注意が必要である。

3 バルクシリコンの熱伝導計算

我々はバルクシリコンに前述の手法を適用し、その妥当性を検証した。始めに、高温(500 K, 1000 K) におけるFPMD を行い変位と力のデータを蓄えた。そのデータをもとに、3 次から6 次までの非調和IFC を推定した。得られた3 次および4 次IFC の精度は先行研究と遜色がないこと、サンプリング温度への依存性は小さいことを確認した。また、フィッティングの際に5 次、6 次の補正項を考慮しない場合、3 次や4 次のパラメータがやや過小評価される事も明らかにした。得られたIFC をBTE およびNEMD と組み合わせて熱伝導率を評価した。得られた温度依存性を図1 に示す。図からわかるようにBTE による計算結果は実験値と非常に良く一致している。このことから、高温のFPMD の結果をもとに推定したパラメータの精度に関しては満足できると言える。

NEMD では熱流方向に有限サイズLz のシミュレーションセルを用意し、セルの両端に温度差を設けることで非平衡状態が達成される。長時間シミュレーションの後に系が定常状態に達すれば、Fourier の法則から熱伝導率κℓ(Lz) が得られる。ところで、平均自由行程ℓ がLz より大きなフォノンは有限サイズ効果の影響を大きく受けるため、Lz/ℓ が小さい領域では熱伝導率のサイズ依存性が顕著になる。図1におけるNEMD の結果は、Lz = 138, 207, 276 nm におけるκℓ(Lz) から外挿によってκℓ(∞) を見積もった値を示している。バルクシリコンの平均自由行程は非常に長く(数μm)、そのためLz/ℓ が十分に大きくとれないことが過小評価の一因として考えられる。

4 I 型クラスレートにおけるフォノン特性の第一原理計算

I 型クラスレート化合物は次世代の熱電材料として注目を浴びている物質であり、ユニットセルに54原子を含む複雑な構造を持っている。構造は12 面体および14 面体を構成するケージ原子とそれら多面体で構成される籠に内包されるゲスト原子から構成される。I 型クラスレートの熱伝導率は非常に低く、籠内のゲスト原子の稼働領域と格子熱伝導率に相関が見られることが知られている。ところが、ゲスト原子が熱伝導率に与える影響については合意が得られていない。BTE に基づく簡単な解析によると、熱伝導率κ は比熱C、群速度v および緩和時間τ を用いてκ = 1/3Cν2τ と書ける。I 型クラスレート化合物の低い熱伝導率は、群速度の低下によってもたらされるのか、あるいは緩和時間の減少によるものなのかを理解することは非常に重要である。

我々は、本研究の重要な課題としてI 型クラスレート化合物の格子熱伝導率を第一原理的に見積もった。Ba8Ga(16)Ge(30) (BGG) を主な計算対象とし、実験で観測される対称性を仮定してIFC の計算を行った。さらに、シリコンの場合と同様にBTE およびNEMD による熱伝導率の見積もりを行った。その結果T =100 K における格子熱伝導率としてBTE とNEMD でそれぞれ0.82 および1.4±0.2 W/mK という値を得た。現実の構造の不規則性を考慮すると、これらの値は実験値の1.1-1.9 W/mK と非常によい一致を示していると言える。また、熱伝導率の低いクラスレートではLz/ℓ を十分大きくとれるのでNEMD の精度が期待でき、さらには、ユニットセルが大きいのでBTE よりも計算効率が高くなることを確認した。

ゲスト原子の効果を考慮するため、我々はBGG 系の第一原理計算の結果を用い、Ba ゲストを除いたGa(16)Ge(30) 系のALM を構築した。このような取り扱いが許されるのは、ケージ原子間の相互作用に比べてケージ{ゲスト間相互作用が非常に弱いためである。ゲスト無しのモデルに対しBTE に基づく解析を行った結果、フォノンの緩和時間τ がBGG 系と比較して一桁以上上昇することが明らかになった(図2)。低エネルギーのフォノンはゲスト原子の影響をより強く受けていることが図からわかる。一方、ゲスト原子の有無は群速度に大きな変化もたらさないことも明らかになった。すなわち、I 型クラスレートの低い熱伝導率は、ケージ{ゲスト間相互作用によりフォノン群速度よりもむしろフォノン緩和時間が低下することによってもたらされるといえる。さらに我々はBa 原子のポテンシャル局面がdouble-well構造を持つBa8Ga(16)Sn(30) に対しても振動解析を行い、BGG との違いについて解析を行った。

5 まとめと展望

本研究では、第一原理に基づく非調和格子モデルの構築およびBoltzmann 輸送方程式あるいは非平衡分子動力学法との組み合わせによる格子熱伝導率計算の枠組みを開発した。本手法はバルクシリコンの様な単純な系からI 型クラスレート化合物のような複雑な系まで適用可能であり、二桁も値の異なる熱伝導率を定量的に計算できることを確認した。また、第一原理計算に基づく数値実験から、I 型クラスレートではケージ{ゲスト間相互作用によるフォノン緩和時間の大幅な減少が熱伝導率低下の重要な要素であるという示唆が得られた。本手法はバルク系のみでなく表面系や界面系にも適用可能であるため、例えばナノワイヤにおいて表面構造が熱伝導率に与える影響や、特にデバイス接続面において重要な界面熱抵抗などを非経験的に取り扱うことが可能になると期待できる。

図1 バルクシリコンにおける格子熱伝導率の温度依存性。計算(BTE:実線, NEMD:四角) および測定値(三角) の比較。

図2 B8Ga(16)Ge(30) (filled circle) およびGa(16)Ge(30) (open circle) におけるフォノン緩和時間の振動数依存性。T = 100 K.

格子熱伝導率は電子輸送係数と共に、(1)集積回路などに用いる極微材料の冷却効率や(2)熱電変換材料の性能係数ZT等に影響を及ぼす重要な基本物理量である。電子輸送係数の非経験的予測へ向けての試みは進展しているものの、格子熱伝導率に関しては、これまで比較的簡単な系に対する研究に限られてきた。包括的な計算を、格子熱伝導率に関しても行えることは重要である。本学位請求論文(以下、本博士論文と記す)では、格子熱伝導率の第一原理計算に関する二つの手法を提唱し、それらを格子熱伝導率が大きく異なる固体結晶(シリコンと熱電材料I型クラストレート)に適用することにより、定量的な計算が可能であることについて述べられている。また、それぞれの計算手法の特徴が具体的に示されている。

本博士論文は英語で書かれており、8つの章より構成されている。第1章から第4章までは導入部である。まず第1章では、熱伝導率の説明とこれまでの理論予測の事例に関して紹介されている。また、先行研究に対する本博士論文の位置づけについてされている。第2章は、本博士論文で用いる密度汎関数理論や原子間力の計算に関する紹介に充でられている。第3章では、原子間力にもとついてフォノンの調和項および非調和項を計算するための手法が、結晶の対称性の取り扱いを含めて説明されている。第4章は、フォノンの基本的性質についての紹介に充てられている。

博士論文の主要な結果は、第5章から第7章にまとめられている。

第5章では、フォノンの第一原理計算結果を用いて格子熱伝導率を求めるための、二つの計算手法について述べられている。一つ目は、調和フォノン系を非摂動系とし、非調和項によるフォノン散乱の影響を、ボルツマン方程式の緩和時間近似に基づき求める方法である。計算された緩和時間とフォノンの群速度から格子熱伝導率を求めることができる。二つ目は、温度差を与えたフォノン系に対して非平衡分子動力学計算を行い、熱流の平均値を求める方法である。熱流と温度差の比から格子熱伝導率を求めることができる。この二つの方法の長所、短所についての議論もされている。

第6章では、シリコン結晶に対して行った結果について述べている。まず、フォノンの調和項および非調和項を第一原理分子動力学計算から求める方法について紹介し、得られた結果が実験値をよく再現している事を確認している。また、フォノン系を3から6次までの非調和項で近似するときの精度について詳細な確認を行っている。さて、計算結果であるが、ボルツマン方程式において3次までの非調和項を用いた場合、100Kから800Kまでの広い温度範囲において、格子熱伝導率の実験値を良く再現する事が示されている。すなわち、この系に対する定量計算を可能にするという当初の目標が達成されている。次に、非平衡分子動力学計算結果について述べてある。計算結果のサイズ依存性が大きいので、データを補外して結果を得ているが、その結果は系統的に実験結果をやや下回っていることが示されている。その原因は、補外の精度が不十分であるとしている。この結果は、シリコンのような散乱長の大きな系での非平衡分子動力学計算には超大規模な計算が必要となり、あまり効率が良くないことを示唆している。

第7章では、I型クラストレートに適用した結果について述べている。非平衡分子動力学計算を適用した場合、この系が小さな散乱長を持つゆえ、サイズ依存性が弱く、しかも格子熱伝導率の実験値を良く再現できることが示されている。また、低温領域では、ボルツマン方程式を用いた計算も実験値をよく再現できることが示されている。さらに、この結晶にBa原子をゲスト原子として導入した時の格子熱伝導率の計算を行い、それが著しく小さくなる事を見出した。この変化の原因を突き止めており、その理由が、フォノン群速度よりも、フォノン緩和時間に対して、ゲスト原子による影響が大きいことにあるとしている。

最後の第8章では、得られた結果がまとめられた。

以上、各章の紹介と共に本博士論文で得られた知見を解説した、本博士論文は、格子熱伝導度に関する先駆的理論研究として意義あるものと認められる。この物理量を第一原理的に計算する手法を提案し、格子熱伝導度が大きいシリコンと、それが小さいクラストレートに対して定量的に実験を再現しただけでなくごゲスト原子の影響についての新しい知見を得ることにも成功している。今後、同手法を様々な結晶に適用してより深い理解を得ることは今後の重要な課題であるが、本博士論文はその理論研究の端緒となるものと期待される。以上の評価により、審査員全員が学位論文として十分なレベルにあり、博士(理学)の学位を授与できると判断した。

なお、本博士論文の第3章と第6章の内容は、共同研究者との共同研究を含むが、それはPhysical Review B誌で公表されている。本博士論文提出者が主体となって計算および結果の解釈を行ったものであり、本博士論文提出者の寄与が十分であると判断される。またこの件に関して、共同研究者の常行真司氏、合田業弘氏から同意承諾書が提出されている。