東京大学微細加工プラットフォーム
超微細リソグラフィー・ナノ計測拠点『マルチフィジックスシミュレーションセミナー I』

半導体デバイスシミュレーションができるまで
2017年7月26日（水）東京大学 VDEC

池野理門

内容
半導体デバイスシミュレーションの基礎
シミュレーションの目的
半導体デバイスの支配方程式
物理モデル（キャリア移動度、生成消滅項）
半導体デバイスシミュレーションの数値解析手法
様々な離散化の手法
有限体積法による支配方程式の離散化
線形化（ニュートン法）と連立方程式の行列解法
大規模疎行列の高速解法
半導体デバイスシミュレーションの展望
半導体技術を取り巻く環境の変化
半導体デバイスシミュレーションの進化
2

半導体デバイスシミュレーションの基礎

3

シミュレーションとはモデリングと数値解析に基づく物理現象の予測・分析
現象の原理の理解と実例への数値解析手法の適用
対象
固体力学：構造、応力、変形
流体力学：流動、圧力、粘性、圧縮・非圧縮
熱：伝熱
光学：波動、光線
電気、電磁気支配方程式
ラプラス方程式（ポアソン方程式）
波動方程式
マックスウェル方程式4

半導体デバイスシミュレーションの目的行うこととその効果
半導体内部の電気伝導現象に関する支配方程式を数値計算手法によって解く
半導体内部の電位分布、キャリア分布を再現し、デバイスの電気的特性を製造前に理解するデバイス製造を行わずにデバイスの特性を予測する
デバイス構造や製造プロセスの設計、改善を効率化する

FinFET の構造と不純物濃度、キャリア密度、電流電圧特性
https://www.silvaco.co.jp/products/vwf/atlas/device3d/device3d_br.html

5半導体デバイスシミュレーションのための物理
半導体の状態を記述する物理量
ポテンシャル（電位）分布未知数とする（これらを求める）
キャリア分布（電子密度 、正孔密度 ）固定電荷（不純物）分布（ドナー濃度 、アクセプタ濃度 ）キャリアの移動や生成・消滅の様子（電流）pn接合（熱平衡）
ポテンシャルp型

空乏層 n型

pn接合（正バイアス）
p型n型
電子

電子

キャリア生成・消滅
直接
間接 キャリア生成
再結合 再結合 （熱、光等）

フェルミレベル

正孔==正孔(())≅≅

キャリアの移動：ドリフト（電位勾配）と拡散（密度勾配）の合計
6⋅=からのズレを補正する方向

半導体デバイスシミュレーションの支配方程式[グラディエント / 勾配][ダイバージェンス/発散・湧き出し]⋅= div =+!

#+⋅ -⋅
00+−−&ドリフト⋅ 92 − /)発散= −/ ⋅ ( −+⋅ 12 ⋅⋅ 15 ⋅⋅⋅, ) = grad) =) ) ),,! $ &キャリア移動の力学
92 = −/ ⋅ 12 ⋅ ⋅+/⋅⋅ 6 = 7 （Maxwell 方程式）

6 = -88=−7=/⋅( −+$

電位と電荷の関係%−+増減−2⋅0拡散生成消滅項)= GR
7（電子電流）= −/ ⋅ GR= GR5⋅

2⋅[ポアソン方程式][電子電流連続式][正孔電流連続式]

半導体デバイスシミュレーションと物理モデル（1）キャリア移動度 12 , 12
電界に対するキャリアの移動のし易さを表すパラメータ
移動度はキャリアに対する種々の「散乱」要因から決まる
ある「散乱」 : が緩和時間 ;: で起こる時、その散乱によって決まる移動度は 1=/;:/= （= はキャリアの有効質量）

キャリア速度
散乱

キャリア速度
散乱

平均速度

加速@=/0/=⋅A

時間@B=/;:/=⋅A
移動度;:

時間

