複素周波数領域有限差分法によるプラズモニック導波路の設計及び特性検証

複素周波数領域有限差分法によるプラズモニック導波路の設計及び

特性検証

呉

迪

†a)

_浜島

_功

†

_{井上修一郎}

††

_{大貫進一郎}

†b)

Design and Characteristic Verification of Plasmonic Waveguides

Using Finite-Difference Complex-Frequency-Domain Method

Di WU

†a)

, Kou HAMASHIMA

†

, Shuichiro INOUE

††

, and Shinichiro OHNUKI

†b)

あらまし 近年，表面プラズモンを用いて，光エネルギーを局所領域に集中して伝搬させるプラズモニック導 波路が注目を集めている．本論文では，著者らが近年開発した，導波路構造の電磁界解析に適した複素周波数領 域有限差分法により，誘電体装荷グレーティング構造を用いたプラズモニック導波路，量子ウォークの検証を目 的とした金属ストライプ型プラズモニック導波路アレーの解析を行い，その励振及び伝搬特性を検証する． キーワード 複素周波数領域有限差分法，高速逆ラプラス変換，プラズモニック導波路，グレーティング構造， 金属ストライプ，量子ウォーク

1.

ま え が き

近年，集積回路の更なる高速化，小型化及び省電力

化を目指して光回路の研究が精力的に進められてい

る

[1], [2]

．従来の光回路に用いられる光導波路は誘電

体材料で作られるものが一般的であるが，光閉じ込

め効果には限界が存在する．そのため，光の回折限界

を超えた小型化には困難が生じた．そこで，金属ナ

ノ構造体における表面プラズモンを利用したプラズ

モニック導波路が検討されている

[3]

．表面プラズモ

ンは光と結合した自由電子の集団振動であり，金属ナ

ノ構造体と誘電体の界面を光と同じ周波数で伝搬す

る．このため，回折限界を超えた小型化，高速な信号

処理が実現できる．本論文では，著者らが近年開発し

た，導波路構造の電磁界解析に適した複素周波数領域

有限差分法

(Finite-Difference

Complex-Frequency-†_{日本大学理工学部電気工学科，東京都}

Department of Electrical Engineering, College of Science and Technology, Nihon University, 1–8–14 Kanda-Surugadai, Chiyoda-ku, Tokyo, 101–8308 Japan

††_{日本大学大学院理工学研究科量子理工学専攻，東京都}

Graduate School of Quantum Science and Technology, Nihon University, 1–8–14 Kanda-Surugadai, Chiyoda-ku, Tokyo, 101–8308 Japan a) E-mail: [email protected] b) E-mail: [email protected]

Domain: FDCFD) [4], [5]

を用いて高効率かつ高精度

にプラズモニック導波路を設計し，その励振及び伝搬

特性を検証する．

プラズモニック導波路の設計及びその特性解析に必須

となる電磁界数値解析は，これまで多くの手法が開発・

検討されている

[6]

∼

[13]

．中でも，時間領域有限差分法

(Finite-Difference Time-Domain: FDTD) [9]

∼

[12]

は商用シミュレータにも搭載され，電磁界時間応答解

析の代表的な手法として，近年最も広く利用される手

法の一つである．しかしながら，空間刻みに依存する

安定条件

[10]

より，プラズモニック導波路のように，

ナノメートルの断面構造に対し，長さがミリメートル

からセンチメートルとなるマルチスケールの構造を解

析する際は，計算コストが膨大になる．無条件安定な

陰的

FDTD

法

[14], [15]

も開発されているが，計算速

度の向上を実現できる一方で，時間間隔を大きくする

につれて計算精度に劣化が生ずる．

周波数領域差分法

(Finite-Difference

Frequency-Domain: FDFD)

は薄膜

[16]

や導波路構造

[17]

など

に対して周波数応答を効率的に求められる手法として

知られる．解析空間は

Yee

格子

[9]

で離散化しているた

め，解析モデルの作成が容易であり，他の差分法とのモ

デル共有が可能などの利点がある

[5]

．また，

subpixel

smoothing

などの手法を用いることによって，曲面構

造を高精度にモデル化でき

[18], [19]

，定式化の改良に

よる異方性

[20]

や非線形媒質

[21]

の解析も比較的容

易に実現できる．著者らは周波数領域差分法を複素周

波数領域

(Complex-Frequency-Domain: CFD)

に拡

張し，高速逆ラプラス変換法

