計算数理工学論文集 Vol. 10 (2010 年 12 月), 論文 No. 03-101210

JASCOME

局所ゲージ変換を用いた磁場中電子の新たな有限要素法
NOVEL FEM WITH PEIERLS PHASE FOR ELECTRON WAVES IN MAGNETIC FIELDS
植田毅 1) 、宮川悠 2)

Tsuyoshi Ueta and Yuu Miyagawa
1) 千葉大学総合メディア基盤センター(〒 263-8522

千葉市稲毛区弥生町 1-33,

E-mail: ueta@faculty.chiba-u.jp)

2) 千葉大学大学院自然科学研究科(〒 263-8522

千葉市稲毛区弥生町 1-33)

The ﬁnite element method has already been extended for analyses of transport propertiesof the electron wave of two-dimensional electron systems in magnetic ﬁelds. Not a few
proposals of the new formulation or improvement have been reported, whereas thereexists almost no report of the application to the analysis of a realistic system. In thepresent paper, it is pointed out that analyses with FEM of conventional formulation do
not give accurate results in case of a large system or of a strong magnetic ﬁeld and thatthe selected gauge inﬂuences numerical results greatly. Furthermore, the conceptually
new ﬁnite element method which solves the diﬃculty fundamentally is proposed. The
formulation in which matrix elements in the absence of a magnetic ﬁeld are multipliedby the Peierls phase is very simple. In order to show its predominancy, the numericalexamples of the proposed formulation are presented.
Key Words : Two-dimensional Electron Systems, Magnetic Fields, Magnetic Translation
Operator, Finite Element Method, Peierls Phase

1. はじめに

本論文では、磁場中の電子波解析に有限要素法が用いら絶縁体-半導体ヘテロ界面に形成される２次元電子系におれない原因を指摘し、それを解消する全く新しい定式化を行いて、デバイスの構造のサイズを超えるバリスティックな伝う。文献 (11) の計算例と同じ系の数値計算結果を示し、本論導が実現されている。このような電子系を用いたデバイスの文で提案する手法の簡便性、優位性を示す。

設計においては複雑な形状の系内での電子波の振舞いを厳密に解析することが不可欠である。また、外界からの制御パラメータの一つとして磁場が用いられるため、磁場中での電子波の解析をする必要がある (1) 。

