数値計算による半導体物性解析手法について-量子力学の基礎清水 共∗ , 森

優策† , 内海 太禄‡Analysis of Semiconductor Physics by use of Numerical Simulation- a Basic Problem of Quantum Mechanics Tomo SHIMIZU, Yusaku MORI, Taro UTSUMI

Abstract
This paper deals with a method of numerical simulation in consideration of a quantum
eﬀect at nano-devices. With the problem of the one dimensional infinite square well,a calculation method to evaluate the electron distribution is suggested. Schrödinger
equation and Poisson equation in x-direction perpendicular to the well potential areanalyzed. CLAPACK routines are used in the simulation.

Keywords: numerical simulation，Poisson equation, Schrödinger equation, diﬀerence method

はじめに301

現在，身の回りの生活空間にはたくさんの電子機器があふれている 1) 。これらの電子機器はそのほとんどに半導体が組み込まれている。このように半導体産業は社会システムを支える重要な基幹産業である。半導体産業は，1947 年の点接触型トランジスタの実験に始まり，真空管に変わる固体デバイスであるトランジスタが登場して，1960 年代にMulti Gate

0.6

20

FD SOI
0.4

150.2

1002010IC(Intefrated Circuit) が開発されて多くの電子機器が世界中に普及した。

250.8

2012

2014

2016

2018

2020

2022

2024

MPU Physical Gate Length (nm)Extended Planar Bulk
EOT: Equivalent Oxide Thickness (nm)152026

Year of production

図 1:

ITRS2011 による微細化トレンド (Highperformance Logic Technology)

IC の発展は，搭載するデバイス (トランジスタなどの素子) の加工寸法を小さくする「微細化」技術の進歩に対応する。デバイスの微細化は，様々な恩量・薄膜化，低消費電力化，低価格化，信頼性向上恵を与える。IC の高性能化，高集積化，小型・軽など，様々なメリットを与え半導体デバイスの中心的役割を担ってきた MOSFET を例にすると比例縮∗ 香川高等専門学校詫間キャンパス 電子システム工学科

小則に従った単純なデバイス構造の微細化の追求が† 香川高等専門学校専攻科 電子情報通信工学専攻

半導体産業の発展に直結した。

‡ 長岡技術科学大学 電気電子情報工学課程

97

独立行政法人国立高等専門学校機構香川高等専門学校研究紀要 6(2015)

Source伴い，デバイス特性を解析するため利用される従来

Gate

Gate

Sourceの比較的容易な解析式は不十分である。特に，デバ

Drain

Drain

イスの微細化に伴うキャリアすなわち電子の扱いが重要となる。電子の波動性を扱う量子力学的効果はSiO2 BOXとして重要なものに，干渉，トンネリング，エネル

Si substrate

Si substrate

FDSOI MOSFET

Bulk MOSFET

ギーの量子化の三つの効果があげられる。

(Fully Depleted Silicon on Insulator)(A)

デバイス特性の解析方法の検討が必要である一方(B)で，個人使用のパソコンの性能向上に伴い，一昔前
Drain

Gate
Sourceのスーパーコンピュータ並の非常に高性能な数値計

Gate

Drain

算環境が安価で入手可能である。これは，大きな研Source

究所や大学に所属しなくても学生が非常に安価な手
SiO2 BOXGate
DG FET(Double Gate)(C)段で高度な数値計算環境を利用可能であることを意

Si substrate

味する。電子の波動性を議論する上で重要な量子力

Tri Gate (Fin) FET

学の基本方程式である波動方程式は，Schrödinger(D)

方程式であり二次の偏微分方程式である。
本報告書では，極微細な半導体デバイスの特性解

図 2: 電界効果型トランジスタの 3 次元構造の比較：

析において，極小領域に対する電子の取り扱い，す(A)Bulk MOSFET, (B)FDSOI MOSFET, (C)DG
FET, (D) Tri-Gate FET

なわち，量子力学的効果の導入の重要性を鑑みて，その基礎となる１次元の井戸型ポテンシャル問題を数値計算により評価する。これは，図 2(C) で紹介した DG FET の gate 幅方向の電子閉じ込めを考

図 1 は，ITRS2011 が公開している半導体デバイ慮したデバイスモデルを想定したデバイス特性解析

スの物理ゲート長とゲート酸化膜厚のトレンドであの基礎となる。

る 2) 。前述したように，デバイスの微細化の結果としてデバイス寸法が nm オーダーへ突入し，近年，一2

般社会においても話題とされる”ナノテクノロジー”

数値計算ライブラリ導入が半導体技術の主戦場となる。酸化膜厚に注目すれ工学分野において，代数的な方法で解を得ることば，EOT(Equivalent Oxide Thickness) で 1nm 以が困難な問題を数値を用いて近似的に解く手法がコ下が主流となり，この厚さは原子数層分に対応する。