(Fast Inverse Laplace

Transform: FILT)

と併用する

FDCFD-FILT

法を近

年開発した

[4]

．本手法は電磁界の定常応答と過渡応

答のいずれも効率良く解析することが可能である．更

に，任意時刻における応答波形を独立に計算できるた

め，計算精度が時間間隔に依存しない特徴や並列計算

に適している利点がある．

本論文では

FDCFD-FILT

法を用いた電磁デバイ

スの設計及びその検証例として，誘電体装荷グレー

ティング

(Dielectric Loaded Gratings: DLGs)

を用

いたプラズモニック導波路及び，量子ウォーク実現に

向けた金属ストライプ型表面プラズモニック導波路ア

レー

(Plasmonic Waveguide Array: PWA)

を紹介す

る．

DLGs

を用いたプラズモニック導波路は光ファイ

バ間の通信を効率的に行えるデバイスとして提案され

ている

[22]

．構造が単純なため，製造は比較的容易で

ある．ここでは，

FDCFD-FILT

法を用いて，表面プ

ラズモンの励振及び伝搬特性を検証する．量子ウォー

ク実現に向けた金属ストライプ型

PWA

では，電磁

界分布を古典的に計算し，量子ウォークへの実用性を

検証する．量子ウォークは量子力学的な粒子ランダム

ウォークであり，その原理検証は量子情報処理への応

用を目指して

2010

年頃から世界各国で精力的に行わ

れている

[23]

∼

[25]

．

2.

解 析 手 法

2. 1

複素周波数領域有限差分法

複素周波数領域有限差分法

(FDCFD)

では，複素周

波数領域

(s :=

σ + jω)

での

Maxwell

方程式を扱う．

入射界と散乱界を分離した

Maxwell

方程式は次式で

表現できる．



∇ ×

sμ

−sε

∇ ×



E

s

H

s



=



s(μ − μ

0

)

H

i

s(ε − ε

0

)

E

i



(1)

ここで，時間因子は

e

st

，

E

i

と

H

i

は入射電界及び

入射磁界，

E

s

と

H

s

は線形方程式の未知ベクトルで

ある散乱電磁界を表す．また，

ε

と

μ

はそれぞれ散乱

体の誘電率と透磁率である．

導波路に用いる金属の周波数分散は，次式で示す複

素周波数領域の

Lorentz-Drude

モデル

[26]

より，複

素比誘電率を

ε

r

(:=

ε/ε

0

)

考慮する．

ε

r

(

s) = 1 +

f

0

ω

2 p

s

2

₊

sΓ

0

+

N



j=1

f

j

ω

p2

s

2

₊

sΓ

_j

₊

ω

2 j

, (2)

ここで，

ω

p

はプラズマ周波数であり，

ω

j

，

f

j

，

Γ

j

は

それぞれ，各モードにおける振動子の角周波数，強度，

ダンピング係数を表す．

式

(1)

を

Yee

格子で差分近似することにより，以下

の線形方程式が得られる

[4]

．

⎡

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎣

sμ

−D

z

D

y

sμ

D

z

−D

x

sμ −D

y

D

x

−D

z

D

y

−sε

D

z

−D

x

−sε

−D

y

D

x

−sε

⎤

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎦

⎡

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎣

H

s x

H

s y

H

s z

E

s x

E

s y

E

s z

⎤

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎦

=

⎡

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎣

E

i x

E

i y

E

i z

H

i x

H

i y

H

i z

⎤

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎦

(3)

ここで，

E

i

=

s(ε − ε

0

)

E

i

，

H

i

=

s(μ − μ

0

)

H

i

．

また，

D

x

，

D

y

，

D

z

は以下のような

x

，

y

，

z

方向の

差分作用素を示す．

D

x

E

ws

=

E

s w

|

i+1,j,k

− E

ws

|

i,j,k

Δ

x

,

D

y

E

ws

=

E

s w

|

i,j+1,k

− E

ws

|

i,j,k

Δ

y

,

D

z

E

ws

=

E

s w

|

i,j,k+1

− E

ws

|

i,j,k

Δ

z

,

D

x

H

ws

=

H

s w

|

i+0.5,j,k

− H

ws

|

i−0.5j,k

Δ

x

,

D

y

H

ws

=

H

s w

|

i,j+0.5,k

− H

ws

|

i,j−0.5,k

Δ

y

,

D

z

H

ws

=

H

s w

|

i,j,k+0.5

− H

ws

|

i,j,k−0.5