2. これまでの定式化とその問題点ベクトルポテンシャルが A(r) であるとき、電子波はシュレディンガー方程式}
1 {(−ih̄∇ − qA)2 + V (r) ψ = Eψ2m

これまで、磁場中の電子波の解析にモードマッチング法(2)、差分法、強束縛近似 (1) 、有限要素法 (3, 4) 、境界要素法(1)などが用いられてきた。複雑な形状の系の解析には境に従う。m、q 、E はそれぞれ電子の有効質量、電荷、エネル

界要素法が最も便利であり、磁場が印加された場合にも拡張ギーである。h̄、i はそれぞれプランク定数 h を 2π で割ったされ、多くの成果をあげてきた (6, 7) 。

定数および虚数単位である。

(4, 5)

しかしながら、系内でポテンシャルの変化がある場合やス２次元電子系を x-y 平面にとり、強さ B の一様磁場が z 軸ピン軌道相互作用、ラシュバ効果 (8) がある場合には境界要正方向に印加されているとすると、ベクトルポテンシャルは素法による解析は困難であり、有限要素法も磁場中の電子波A = B (−(1 − s)y, sx, 0)（s は任意の実数）と選べる。これのために拡張され、(9, 10) 改良されてきた (11) 。磁場中電子までの定式化では波の有限要素法の定式化と簡単な計算例の報告は多いが、それを用いた物理現象の解析 (9) 、現実的デバイスの特性の解析の報告は少ない。

2010 年 9 月 30 日受付，2010 年 11 月 10 日受理{}h̄2 2qBh̄∂ψ∂ψ∇ ψ+i−(1 − s)y+ sx2mm∂x∂y2 2 {}q B+(1 − s)2 y 2 + s2 x2 ψ = {E − V (r)} ψ (2)2m−と展開して、左辺の第２、３項を磁場が印加されることにより付加されるポテンシャルとして扱う。すなわち、局所的にy波数の変化する電磁波と同様な定式化を行う。
解析領域を小さな領域（有限要素）に分割し、j 番目の節d2点の座標を r j とする。j 番目の節点でのみ値が 1、他の全ての節点では値が 0、つまり{Nj (r k ) =xO1

j=k0

j ̸= k(3)-d2-ddを満たし、節点 j を囲む節点との間を線形につなぐ r の連続実関数を j 番目の節点における形状関数（補間関数）Nj (r)Fig. 1とする。

The rectangular quantum dot investigated in ref.(11)形状関数 Nj (r) を用いて領域内の波動関数を

ψ(r) =∑Nj (r)ψj(4)j(11) の Table I に示し、「要素数を増やすとゲージによらずと展開し、これを試行関数とする。このとき、汎関数

L = ψ † (K − M) ψ − ψ † Q同じ値に収束する」と述べている。しかし、その Table I を(5)を ψ † で変分することにより、{ψi } についての連立方程式(K − M) ψ = Q

いて解析している。エネルギー固有値の低い３つについて、√kd = 2mEd/h̄ で定義される無次元化した波数の値を文献見ると、要素数 1250 の場合においても、固有波数の小数点第 3 桁目はゲージのよって異なっている。
また、ベクトルポテンシャル依存性が明瞭になる計算例として、Fig.2 に示す x 方向に幅 d、y 方向 8d の長方形領域の中(6)央、左右に幅 d のエミッタ (x < −d/2) とコレクタ (x > d/2)を得る。ここで、ψ を各節点の波動関数 ψi を第 i 成分とすを取り付けた十字構造を考える。エミッタ、コレクタ部分のる列ベクトルとして定義した。また、行列 K、M およびベクyトル Q は以下のように定義される。
()∫∂Nj ∂Nk∂Nj ∂Nk
Kjk =+dr(7)∂y ∂yvjk ∂x ∂x∫
Mjk = ENj Nk drvjk{∫∂NkqBh̄ ∑−(1 − s)yk
Nj Nl−idrm∂xvjkl}∫∂Nk+sxk
Nj Nldr∂yvjk]∑[}
q2 B 2 {−V (r l ) +(1 − s)2 yl2 + s2 x2l2ml∫
×
Nj Nk Nl dr(8)vjk∫∑
Qj =q(r k )Nj Nk dS(9)SNk4dEmitterCollectorxO-d2d2-4dここで、vjk は節点 j と節点 k を共に含む要素内での積分で∫あることを示す。また、 S は境界上で波動関数の法線方向N微分の値が与えたれた部分の面積積分であり、q(r) は境界上

Fig. 2

The geometry of a cross-shaped conductor.

で与えられた境界面外向き法線方向微分の値である。
この定式化を用いると、単純な差分法を用いているときから指摘されていたことではあるが、行列 M の最後の項は系が大きな場合、磁場が強い場合に非常に大きくなる可能性がある。その振舞いは s、すなわちベクトルポテンシャルのゲージの取り方に依存する。
文献 (11) では、横幅 2d、縦 d の長方形量子ドット (Fig.1)
について、境界上で波動関数が 0 となる境界条件の固有値問
題を計算例として、4 種のベクトルポテンシャルの場合につ境界条件の取り扱いは文献 (9) の II-B 節の取り扱いに従う。
ただし、ここでは計算を簡便にするために、エミッタ、コレクタ部分には磁場がかかっていないものとし、これらの部分のベクトルポテンシャルはそれぞれ磁場がかかっている部分と連続になるようなゲージにとる。
磁場がかかっている部分のゲージを A = (−By, 0, 0) として、エミッタから kd = 10 の基本モード波が入射した場合の透過率 T 、反射率 R を磁場 B̃ ≡ qBd2 /h̄ = 0 ∼ 5 に対して、接点数 Ns を変化させて計算した。確率保存則から T + R = 1られる。

を満たさなければならないが、計算結果では Ns ≥ 2000 にそこで、ベクトルポテンシャルを A = B/2(−y, x, 0)、A =対しては、|T + R − 1| の誤差は 0.02%以下になっている。しB(0, x, 0) と し た 場 合 の 計 算 結 果 を Fig.5, 6 に 示 す。A =B(0, x, 0) の場合には、このような系でも計算精度の良い境

1界要素法で計算した透過率および波動関数が非常によく一致

B=0
B=1
B=2
B=3
B=4
B=50.80.6

している。また、Fig.6 (b) では位相は |y| が大きくなろうとも激しく振動していない。これらの結果から、原点から離れたところでの位相の激しい変化が収束を悪くしていることがT分かる。 実際、A = (0, Bx, 0) から A = (−By, 0, 0) にゲー

0.4(a)

(b)

0.2
0
1000

2000

3000

4000

5000

6000

7000

NsFig. 3

The convergence characteristic of the transmission

probability.
かしながら、Fig.3 に示すように、B̃ = 0 の場合、透過率 Tは Ns = 3000 で一定値に収束しているが、磁場が大きくなると Ns を大きくしても透過率 T はほとんど収束していないことが分かる。Ns = 3000 のとき、接点間の距離が 1/8 波長

Fig. 5になっており、B̃ = 0 の場合の収束は理解できる。磁場がかphase (b) for kd = 10 and B̃ = 20 for A = B/2(−y, x, 0).

Density plots of the probability density (a) and the

かっている場合の収束の悪さの原因を見るために、B̃ = 20の場合の確率密度 |ψ|2 と波動関数の位相 Im (ln ψ) を Fig.4(a)(b)に示す。この場合、電子のサイクロトロン半径 rc とエミッタ(a)(b)42πFig. 60Fig. 40

Density plots of the probability density (a) and thephase (b) for kd = 10 and B̃ = 20 for A = (0, Bx, 0).

Density plots of the probability density (a) and thephase (b) of a wave function for kd = 10 and B̃ = 20 in
A = (−By, 0, 0), which is an insuﬃcient convergence case.

ジ変換した場合、波動関数に exp (−iqBxy/h̄) の位相因子がかかる。したがって、x 方向に y と B に比例した位相変化が加わる。その意味で、この系の解析に A = (−By, 0, 0) は最の幅 d の比は rc /d = kd/B̃ = 0.5 であり、電子の古典軌道かも不適なゲージであったと言える。この因子は |ψ|2 には全くらすると、確率密度は系全体に広がるものと考えられるが、現れないため、確率密度を確認しただけではその振動に気がFig.4 (a) では系の端まで広がっていないことが分かる。まつかない。

た、|ψ| は |y| が大きくなっても激しい変化はしないが、位系の形が特定の方向に広がりがある場合にはその方向に相は |y| が大きくなると激しく変化していることがわかる。

依存性がないゲージを取れば解決できる。しかし、系が円形

この場合、|y| が大きくなったときに激しく変化する位相をなど等方的で大きな場合には大局的ゲージ変換では対処でき正しく表現できていないことが原因で収束していないと考えない。そこで、本論文ではこの問題を根本的に解決する新し2

い概念の定式化を提案し、計算例によりその有効性を示す。

は、点 r j におけるベクトルポテンシャルが 0 になるように局所ゲージ変換した波動関数をつないで系全体の波動関数を3. 磁気的並進演算子を用いた定式化
ベクトルポテンシャルを対称ゲージ A = B/2(−y, x, 0) として、原点近傍でシュレディンガー方程式 (1) を考える。原点近傍では A ∼ 0 であるから、原点近傍での波動関数は磁場がかかっていない場合のものに非常に近いと考えられる。
ここで、磁気的並進演算子()iT (R) ≡ exp
R · (p + qA) ,h̄表したことになっている。形状関数に付く位相因子は局所的ゲージから元のゲージに戻すための因子である。
簡単のためにノイマン境界条件の部分を省略すると汎関数は次のように与えられる。
∫ [1L ={(p − qA) ψ}∗ · {(p − qA) ψ}
2m
V− ψ ∗ (E − V (r)) ψ] dr,(10)を定義すると、この演算子は第一項[(p − qA)2 , T (R)] = 0∫ {}∗q1 ∑ ∗ei 2h̄ B ·(r j ×r ) pNj (r)ψj2mVj,k}{ q· ei 2h̄ B ·(r k ×r ) pNk (r) drψk∫qh̄2 ∑ ∗−ψjei 2h̄ B ·{(r k −r j )×r }
2mvjk

=(11)

および=q

T (R)V (r)ϕ(r)e−i 2h̄ B ·(R×r ) V (r + R)ϕ(r + R)==

V (r + R)T (R)ϕ(r)j,k× (∇Nj (r)) · (∇Nk (r)) drψk ,}∗∑ ∗ ∫ { i q B ·(r ×r )jE
ψje 2h̄Nj (r)(12)
第二項=なる関係を満たす。ただし、これらの関係式は対称ゲージのj,kときのみ成り立つ。ここで、[ , ] は [A, B] ≡ AB − BA で定==ET (R)ψ(r)∑=j,kとなる。また、∑=}

qi 2h̄ B ·(r k ×r )Vl ψj∗∫lNk (r) drψkei 2h̄ B ·{(r k −r j )×r }qvjk×Nj (r)Nl (r)Nk (r)drψk .

(14)であり、ψ(r + R) は r = −R においてベクトルポテンシャqV{j,k,lqT (R)ψ(r) = e−i 2h̄ B ·(R×r ) ψ(r + R)

ei 2h̄ B ·{(r k −r j )×r }×Nj (r)Nk (r)drψk ,()∫ {}∗ ∑qi 2h̄ B ·(r j ×r )∗eψjNj (r)Vl Nl (r)× e(13)}Nk (r) drψkvjkj,k第三項

qi 2h̄ B ·(r k ×r )× e∑ ∗∫Eψj

義される交換子である。
式 (1) の両辺に T (R) を演算すると[]
1T (R)(p − qA)2 + V (r) ψ(r)
2m
[]
1
=
(p − qA)2 + V (r + R) T (R)ψ(r)
2m

V{積分∫

dr は r j を r k を２頂点にもつ有限要素内の体積vjk

ルが 0 になるようなゲージ、即ち、A = 21 B × (r + R) の場積分であるから、r = r k + ξ により定義される新しい積分

合の r = −R 近傍の波動関数を表す。

変数 ξ を導入すると、行列要素 Kjk は∫qh̄2Kjk = −ei 2h̄ B ·{(r k −r j )×(r k +ξ)}
2m vjk

ここで、形状関数 Nj (r) は点 r j の近傍でのみ値を持つから、ψj Nj (r) は ψ(r − r j ) をよく近似するから、系全体の波動関数を

ψ(r) =∑qi 2h̄ B ·(r j ×r )

ψj e

Nj (r)(15)=

qi 2h̄

ej∫==≈を得る。従って、波動関数を式 (15) のように展開すると、近似的にベクトルポテンシャルを考える必要がなくなる。これ

qvjk

この波動関数にベクトルポテンシャルを含む運動エネル1(p − qA(r))2 ψ(r)2m∑q1
(p − qA(r))2
ψj ei 2h̄ B ·(r j ×r ) Nj (r)2mj∑
qi 2h̄ B ·(r j ×r ) 1
ψj e{p − q (A(r) − A(r j ))}2 Nj (r)2mj∑q1 2
ψj ei 2h̄ B ·(r j ×r )p Nj (r)(16)2m
jei 2h̄ B ·{(r k −r j )×ξ} (∇Nj (ξ)) · (∇Nk (ξ)) dξ×と表す。
ギー演算子を演算すると

× (∇Nj (ξ)) · (∇Nk (ξ)) dξ()
2B ·(r k ×r j ) − h̄
2mqi 2h̄()∫2B ·(r k ×r j ) − h̄(∇Nj ) · (∇Nk ) dξ
2mvjk

≈e=ei 2h̄ B ·(r k ×r j ) K0jkq(17)と表せる。ここで、磁場が印加されていない場合の行列要素を K0jk と表した。２行目から３行目への変形において、有限要素を貫く磁束は磁束量子 h/q に比べ十分小さいとして積分内の位相因子を 1 と近似した。

Table 1

同様に、∫
Mjkei 2h̄ B ·{(r k −r j )×r }q≡

Variations in 3 lowest eigen wavenumbers kd for 3

types of gauge.

vjk

gauge

k0 dref.(11) TABLE I

15.079220411952475

5.070

25.079220411952162

5.070

35.079220411952306

5.071

gauge

k1 dref.(11) TABLE I

15.2358151469960035

5.224

磁場中の電子に対する有限要素法は近似的に磁場が印加25.235815146996086

5.223

されていない場合の行列に所謂 Peierls 位相因子をかけるだ

3

5.235815146996177

5.228

けで定式化できることが分かる。これは文献 (12) で Peierls

gauge

k2 dref.(11) TABLE I

が導出した磁場中電子の強束縛近似した Schrödinger 方程式

15.5357147392872115

5.520（差分近似と同等 (1) ）のアナロジーとなっている。

2

5.535714739287768

5.520

3

5.535714739287288

5.528

× {ENj (r)Nk (r)()}∑−Nj (r)V (r k )Nl (r) Nk (r) drlqei 2h̄ B ·(r k ×r j ) M0jk

≈

(18)

を得る。

ここで、式 (17)、(18) の位相因子は A = B/2(−y, x, 0) を用いて

[ ∫ rj]qqei 2h̄ B ·(r k ×r j ) = exp i ↗ A(r) · drh̄ r k(19)と書ける。ここで、積分記号の矢印は点 r k から点 r j へ直線で結ぶ経路に沿って積分することを意味する (6) 。右辺の表現を用いると、rotA = (0, 0, B) を満たす任意のゲージに対して=

Mjk=の (d))、B(−y, 0, 0) (gauge 2、(a))、B(0, x, 0) (gauge 3、(c))
について計算を行った。N = 20 の場合（要素数 1600）に、低い方から３つのエネルギー準位に対応する無次元化した波数 k0 d, k1 d, k2 d の値を文献 (11) の要素数が 1250 の場合の結果と比較したものを Table.1 に示す。

[Kjk]
∫q rjexp i ↗ A(r) · dr K0jkh̄ r k]
[ ∫ rjqexp i ↗ A(r) · dr M0jkh̄ r kて、３種類のゲージ A = B(−y, x, 0)/2 (gauge 1、文献 (11)(20)

Table.1 より、本論文で提案した手法では、異なるゲージを用いて計算したとしても、計算結果は小数点第 13 位まで(21)（計算精度で）一致している。これは他の分割数に対しても同様である。

が成り立つ。

次に、Fig.7 に gauge 1 の場合の最低 3 準位の波数の要素有限要素内の体積積分中の位相因子を 1 とする近似の影響

数依存性を示す。要素数に対して速やかに収束していることで数値解の精度が落ちるが、有限要素を小さくすればその近似精度はシステマティックに上がる。この近似が重大な影響6.0を及ぼす場合は、有限要素内の積分において、位相因子をポ5.8

テンシャル同様に接点間で形状関数を用いて補間し、積分を5.6kd実行すればよい。その場合の行列 K、M は()qqh̄2 ∑ i 2h̄
Kjk = ei 2h̄ B ·(r k ×r j ) −e B ·{(r k −r j )×r l }2ml∫×(∇Nj (ξ)) · (∇Nk (ξ)) Nl (ξ)dξ
(22)

5.4
5.2
5.0vjkMjk=∑ i q B · (r −r )×rqn}
ei 2h̄ B ·(r k ×r j )e 2h̄ { k j1000

2000

3000

4000

5000

6000

Number of Elementsn∫

Fig. 7{ENj (r)Nk (r)Nn (r)×vjk−0∑

number of elements}

V (r k )Nj (r)Nl (r)Nk (r)Nn (r) dr

Variations in 3 lowest eigen wavenumbers kd for the(23)lと与えられる。

4. 数値計算例が分かる。残る誤差は積分内の位相因子の近似の影響と考えられる。
磁場の強い場合の例として、文献 (11) の Fig.4 と同じ B̃ =100 における 9 番目のエネルギー準位（kd = 10.0247）でのここでは、式 (17), (18) を用いた計算例を示す。提案手法

確率密度のようすを Fig.8 に示す。近似的な式 (20), (21) をの近似精度、ゲージ依存性を確認するために、文献 (11) にお用いた解析であるにもかかわらず、磁場が強い場合にもよくいて用いられた、Fig.1 の長方形量子ドット内の固有値問題一致している。

を解析する。長さ d を N 等分するように有限要素に分割し次に、系の異方性が大きな場合の例として、Fig.2 から左より、本手法はどのようなゲージをとっても計算精度内で全d2く同じ解を与えることを示した。本手法では非常に簡便に磁場の効果を取り入れることができ、実用的であり、また、系0の形状に関係なく、伝導問題においてはリードの処理が容易であるゲージを選ぶなど優位性があることが示された。今後、本手法を用いた様々な分野での磁場中電子波の有限要素-d2-dFig. 80d解析が期待される。

参考文献

Density plot of the probability density of the 9thenergy level (kd = 10.0247) for B̃ = 100.

(1) Supriyo Datta :Electronic Transport in MesoscopicSystems, Cambridge University Press, (1997).
右の導波路を取り去った縦 8d、横 d の長方形量子ドットの固(2) R. L. Schult, H. W. Wyld and D. G. Ravenhall : Quan-有値問題の解析例を示す。この場合にも計算結果は計算精度

tum Hall Eﬀect and General Narrow-Wire Circuits,の範囲で全くゲージには依存しないが、Fig.4 においては最Phys. Rev. B, 41 (1990), pp. 12760–12780.

も不利と考えられたゲージ A = (−By, 0, 0) のときの結果を(3) Craig S. Lent and David J. Kirkner :The quantum示す。B̃ = 20 のときに、伝導問題で用いている kd = 10 にtransmitting boundary method, J. Appl. Phys., 67

最も近い２つのエネルギー固有状態 kd = 10.0159, 9.95098(1990), pp. 6353–6359.の確率密度を Fig.9 に示す。確率密度は長方形内全体に広が(a)
4d(b)
4d

(4) L. Ramdas Ram-Mohan : Finite Element and Boundary
Element Applications in Quantum Mechanics, Oxford
University Press, (2002).
(5) H. R. Frohne, M. J. Mc Lennan and S. Datta : An Eﬃcient Method for the Analysis of Electron Waveguides,
J. Appl. Phys., 66 (1989), pp. 2699–2705.
(6) Tsuyoshi Ueta : Green’s Function of a Charged Particle
in Magnetic Fields, J. Phys. Soc. Jpn., 61 (1992),
pp. 4314–4324.
(7) Tsuyoshi Ueta : Boundary Element Method for Elec-00tron Waves in Uniform Magnetic Fields, Engineering
Analysis with Boundary Elements, 17 (1996), pp. 69–74.
(8) Y. A. Bychkov, E. I. Rashba : Oscillatory eﬀects and themagnetic susceptibility of carriers in inversion layers, J.
Phys. C: Solid State Phys., 17 (1984), pp. 6039–6045.
(9) Yongjiang Wang, Jian Wang and Hong Guo : Magne--4d-4d-d2 d2-d2 d2Fig. 9 Density Plots of the probability density for two eigen
states kd = 10.0159 (a) and kd = 9.95098 (b) within the8d × d rectangle.
り、Fig.4 のようになることはない。伝導問題では kd = 10であるが、Fig.6 の場合には kd = 9.95098 の状態に対応した状態を励起していることが分かる。

toconductance of a stadium-shaped quantum dot: Aﬁnite-element-method approach,Phys. Rev. B,49(1994), pp. 1928–1934.
(10) Manhua Leng and Craig S. Lent :Quantum trans-mitting boundary method in a magnetic ﬁeld, J. Appl.
Phys., 76 (1994), pp. 2240–2249.
(11) Koichi Hirayama, Yasuhiro Honma, Yoshio Hayashi
and Masanori Koshiba : A Novel Finite-Element Formulation for the Analysis of the Energy Levels of a5. むすび
磁気的並進演算子とベクトルポテンシャルを含むゲージ不変な運動エネルギー演算子との交換関係を考えることによQuantum Cavity in a Magnetic Field, IEEE PHOTONICS TECHNOLOGY LETTERS, 10 (1998), pp. 1359–1361.

り、形状関数がピークを持つ節点でベクトルポテンシャルが(12) R. E. Peierls : Zur Theorie des Diamagnetismus von

0 になるような局所的なゲージ変換を行った形状関数を用い

Leitungselektronen, Z. Phys., 80 (1933), pp. 763–791.て波動関数を展開すると、位相因子だけで磁場の効果を取り込む有限要素法を定式化した。固有値問題に適用することに(13) E. Brown:

Bloch Electrons in a Uniform MagneticField, Phys. Rev., 133 (1964), pp. A1038–A1044.