移動度と拡散係数の間には「アインシュタインの関係式」が成り立っている。（熱平衡状態から導出）/12 =>?
82

半導体中の散乱要因と移動度モデル
フォノン散乱熱振動による結晶格子の周期性の乱れ
不純物散乱（イオン化、中性）
固定電荷との電気的干渉、結晶格子点外の原子との干渉
表面散乱
MOS界面などでの縦電界による界面との衝突
キャリア間散乱キャリア密度が高い場所での互いの干渉
速度飽和
高電界中でキャリア速度が電界に比例しなくなる際の補正
複数の散乱要因が存在する場合の移動度
原理的には「マーティセンの法則」（Matthiessen’s rule）に従うが、実用的には様々な複合モデルが提案されている
1/1Total = 1/1A+ 1/1B+ 1/1C+⋯
9

半導体デバイスシミュレーションと物理モデル（2）キャリアの生成消滅項
熱や光によるキャリアの発生と再結合による消滅の差
基本的には熱平衡状態からのズレを戻す方向に働く
Shockley-Read-Hall（SRH）再結合
バンドギャップ中のトラップ準位を介した間接遷移−+ ;2 ( +L MNO =;5 +Auger 再結合
直接遷移による再結合
L P=(−)(Q2 ⋅)+ Q5 ⋅ )インパクト・イオン化（Impact Ionization）
高エネルギーのキャリアの衝突による電子正孔対の励起
R SS =

1T ⋅ 9 + T5 ⋅ 95/ 2 2

モデル式中に電流値を含むので、定式化が複雑になる
10

半導体デバイスシミュレーションの数値解析手法

11

シミュレーション（数値解析）実現までの流れ
支配方程式（偏微分方程式）

空間メッシュ空間の離散化• 差分法• 有限要素法• 有限体積法：

時間の離散化• 陽解法• 陰解法

連立方程式（非線形）
線形化（Newton法・Tailor展開）
行列解法

連立一次方程式（大規模疎行列）
12

時間ステップさまざまな離散化手法（1）
有限差分法（Finite Differential Method; FDM） による離散化
Y(!)
U–1Uℎ

U+1

(U–1, ]+1)
U+2

BY
Y [\ − Y
1階微分≈
B! !=!Uℎℎ

Y [ − 2Y [\ + Y
2階微分 B Y≈
B! !=!Uℎ(U–1, ]–1)

(U+1, ]+1)

(U, ])(U+1, ]–1)

• 各格子点に関する連立（一次）方程式を作る• 基本的には構造格子を対象にするので、複雑な形状の精度に課題

13さまざまな離散化手法（2）
有限要素法（Finite Element Method; FEM）による離散化
接点要素• 要素内部の状態を接点の未知数値を変数とする補間関数で表す• 各要素に対する方程式の係数を組み合わせて全体の方程式を作る• 要素形状が任意なので、複雑な形状にも対応できる

14

さまざまな離散化手法（3）有限体積法（Finite Volume Method; FVM）による離散化
コントロールボリューム接点• 支配方程式をコントロールボリューム内で積分し、それを離散化する• 有限要素法と比較して簡便な定式化ができる• コントロールボリュームの流れの出入りを定式化する• 要素形状が任意なので、複雑な形状にも対応できる

15

離散化手法の比較
支配方程式（偏微分方程式）
有限差分法
差分近似

連立一次方程式

有限要素法

有限体積法

要素内部の近似・積分

コントロールボリューム積分

重み付き残差式

支配方程式（積分形）

離散化近似＋線形化

離散化近似＋線形化

連立一次方程式

連立一次方程式• 構造化格子• 物理量は非保存• 物理量は保存• 物理量は非保存• 複雑な形状に対応• 複雑な形状に対応• 複雑な形状に課題• 計算時間がかかる• 高次の精度に課題

16