Δ

z

,

式

(3)

に 示 す 線 形 方 程 式 を 解 け ば 未 知 電 磁 界 で

ある

E

s

と

H

s

は求まるが，

|E

s

|/|H

s

| =



μ/ε

の 関 係 が あ る た め 電 界 と 磁 界 の オ ー ダ は 異 な る ．

E

s

₌



_μ/εE

s

_{のような規格化を行うことで数値}

的な条件を改善することは可能であるが，媒質境界面

の処理は複雑になるため，次式に示す波動方程式を考

慮する．

∇ × (−μ

−1

_{∇ × E}

s

₎

_{− s}

2

_εE

s

=

sE

i

−

∇ × μ

−1

H

i

(4

a)

∇ × (ε

−1

_{∇ × H}

s

_{) +}

_s

2

_εμH

s

=

sH

i

+

∇ × ε

−1

E

i

(4

b)

本論文では，金属表面における表面プラズモンを求め

るため，式

(4

a)

による電界

E

s

を評価する．金属の

透磁率

μ

を真空中の透磁率

μ

0

で近似できる場合は，

次式に変形できる．

∇ × (∇ × E

s

₎

_{− s}

2

_εμE

s

₌

_sμE

i

_,

₍₅₎

式

(5)

における未知ベクトルは散乱電界

E

s

のみと

なる．

FDCFD

法で必要な吸収境界条件に，

CPML

(Con-volutional Perfectly Matched Layer)

を適用する場

合は，ナブラ演算子

∇

を

∇

c

に置き換える．すな

わち，

∇

c

× (

∇

c

× E

s

)

− s

2

εμE

s

=

sμE

i

(6

a)

∇

c

=

_c

∂

x

∂x

a

x

+

∂

c

y

∂y

a

y

+

∂

c

z

∂z

a

z

(6

b)

c

w

=

κ

w

+

σ

w

α

w

+

sε

0

, w = x, y, z

(7)

ここで，

a

_x

，

a

_y

，

a

_z

は単位ベクトル，

κ

w

，

σ

w

，

α

w

は

CPML

のパラメータである

[27]

．

式

(5)

を離散化することで次式に示す線形方程式を

得る．

⎡

⎢

⎣

D

2 yz

D

x

D

y

D

x

D

z

D

x

D

y

D

xz2

D

y

D

z

D

x

D

z

D

y

D

z

D

xy2

⎤

⎥

⎦

⎡

⎢

⎣

E

s x

E

s y

E

s z

⎤

⎥

⎦=

⎡

⎢

⎣

sμE

i x

sμE

i y

sμE

i z

⎤

⎥

⎦, (8)

ここで，方程式左辺の係数行列

D

yz2

，

D

x

D

y

，

D

x

D

z

などは

2

階の差分作用素であり，以下の式で表すこと

ができる．

D

2 yz

E

xs

=

1

c

y

|

i,j,k

Δ

y

2



E

s x

|

i,j+1,k

− E

xs

|

i,j,k

c

y

|

i,j+0.5,k

− E

s x

|

i,j,k

− E

xs

|

i,j−1,k

c

y

|

_i,j−0.5,k

+

1

c

z

|

_i,j,k

Δ

z

2



E

s x

|

i,j,k+1

− E

xs

|

i,j,k

c

z

|

i,j,k+0.5

− E

xs

|

i,j,k

− E

xs

|

i,j,k−1

c

z

|

i,j,k−0.5

− s

2

_ε|

i,j,k

μ|

i,j,k

E

xs

|

i,j,k

,

D

x

D

y

E

sx

=

E

s

x

|

i+1,j,k

+

E

xs

|

i,j−1,k

− E

sx

|

i+1,j−1,k

− E

xs

|

i,j,k

c

x(i,j,k)

c

y(i,j+0.5,k)

Δ

xΔy

,

D

x

D

z

E

xs

=

E

s

x

|

i+1,j,k

+

E

xs

|

i,j,k−1

− E

xs

|

i+1,j,k−1

− E

xs

|

i,j,k

c

x(i,j,k)

c

z(i,j,k+0.5)

Δ

xΔz

,

ここで，

Δ

x

，

Δ

y

，

Δ

z

は

Yee

格子各辺の長さ，

i

，

j

，

k

はそれぞれ

x

，

y

，

z

方向の座標を表す．全解析領域

にこの差分近似を適用すると，次の線形方程式に帰着

できる

Ax = b

(9)

ここで

A

は係数行列，

x

は求める未知電界，

b

は入

射電界を表す既知ベクトルである．係数行列

A

は有

限差分スキームにより作られた疎行列となる．この線

形方程式を解くために，直接法や大規模行列に対して

は反復法が利用される

[28]

．

2. 2

高速逆ラプラス変換

著者らが開発した電磁界解析法では，複素周波数領

域で求めた電磁界に高速逆ラプラス変換法

(FILT)

を

併用して時間応答解析を行う．

FILT

法は細野

[29]

に

より開発された，数値的に逆ラプラス変換する手法で

ある．複素平面で有限の特異点を計算することにより，

任意時刻での応答波形を求める．これまでに，線路

解析

[30]

，電磁波散乱解析

[31]

から，地中レーダ

[32]

，

超高速磁気記録用アンテナの設計

[33]

などの応用ま

で，報告例がある．本論文では，複素周波数領域有限

差分法で求めた電磁界を，高速逆ラプラス変換を用い

て時間領域に変換する．開発手法は安定条件の制限が

ないため，任意時刻における応答波形を独立に計算す

ることが可能であり，並列計算を容易に実現できる．

FILT

法では，指数関数を双曲線関数で近似する．

更に，留数定理を用いて，

Bromwich

積分を以下の形

に書き換える．

f(t) ≈ f

ec

(

t, a) = e

a

t

∞



n=1

F

n

,

(10)

ここで，

F

n

= (

−1)

n

Im[

F (s)],

(11

a)

s = a

+

j(n − 0.5)π

t

(11

b)

F (s)

は時間関数

f(t)

の像関数を表す，係数

a

は

FILT

法の近似パラメータであり，おおよそ

10

−a

の精度が

保証される

[29]

．

数値解析を行う際，無限級数を有限で打ち切る必要

がある．級数に対する収束を加速させるために，オイ

ラー変換を適用する．実際の逆変換を行う際は以下の

式を用いる．

図 1 高速逆ラプラス変換の並列アルゴリズム

f

ec

(

t, a) = e

a

t



_k



n=1

F

n

+

1

A

p0 p



q=1

A

pq

F

k+q



,

(12)

ここで，

A

pp

= 1

, A

p0

= 2

p

, A

pq

=

A

pq−1

−

_{q! (p+1−q)!}

(

p+1)!

(13)

k

は無限級数の打ち切り項数，

p

は収束を加速するた

めのオイラー変換の項数である．

図

1

に高速逆ラプラス変換のアルゴリズムを示す．

複素周波数領域で求めた電磁界を式

(12)

に代入する

ことにより，時間領域の近似解である

f

ec

(

t, a)

が求ま

る．また，各観測時間及び級数の各項における計算は

全て独立に行えるため，並列計算を容易に実現できる．

3.

数 値 結 果

3. 1 FDCFD-FILT

の計算精度及び計算速度の

検証

複素周波数領域有限差分（

FDCFD

）法の計算精度

を検証するために，図

2

に示す金円柱による電磁波の

散乱解析を行う．図

3

は金円柱近傍における電界強度

の波長応答である．ここで，円柱の半径

r

は

10nm

，

観測点は円柱表面から

5nm

離れた点とした．

FDCFD

法による解析結果は，数学的に正しい厳密解と

3

桁程

度一致した．また，

FDCFD

法と

FILT

法とを併用し，

電磁界の時間応答解析を行った．図

4

に変調

Gaussian

パルスを入射した際の金属円柱近傍における電界の時

間応答を示す．

FDCFD-FILT

法と

FDTD

法の解析

図 2 金属円柱の解析モデル 図 3 金属円柱近傍における電界強度の波長応答特性 図 4 変調 Gaussian パルスを入射した際の電界時間応答

結果を比較すると，相対誤差は

0.1%

以下となること

を確認した．また，図

5

は

15fs

から

20fs

の応答を拡

大したもので，時間間隔

Δt

を変えて解析結果を比較

した．

FDCFD-FILT

法は安定条件の制限がないため，

Δt

を変化させた二つの結果は計算精度の範囲で完全

に一致している．

次に，

FDCFD-FILT

法の計算時間を検証するため，

表

1

に解析時間を示す．単一計算ノードの場合，時間

間隔が

FDTD

法における安定条件の約

1500

倍であ

る

Δt = 0.3fs

より大きいと，本手法の優位性が確認

できる．また，本手法では完全並列計算が実現できる

ため，並列化効率はほぼ

100%

となり，計算時間は使

用する計算ノード数に反比例することも確認できる．

図 5 変調 Gaussian パルス入射による応答の拡大図 表 1 時間間隔Δt を変化させた際の解析時間 図 6 誘電体装荷グレーティングを用いたプラズモニック 導波路

3. 2 DLGs

を用いたプラズモニック導波路

図

6

に誘電体装荷グレーティング

(Dielectric Loaded

Gratings: DLGs)

を用いたプラズモニック導波路の

モデルを示す．シリカクラッドの上に薄い金薄膜を設

置し，更にその上にシリカの入出力用の

DLGs

や導

波路クラッドを装着する構造である．入力

DLGs

に

Gaussian

ビームを入射すると，入射光は

DLGs

と結

合することにより表面プラズモンを励振し，導波路に

沿って伝搬される．また，励振された表面プラズモン

が出力

DLGs

に到達すると，出力光として放射される．

従来の時間領域解析法で

DLGs

による出力光の放射

を解析する場合，ビーム入射の初期状態から電磁界の

逐次計算を行うが，

FDCFD-FILT

法では，表面プラ

図 7 入力誘電体装荷グレーティング構造 図 8 グレーティング周期λg及び空気層幅Waの反射率 R と結合率 C 特性 (tAu= 20nm)

ズモンが出力

DLGs

に到達した時点など，任意時刻で

の解析が可能となるため，より効率的に特性を検証で

きる．

図

7

は入力

DLGs

構造の拡大図であり，

z

方向に

は一様な断面を有するモデルを想定した．ここで，

λ

g

は

DLGs

周期，

t

Au

は金薄膜の厚み，

W

a

は空気層の

幅を示す．

DLGs

及び金属薄膜上部クラッドの高さ

h

は入射光の半波長に固定し，

DLGs

の長さは

10

μm

と

した．また，入射光として波長

1.55

μm

の

Gaussian

ビームを仮定する．入力

DLGs

の反射率

R

及び結合

率

C

を評価するために，

DLGs

から

+

y

方向に

3

μm

離れた観測面

A (

Σ

A

)

及び

+

x

方向に

1

μm

離れた観

測面

B (

Σ

B

)

を通過するエネルギー比を次式により計

算する．

R =



ΣA

(

E

s x

)

2

· Δx



ΣA

(

E

i x

)

2

· Δx

,

(14

a)

C =



ΣB

(

E

tot y

)

2

· Δy



ΣA

(

E

i x

)

2

· Δx

,

(14

b)

ここで，

E

ytot

=

E

sy

+

E

yi

．

まず，入力

DLGs

による表面プラズモンの励振及び

伝搬効率を評価する．図

8

に

DLGs

周期

λ

g

及び空気

層

W

a

を変化させた際の，

a)

観測面

A

で計算した反

射率

R

，

b)