ンピュータ技術の発展によって注目されている。この

これは，従来から行われているバルク結晶的な物性ような数値解析シミュレーションでは多次元行列演特性を利用したデバイス解析から，デバイス・結晶算や高速フーリエ変換等の数学的処理が必要不可欠を構成す分子や原子の特性を考慮したデバイス解析である。これらの計算処理のアルゴリズムは確立さの重要性と必要性が生じることに対応する。

れた技術である。そのため，多数の計算ライブラリまた，EOT の薄膜化に代表されるように比例縮がさまざまな方法で提供されている。そこで，デバ小則に従ったデバイスサイズの微細化は，従来のイス解析に利用する数値計算では，学生が個々に計MOSFET 構造の物理的限界に直面し，多くの新しいデバイス構造が生み出されている。図 2 は，従来算ライブラリを構築する事を比較的重要性の低い作の半導体デバイスの代表格であるである MOSFETの信頼性が高いオープンソースソフトウェアであるの 3 次元構造と，近年注目される電界効果型トラ

有名な線形代数演算ライブラリの LAPACK(Linearンジスタである FDSOI MOSFET，DG FET, Tri-Algebra PACKage) をここでは採用する。なかでもプラットフォームとして現在主流である Windows業として割り切り，非常に高速に実行可能で，計算3)

Gate FET の 3 次元構造を示している 。図 (C) のDouble Gate 構造に注目すれば，source-drain 間の電流路であるチャネルが二つの gate によって挟まOS で，C 言語で利用可な数値計算ライブラリである CLAPACK を導入して数値計算を実行する 4) 。

れ，デバイス構造の微細化により gate 方向に量子

表 1 は，この報告書で使用する CLAPACK の二井戸が形成されチャネルを形成する電子の閉じ込めつの関数である。dgesv 関数は Poisson 方程式に，が起こる。

dstev 関数は Schrödinger 方程式を数値計算するために利用する。

新しいデバイス構造の登場やデバイスの微細化に98清水共 他 : 数値計算による半導体物性解析手法について 清量子力学の基礎清式として表される。
表 1: CLAPACK 関数
関数名h̄2−2m用途
一般行列の連立一次方程式を計算

dgesv
dstev

実対称三重対角行列の固有値1(∆x)−2ϕ(xj )と固有ベクトルを計算