有限体積法による定式化：微分形 積分形（1）ガウスの（発散）定理に基づく支配方程式の書き換え^`⋅ _ dv = ^ (_ ⋅ a) dsbk: 積分体積:: 積分表面

意味：内部での流れの増減は、表面での出入りを見ればわかる
ポアソン方程式の書き換え⋅ -⋅体積積分発散定理^⋅ -⋅^-⋅mM

d: 法線ベクトルf g= −/ ⋅ ( −j!

+−)f g + hg
if
⋅ hg
=f g +ig
2階微分

B@ = −/ ^−+−B@

⋅ d Bl = −/ ^−+−

B@mm

17

１階微分＋積分有限体積法による定式化：微分形 積分形（2）ガウスの（発散）定理に基づく基本方程式の書き換え（つづき）
電子電流連続式の書き換え
0体積積分
^m0+

⋅ 12 ⋅

B@ + ^m⋅⋅ 12 ⋅−
⋅2⋅−= / ⋅ GR
2⋅2階微分

B@ = ^ GR B@m

発散定理
^m0

B@ + ^M12 ⋅⋅−2⋅⋅ d Bl = ^ GR B@

１階微分＋積分⋅ d Bl = ^ GR B@

１階微分＋積分m

正孔電流連続式の書き換え（結果のみ）^m0

B@ − ^M15 ⋅⋅+185⋅m有限体積法による定式化：離散化（1）ポアソン方程式 （2次元構造格子、均等メッシュの例）
積分形
^M

-⋅

⋅ d Bl = −/ ^m

離散形-n−j!

oj$ + -+-

2−p−oj!

oj$

j! + -j$= −/ o − o +−+−

q−o−j!
j$⋅ j! ⋅ j$q− 2 A + A#B@

( o,整理後An+p+ A#2+19o + /ko−o+o, o)( n,−n, n)=0有限体積法による定式化：離散化（2）電子電流連続式
積分形^m0B@ + ^12 ⋅M⋅−2⋅⋅ d Bl = ^ GR B@m

とりあえず時間変化は考えない（定常解）

離散形

s12

n+2n+o整理後
vw 2−2⋅n−on−

j!

j$ + 12p+o

⋅

p−o

j$
2j!
2+ o
2− oq+ o
q− o+12⋅
j! + 12
⋅
j!t
2
j$
2
j$
n− o
p− o
2− oq− oj$ + 2
j$ + 2
j! + 2
j!
j!
j!
j$
j$
= GR ⋅ Δ! ⋅ Δ$

2−o+ w2 p + o
p− o2+w#2 2 + o
2 − o + w#q+ o222n+ o +p+ o + #2+ o + #q+o

20q−oo x− k ⋅ GR=0ここまでの定式化結果のおさらい（その1）離散化まで済んだ現在の問題：
未知数 3 、式数 3 の連立方程式（ ：格子点数）
各式は一般に非線形（2次以上の項、指数、分数を含む）
各式は中心および近隣の格子点上の未知数( , , )の関数

z2 = ( 2 ,2, 2)

zp = ( p ,zq = ( q ,zo = ( o ,o, o)

点~点の数分点•

ozn = ( n ,o)n{n| (zo , zn , zp , zq , z2 )=0o)}o{n (zo , zn , zp , zq , z2 )=0点p, p)q, q))5o q (zo , zn , zp , zq , z2 )=0計3 本（ =格子点数）の方程式が得られる

n, n)

21

半導体特有の離散化テクニック
Sharfetter-Gummel 法
キャリア密度は指数的に変化するため、線形の補完（差分）による平均や勾配は精度が悪くなる可能性がある
電流連続式の離散化に指数関数を導入して改善を図る
電子電流92= 12 ⋅−/

B
⋅
−
B!

B2⋅
B!

92−/ ⋅ 2
12 Bn− o@≡
⋅
=‚⋅
j!
2 B!
[とりあえず定数とする]
Q€ ≡

線形差分
92
= 12−/

n+

92o

⋅

n−j!

o−

2n−j!

o

中点（積分表面）では電界は線形補間で良いが電子密度は不正確•
Bから Q€ = @ ⋅ −を解いてSharfetter-Gummel 法B!
…( †ƒ⋅‡) −92−1n− oo⋅
n= 12⋅! = o + ƒ⋅„( n − o)…( †‡) − 1j!
−/−122

【参考】 Sharfetter-Gummel 法の導出! = Q(!) ⋅と仮定する （指数項と変動分に分ける）
B
これを Q€ = @ ⋅ −に代入すると
B!
Bƒ−Q ! ⋅ ƒ + Q(!) ⋅ @ ⋅ ƒQ€ = @ ⋅ Q !
から •B!
B− Q ! = Q€ ⋅ ƒ が得られる
B!
Q€Q€
ƒ
! = + Q\ ⋅+ Q\ すなわちこれを解いて Q ! = ⋅@@
境界条件

以上より

9ƒ

Q€
+ Q\o =@Q€
+ Q\ ⋅n =@! =o+ƒ⋅„から、定数は−1( n−ƒ⋅„ − 1ƒ23o)

ƒn− oQ\ = ƒ⋅„−1ƒ⋅„ −o⋅nQ€ = @ ⋅ƒ⋅„ − 1

過渡解析のための時間離散化の方法
時間変化の解析の方法は以下の二通りが考えられる^

m0

B@ + ^① 陽解法：
00+Δ0M12 ⋅⋅−2⋅• は1つ前の時刻の値（!
時間o−!

j0′o

⋅ d Bl = ^ GR B@m’•など）から直接求める⋅ j! ⋅ j$ + s12′n + ′o ′n − ′o⋅j$ + ⋯2j!

• 初期状態が求まればそれ以降は高速に解析が進む• 安定性と精度のために時間刻みを小さくする② 陰解法： 全ての値を各時刻でその都度求め直す
00+Δ0!
時間o−!

j0′o⋅ j! ⋅ j$ + s12n+2o⋅n−j!

oj$ + ⋯

• 各時刻において未知数の求解が必要• 比較的大きな時間刻みでも結果が安定24ここまでの定式化結果のおさらい（その1）[再掲]
離散化まで済んだ現在の問題：
未知数 3 、式数 3 の連立方程式（ ：格子点数）
各式は一般に非線形（2次以上の項、指数、分数を含む）
各式は中心および近隣の格子点上の未知数( , , )の関数

z2 = ( 2 ,2, 2)

zp = ( p ,zq = ( q ,zo = ( o ,o, o)

点~点の数分点•

ozn = ( n ,o)n{n| (zo , zn , zp , zq , z2 )=0o)}o{n (zo , zn , zp , zq , z2 )=0点p, p)q, q))5o q (zo , zn , zp , zq , z2 )=0

計3 本（ =格子点数）の方程式が得られるn, n)25

非線形連立方程式の求解
非線形方程式の Newton 法による解法（1元の場合）• 以下の ) ! に1次のTailor展開を施した下記の方程式を解く。
B) ∗∗) ! +! ⋅ Š! = 0B!
• 求めた Š! を用いて !* を更新する。 (!* + Š! → !*)• これを Š! が十分小さくなるまで繰り返し、得られた !* を解 ! とする。
26

連立方程式の線形化
多元方程式の線形化
多変数方程式のTailor展開は一般に以下のようになる) g∗) ∗) ∗) ∗+! ⋅ Š!\ +! ⋅ Š! +! ⋅ Š!Ž + ⋯ = 0!\!
!Ž

複数の式を連立する
以下のような行列形式に書ける式の数と未知数の数がともに 3 なので、正方行列になる)\ ∗)\ ∗)\Jacobi 行列∗!
!\)!∗!\⋮)Ž• ∗!
!\!
!
)!∗!
⋮)Ž• ∗!
!

⋯⋯⋱⋯!
!Ž•)
!∗!Ž•⋮)Ž• ∗!
!Ž•27Š!\Š!
⋮Š!Ž•)\ g∗) g∗=−⋮)Ž• g∗· Š! = “ここまでの定式化結果のおさらい（その2）線形化が済んだ現在の問題：
非線形連立方程式を線形化（Tailor展開）しNewton法で解く大規模（3 ⨯3 ）連立一次方程式を繰り返し解く
スタート（初期解 !* 設定）! = !* における行列 と右辺“を設定（連立方程式 · Š! = “ を立式）
Š! を求める（Š! = ⁻¹·“ を計算）

： 3 ⨯3 のJacobi 行列
Š!, “：サイズ 3 のベクトル3 ⨯3 行列の逆行列の計算または同等の計算量が必要!* を更新（!* + Š! → !*）
No

収束した?
（Š! が小さい?）

Yes

28

終了（!* を解 ! とする）

連立方程式の求解（行列解法）
大規模な連立一次方程式を解く場合の問題
直接法（真面目な求解）はメモリ使用量や計算時間が膨大
メモリ使用量： O( ²) 、計算量：O( ³)
直接法の例： Gauss消去法
š\,\
š ,\
⋮
š2,\

š\,
š ,
⋮
š2,1 š′\,
01
⋮
⋮
00

⋯ š\,2⋯ š ,2⋱⋮⋯ š2,2!\!
⋮!2=① 前進消去

⋯ š′\,2⋯ š′ ,2⋱ ⋮⋯ 1!\!
⋮!2=›\›⋮›2② 後退代入!2 = ›′2!2 \ = ›′2 \ − š′2 \,2 !2›′\›′⋮›′2⋮!\ = ›′\ − (šœ\, ! + ⋯ + šœ\,2 !2 )29

LU分解と直接解法

行列 を = •·ž の形の積に分解しておく•: 下三角行列、ž: 上三角行列
š\,\
š ,\
⋮
š2,\š\,š ,⋮š2,⋯ š\,2⋯ š ,2⋱⋮⋯ š2,2=1

Ÿ ,\⋮Ÿ2,\01
⋮

⋯ 0
⋯ 0
⋱ ⋮
⋯ 1

Ÿ2,Y\,\0
⋮
0Y\,
Y ,
⋮
0= •·ž の形に分解できているとして・・·! = “•·ž! = “•· = “ž·! == •-1·“! = ž-1·⋯ Y\,2⋯ Y ,2⋱⋮⋯ Y2,2① 前進代入で計算可能② 後退代入で計算可能

LU分解のメリットLU分解自体の計算量を考えるとトータルの計算量は不変
一つの を複数の“に対して解く場合は分解が一度で済む
30

反復解法による行列計算（1）
直接解法はメモリ使用量や計算量が多い
大規模行列では反復解法による収束計算が一般的
反復解法の種類
定常法： 同じ行列、ベクトルを使用した漸化式・遅い
非定常法：残差ベクトルの基底の直交性を利用・速い定常法の例（Jacobi法）·g = “ の代わりに g(>+1) = ¢·g(>) + £ を繰り返し解く
どんな ¢ と £ を使えば良い?
!

( [\)

1
=
› − š ,\ !\ − ⋯ − š ,2 !2š,U 以外その心は・・・

š ,\ !\ + ⋯ + š , !

( [\)+ ⋯ + š ,2 !2 = ›を !U(>+1) について（それ以外は既知として）解いている
31

反復解法による行列計算（2）
より積極的に真の解を目指すような解の修正を考える
下記漸化式で修正方向ベクトル ¤U と係数 TU を求める
g( [\) = g( ) + T2 \⋅ ¤ ⟶ g = g€ + ¦g( )T\ ⋅ ¤\§€T ⋅¤g€T€ ⋅ ¤€g(\)¤U, TU の求め方については、様々な方法が提案されている
共役勾配法/CG法 （Conjugate Gradient Method）
BCG法 （Bi-Conjugate Gradient Method）非対称行列
BiCGSTAB法 （Bi-Conjugate Gradient STABilized）
32

行列計算のテクニック
デバイスシミュレーションで扱うJacobi行列をよく見てみると1本の式に現れる未知数の数は 3 よりもずっと少ない
行列中の殆どの要素は 0 （=疎行列）
非ゼロの行列要素は対角に対して対称に現れる構造格子（直交格子）の場合は対角要素に並行した非ゼロ要素列が現れるバンド行列になる
疎行列

バンド行列nnwoewoess33

大規模疎行列の行列解法
大規模疎行列（large sparse matrix）ならではの解法がある行列の（メモリへの）格納
非ゼロ要素は無視（不完全）LU分解を利用した反復解法
行列の分割と並列化（行列の入れ替え）「疎行列らしさ」が失われてしまう例
疎行列（バンド行列）•（真面目な）
LU分解

žx

フィルイン（Fill-in）元の行列では０だった所に数値が現れる
34

不完全LU分解による前処理と反復解法
大規模疎行列は「手を抜いて」解くことができる
反復解法で収束させれば良い
厳密な逆行列は要らない
LU 分解でフィルインを手抜く＝不完全LU分解
疎行列
•ž
不完全
LU分解

x

フィルイン
レベル０·! = “ の代わりに (•ž) -1 ·! = (•ž) -1“ を解く(•ž) -1 = ¨+(•ž) -1© は ¨ に近いと考えられる（© ≡ – •·ž：誤差行列）35

レベル0ではフィルインは考えない

大規模行列解法の高速化・並列化
最大の課題である不完全ILU分解の並列処理による高速化
対象の行列をいくつかのブロックに分割して並列処理する
4x4 のブロックを4台のPEに分散

PE#1
PE#2

デバイス構造同時にスタート

順番にスタート

PE#0

領域分割と行（接点）の順序入れ替え

PE#3• 処理を開始するまでの• デバイスの構造（形状）を反映した分割

待機時間が生じる• PE間の通信量が多い• 領域間の境界面も独立した領域とする• 並列度向上＆通信量削減

36

半導体デバイスシミュレーションの数値解析のまとめ
支配方程式（ポアソン方程式、電子・正孔電流連続式）
空間の離散化 ・・・ 有限体積法
時間の離散化 ・・・ 陰解法• キャリア移動度• 生成消滅項• 積分形への変換（Gaussの定理）• Sharfetter-Gummel法• Tailor 展開（一次）

線形化• Newton法による連立一次方程式（大規模疎行列）

反復収束求解• 反復解法（共役勾配法）• 不完全LU分解

高速行列解法• 領域分割・並列化

解析結果（未知数分布 デバイス特性）
37• デバイス性能予測• プロセス評価・改善

半導体デバイスシミュレーションの展望38

半導体デバイスシミュレーションの現状と展望
半導体技術をとりまく環境の変化
微細化の減速・限界と多様なブレークスルーの方向性
構造：FinFET、３D集積、微細形状
材料：SiC、SiGe、III-V、HBT、high/low-k
アプリケーション領域の拡大
Flash メモリ、相変化・抵抗メモリ
パワーデバイス
イメージセンサ、太陽電池
新たな課題
高信頼性（経年劣化、ソフトエラー）
生産性（バラつき、歩留まり、コスト）
39

半導体デバイスシミュレーションの現状と展望
半導体デバイスシミュレーションの進化
新たな物理現象とモデルへの対応（マルチフィジックス）
熱、ストレス、光子吸収・放射、電磁界高エネルギーキャリア輸送、キャリア量子化、トラップ
化合物材料、ヘテロ接合、結晶方位
複雑な解析対象への対応
静電破壊、電荷漏洩
ソフトエラー、ラッチアップ、熱暴走
多様な構造への対応
大規模3次元解析、回路-デバイス連成解析
メッシュ生成、高速ソルバー、並列化
40ありがとうございました

41