入射光と入力

DLGs

との結合率

C

を示す．

ここで，金薄膜の厚みは

20nm

と固定した．

DLGs

の

反射率は

0.30

で最小となり，このときの結合率は

0.15

図 9 グレーティング周期λg及び空気層幅Waの反射率 R と結合率 C 特性 (tAu= 100nm) 図 10 金薄膜厚みtAuを変化させた際のクラッド内に伝 搬するエネルギーの伝搬損失L 図 11 プラズモニック導波路近傍における電界強度分布

となる．

図

9

に 金 薄 膜 の 厚 み

t

Au

を

100nm

に 変 化 し た

際の解析結果を示す．観測面

A

における反射率は，

λ

g

= 1

.20nm

，

W

a

= 0

.48nm

で

R = 0.17

と最小に

なる．このときの結合率は

0.26

となる．

図

10

に金薄膜の厚み

t

Au

を

20nm

から

100nm

に変

化させた際の，クラッド内を伝搬するエネルギーの伝

搬損失

L

を示す．厚みに対する

DLGs

のパラメータは

結合率が最大となるように設定した．伝搬距離

100

μm

に対するエネルギーの伝搬損失は，

t

Au

= 100nm

の

場合

−2.5dB

程度となることが確認できる．

図

11

に金薄膜の厚み

