(57)【要約】      （修正有）【課題】従来に比べて少ない計算量で電子−フォノン間相互作用を考慮しつつデバイスの電気的特性を原子論的に算出するデバイス・シミュレータを提供する。【解決手段】デバイスの構造に関するデータを格納する記憶部２と、デバイスの構造に関するデータを読み出して、デバイスの構造に応じた電子のバンド構造とフォノンの分散関係とを算出するモデル生成部１０１と、電子のバンド構造およびフォノンの分散関係に基づいて、波数依存性を有さない仮想的なフォノンである仮想フォノンと、仮想フォノンの散乱強度を決定する仮想フォノンモード決定部１０２と、モデル生成部が算出した電子のバンド構造と仮想フォノンモード決定部が決定した仮想フォノンおよび仮想フォノンの散乱強度とに基づいて、デバイスの電気的特性をシミュレートするシミュレータ部１０３と、を有する。【選択図】図２
2016024604.tif 000002
【特許請求の範囲】
【請求項１】
電子−フォノン間相互作用に影響されるデバイスの電気的特性をシミュレートするデバイス・シミュレータであって、
シミュレートされるデバイスの構造に関するデータを格納する記憶部と、
前記記憶部から前記デバイスの構造に関するデータを読み出して、前記デバイスの構造に応じた電子のバンド構造とフォノンの分散関係とを算出するモデル生成部と、
前記モデル生成部が算出した前記電子のバンド構造および前記フォノンの分散関係に基づいて、波数依存性を有さない仮想的なフォノンである仮想フォノンと、前記仮想フォノンの散乱強度を決定する仮想フォノンモード決定部と、
前記モデル生成部が算出した前記電子のバンド構造と前記仮想フォノンモード決定部が決定した前記仮想フォノンおよび前記仮想フォノンの散乱強度とに基づいて、前記デバイスの電気的特性をシミュレートするシミュレータ部と、を有するデバイス・シミュレータ。
【請求項２】
前記仮想フォノンは、複数のエネルギーレベルを有し、当該複数のエネルギーレベルの隣接するエネルギーレベル間のエネルギー差はいずれも等しい、ことを特徴とする請求項１に記載のデバイス・シミュレータ。
【請求項３】
前記仮想フォノンモード決定部は、前記電子のバンド構造および前記フォノンの分散関係に基づいて第一原理計算により求めた散乱時間のエネルギー依存性および電子移動度の電子密度依存性を用いて前記仮想フォノンおよび前記仮想フォノンの散乱強度を決定する、請求項１または２に記載のデバイス・シミュレータ。
【請求項４】
コンピュータを用いて電子−フォノン間相互作用に影響されるデバイスの電気的特性をシミュレートするデバイス・シミュレーションの方法であって、
演算部が、記憶部に格納されたシミュレートされるデバイスの構造に関するデータを読み出して、前記デバイスの構造に応じた電子のバンド構造とフォノンの分散関係とを算出するモデル生成ステップと、
前記演算部が、前記モデル生成ステップにおいて算出した前記電子のバンド構造および前記フォノンの分散関係に基づいて、波数依存性を有さない仮想的なフォノンである仮想フォノンと、前記仮想フォノンの散乱強度を決定する仮想フォノンモード決定ステップと、
前記演算部が、前記モデル生成ステップにおいて算出した前記電子のバンド構造と前記仮想フォノンモード決定ステップにおいて決定した前記仮想フォノンおよび前記仮想フォノンの散乱強度とに基づいて、前記デバイスの電気的特性をシミュレートするステップと、を有するデバイス・シミュレーションの方法。
【請求項５】
コンピュータに電子−フォノン間相互作用に影響されるデバイスの電気的特性をシミュレートさせるデバイス・シミュレーション・プログラムであって、
前記デバイス・シミュレーション・プログラムは、前記コンピュータに、
シミュレートされるデバイスの構造に応じた電子のバンド構造とフォノンの分散関係とを算出するモデル生成ステップと、
前記モデル生成ステップにおいて算出した前記電子のバンド構造および前記フォノンの分散関係に基づいて、波数依存性を有さない仮想的なフォノンである仮想フォノンと、前記仮想フォノンの散乱強度を決定する仮想フォノンモード決定ステップと、
前記モデル生成ステップにおいて算出した前記電子のバンド構造と前記仮想フォノンモード決定ステップにおいて決定した前記仮想フォノンおよび前記仮想フォノンの散乱強度とに基づいて、前記デバイスの電気的特性をシミュレートするステップと、
を実行させることを特徴とするデバイス・シミュレーション・プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、デバイス・シミュレータに関し、特に、電子−フォノン間の相互作用を考慮してデバイスの電気的特性をシミュレートするデバイス・シミュレータに関する。
【背景技術】
【０００２】
デバイス設計においてシミュレーションによりデバイスの性能を見極めることは極めて重要である。例えば、近年、プレーナ型のトランジスタにとってかわる新しいトランジスタとしてシリコン・ナノワイヤ型（以下、ＳｉＮＷ型）のトランジスタが注目されているが、そのような新しいデバイスを設計する場合、その性能を決定するデバイス特性をコンピュータ上でシミュレートすることは極めて重要である。
【０００３】
特許文献１には、ゲート電極が作るポテンシャルを、一方向に無限長であり他の方向においてゲート幅を有する面電荷を周期的に与えることにより、ナノデバイスにおけるソース・ドレイン間電流電圧特性を計算する装置および方法が記載されている。
【０００４】
非特許文献１では、原子論的な方法によりあらゆる全ての電子−フォノン間相互作用を考慮してシリコン・ナノワイヤの量子輸送特性を求める手法が述べられている。
【先行技術文献】
【特許文献】
【０００５】
【特許文献１】特開２０１４−０１７３８２号公報
【非特許文献】
【０００６】
【非特許文献１】Wenxing Zhang, et al., Atomistic modeling of electron-phonon coupling and transport properties in n-type [110] silicon nanowires（アトミスティック・モデリング・オブ・エレクトロン−フォノン・カップリング・アンド・トランスポート・プロパティーズ・イン・ｎタイプ［１１０］シリコン・ナノワイヤーズ）, Physical Review B, (2010), 82, 115319
【発明の概要】
【発明が解決しようとする課題】
【０００７】
例えばシリコン・ナノワイヤの輸送特性では、電子−フォノン間相互作用が大きな影響を与える。しかしながら、輸送特性シミュレーションにおいて電子−フォノン間相互作用を考慮する場合、その計算には膨大な時間を要する。
【０００８】
本発明は、従来に比べて少ない計算量で電子−フォノン間相互作用を考慮しつつデバイスの電気的特性を原子論的に算出することができるデバイス・シミュレータを提案する。
【課題を解決するための手段】
【０００９】
本発明の一態様は、
電子−フォノン間相互作用に影響されるデバイスの電気的特性をシミュレートするデバイス・シミュレータであって、
シミュレートされるデバイスの構造に関するデータを格納する記憶部と、
記憶部からデバイスの構造に関するデータを読み出して、デバイスの構造に応じた電子のバンド構造とフォノンの分散関係とを算出するモデル生成部と、
モデル生成部が算出した電子のバンド構造およびフォノンの分散関係に基づいて、波数依存性を有さない仮想的なフォノンである仮想フォノンと、仮想フォノンの散乱強度を決定する仮想フォノンモード決定部と、
モデル生成部が算出した電子のバンド構造と仮想フォノンモード決定部が決定した仮想フォノンおよび仮想フォノンの散乱強度とに基づいて、デバイスの電気的特性をシミュレートするシミュレータ部と、を有するデバイス・シミュレータである。
【発明の効果】
【００１０】
本発明は、少ない計算量で電子−フォノン間相互作用を考慮してデバイスの電気的特性を原子論的に算出することができる。そのため、原子論的なデバイス・シミュレーションにおいて従来に比べ計算コストが低減されたデバイス・シミュレータが供される。
【図面の簡単な説明】
【００１１】
【図１】本発明の実施形態によるデバイス・シミュレータの構成を示すブロック図
【図２】同デバイス・シミュレータの機能的構成を示すブロック図
【図３】同デバイス・シミュレータの動作の流れを示すフローチャート
【図４】シミュレートするデバイスの構造の概略を示す模式図
【図５】原子軌道を説明する図
【図６】＜１００＞シリコン・ナノワイヤの構造を説明する図
【図７】第一原理計算により求めた電子のバンド構造を示す図
【図８】フォノン振動を説明する図
【図９】第一原理計算により求めたフォノンの分散関係を示す図
【図１０Ａ】第一原理計算により求めた散乱時間の電子エネルギー依存性と、仮想フォノン（エネルギー間隔２０ｍｅＶ）より求めた散乱時間の電子エネルギー依存性とを示すグラフ
【図１０Ｂ】仮想フォノン（エネルギー間隔２０ｍｅＶ）の各モードを示す図
【図１０Ｃ】第一原理計算により求めた電子移動度の電子密度依存性と、仮想フォノン（エネルギー間隔２０ｍｅＶ）より求めた電子移動度の電子密度依存性とを示すグラフ
【図１１Ａ】第一原理計算により求めた散乱時間の電子エネルギー依存性と、仮想フォノン（エネルギー間隔１０ｍｅＶ）より求めた散乱時間の電子エネルギー依存性とを示すグラフ
【図１１Ｂ】仮想フォノン（エネルギー間隔１０ｍｅＶ）の各モードを示す図
【図１１Ｃ】第一原理計算により求めた電子移動度の電子密度依存性と、仮想フォノン（エネルギー間隔１０ｍｅＶ）より求めた電子移動度の電子密度依存性とを示すグラフ
【図１２Ａ】第一原理計算により求めた散乱時間の電子エネルギー依存性と、仮想フォノン（エネルギー間隔５ｍｅＶ）より求めた散乱時間の電子エネルギー依存性とを示すグラフ
【図１２Ｂ】仮想フォノン（エネルギー間隔５ｍｅＶ）の各モードを示す図
【図１２Ｃ】第一原理計算により求めた電子移動度の電子密度依存性と、仮想フォノン（エネルギー間隔５ｍｅＶ）より求めた電子移動度の電子密度依存性とを示すグラフ
【図１３】デバイス・シミュレーションの概要を示す模式図
【図１４】非平衡グリーン関数法の概要を示す模式図
【図１５】デバイス・シミュレーションにおける定常状態を説明する模式図
【図１６】デバイス・シミュレーションの概要を示す模式図
【図１７】グリーン関数の構成を説明する模式図
【図１８】デバイス・シミュレーションにより得られたデバイスの電流電圧特性のグラフ
【発明を実施するための形態】
【００１２】
以下、本発明の実施形態として、非弾性散乱過程である電子−フォノン散乱過程を考慮した原子論的デバイス・シミュレータを説明する。以下の説明では、シミュレーションの例として、断面形状が２．２ｎｍ×２．２ｎｍの正方形である＜１００＞シリコン・ナノワイヤの量子輸送シミュレーションを説明する。当業者であれば、以下の説明から、本発明の実施形態の構成、作用および有利な効果について、明確に理解することができるであろう。なお、本実施形態によるデバイス・シミュレータでは、言うまでも無く、上記形状のシリコン・ナノワイヤ以外のデバイスについてのシミュレーションも可能である。
【００１３】
図１は、本実施形態によるデバイス・シミュレータ１０の構成を示すブロック図である。デバイス・シミュレータ１０は、データの入力を受け付ける入力部３、様々なデータおよびプログラム等を記憶する記憶部２、記憶部２に格納されたプログラムを実行することにより入力部３が受け付けたデータおよび記憶部２に記憶されたデータを読み出して様々な演算を行い演算結果を出力する演算部１、演算結果等のデータを外部に出力する出力部４を有する。
【００１４】
入力部３は、外部の演算装置に対するデータ・インタフェース、および、ユーザからのデータ入力を受け付けるユーザ・インタフェースを含み、ユーザ・インタフェースは、マウス、キーボードといった入力装置を含む。
【００１５】
記憶部２は、各種データおよび演算部１が実行するプログラム等を格納する。記憶部２は、例えば、ＲＯＭ（Read-Only memory）およびＲＡＭ（Random Access Memory）である。記憶部２は、本発明の実施形態によるデバイス・シミュレーション方法をコンピュータで実現するためのデバイス・シミュレーション・プログラム２０４（図２）を格納している。併せて記憶部２には、後述するように、シミュレートされるデバイスの構造に関するデータ（デバイス構造データ２００（図２））、デバイスの電子のバンド構造やフォノンの分散関係に関するデータ（モデルデータ２０１（図２））、後で詳しく説明する仮想フォノンの精度を評価するためリファレンス・データ（評価用データ２０２（図２））、および、決定された仮想フォノンのエネルギーレベルや散乱強度に関するデータ（仮想フォノン・データ２０３（図２））を記憶している。本デバイス・シミュレータ１０は、これらのデータを用いて所定の構造を有するデバイスについての電子−フォノン相互作用を考慮した原子論的デバイス・シミュレーションを実行する。
【００１６】
演算部１は、所謂、中央処理装置（ＣＰＵ、Central Processing Unit）である。演算部１は、装置１０の全体の動作を制御する制御部としても動作する。なお、演算部１の構成は、ＣＰＵおよびその周辺の補助回路に限定されない。演算部１は、所定の演算処理に特化したＧＰＵ（Graphics Processing Unit）等でもよい。あるいは、演算部１は、ＡＳＩＣ（Application Specific Integrated Circuit）や、ＦＰＧＡ（Field-Programmable Gate Array）のようなプログラマブルロジックデバイスや、マイクロコントローラ、といったプロセッサとして実装可能である。演算部１は、ＣＰＵ等の上記要素を複数組み合わせて構成されてもよく、かつ演算部１の構成は上述の例に限定されるものではない。演算部１は、記憶部２に格納されたデバイス・シミュレーション・プログラム２０４（図２）を実行することにより、本発明の実施形態によるデバイス・シミュレーション方法を実行することができる。また、デバイス・シミュレーション・プログラム２０４（図２）は、フレキシブル・ディスク、光学ディスク、フラッシュメモリ等に記録されてよく、また、インターネット等のネットワークを介して送信されてよい。
【００１７】
出力部４は、演算部１が行ったデバイス・シミュレーションの結果を出力する。出力部４は、例えば、モニタディスプレイ、あるいは、プリンタである。あるいは、出力部４は、ハードディスクやフラッシュメモリといった外部記憶装置や、ネットワークを介して接続された他のコンピュータとの接続を実現するためのインタフェース機構を備えてもよい。
【００１８】
図２を参照し、演算部１が、記憶部２に格納されたデバイス・シミュレーション・プログラム２０４を実行することにより、演算部１により実現されるデバイス・シミュレータ１０の機能的構成について説明する。
【００１９】
演算部１がデバイス・シミュレーション・プログラム２０４を実行することにより実現されるモデル生成部１０１は、入力部３から入力されて記憶部２に記憶されているシミュレートするべきデバイスの構造に関するデータ（デバイス構造データ２００）を読み出し、当該データに基づいて、デバイスの電子のバンド構造、および、フォノンの分散関係等を原子論に基づく第一原理的な演算により導出する。
【００２０】
電子のバンド構造の算出は、強結合近似法や等価モデル（ＥＭモデル）といった既知の原子論的モデルに従えばよい。フォノンの分散関係の算出は、キーティング（Keating）原子価力場法といった既知の原子論的モデルに従えばよい。算出された電子のバンド構造やフォノンの分散関係は、モデルデータ２０１として記憶部２に格納される。
【００２１】
モデル生成部１０１は、さらに、電子のバンド構造やフォノンの分散関係に加え、シミュレートするべきデバイスの、散乱時間のエネルギー依存性と、電子移動度の電子密度依存性とを原子論的に計算し、これらのデータを評価用データ２０２として記憶部２に格納する。
【００２２】
次に、演算部１がデバイス・シミュレーション・プログラム２０４を実行することにより実現される仮想フォノンモード決定部１０２は、特定のエネルギーレベルを有し波数依存性を有さない仮想的なフォノン（仮想フォノン）、および、各仮想フォノンの散乱強度を、評価用データ２０２（上述の散乱時間のエネルギー依存性と電子移動度の電子密度依存性）を用いて、決定する。
【００２３】
すなわち、仮想フォノンモード決定部１０２は、評価用データ２０２として記憶されている第一原理計算による散乱時間のエネルギー依存性を用いて、変分法により、仮想フォノンの散乱強度を決定する。そして仮想フォノンモード決定部１０２は、決定された仮想フォノンの散乱強度から電子移動度の電子密度依存性を算出し、評価用データ２０２として記憶されている、先に求めた第一原理計算による電子移動度の電子密度依存性と比較し、仮想フォノンの適合性を評価する。適合性評価において所定レベル以上の高評価を得た仮想フォノンは、後のデバイス・シミュレーションで使用される。そのため、高評価を得た当該仮想フォノンのエネルギーレベルおよび散乱強度は、仮想フォノン・データ２０３として記憶部２に格納される。
【００２４】
演算部１がデバイス・シミュレーション・プログラム２０４を実行することにより実現されるシミュレータ部１０３は、仮想フォノンモード決定部１０２が決定した仮想フォノンおよびその散乱強度を用いて、デバイス・シミュレーションを行い、デバイスの特性（電流電圧特性等）を算出する。
【００２５】
次に、図３〜図１８を参照し、デバイス・シミュレータ１０の演算部１がデバイス・シミュレーション・プログラム２０４を実行したときにデバイス・シミュレータ１０が行う動作について説明する。
【００２６】
図３は、演算部１が実行する処理の流れを示すフローチャートである。
【００２７】
ステップＳ１において、演算部１は、モデル生成部１０１（図２）として動作することにより、シミュレートするべきデバイス（デバイス３００（図４））の構造に関するデータ（デバイス構造データ２００（図２））を記憶部２から読み出し、当該データに基づいて、デバイスの電子のバンド構造を計算する。なお、デバイス構造データ２００（図２）は、ユーザが入力部３を介して、シミュレータ１０に入力してもよいし、入力部３を介して接続される外部の機器から、データ通信によりシミュレータ１０に入力されてもよい。
【００２８】
図４は、シミュレートするべきデバイスの構造の例を示した模式図である。図４（ａ）は、デバイス３００の立面図であり、図４（ｂ）は、デバイス３００の側面図である。ここではデバイス３００は、断面形状が、２．２ナノメートル×２．２ナノメートルの正方形形状を有する＜１００＞シリコン・ナノワイヤである。デバイス３００のシリコン・ナノワイヤは、長手方向に３０ナノメートルの長さを有し、ソース電極４０１側１０ナノメートル長の領域３０１は、ｎ型シリコンであり、絶縁膜４０４および金属４０５を介してゲート４０２と対向する１０ナノメートル長の中央領域３０２は、真性シリコンであり、ドレイン電極４０３側１０ナノメートル長の領域３０３は、ｎ型シリコンであるとする。ここでのデバイス３００は、ゲート長１０ナノメートル、酸化膜厚１ナノメートルのゲート・オール・アラウンド（ＧＡＡ）型シリコン・ナノワイヤ・トランジスタである。
【００２９】
ここで、図５〜図７を参照し、デバイス３００の電子のバンド構造の決定の方法について説明する。
【００３０】
モデル生成部１０１（演算部１）（図２）は、強結合近似法を用いて、デバイス３００のシリコン・ナノワイヤの電子のバンド構造を算出する。モデル生成部１０１（演算部１）は、算出した電子のバンド構造を、モデルデータ２０１として記憶部２に格納する。
【００３１】
図５は、デバイス３００の原子５０１および原子５０１の原子軌道５０２の模式図である。モデル生成部１０１（演算部１）（図２）は、強結合近似法に従い、各原子５０１に属する原子軌道５０２を最小単位とし、最近接原子軌道間の相互作用までを考慮してハミルトニアンを決定する。すなわち、モデル生成部１０１（演算部１）は、強結合近似法に従い、各原子５０１に属する原子軌道５０２の原子軌道関数
【数１】
2016024604.tif 000003
を基底として、電子のハミルトニアンＨｅを行列形式で決定し、当該ハミルトニアン行列の固有値を算出することにより、電子のバンド構造を算出する。ここでの原子軌道関数φｎ，λＴＢの添え字ｎは、原子の番号（ｎ：ｎ＝１，２，３，．．．，ＮＡ）（ＮＡは原子数。）であり、添え字λは、原子軌道の種類（ｓ軌道、ｐｘ軌道、など）を表す。
【００３２】
原子の数ＮＡとは、考慮する系に含まれる原子の総数である。ここで、図６は、シリコン・ナノワイヤの構造を示す図である。図６から判るように、＜１００＞シリコン・ナノワイヤは、長手方向（図の左右方向）にＡ、Ｂ、Ｃ、Ｄの４種類の原子層が周期的に繰り返される周期性を有する。そのため、モデル生成部１０１（演算部１）（図２）がするバンド計算においては、これら４つの原子層から構成される繰り返しの基本単位（単位胞）（図６のａ）に含まれる原子５０１のみを考慮すれば足りる。断面形状２．２ｎｍ×２．２ｎｍの正方形の＜１００＞シリコン・ナノワイヤ（図４）の場合は単位胞に含まれる原子数ＮＡは、ＮＡ＝１４５でよい。
【００３３】
一方、各原子について、１０種類の軌道（ｓ軌道１種類、ｐ軌道３種類、ｄ軌道５種類、ｓ＊軌道１種類の１０種の軌道）を考慮すればシリコンの伝導帯をよく記述できることが知られている。したがって、モデル生成部１０１（演算部１）（図２）がデバイス３００の電子のバンド構造を強結合近似法により計算する場合、必要な基底関数の数ＭＴＢは、ＮＡ×１０＝１４５０となる。そのため、モデル生成部１０１（演算部１）（図２）が決定するハミルトニアン行列のサイズは、ＭＴＢ×ＭＴＢ＝１４５０×１４５０となる。モデル生成部１０１（演算部１）（図２）は、このハミルトニアン行列の固有方程式を生成して固有値問題を解くことによって、電子のバンド構造を算出する。図７は、モデル生成部１０１（演算部１）（図２）がそのようにして強結合近似法により算出した電子のバンド構造を示す図である。算出したバンド構造は、モデルデータ２０１として、記憶部２０１に格納される。
【００３４】
なお、強結合近似法に代えて、等価モデル（ＥＭモデル）を用いて、デバイス３００（図４）の電子のバンド構造を算出し、モデルデータ２０１を生成してもよい。（等価モデルの詳細については、G. Mil'nikov, N. Mori, and Y. Kamakura, "Equivalent transport models in atomistic quantum wires", Physical Review B, 85, 035317 (2012)を参照。）等価モデルでは、強結合近似法との比較において極めて少数の基底関数φｊＥＭ（ｊ：ｊ＝１，２，．．．，ＭＥＭ；ＭＥＭ＜＜ＭＴＢ）を構築することにより電子状態を記述する。そのため、電子のバンド構造の計算にかかる計算コストならびにデバイス・シミュレーションに要する計算コストの削減の観点で、等価モデルは、強結合近似法よりも有利である。図７にＭＥＭ＝４４として等価モデルに従ってモデル生成部１０１（演算部１）（図２）が算出した電子のバンド構造を示す。ＭＥＭ＝４４（＜＜ＭＴＢ＝１４５０）の等価モデルにより計算されたデバイス３００（図４）の電子のバンド構造は、とくにエネルギー範囲０＜Ｅ−Ｅｃ＜０．５ｅＶの範囲において、強結合近似法による結果をよく再現している（よく一致している）ことがわかる。このように強結合近似法に代えて等価モデルによりバンド構造を決定する場合でも、決定されたバンド構造はモデルデータ２０１として、記憶部２０１に格納される。
【００３５】
図３に戻り、ステップＳ２において、モデル生成部１０１（演算部１）は、記憶部２に格納されているデバイス構造データ２００（図２）に基づいて、シミュレートされるデバイス３００（図４）におけるフォノンの分散関係を算出する。本ステップでは、モデル生成部１０１（演算部１）は、フォノンの分散関係を、キーティング（Keating）原子価力場法を用いて算出する。モデル生成部１０１がこのようにして導出した分散関係は、モデルデータ２０１（図２）として記憶部２に格納される。ただし、本発明の実施形態では、フォノンの分散関係の導出の手法について一切の限定はなく、ここで説明するキーティング原子価力場法は、その一例に過ぎない。
【００３６】
キーティング原子価力場法においては、図８に示されるように原子５０１の間をバネ５０３で結んだモデルを用い、原子系についての古典力学的な運動方程式を解くことにより、基準モード（全ての原子５０１が同じ振動数で同位相で振動するモード）を求める。なお、図８には、原子間の結合方向に沿ったバネ５０３のみが示されているが、計算では、結合角方向に関するバネも考慮する。
【００３７】
図９は、そのようにしてモデル生成部１０１が算出したデバイス３００のフォノンの分散関係のプロットである。先に述べたようにデバイス３００は、断面形状が２．２ナノメートル×２．２ナノメートルの正方形形状の＜１００＞シリコン・ナノワイヤであり、単位胞ａ（図６）に含まれる原子の数は、原子数ＮＡ＝１４５個であり、かつ各原子が、ｘ，ｙ，ｚ方向の３つの自由度を有するため、図９に示すように、フォノンの基準モードは、フォノン波数ベクトルｑに対し、３ＮＡ＝４３５個のモード（ブランチ）を有する。（実際には、計算に際しシリコン・ナノワイヤの表面を水素原子で終端しているため、ブランチの数は４７１個になる。）
【００３８】
図３に戻り、次にモデル生成部１０１（演算部１）（図２）は、ステップＳ３において、ステップＳ１およびステップＳ２で算出した結果（モデルデータ２０１）を用いて、つまり、電子のバンド構造およびフォノンの分散関係に基づいて、デバイス３００における、散乱時間の電子エネルギー依存性、および、電子移動度の電子密度依存性を原子論に基づいて第一原理的に計算する。すなわち、本ステップにおいては、モデル生成部１０１（演算部１）（図２）は、図９に示した全てのフォノンモードのブランチを考慮して、散乱時間の電子エネルギー依存性、および、電子移動度の電子密度依存性を算出する。図１０Ａ、図１１Ａ、および、図１２Ａのプロット６０１は、モデル生成部１０１（演算部１）（図２）がそのようにして算出した散乱時間の電子エネルギー依存性のプロットであり、図１０Ｃ、図１１Ｃ、および、図１２Ｃのプロット６０５は、モデル生成部１０１（演算部１）（図２）がそのようにして算出した電子移動度の電子密度依存性のプロットである。このようにしてモデル生成部１０１（演算部１）（図２）が第一原理的に求めた散乱時間の電子エネルギー依存性、および、電子移動度の電子密度依存性は、評価用データ２０２として、記憶部２に格納される。
【００３９】
次に、ステップＳ４からステップＳ８までの一連の処理において、演算部１は、仮想フォノンモード決定部１０２（図２）として動作して、特定のエネルギーレベルを有しかつ波数依存性を有さない仮想的なフォノン（仮想フォノン）、および、各仮想フォノンの散乱強度を決定する。
【００４０】
ここで、本発明の実施形態において導入される仮想フォノンについて説明する。
【００４１】
始状態ｉにいる電子が、ブランチν、波数ベクトルｑのフォノンを吸収・放出し、終状態ｊとなる電子−フォノン間相互作用過程を、Ｔｉｊνｑａｊ†ａｉ（ｂνｑ＋ｂν，−ｑ†）（ここで、Ｔｉｊνｑは、相互作用強度であり、ａｊ†（ａｉ）は電子の生成（消滅）演算子であり、ｂν，−ｑ†（ｂνｑ）は、フォノンの生成（消滅）演算子である。）とすれば、原子論に基づいて第一原理的に電子の散乱時間（散乱確率）や、移動度や、デバイス特性を計算しようとする場合、フォノンモードに関する多くの項の和、すなわち、
【数２】
2016024604.tif 000004
という形を持つ式を計算する必要がある。デバイス３００（図４）の場合、先に述べたように４７１（＝ＮＢ）個のブランチがあり、図９に示した分散関係を精度よく表現するために、波数ベクトル方向にも十分に細かくメッシュを切る（典型的には、波数ベクトル方向に、６００点程度のメッシュを構成する）場合、上記Σνｑは、２８２，６００項程度の和となる。この計算が、デバイス・シミュレーションにおいて多大な計算コストを要求する。
【００４２】
そこで、本発明の実施形態においては、演算部１（図２）は、式（１）を、
【数３】
2016024604.tif 000005
と置き換えて計算を行う（δは、クロネッカーのデルタを表す）。こうすることで、デバイス・シミュレーションに要する計算量を劇的に減少させることができる。すなわち、ブランチの指数νおよび波数ベクトルｑに依存するフォノンｈ−ων（ｑ）（ｈ−は、プランク定数ｈを２πで除した換算プランク定数を表す。）を、等エネルギー間隔で配置される（すなわち、ｈ−ωμ＝（μ−１）ｈ−Δω（μ）（ｈ−Δωは、仮想フォノンのエネルギー間隔、μは、μ＝１，２，３，．．．，ＮＥ；ＮＥ＜＜ＮＢ））波数ベクトルに依存しない少数個の仮想的なフォノンに置き換えて計算を行う。（なお、仮想フォノンのエネルギー間隔ｈ−Δωは、必ずしも等間隔でなくともよい。）
【００４３】
仮想フォノンの散乱強度Ｔ−μ（Ｔ−は、式（２）におけるＴの上にバーを付した文字記号を表す。）は、第一原理的に決定することができないため、本発明の実施形態においては、演算部１（図２）は、ステップＳ４〜ステップＳ８の処理を行うことにより、仮想フォノンの散乱強度Ｔ−μを決定する。なお、図１０Ｂ、図１１Ｂ、および、図１２Ｂのプロット６０４ａ，ｂ，ｃは、そのような仮想フォノンの例である。以下にステップＳ４〜ステップＳ８での処理を説明する。
【００４４】
図３のステップＳ４において、仮想フォノンモード決定部１０２（演算部１）（図２）は、エネルギー間隔ｈ−Δωを２０ｍｅＶとして、等間隔に並んだ、波数に依存しない仮想フォノンのモードを設定する。ここでの、エネルギー間隔ｈ−Δω＝２０ｍｅＶは、第一原理的に求めたフォノンのエネルギー（図９）のおよそ最大レベルの値（６０ｍｅＶ）を予め定めておいた仮想フォノン個数初期値（ここでは３）で除して求めた値である。図１０Ｂは、このようにして求めた仮想フォノンのモード６０４ａをプロットした図である。
【００４５】
図３のステップＳ５において、仮想フォノンモード決定部１０２（演算部１）（図２）は、仮想フォノンの散乱強度のパラメータＴ−μに適当な初期値を与え、当該散乱強度のパラメータＴ−μの初期値に基づいて、散乱時間の電子エネルギー依存性を計算する。そして、仮想フォノンモード決定部１０２は、仮想フォノンの散乱強度のパラメータＴ−μに基づく散乱時間の電子エネルギー依存性が、ステップＳ３において第一原理的に求めた散乱時間の電子エネルギー依存性（図１０Ａ等のプロット６０１）に可能な限り一致するように、変分法を用いて仮想フォノンの散乱強度のパラメータＴ−μを調整する。図１０Ａのプロット６０２ａは、仮想フォノンモード決定部１０２がそのようにして求めた、仮想フォノンの散乱強度のパラメータＴ−μに基づく散乱時間の電子エネルギー依存性であり、同図の差し込み図中のチャート６０３ａは、そのときの仮想フォノンの散乱強度のパラメータＴ−μの値を示している。
【００４６】
次に、図３のステップＳ６において、仮想フォノンモード決定部１０２（演算部１）（図２）は、ステップＳ４で設定した仮想フォノン、および、ステップＳ５で求めた仮想フォノンの散乱強度のパラメータＴ−μに基づいて電子移動度の電子密度依存性を算出する。図１０Ｃのプロット６０６ａは、そのようにして算出された電子移動度の電子密度依存性であり、もう一方のプロット６０５は、ステップＳ３において第一原理的に求めた電子移動度の電子密度依存性である。
【００４７】
次に、図３のステップＳ７において、仮想フォノンモード決定部１０２（演算部１）（図２）は、ステップＳ６の結果を用いてステップＳ４において設定した仮想フォノンを評価する。例えば、仮想フォノンモード決定部１０２は、ステップＳ３で第一原理的に求めた電子移動度（すなわちプロット６０５（図１０Ｃ））と、ステップＳ４で設定した仮想フォノンに基づいて求めた電子移動度（すなわちプロット６０６ａ（図１０Ｃ））との差分をとり、当該差分の絶対値を電子密度軸（図１０Ｃ）に沿って数値的に積分することにより、両者の乖離の程度を定量化する。そして、仮想フォノンモード決定部１０２は、そのようにしてもとめた両者の乖離の程度を、予め定めた閾値と比較することにより、ステップＳ４で設定した仮想フォノンの適合度（ステップＳ３で求めた電子移動度との一致の程度）を評価する。仮想フォノンモード決定部１０２（演算部１）（図２）が、両者は閾値以上にかけ離れている、と判断した場合（ステップＳ７における「ＮＯ」）、処理は、ステップＳ８へ進む。仮想フォノンモード決定部１０２（図２）が、両者はかけ離れていない、と判断した場合（ステップＳ７における「ＹＥＳ」）、処理は、ステップＳ９へ進む。図１０Ａ，図１０Ｂ，図１０Ｃに示した仮想フォノンの例（エネルギー間隔：２０ｍｅＶ、仮想フォノン数３（初期値）の例）では、仮想フォノンモード決定部１０２（図２）は、両者が閾値を上回る程度に乖離していると判断するため、処理は、ステップＳ８へ進む。
【００４８】
図３に戻って、仮想フォノン決定部１０２は、ステップＳ８において、エネルギー間隔をｈ−Δωを狭く設定し直して、仮想フォノンの数を増やす。ここでは、エネルギー間隔を２０ｍｅＶから１０ｍｅＶに変更し、仮想フォノン数を３から６に増やす。
【００４９】
上記ステップＳ８から戻ってきたステップＳ４においては、ステップＳ８で定めたエネルギー間隔に従って、仮想フォノン決定部１０２が、改めて仮想フォノンを設定する。
【００５０】
つまり、２回目のステップＳ４において、仮想フォノンモード決定部１０２（演算部１）（図２）は、エネルギー間隔ｈ−Δωを１０ｍｅＶとして、等間隔に並ぶ、波数に依存しない仮想フォノンのモードを設定する。図１１Ｂは、このようにして求めた仮想フォノンのモード６０４ｂをプロットした図である。
【００５１】
図３のステップＳ５において、仮想フォノンモード決定部１０２（演算部１）（図２）は、１回目のステップＳ５と同様に、仮想フォノンの散乱強度のパラメータＴ−μに適当な初期値を与え、当該散乱強度のパラメータＴ−μの初期値に基づいて、散乱時間の電子エネルギー依存性を計算する。そして、仮想フォノンモード決定部１０２は、仮想フォノンの散乱強度のパラメータＴ−μに基づく散乱時間の電子エネルギー依存性が、ステップＳ３において第一原理的に求めた散乱時間の電子エネルギー依存性（図１１Ａ等のプロット６０１）に可能な限り一致するように、変分法を用いて仮想フォノンの散乱強度のパラメータＴ−μを調整する。図１１Ａのプロット６０２ｂは、仮想フォノンモード決定部１０２がそのようにして求めた、仮想フォノンの散乱強度のパラメータＴ−μに基づく散乱時間の電子エネルギー依存性であり、同図の差し込み図中のチャート６０３ｂは、そのときの仮想フォノンの散乱強度のパラメータＴ−μの値を示している。
【００５２】
次に、図３のステップＳ６において、仮想フォノンモード決定部１０２（演算部１）（図２）は、直前のステップＳ４で設定した仮想フォノン、および、直前のステップＳ５で求めた仮想フォノンの散乱強度のパラメータＴ−μに基づいて電子移動度の電子密度依存性を算出する。図１１Ｃのプロット６０６ｂは、そのようにして算出された電子移動度の電子密度依存性であり、もう一方のプロット６０５は、ステップＳ３において第一原理的に求めた電子移動度の電子密度依存性である。
【００５３】
次に、図３のステップＳ７において、仮想フォノンモード決定部１０２（演算部１）（図２）は、直前のステップＳ６の結果を用いて直前のステップＳ４において設定した仮想フォノンを評価する。図１１Ａ，図１１Ｂ，図１１Ｃに示した仮想フォノンの例（エネルギー間隔：１０ｍｅＶ、仮想フォノン数６の例）では、仮想フォノンモード決定部１０２は、両者が閾値を上回る程度に乖離していると判断するため、処理は、ステップＳ８へ進む。
【００５４】
図３に戻って、仮想フォノン決定部１０２はステップＳ８において、エネルギー間隔をｈ−Δωを狭く設定し直して、仮想フォノンの数を増やす。ここでは、エネルギー間隔を１０ｍｅＶから５ｍｅＶに変更し、仮想フォノン数を６から１３に増やす。
【００５５】
上記ステップＳ８から戻ってきたステップＳ４においては、直前のステップＳ８で定めたエネルギー間隔に従って、仮想フォノン決定部１０２（演算部１）（図２）が、改めて仮想フォノンを設定する。
【００５６】
３回目のステップＳ４において、仮想フォノンモード決定部１０２（演算部１）（図２）は、エネルギー間隔ｈ−Δωを５ｍｅＶとして、等間隔に並ぶ、波数に依存しない仮想フォノンのモードを設定する。図１２Ｂは、このようにして求めた仮想フォノンのモード６０４ｃをプロットした図である。
【００５７】
図３のステップＳ５において、仮想フォノンモード決定部１０２（演算部１）（図２）は、１回目および２回目のステップＳ５と同様に、仮想フォノンの散乱強度のパラメータＴ−μに適当な初期値を与え、当該散乱強度のパラメータＴ−μの初期値に基づいて、散乱時間の電子エネルギー依存性を計算する。そして、仮想フォノンモード決定部１０２は、仮想フォノンの散乱強度のパラメータＴ−μに基づく散乱時間の電子エネルギー依存性が、ステップＳ３にて第一原理的に求めた散乱時間の電子エネルギー依存性（図１２Ａ等のプロット６０１）に可能な限り一致するように、変分法を用いて仮想フォノンの散乱強度のパラメータＴ−μを調整する。図１２Ａのプロット６０２ｃは、仮想フォノンモード決定部１０２がそのようにして求めた、仮想フォノンの散乱強度のパラメータＴ−μに基づく散乱時間の電子エネルギー依存性であり、同図の差し込み図中のチャート６０３ｃは、そのときの仮想フォノンの散乱強度のパラメータＴ−μの値を示している。
【００５８】
次に、図３のステップＳ６において、仮想フォノンモード決定部１０２（演算部１）（図２）は、直前のステップＳ４で設定した仮想フォノン、および、直前のステップＳ５で求めた仮想フォノンの散乱強度のパラメータＴ−μに基づいて電子移動度の電子密度依存性を算出する。図１２Ｃのプロット６０６ｃは、そのようにして算出された電子移動度の電子密度依存性であり、もう一方のプロット６０５は、ステップＳ３において第一原理的に求めた電子移動度の電子密度依存性である。
【００５９】
次に、図３のステップＳ７において、仮想フォノンモード決定部１０２（演算部１）（図２）は、直前のステップＳ６の結果を用いて直前のステップＳ４において設定した仮想フォノンを評価する。図１２Ａ，図１２Ｂ，図１２Ｃに示した仮想フォノンの例（エネルギー間隔：５ｍｅＶ、仮想フォノン数１３の例）では、仮想フォノンモード決定部１０２（図２）は、両者が閾値を超えない程度によく一致していると判断し、処理は、ステップＳ９へ進む。
【００６０】
仮想フォノンモード決定部１０２（図２）は、図３のステップＳ９において、仮想フォノンモードおよび散乱強度のパラメータを、仮想フォノン・データ２０３として記憶部２に格納する。
【００６１】
仮想フォノン・データ２０３（図２）は、デバイス３００について求めた仮想フォノンに関するデータ（仮想フォノンのエネルギーレベルや散乱強度）を格納する。なお、仮想フォノン・データ２０３は、ソース−ドレイン間距離（デバイス３００の長手方向寸法）や、ゲートの構造がデバイス３００とは異なるデバイスについても、使用することができる。したがって、そのようなデバイスについてデバイス・シミュレーションを行う場合には、デバイス３００について求めた仮想フォノン・データ２０３を使用することができる。よって、そのような場合には、ステップＳ１からステップＳ９までの処理は、スキップしてよい。
【００６２】
次に、図３のステップＳ１０において、演算部１は、シミュレータ部１０３（図２）として動作して、デバイス３００についてのデバイス・シミュレーションを行い、デバイス３００の電流電圧特性を出力する。
【００６３】
本発明の実施形態においてシミュレータ部１０３（演算部１）（図２）は、ステップＳ７で評価用データとよく一致していると判断された、ステップＳ４で設定した仮想フォノンおよびステップＳ５で求めた仮想フォノンの散乱強度パラメータＴ−μ（仮想フォノン・データ２０３として記憶部２に格納されたデータ）を用いてデバイス３００の特性（電流電圧特性）を求める。つまり、シミュレータ部１０３は、仮想フォノンとその散乱強度を用いてデバイス３００における電子−フォノン間相互作用を考慮して、デバイス３００の電気的特性（電流電圧特性）を求める。そのため、従来の手法で電子−フォノン間相互作用を考慮してデバイスの電気的特性を計算する場合に比べ、計算コストが低減される。
【００６４】
ここで、シミュレータ部１０３が実行するデバイス・シミュレーションについて説明する。
【００６５】
一般に、定常状態のデバイス・シミュレーションでは、輸送方程式とポアソンの式を連立させて自己無撞着的に解くことにより、注目するデバイスの電流電圧特性を得ている（図１３参照）。
【００６６】
デバイス内のポテンシャル分布φ（ｘ）（ｘはデバイス内の位置を指すベクトル量。）が決定されれば、輸送方程式を解くことにより、デバイス内の電子密度分布ｎ（ｘ）（ｘはデバイス内の位置を指すベクトル量。）や電流が得られる（φ→ｎ）。しかしながらポテンシャル分布φ（ｘ）は、ポアソンの式を通して電子密度分布ｎ（ｘ）に依存する（ｎ→φ）ため、デバイス内のポテンシャル分布φ（ｘ）と電子密度分布ｎ（ｘ）とを決定するために、輸送方程式とポアソンの式を連立させて解く。
【００６７】
通例、バイアスに関する条件のもとで、輸送方程式とポアソンの式とについて自己無撞着計算を行うことにより、電子密度分布や電流値が得られ、それらを纏めることによりデバイスの電流電圧特性が得られる。
【００６８】
これまで、上記輸送方程式としては、ボルツマン輸送方程式がもっぱら用いられてきたが、図４に示したデバイス３００のような極微細デバイスでは、量子力学的な効果を無視することができないため、シミュレータ部１０３（演算部１）（図２）は、ステップＳ１０において、非平衡グリーン関数（ＮＥＧＦ）法を用いた量子効果を考慮したデバイス・シミュレーションを行う。
【００６９】
次にＮＥＧＦ法について説明する。
【００７０】
ＮＥＧＦ法では、図１４に示したように、ダイソンの式を解いてグリーン関数Ｇｒを求め、ケルディッシュの式を解いて相関関数Ｇ＜を求める。ダイソンの式は、本質的には量子力学的な運動方程式である。ダイソンの式を解いてグリーン関数Ｇｒを求めることにより、デバイス内の電子の状態密度がわかる。相関関数Ｇ＜は、デバイス内の電子密度を表しており、ケルディッシュの式により、各状態への電子の出入りをカウントすることにより、相関関数Ｇ＜が決定される。
【００７１】
デバイス・シミュレーションにおいては、シミュレータ部１０３（演算部１）が、ダイソンの式やケルディッシュの式とポアソンの式とを自己無撞着的に解き、デバイス内のポテンシャル分布および電子密度分布を算出する。
【００７２】
なお、シミュレータ部１０３がするデバイス・シミュレーションにおいては、デバイス３００（図４）に電源が接続され（図１５）、ソース電極から電子が注入され続け、ドレインからは電子が流出し続けて、結果として定常に達した状態を計算する（図１６）。デバイス内の状態は、上記したように輸送方程式を解くことにより求め、デバイス３００と電極との接続部分は境界条件として取り込む。このとき、ソース電極、ドレイン電極、および、フォノンは、図１６に示したように熱平衡状態に保たれていることを仮定する。
【００７３】
図１７は、グリーン関数Ｇｒを模式的に説明する図である。ステップＳ１で先に求めた電子のバンド構造はデバイス領域のハミルトニアンとしてグリーン関数Ｇｒの計算に用いられる。また、電子−フォノン間相互作用は、自己エネルギーΣｓｃａｔとしてグリーン関数Ｇｒに導入される。
【００７４】
シミュレータ部１０３は、仮想フォノンおよびその散乱強度のパラメータＴ−μに基づいて自己エネルギーを求める。すなわち、シミュレータ部１０３は、
【数４】
2016024604.tif 000006
を計算することにより自己エネルギーを算出する。
【００７５】
従来においては、この自己エネルギーの算出において、
【数５】
2016024604.tif 000007
という計算が必要であった。（ここで、＜ｊ｜・・・｜ｉ＞等の項は、電子状態ｉ，ｊ等の波動関数に関する空間積分を表している。）図１４に示したように、この自己エネルギーの計算は、各自己無撞着ループ（図１６参照。）内で行う必要がある。デバイス３００（図４）のような２．２ナノメートル×２．２ナノメートルのシリコン・ナノワイヤ・トランジスタの場合、この従来の方法では、あるエネルギー点におけるΣｓｃａｔの計算に一般的なパーソナル・コンピュータで１日以上の処理時間を要する。これに対し、本発明の実施形態によるデバイス・シミュレータにおいては、シミュレータ部１０３（演算部１）（図２）は、式（４）の代わりに式（３）を用いるため、モード数が劇的に減少し、デバイス・シミュレーションの高速化が実現される。デバイス３００についてのデバイス・シミュレーションの場合、式（４）を用いる従来の手法で一般的なＰＣで一日以上を要した計算が、本発明の実施形態においては、同じ一般的なＰＣを用いて１分以内で完了する。
【００７６】
図１８は、そのようにしてシミュレータ部１０３が求めたデバイス３００の電流電圧特性である。図１８の結果を出力するために、一般的なＰＣでバイアスポイント１点あたり、およそ１時間程度の時間を要した。このようにシミュレータ部１０３（図２）は、デバイス・シミュレーションにおいて、自己エネルギーΣｓｃａｔの計算を、仮想フォノンおよびその散乱強度パラメータＴ−μを用いて実行する（式（３））。そのため、従来の手法（式（４））に比べ、計算コストが劇的に低減される。
【００７７】
本発明の実施形態によるデバイス・シミュレータ１０は、シミュレーション対象のデバイスについて、仮想フォノンを設定し、設定した仮想フォノンおよびその散乱強度パラメータを用いて、電子−フォノン間相互作用（自己エネルギー）を算出する。そうすることで、従来の第一原理的な計算に比べて、劇的に計算量が低減されている。なお、デバイス・シミュレータ１０は、上で説明した、第一原理計算的な電子のバンド構造およびフォノンの分散関係（モデルデータ２０１）の算出において使用する計算アルゴリズム、散乱時間の電子エネルギー依存性および電子移動度の電子密度依存性（評価用データ２０２）の算出において使用する計算アルゴリズム、仮想フォノンの散乱強度パラメータの算出において使用する計算アルゴリズム、デバイス３００の電流電圧特性のシミュレートにおいて使用する種々の計算アルゴリズムを、デバイス・シミュレーション・プログラム２０４内に保持している。
【００７８】
さらに、本実施形態のデバイス・シミュレータ１０によれば、
１．電子状態ｉ，ｊの波動関数に関する空間積分が解析的に行える（数値計算が不要になる）。
２．フォノンモードの波数ｑ依存性がなくなるため、グリーン関数の対角ブロックのみを計算すればよくなる（波数ｑ依存性がある場合には、全ブロックの計算が必要）。
３．フォノン・エネルギーの間隔が等間隔であるため、装置内に保持しておく必要があるグリーン関数の数が低減され（外部記憶装置に保持させる必要がなくなる）。
といった効果があり、これらの効果もまた、計算の高速化に資するものである。
【符号の説明】
【００７９】
１  ・・・  演算部
２  ・・・  記憶部
３  ・・・  入力部
４  ・・・  出力部
１０  ・・・  デバイス・シミュレータ
１０１  ・・・  モデル生成部
１０２  ・・・  仮想フォノンモード決定部
１０３  ・・・  シミュレータ部
２０１  ・・・  モデルデータ
２０２  ・・・  評価用データ
２０３  ・・・  仮想フォノン・データ
【図１】
2016024604.tif 000008
【図２】
2016024604.tif 000009
【図３】
2016024604.tif 000010
【図４】
2016024604.tif 000011
【図５】
2016024604.tif 000012
【図６】
2016024604.tif 000013
【図７】
2016024604.tif 000014
【図８】
2016024604.tif 000015
【図９】
2016024604.tif 000016
【図１０Ａ】
2016024604.tif 000017
【図１０Ｂ】
2016024604.tif 000018
【図１０Ｃ】
2016024604.tif 000019
【図１１Ａ】
2016024604.tif 000020
【図１１Ｂ】
2016024604.tif 000021
【図１１Ｃ】
2016024604.tif 000022
【図１２Ａ】
2016024604.tif 000023
【図１２Ｂ】
2016024604.tif 000024
【図１２Ｃ】
2016024604.tif 000025
【図１３】
2016024604.tif 000026
【図１４】
2016024604.tif 000027
【図１５】
2016024604.tif 000028
【図１６】
2016024604.tif 000029
【図１７】
2016024604.tif 000030
【図１８】
2016024604.tif 000031