はじめに

近年，薄膜積層技術の発展に伴い厚さ数 Aから数百μmまでの薄膜製作が可能となり，主に集積回路の基板として使用される．このような薄膜デバイスの運用は発熱制御や放熱性能によって左右されるため薄膜の熱伝導伝熱解析は極めて重要である．最近，MEMS分野の熱伝導実験の結果からフーリエ法則(Fourier's Law)の適用に限界があることやエピタキシャルな(Epitaxial)接合面でも大きな熱抵抗が存在することが議論され，異種物質の多層積層によって作られたSuperlatticeによる優れた断熱物質の開発も試みられている．

現在薄膜の熱伝導においては薄膜の厚さによる熱伝導率の劣化と境界面で起こる熱抵抗とが相当注目されているが，そのメカニズム(Mechanism)に関しては未だに明確な説明はなされていないのが事実である．本研究では有限のアルゴン固体を用いた古典分子動力学法から薄膜のサイズ効果による熱伝導率の劣化と異種物質の境界面で起こる熱抵抗のメカニズムを探る解析を行った．

本研究の結果によると薄膜の場合，その厚さはプォノンの平均自由行程(Mean Free Path, MFP)を決定する支配的な因子になって，熱伝導率の劣化はMFPの減少に起因することが確認された．更に分子動力学のシミュレーション結果はバルク状態のプォノンMFPの評価を与え，厚さによる熱伝導率の依存性が定量的に計算できることを示した．そして境界面で起こる界面熱抵抗(Thermal Boundary Resistance, TBR)は巨視系の界面で起こるエネルギー反射の説明のために使われている音響インピーダンスの不一致モデル(Acoustic Impedance Mismatch Model, AIMM)と同じように界面でのエネルギーの反射によるものとして仮定し，TBRによる熱エネルギー伝達の低下が見事に見積もられる新しいモデルを開発した．

計算方法

薄膜のサイズ効果による熱伝導率の変化は図1のようにアルゴン分子をfcc <111>面として並べた違う長さを持つ固体系のMD計算から調べた．図の最下の3層と最上3層を固定分子層とし，隣接する3層は速度スケーリング(Velocity Scaling)によって温度差を与える温度制御層である．下の温度制御層は高温に維持され，上の温度制御層は低温に維持されるので熱エネルギーの流れは下から上の方向である．

熱伝達方向に垂直方向であるx-y面は周期境界条件(Periodic Boundary Condition, PBC)にしたが，この条件はまさに計算系が実在の薄膜のような形状になっていることを保証する．境界面のTBRは図1の下半部にはアルゴン分子を配置し，上半部は他の物性を仮定した分子を配置することによって境界面で起こる温度ジャンプと熱流束の低下を調べた．表1は本研究で用いたL-J(12-6)ポテンシャルのパラメータで，運動方程式の積分には誤差がΔt3に比例する速度Verlet法を用いた．

固体系の場合，気体系あるいは液体系とは違って系の内部圧力による熱伝導率がかなり影響を受けると予想されるので，予め内部圧力が0になるような分子間距離を決定する必要がある．固体アルゴンの場合，系の平均温度による内部圧力0に対応する分子間距離としてアルゴン直径との比を単位にした式(1)を提案する．

本研究で計算した系の平均温度は主に10 Kと40 Kで，これらの温度に対する分子間距離は各々1.0965σと1.1115σである．系の温度勾配は両端をそれぞれ設定された高温と低温に維持することによって作られた．系が十分安定されると与えたエネルギーと奪ったエネルギーは同じになるはずで，結局系は非平衡定常状態(Non Equilibrium Steady State, NESS)になる．

図2は熱伝導率のサイズ効果を求めるため平均温度が40 K，その大きさが68Ax59Ax55A (分子配置18x18x18に相当)である系のMD計算結果の例である．一つ区間の計算時間は200 ps (Δt=1 fs)で，温度制御は20回の計算当たりに1回施した．図の(b)から分かるように系に与えたエネルギーと系から奪ったエネルギーが同じでNESSになっていることが確認できる．

シミュレーション結果と検討

内部応力の熱伝導率への影響

自由電子(Free Electron)による熱伝導への寄与が大きい金属と異なり，アルゴンのような電気絶縁体の熱伝導は全て格子振動に起因するので，固体が圧縮と引張状態での熱伝導率はゼロ応力の場合と異なると予測される．この影響を調べるために分子間距離を変化させて圧縮，ゼロ応力，引張状態で熱伝導率を測定した．図3は系の平均温度が10K，40Kの場合に系の内部圧力変化に対する熱伝導率の挙動を示すものである．系が圧縮され，分子間距離が近付くほど熱伝導率が増加することと，その傾向は系が低温になるほど著しいことを示している．この結果は分子間距離が振動エネルギーの伝達に支配的な因子であることを意味し，従来の実験やシミュレーション結果と一致する．

サイズ効果による熱伝導率の劣化

薄膜の厚さがμmのオーダ以下になると，既に述べたようにマクロ固体系の伝熱解析に使われているフーリエ法則からの熱伝導率と比べて相当低い値になることはよく知られている．本研究では薄膜のサイズ効果による熱伝導率の劣化の様子を調べるため，熱伝達方向への分子配置を6層から54層まで6層ずつ増やしながらMD計算を行った．熱伝達方向に垂直な面の分子配置は全て18x18で，両端の温度差は系の長さによって4 Kから6 Kまで取ったが，面積の大きさと両端の温度差は系の熱伝導率に影響を与えないことを確認した．

計算は系の平均温度が10 Kと40 Kの二つケースに対して行った．図4は系の長さが100A以下のものは除いて，熱伝導率の逆数と系の長さの逆数との結果を示すものである．系が長くなるほど熱伝導率が増えて行くことを示している．図には比較のためバルク状態での固体アルゴンの実験値も一緒に載せたが，確かに系の長さが或る程度以上になると熱伝導率がバルク値と同じになることを示唆している．データは3.1節と同じように一つの計算区間は200 psで，10回の計算を施した平均値である．

図4に表した近似式，そして熱伝導率とフォノンMFPとの関係である式(1)を用いることによって薄膜の厚さによるフォノンMFPの変化が計算できる．

上式によると熱伝導率はフォノンのMFPに比例することが分かる．従ってフォノンの衝突が頻繁に起こってフォノンのMFPが短くなると熱伝導率は減少し，逆に衝突の頻度が低くなってMFPが長くなると熱伝導率は増加する．フォノンMFPは様々なフォノンの散乱(Scattering)によって決定されるが，本研究のように対象系が完全結晶(Perfect Crystal)の場合には散乱は系の境界とフォノンとの衝突，そしてフォノン同士の衝突によって起こる．熱エネルギーを運ぶ個々のフォノンの伝搬方向は気体分子のようにあらゆる方向を向いて行く．

2節で説明したように熱伝導率方向であるz軸の方向には長さが大変短いことと，それと垂直なx-y面は周期境界条件によって無限に広がっていることを考えると，系の長さである薄膜の厚さがフォノンのMFPを決める支配的な因子であることを意味する．もう一つの散乱であるフォノン同士の衝突は系内に生成されたフォノンの数による．これは系の温度が高温になるほどフォノンの数が増加するので，Tに比例してフォノンMFPは減少するはずである．それで上記の散乱過程を同時に考慮するとフォノンのMFPは次のように表現できる．

上式でlbulkは系と同じ物質のバルク状態でのフォノンのMFPを，そしてlsysは系の伝熱方向への厚さを表す．式(3)から薄膜の厚さが増加するほど右辺の第2項は0に近付いて，lpは結局バルク状態のフォノンMFPであるlbulkと同じになることが分かる．

しかし分子層が数十から数百層までの超薄膜の場合には1/lsysと比べて1/lbulkは無視できるほど小さくなってlpはほぼ薄膜の厚さlsysと同じになると予想される．即ち式(3)の関係を考慮すると1/λと1/lsysは図4のように線形関係になるはずである．式(2)と(3)から熱伝導率の逆数を求めると式

上式右辺の2番目項に含まっているlbulkは系の平均温度に依存する常数であるので，結局この項は一つの常数として表現できる．

バルク状態ではlsysは無限大として考えられるので上式右辺にある第1項は0になるので，lbulkは結局図4の近似式の傾きAと切片Bとの比(lbulk=A/B)になる．図4の近似式からバルク状態でのフォノンMFPを計算すると平均温度10 Kの場合にはlbulk @ 10 K=122 A，平均温度40 Kの場合にはlbulk @ 40 K=34 Aである．これらの長さをアルゴン分子の直径と比べると各々36倍と10倍くらいに相当する．このように求められたlbulkと薄膜の厚さ，lsysを式(4)に代入することによって厚さによる熱伝導率の変化が計算できる．

図5は上式から計算した薄膜熱伝導率と図4に示したMD計算による熱伝導率を比較したもので，互いによく一致している．

界面熱抵抗の新しいモデル(ERM)

マクロ系の接触面で起こる熱抵抗に関する研究は相当行っており，表面の加工精度を上げるなどの接触面積を増加する方法を用いて或る程度減少させることがよく知られている．しかしミクロ系の場合には熱抵抗のメカニズムに関して現在まで明確な説明が与えていない．最近の研究によると異種分子によるSupperlatticeのような接触面間の距離が分子間距離くらいで，更に境界において分子同士がピタキシャルに配置された場合にも境界面で熱抵抗が存在することを松本らが報告したが，TBRの予測まではできなかった．

松本らは参照系と計算系との熱流束比をマクロ系の境界で起こる波の反射現象と同一なものとして捕まえた．即ち参照系と計算系との熱流束比がAIMMによるエネルギー反射率(Energy Reflection Coefficient, ERC)と同じになると仮定したが，定量的に満足すべき結果は得られなかった．しかし松本らの試みはTBRの正確な予想はできなかったが，従来のモデルと比べて注目する価値はあると思われる．

現在までTBRの予測にはAcoustic Mismatch Model(AMM)とDiffuse Mismatch Model(DMM)と呼ばれるLittleとSwartzによって開発されたモデルが使われているが，実験値との比べからTBRの正確な予測ができないことがよく知られている．更にこれらモデルは大変使い難い形であるので，工学設計とか評価の目的で簡単に応用できない欠点も持っている．これに反して松本らの考え方はかなり理解しやすくて，式も簡単である長所もある．松本らのAIMMによるERC計算はAMMとDMMに比べて改善された結果を示したが，正確な予測ができなかった．そのことはマクロ系に使うため作られた式をミクロ系にそのまま適用したことに原因があると思われる．従ってマクロ系に使われるAIMMをミクロ系の条件を考慮して修正すると分子レベルでのTBRを定量的に予測できるモデルの開発が期待される．

AIMMは境界面では(1)両側に作用する力が同じ，(2)両側の変位(Displacement)が同じでなければならない仮定から出発する．巨視系に対しては上記の二つの境界条件は確かである．しかし分子レベルのスケルから考えると，一つの分子がある方向に運動すると回りの分子同士もその影響を受けて同じ方向に動くはずであるが，同じ変位で動くことは不自然である．特に境界面の場合には両側分子の質量が違うので変位の差があるはずである．しかし分子間の相互作用力は2体ポテンシャル(Two Body Potential)を使ったので同じになることは容易に分かる．従って境界条件(1)は巨視系の場合と同一で，境界条件(2)は両側の変位比は質量比或はポテンシャル比と同じであると変更されなければならない．

上記の修正された境界条件のもとでAIMM理論の適用からERCを計算すると式(6)として与えられる．

図8の(a)は系の半分がポテンシャルは同じで，質量だけ違う場合と質量は同じであるがポテンシャルだけが違う場合に式(6)から求めたERCである．図の(b)はポテンシャルと質量とも違う場合のERCの予想値を示すことである．

図8から分かるように本研究で開発されたERMによるERCの予測はMD計算結果と比べると相当正確な値を与える．一番大きいずれは予想値から約10％以内である．しかし式(6)に関して述べなければならないことは修正した境界条件から理論的に求められる理論式はαとβの割り算ではなく，掛け算になるように思われる．しかし振動理論にもとついた振幅比の計算によると，質量とポテンシャルとも違う場合の振幅比は質量だけ或はポテンシャルだけが違う場合の振幅比より小さくなるので，αとβの割り算によってERCを計算しなければならないことが確認された．

本論文は「分子動力学法による薄膜の熱伝導に関する研究」と題し，ナノ・マイクロスケール固体薄膜の膜圧方向の熱伝導特性を分子動力学法シミュレーションを用いて検討し，薄膜厚さ減少に伴う熱伝導率低下の現象と異種物質間の境界での界面熱抵抗現象のメカニズムについて議論し，界面熱抵抗に関する定量的予測の可能な現象論的モデルを提案したものであり，論文は全4章よりなっている．

第1章は，「序論」であり，本研究と関連して，新素材として期待されているスーパーラティス(super-lattice)と呼ばれる多層構造薄膜の生成と応用などについて述べるともに，薄膜の熱伝導率および固体接合面の界面熱抵抗の予測の必要性とこれらに関する巨視系と微視系に関する既存の研究，従来の界面熱抵抗のモデルなどの研究をレビューし，本論文の研究目的について述べている．

第2章は，「計算方法」であり，本研究で用いた分子動力学法シミュレーションの概要と実際に計算で用いた計算系の構成，系内部に温度勾配を生成させる方法や熱伝導率や界面熱抵抗を評価するための手法について述べている．なお，計算システムの簡単化と一般化のために，原子間ポテンシャルとしてはLennard-Jones (12,6)を採用している．

第3章は，「計算結果」であり，本手法を用いることによって，計算系の境界条件や温度制御間隔，温度勾配などの計算条件によらずに一定の熱伝導率の評価が可能であることを示して，分子動力学法シミュレーションの妥当性を示した後に，下記の2つの計算結果に基づく考察を行っている．

系の内部応力が熱伝導率に及ぼす影響について検討し，系の内部圧力によって，平均分子間間隔が変化し，ポテンシャルの非調和効果によって有効なバネ定数が変化することで熱伝導率が変化すると解釈している．また，系の温度による薄膜の熱伝導率の変化を分子動力学法シミュレーションによって求めた．高温部では熱伝導率が絶対温度に反比例する結果が得られ，フォノンのウムクラップ散乱モデルで説明がつくことがわかった．一方で，低温部分では古典分子動力学法での計算では量子力学的な比熱の減少は計算されないために温度の低下とともに一定の熱伝導率に収束する様相が計算された．また，高温部では，アルゴンの熱伝導率の実験結果ともおおよその一致を示しており，分子動力学法が熱伝導の計算に有用であることを確認している．さらに，薄膜の厚さを変えた場合の熱伝導率のシミュレーション結果から厚さ方向のサイズ効果による熱伝導率の変化を確認し，逆に，熱伝導率の計算結果からフォノンの平均自由行程を求めて整理している．また，薄膜の厚さによる熱伝導率の変化を定量的に予測できるモデルを提案している．

異種固体物質の界面における界面熱抵抗について，仮想的に質量のみが異なる物質間の界面およびLennard-Jonesポテンシャルのエネルギーのみが異なる物質間の界面について分子動力学法シミュレーションを行い，一定の熱流束がある場合にこれらの界面で温度のジャンプが現れることから界面熱抵抗を計算している．さらに，マクロな波動伝搬モデル（音響伝搬モデル）の界面での境界条件を分子スケールでの挙動に対応させて変化させることで，界面熱抵抗を定量的に求めることのできる現象論モデルを提案している．このモデルによって，任意の質量比およびポテンシャルエネルギーの比の物質間の界面熱抵抗を計算できる．また，質量比とエネルギー比が一定の条件を満たした場合には全く界面熱抵抗がなくなる場合もあることを示した．

第4章は「結論」であり，上記の研究結果をまとめたものである．

以上を要するに，本論文は分子動力学法を用いて薄膜の熱伝導率および異種固体物質接合界面における界面熱抵抗についてシミュレーションを行い，その結果に基づき，明瞭な現象論モデルを提案したものであり，マイクロ・ナノスケールの熱伝導機構に関する重要な知見を与えており，分子熱工学の発展に寄与するものと考えられる．

よって本論文は博士（工学）の学位請求論文として合格と認められる．