t

Au

= 100nm

，

DLGs

周期

λ

g

= 1

.20nm

，空気層

W

a

= 0

.48nm

の際の

DLGs

近傍における電界分布を示す．金薄膜表面のみに電界

が強く集中しながら，

+

x

方向への伝搬が確認できる．

図 12 表面プラズモンが出力 DLGs に到達した際の電界 分布 図 13 導波路上空における電界強度分布

次に，出力

DLGs

を設置した場合の解析結果を示

す．入力側に

50

μm

離れた点から，長さ

10

μm

の出力

DLGs

を設置した．図

12

に伝搬された表面プラズモ

ンが出力

DLGs

に到達した時点における電界の空間分

布を示す．出力

DLGs

により，クラッド内に閉じ込め

られた光が導波路上方に放出されることを確認できる．

図

13

は金薄膜表面から

x

軸に対して平行に

2

∼

3

μm

離れた面での電界強度分布である．ここで，

−10μm

か

ら

0

μm

に入力

DLGs

，

50

μm

から

60

μm

に出力

DLGs

を設置した．入出力

DLGs

が設置された部分でのみ電

界は強くなる．特に出力部から放射される電界強度は

約

0.5V/m

に対し，クラッド部の電界強度は

0.05V/m

以下となる．また，次式により定義するエネルギーの