3(2 {ϕ(xj+1 ) + ϕ(xj−1 )ª)+V (xj )ϕ(xj ) = εϕ(xj )上式を α, βj , γ, E を導入して書き直すと次式となる。

Schrödinger 方程式の差分化

αϕ(xj+1 ) + βj ϕ(xj ) + γϕ(xj−1 ) = Eϕ(xj )量子力学の基本方程式である Schrödinger 方程式(5)に注目して，その数値シミュレーション方法を紹介する。電子が de Broglie 波長程度の領域に閉じ込

ここで,められたときの量子力学の基本問題である 1 次元(x 軸方向) に対して，井戸型ポテンシャル問題を考える。質量 m[kg] の 1 つの電子が，幅 L[m] の空間に以下に示すポテンシャル V (x) によって閉じ込めαられている。ここで，井戸幅方向である x 軸の原点(x = 0) を井戸端として，もう一方の端を x = L とする。
(0 0≤x≤LV (x) =(1)∞ 0 > x, L < xSchrödinger 方程式は次式である。
¾½h̄2 d2+V(x)ϕ(x) = εϕ(x)−2m dx2(4)1=γ=βj=−(E=−2mεh̄22(∆x))
2m+ 2 V (xj )h̄(∆x)22である。式 (1) の境界条件により，電子は井戸の外で存在できない。すなわち，波動関数 ϕ(x) は井戸端で零となる。

(2)ϕ(x0 ) = ϕ(xN +1 ) = 0

ここで，V (x) はポテンシャル，h̄ は Planck 定数，εは固有エネルギー，ϕ(x) は波動関数である。

(6)

式 (5) と式 (6) の関係から，井戸内の j = 1, 2, ..., N

上記の 1 次元 Schrödinger 方程式は 2 階の微分方の各 xj 点で Schrödinger 方程式を立てると次のよ程式である。よって，数値計算するために微分方程うな N 個の方程式ができる。

式を差分化する。計算に使用する x 方向の分割（メッシュ）間隔を ∆x として離散化することで，井戸内を N 点に分割する。井戸端 x = 0 を x0 として，もαϕ(x2 ) + β1 ϕ(x1 )=Eϕ(x1 )う一方の井戸端 x = L を xN +1 とする。波動関数αϕ(x3 ) + β2 ϕ(x2 ) + γϕ(x1 )=Eϕ(x2 )ϕ(x) は ϕ(xj ) として離散値データとして扱われる。
ここで，j = 0, 1, 2, ..., N, N + 1 である。二回微分αϕ(x4 ) + β3 ϕ(x3 ) + γϕ(x2 )..
.

=Eϕ(x3 )

=Eϕ(xN −2 )+γϕ(xN −2 )=Eϕ(xN −1 )βN ϕ(xN ) + γϕ(xN −1 )=Eϕ(xN )の波動関数は，下記のように差分化される。
¾½dd2dϕ(x)ϕ(x)=dx2dx dx½¾d ϕ(xj+1 ) − ϕ(xj )=dx∆x==αϕ(xN −1 ) + βN −2 ϕ(xN −2 )+γϕ(xN −3 )αϕ(xN ) + βN −1 ϕ(xN −1 )ϕ(xj+1 )−ϕ(xj )ϕ(xj )−ϕ(xj−1 )−∆x∆x∆x12 {ϕ(xj+1 )(∆x)+ϕ(xj−1 ) − 2ϕ(xj )}

(7)
(3)

ここで，式 (7) は固有値問題であるため，量子数と式 (2) を式 (3) によって差分方程式に変換すること

して添字 i を用いて波動関数 ϕ(x) と固有エネルギー

で，位置 x = xj における Schrödinger 方程式は下E をそれぞれ ϕi (x) と固有エネルギー Ei と表記す99

独立行政法人国立高等専門学校機構香川高等専門学校研究紀要 6(2015)

る。上式を行列形式で記述すると次式となる。
⎡⎤
β1 α0 ···0⎢⎥⎢ γ β2α0···0 ⎥⎢⎥⎢ 0γβ3α0···0 ⎥⎢⎥⎢ .
.. ⎥..
..
..
..
..
⎢ .
⎥.
.
.
.
.
. ⎥⎢ .
⎢⎥⎢ 0 ···α0 ⎥0γ βN −2⎢⎥
⎢
⎥
···
0γ
βN −1 α ⎦⎣ 00···0γβN⎡⎤⎤⎡
ϕi (x1 )ϕi (x1 )⎢⎥⎥⎢⎢ ϕi (x2 ) ⎥⎢ ϕi (x2 ) ⎥⎢⎥⎥⎢⎢ ϕ (x ) ⎥⎢ ϕ (x ) ⎥i 3i 3⎢⎥
⎥⎢⎢⎥⎥
⎢
..
..
⎥ = Ei ⎢⎥ (8)×⎢.
.
⎢⎥
⎥⎢⎢⎥
⎥
⎢
⎢ ϕi (xN −2 ) ⎥⎢ ϕi (xN −2 ) ⎥⎢⎥⎥⎢⎢⎥⎥⎢⎣ ϕi (xN −1 ) ⎦⎣ ϕi (xN −1 ) ⎦ϕi (xN )ϕi (xN )0.25

n=2n=3/ 10⁹

0.15|ψn|²0.20

0.10

0.05

0.00
0246810x [nm]図 3: 井戸幅 10[nm] の井戸型ポテンシャルに対する波動関数の分布
1.00
n=1accumulation probabilityn=2

上式を次式の様に表記すれば，行列 A の固有値問題である。

Aϕ~i = Ei ϕ~in=1(9)n=30.80

0.60

0.40

0.20

0.00
024ここで，行列 A が式 (8) から分かるように三重対角

6810x [nm]行列となる。よって，前節の表 1 で紹介した dstev

図 4: 井戸幅 10[nm] の井戸型ポテンシャルに対す関数を使用して Schrödinger 方程式を数値計算する。

る各状態の累積存在確率分布
計算コスト（時間）の増加に直接関係するため，計表 2: 固有エネルギーの計算結果
量子数

状態n基底

1234

第 1 励起第 2 励起第 3 励起算精度には多方面からの検討が必要となる。

固有エネルギー

図 3 は，表 2 の結果と伴に数値計算から得られた

En [meV ]波動関数 ϕn (x) に対して，下式により波動関数が規3.8115.2234.2460.87

格化された n = 1, 2, 3 に対応する基底状態と励起状態の存在確率密度分布である。
Z L
NX|ϕn (x)|2 dx =|ϕn (xj )|2 ∆x0j=1=

表 2 は，井戸幅 L = 10[nm] として数値計算にれているため各状態の電子が井戸内にのみ存在するシャルに対して，解析的に Schrödinger 方程式を解ことを現している。この波動関数も固有エネルギーけば次式として固有エネルギーが得られる。と同様に，井戸型ポテンシャル問題の解析値が良く

2 2h̄ π 2n2mL2(n = 1, 2, ...)

(11)

式 (11) は，井戸が無限高のポテンシャル障壁に囲まより得たデータである。無限障壁高の井戸型ポテン

En =1知られるている。計算結果は解析値と良い一致を示(10)している。

計算された固有エネルギー値は式 (10) の解析値と図 4 は，存在確率密度分布から求めた累積存在確率比較して 1% 程度の誤差が発生している。これは，Pn (x) である。井戸端 x = 0 を基準として Pn (0) =0，井戸端 x = L で Pn (L) = 1 である。

今回の数値計算において使用した１次元方向の少ない分割数 N = 200 が原因である。分割数と誤差の関係は別途計算しているが，本質と異なるため今回

Pn (xj ) =jX
y=1

詳細については割愛する。N が増すにつれて計算の精度が向上することを確認しているが，N の増大は100|ϕn (xy )|2 ∆x(n = 1, 2, ...)
(12)

清水共 他 : 数値計算による半導体物性解析手法について 清量子力学の基礎清Pn (xj ) 分布により，井戸端 (半導体デバイスでは界面) 近傍における電子分布が推測できる。

250

200

ネルギーが低い状態から順に状態を占有する。

Schrödinger 方程式の計算結果から言えば，表 2 で与えられる基底状態から順に電子が占有される。実際の電子の状態占有については，ここでは扱わない。

potential [meV]

電 子 は ，Pauli exclusion principle に 従 い エ15010050

図 4 から，少なくとも井戸型ポテンシャルにおい0ては，電子状態は基底状態から電子のエネルギーが046810x [nm]増加するに連れて井戸端に向かって電子の存在確率

図 5: 井戸内の一様電荷分布によるポテンシャル分布が増加する事が分かる。また，電子は井戸端の極近傍に向かいその存在確率が急激に減少して，ある程ンシャルエネルギーと比較すると，Poisson 方程式度の電子がほとんど存在しない領域が形成されるこによる電子自身のポテンシャルを考慮すれば電子のとが分かる。すなわち，井戸内の電子濃度の変化が固有エネルギーと波動関数が大きく影響を受ける事井戸端近傍における電子分布に大きく影響することが推測される。例えば，電子は自身のポテンシャルが推測される。の影響を受けより界面側へ分布がシフトすることがここでの井戸型ポテンシャル問題は無限高のポテ推測できる。すなわち，実際に電子自身によるポテンシャル障壁を仮定しているが，実際問題としてデンシャルの影響を考慮するためには，Schrödinger

バイスのモデル化には有限高のポテンシャル障壁問方程式と Poisson 方程式を繰り返し計算して，パラ題を扱うことが必要である。そのため，波動関数はメータの変化が収束する様に計算する自己無撞着ポテンシャル障壁内へ侵入する。よって，電子分布(self-consistent) 法を導入することが必要となる。

はポテンシャル障壁側へシフトする。これについては，別の機会に研究紀要で紹介する予定である。

425Poisson 方程式

おわりに半導体デバイスの微細化により，教育現場においても量子効果を取り入れたデバイス特性解析の必要

x 軸方向に対する１次元の Poisson 方程式は次式である。

性が増している現状を踏まえて，Schrödinger 方程式と Poisson 方程式を例に数値シミュレーションの方法を報告する。数値計算においては，線形代数演2

ρ(x)dV (x) = −dx2εa(13)

算ライブラリ LAPACK を採用して利便性の向上を図る。極微細な半導体デバイスの簡易モデル化によ

ここで，V (x) は電位分布，ρ(x) は空間電荷密度分り，この数値計算から量子効果の影響が推測できる

布，εa は誘電率である。前節の Schrödinger 方程式

ことが確認できる。と同様の手順により，方程式を差分化して，x 軸方向に対して離散データを扱うことで，V (xj ) に対して j = 1, 2, ..., N の N 個の方程式を立てることが可参考文献

能となり，連立一次方程式を得る。第 2 節の表 1 で1) ”IC ガイドブック未来を創る半導体”, 社団法人
電子情報技術産業協会, 2012

紹介した dgesv 関数を利用することで，数値計算できる。
図 5 は，前節で使用した井戸幅 L = 10[nm] の１

2) International Technology
Semiconductors (ITRS),(http://public.itrs.net/)次元井戸型ポテンシャルを想定して，電子が 1.2 ×1020 [cm−3 ] とかなり高濃度で井戸内に一様に分布していると仮定した場合に，Poisson 方程式を上記の関数により数値計算して得たポテンシャル分布でRoadmap for2011 Edition3) 平本俊郎 編著, ”集積ナノデバイス”, 丸善株式会社, 2009ある。高濃度の電子分布により，約 230[meV ] のポテンシャル障壁が井戸中央に発生する。

4) Netlib Repository at UTK and ORNL
(http://www.netlib.org/)

前節で計算した表 2 の固有エネルギー値から分かるように，単純に Poisson 方程式から得られたポテ101