入出力比

O

は

12.8%

であることを確認した．

O =



ΣC

(

E

tot x

)

2

· Δx



ΣA

(

E

i x

)

2

· Δx

,

(15)

ここで，

Σ

C

は出力

DLGs

から

+

y

方向に

3

μm

離れ

た観測面

C

である．

3. 3

金属ストライプ型

PWA

量子ウォークは量子情報処理の要素技術であり，近

年盛んに研究されている

[23]

∼

[25]

．図

14

に示す金属

ストライプ型表面プラズモニック導波路アレー

(Plas-monic Waveguide Array: PWA)

は金属ナノ構造をも

図 14 金属ストライプ型 PWA 図 15 金属ストライプ間における表面プラズモンの遷移 図 16 PWAにおける表面プラズモンの遷移

期待されている．ここでは，光ファイバの端部よりス

トライプに光を入射する際の表面プラズモンの遷移及

び，電界分布を解析する．本論文のモデルは

y

軸方向

に一様で，

x

軸方向の金属ストライプの厚みは

20nm

とする．

図

15

に間隔

1.5

μm

をもつ二枚のストライプ間にお

ける電界の強度分布を示す．金属ストライプ

A

に直接

光を入射し，表面プラズモンを励振及び伝搬させる．

図より，プラズモンは

2

枚のストライプ間の遷移を繰

り返しながら伝搬することが確認できる．

次に，量子ウォークに応用を目的とした

31

枚のス

トライプに対して伝搬特性を検証する．中心のストラ

イプに光を入射し，その上下に各

15

枚のストライプ

を設置した際の

PWA

近傍における電界分布を図

16

に示す．中心のストライプに表面プラズモンが励振さ

れ，その後，上下のストライプにプラズモンが遷移し

ながら伝搬する．

図

17

に各観測面における最大値で規格化した電界

の強度分布を示す．表面プラズモンの強度が最大とな

る点は伝搬及び遷移するにつれて中央から

PWA

の両

端に向かって線型的に広がり，入射位置から離れたス

図 17 観測面における電界強度分布

トライプに電界が最も集中するような結果が得られた．

電磁界解析で計算された電界強度は，単一光子が観測

される確率分布に対応する．古典的なランダムウォー

クでは，入力を中心として

Gauss

分布的に粒子の出力

が広がるのと異なり，量子ウォークでは光子が線型的

に広がり，二つの端点に局在する

[34]

．今回の解析結

果はこの量子ウォークの確率分布特性を明確に示して

おり，

PWA

を量子ウォークに応用する際の有用性が

検証できる．

4.

む す び

本論文では，著者らが近年開発した

FDCFD-FILT

法のアルゴリズムを紹介し，電磁界の散乱解析を行う

ことにより本手法の計算精度及び計算速度を検討した．

複素周波数領域の有限差分法（

FDCFD

）と高速逆ラ

プラス変換（

FILT

）を併用することにより任意時刻

における時間応答を独立に計算できる．また，

Yee

セ

ルを用いて差分近似しているため，解析モデルの作成

が容易で，他の差分法とのモデル共有が可能である．

次に，本開発法を用いて誘電体装荷グレーティング

を利用したプラズモニック導波路の解析を行い，表

面プラズモン伝搬効率とグレーティング周期長さ，空

気層及び金薄膜の厚みとの関係を明らかにした．設

計した構造は電界強度

1V/m

をもつ入射光に対して，

50

μm

離れた出力

DLGs

より

0.4V/m

程度の電界強

度が得られた．また，複数金属ストライプを用いたプ

ラズモニック導波路アレーのシミュレーションを行い，

量子ウォーク応用への実用性を示した．

謝辞 本研究を進めるに当たり，東京都立産業技術

センターの山口隆志氏及び本学大学院生の大西崚平氏

から多大な助言を賜り，厚く感謝を申し上げます．ま

た，本研究の一部は，科研費基盤

(C) 17K06401

，及

び日本大学理工学部プロジェクト研究助成金の援助を

受けて行われた．

文

献

[1] L.A. Coldren, S.W. Corzine, and M.L. Mashanovitch, Diode lasers and photonic integrated circuits, John Wiley & Sons, NY, USA, 2012.

[2] S.J. McNab, N. Moll, and Y.A. Vlasov, “Ultra-low loss photonic integrated circuit with membrane-type photonic crystal waveguides,” Opt. Express, vol.11, no.22, pp.2927–2939, 2003.

[3] S. Bozhevolnyi, Plasmonic nano guides and circuits, Jenny Stanford Publishing, NY, USA, 2009. [4] D. Wu, R. Ohnishi, R. Uemura, T. Yamaguchi,

and S. Ohnuki, “Finite-difference complex-frequency-domain method for optical and plasmonic analy-ses,” IEEE Photonics Technol. Lett., vol.30, no.11, pp.1024–1027, 2018.

[5] S. Ohnuki, R. Ohnishi, D. Wu, and T. Yamaguchi, “Time-division parallel FDTD algorithm,” IEEE Photonics Technol. Lett., vol.30, no.24, pp.2143– 2146, 2018.

[6] 電気学会編，計算電磁気学，培風館，東京，2003. [7] J.M. Jin, The finite element method in

electromag-netics, John Wiley & Sons, NY, USA, 2015. [8] M. Koshiba, Optical waveguide theory by the

fi-nite element method (Advances in opto-electronics), Kluwer Academic Pub, Berlin, Germany.

[9] K. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE Trans. Antennas Propag., vol.14, no.3, pp.302–307, May 1966.

[10] A. Taflove and S.C. Hagness, Computational elec-trodynamics, 2nd ed., Artech House, Norwood, MA, USA, 1995.

[11] 宇野 亨，FDTD 法による電磁界及びアンテナ解析，コ ロナ社，東京，1998.

[12] 橋本 修，実践 FDTD 時間領域差分法，森北出版，東京， 2006.

[13] W.C. Chew, J.M. Jin, E. Michielssen, and J. Song, Fast and efficient algorithms in computational, Artech House, MA, USA, 2001.

[14] T. Namiki, “A new FDTD algorithm based on alternating-direction implicit method,” IEEE Trans. Microw. Theory Tech., vol.47, no.10, pp.2003–2007, 1999.

[15] J. Shibayama, M. Muraki, J. Yamauchi, and H. Nakano, “Efficient implicit FDTD algorithm based on locally one-dimensional scheme,” Electron. Lett., vol.41, no.19, pp.1046–1047, 2005.

[16] W.E.I. Sha, W.C.H. Choy, and W.C. Chew, “A com-prehensive study for the plasmonic thin-film solar cell with periodic structure,” Opt. Express, vol.18, no.6, pp.5993–6007, 2010.

[17] F. Xu, Y. Zhang, W. Hong, K. Wu, and T.J. Cui, “Finite-difference frequency-domain algorithm for modeling guided-wave properties of substrate in-tegrated waveguide,” IEEE Trans. Microw. Theory

Tech., vol.51, no.11, pp.2221–2227, Nov. 2003. [18] A. Deinega and I. Valuev, “Subpixel smoothing for

conductive and dispersive media in the finitediffer-ence time-domain method,” Opt. Lett., vol.32, no.23, pp.3429–3431, 2007.

[19] T. Yamaguchi, M. Ebisawa, and S. Ohnuki, “Precise finite difference analysis of Lorentz force acting on metal nanoparticle irradiated with light,” Progress In Electromagnetics Research C, vol.73, pp.81–86, 2017. [20] J.A. Kong, Electromagnetic wave theory, EMW

Pub-lishing, MA, USA, 2008.

[21] N. Karasawa, “Computer simulations of nonlin-ear propagation of an optical pulse using a finite-difference in the frequency-domain method,” IEEE J. Quantum Electron., vol.38, no.6, pp.626–629, 2002. [22] B. Steinberger, A. Hohenau, H. Ditlbacher, A.L.

Stepanov, A. Drezet, F.R. Aussenegg, A. Leitner, and J.R. Krenn, “Dielectric stripes on gold as sur-face plasmon waveguides,” Appl. Phys. Lett., vol.88, no.9, 094104, 2006.

[23] A. Peruzzo, M. Lobino, et al., “Quantum walks of correlated photons,” Science, vol.329, no.5998, pp.1500–1503, 2010.

[24] M. Karski, L. F¨orster, et al., “Quantum walk in posi-tion space with single optically trapped atoms,” Sci-ence, vol.325, no.5937, pp.174–177, 2009.

[25] A. Crespi, R. Osellame, et al., “Anderson localiza-tion of entangled photons in an integrated quantum walk,” Nature Photonics, vol.7, pp.322–328, 2013. [26] A.D. Raki´c, A.B. Djuriˇsi´c, J.M. Elazar, and M.L.

Majewski, “Optical properties of metallic films for vertical-cavity optoelectronic devices,” Appl. Opt., vol.37, no.22, pp.5271–5283, 1998.

[27] J.A. Roden and S.D. Gedney, “Convolution PML (CPML): an efficient FDTD implementation of the CFS-PML for arbitrary media,” Microw. Opt. Tech-nol. Lett., vol.27, no.5, pp.334–339, 2000.

[28] W.H. Press, S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery, “Numerical recipes in C,” 2nd ed., Cambridge University Press, New York, NY, USA: 1992.

[29] T. Hosono, “Numerical inversion of Laplace trans-form and some applications to wave optics,” Radio Sci., vol.16, no.6, pp.1015–1019, Nov.–Dec. 1981. [30] M.S. Mami¸s and M. K¨oksal, “Solution of

eigenprob-lems for state-space transient analysis of transmission lines,” Electric Power System Research, vol.55, no.1, pp.7–14, 2000.

[31] S. Ohnuki, Y. Kitaoka, and T. Takeuchi, “Time-domain solver for 3D electromagnetic problems using the method of moments and the fast inverse laplace transform,” IEICE Trans. Electron., vol.E99–C, no.7, pp.797–800, July 2016.

[32] R. Ozaki and T. Yamasaki, “Analysis of pulse responses from conducting strips with dispersion

medium sandwiched air layer,” IEICE Electronics Express, vol.15, no.6, pp.1–6, 2018.

[33] S. Kishimoto, T. Okada, S. Ohnuki, Y. Ashizawa, and K. Nakagawa, “Efficient analysis of electromag-netic fields for designing nanoscale antennas by using a boundary integral equation method with fast in-verse Laplace transform,” Proc. Electromagn. Res., vol.146, pp.155–165, 2014. [34] 今野紀雄，量子ウォーク，森北出版，東京，2014. （2019 年 5 月 20 日受付，8 月 26 日再受付， 2020年 1 月 15 日公開）

呉

迪 （学生員）

2016日大・理工・電気卒．2018 同大学 院博士前期課程了．同年同大学院博士後期 課程入学，現在に至る．計算電磁気学の研 究に従事．2017 本会電磁界理論研究会学 生優秀発表賞受賞．2018 本会学術奨励賞 受賞．

浜島

功 （学生員）

2018日大・理工・電気卒．同年同大学院 博士前期課程入学，現在に至る．計算電磁 気学の研究に従事．

井上修一郎

1987早大・理工・応物卒．1993 同大学 院博士後期課程了．スタンフォード大学ギ ンツトン研究所ポストドクトラルフェロー． 1996同リサーチアソシエイト．1998 日大・ 原子力研究所・助手．2000 日大・量子科学 研究所・専任講師．2006 助教授を経て同 大学教授．博士 (理学)．量子光学・量子情報・量子プラズモニ クスの研究に従事．

大貫進一郎 （正員）

1991日大・理工・電気卒．2000 同大学院 博士後期課程了．イリノイ大ポストドクト ラル研究員・客員講師を経て，2004 日大・ 理工・電気助手．現在，同大学教授，2011 イリノイ大客員准教授．博士（工学）．計算 電磁気学の研究に従事．2000 年度鹿島学 術振興財団海外派遣研究員．2013 日本磁気学会論文賞，2014 電気学会電気学術振興賞 (進歩賞) 各受賞．