散乱行列法を用いたフォトニック結晶の光学応答解析
岩長 祐伸
東北大学 大学院理学研究科 物理学専攻
フォトニック結晶の光学特性を定量的に解析，評価するために，散乱行列法を用いて，反射，透過，回折の光学スペクトルを SX-7 上で計算した．この稿では周波数領域の計算手法である散乱行列法について概説したあと，計算結果の実例として金属フォトニック結晶の光学スペクトルを中心に紹介する．

1. はじめに1.1. 研究の背景
フォトニック結晶とは，光の波長程度の周期をもつ構造体の総称である．光の究極の制御媒体として注目されてから，20 年ほどの歴史をもつ．近年の実験技術の進歩によって，ミクロンスケールでの光の制御と応用，新奇な物理，光学現象の開拓の舞台として，光科学の分野では非常に魅力的な研究テーマである．
フォトニック結晶には，全方向に十分な大きさを持つ「結晶」と基板に電子線リソグラフィ法などを用いて構造体を作る「スラブ（平板）型」の 2 種類に大別できる．
フォトニック結晶にもっとも特徴的な性質と考えられている「光の閉じ込め」効果の完全な実現には 3 次元結晶を必要とするが，微小サイズの光集積回路用の導波路や面発光レーザー素子としてはスラブ型が有望視されている．実験の面では，加工のしやすさからスラブ型の性質が中心的に研究されてきた[1]．

図 1 フォトニック結晶の概念図．左：光の閉じ込めが可能な woodpile型の 3 次元結晶，右：1 次元フォトニックスラブ．

1.2. 散乱行列法
フォトニック結晶のような構造体の電磁波にたいする応答を計算する手法として，有限差分時間領域法（FDTD 法）が広く用いられている．FDTD 法は Maxwell 方程式を時間差分して解く方法であるから，計算式が自然にベクトル化されており，SX-7のようにベクトル型コンピュータでの運用に適している．実際，本誌にも多くの運用例が報告されている．
今回の共同研究で用いた方法は FDTD 法とは違い，周波数領域での光学スペクトルを得るために，散乱行列を構成する方法を選んだ．電磁波の空間成分を擬平面波展開して扱う．この散乱行列法はその名のとおり，行列演算に帰着できることから，SX-7のベクトル演算機能と並列演算を十分に生かすことが期待できる．
ここでは計算コードの詳細には触れないが，散乱行列法が FDTD 法にくらべて広く認知されていないので，基本的なアルゴリズムについて概説し，散乱行列法の紹介にかえたい．
散乱行列法は一般に知られている転送行列法の一種である．二つの方法は数学的には等価である．その様子を見るために境界条件について考えてみる．簡単のために光はz軸に沿って進む場合を考える（図 2）．境界の接線方向の接続条件は，電場と磁場の接線成分が等しいことである．光の進む向きが ±z方向あることを考慮すると，接続条件はz = z0でEI,x(+) + EI,x(–) = EII,x(+) + EII,x(–)
HI,y(+) + HI,y(–) = HII,y(+) + HII,y(–)
である．真空中にある平面波が +z方向に進む状況をexp(ikz – iωt)と表すと，領域Iの電場EI,,x(+) は振動成分exp(iknI z – iωt)をもつ．同様に電場EI,,x(–) はexp(– iknI z – iωt)に比例する．指数関数の引数に媒質中の屈折率nIが含まれる．

図 2 境界条件と電場成分の組替え．
通常の転送行列法では，組(EI,x(+), EI,x(–))で各層の電場を表現する．この組の並進はと表すことができ，媒体の屈折率nIが複素数で虚部Im(nI)が正のとき（つまり，媒体が光を吸収するとき）にはEI,x(–) はexp[Im(nI)kd]に比例して指数関数的に増大してしまう．
この性質が転送行列法を用いた数値計算を不安定にすることはよく知られている[2]．
一方，散乱行列法では図 2 で示すように組を取り替える．この組替えによる利点は

次式のように並進操作による指数関数的な増大項が表れないことにある：
行列の右下成分において exp 内の符号が転送行列の場合と反対になっている．つまり，発散的な項が現れないという利点がある．
フォトニック結晶のように構造を持った対象でも電磁波の伝播は境界条件と並進の繰り返しで表すことができる．上述の組替えを層の数だけ帰納的に繰り返すことにより，対象からの光学応答を計算できる（図 3）．入射波をu(0)にとると，反対向きの入力d(n+1) = 0 となり，u(n+1) が透過成分，d(0) が反射成分に対応する．層に周期性があるときはu,dは電磁波を擬平面波展開したFourier係数からなるベクトルにとることができる．このとき，u,dは回折光の各次数に対応する成分を自然に含むことになる．

図 3 散乱行列の構成と概念図．左：各層において進行波uと逆向きの進行波dを設定する．右：組替えた組(u(i),d(n+1–i))にたいする散乱行列Sと入力，出力の関係．
以上のように，散乱行列法の特徴は指数関数の発散的な項が現れないアルゴリズムになっており，数値計算の安定性を保証している点にある．散乱行列を構成するときには境界条件と並進を行列演算で表して，帰納的に散乱行列を構成する．したがって，ベクトル演算と並列機能をもつ SX-7 上での運用にとくに適したアルゴリズムになっている．
散乱行列法を用いた電磁波の計算と光学応答に関しては，文献[2]に詳しい．今回の共同研究では文献[3]の計算手順をおおむね踏襲してコード化し，SX-7 上で運用した．これらの文献は同じ対象を扱っており，本質的に同じものである．次節で計算結果を紹介していく．

2. 結果
金属フォトニック結晶は表面プラズモンという金属状態の素励起による効果を顕著に観測することができる対象として注目を浴び，その制御を目指した研究が多くなされている．以下では，今回の共同研究の結果から 2 例を紹介する．

2.1. Lamellar 型金属フォトニック結晶スラブ
図 4 のように x 軸方向にのみ周期性があるとき，単色光の電場ベクトルはE(x,y,z,t) = Σm Em(z)exp(ikmx)exp(–iωt)
km = kx + 2πm/Lx (m = 0, ±1, ±2, ···)のように擬平面波展開できる．Lxは周期である．図 4 の場合，入射光は空気側から来て，電場ベクトルはxz平面内にある．

図 4 金属フォトニック結晶スラブの構造．周期は 500 nm で銀(Ag)と空気(Air)の周期層からなる．図は xz 面での断面図であり，Ag のロッドは紙面に垂直な y 軸方向に無限に長い構造をしている．
つぎに反射，透過スペクトルを示す．図 5 は垂直入射条件( θ = 0º)下での光学応答である．このスペクトルの特徴としては，1.2 eV にスリットを導波路とするモードが存在して，透過率がほぼ 1（効率として 100%），反射が 0 になっている．これは銀の厚さが 300 nm と表皮効果の厚さにくらべて十分に厚いことを考えると，非常に特徴的な結果である．導波路モードの観点からいうと，1.2 eV のモードは TEM モードであり，TE や TM モードのようなカットオフ（周波数の下限）がない．これにより，周期長からかなり長波長域で高い効率の透過が生じている．
また，1.7 eV より高エネルギー側では回折が生じて，導波路モードとは異なる回折モードや表面プラズモンと相関がある反射の減少が現れている．回折スペクトルも同時に計算しているので，その結果をつぎに示す．
図 6 は図 5 と同じ入射光条件下での回折スペクトルを示している．図 6(a)は全回折強度を示しており，図 6(b)は回折光の成分を反射成分（点線）と透過成分（破線）に分解したものである．
図 5 の反射スペクトルと比較すれば，回折光の強度が反射の減少と対応していることがよく分かる．とくに，1.7-2.5 eV の範囲では回折光は透過成分からなり，2.5 eV以上のエネルギー範囲では反射成分が主要である．

金属フォトニック結晶からの回折成分は基板の石英側（透過方向）と空気側（反射方向）に分けることができるが，石英側への 1 次の回折チャンネルが開くのが 1.7 eVであり，空気側への 1 次の回折が可能になるのは 2.5 eV である．この性質を回折成分の分解（図 6(b)）は定量的に示している．

図 5 Lamellar 型金属フォトニック結晶スラブの(a)反射スペクトル，(b)透過スペクトル．入射光は垂直入射で電場ベクトルの偏光は図 4 参照．
回折スペクトルを見ると，1 つの疑問がわく．
「2.7 eV の回折強度のピークは一体なにに起因しているのだろうか？」2.7 eV のエネルギーはいずれの回折次数とも対応していない．
この疑問に答えるために構造パラメータを変えて知見を得ることが有用であると考えられる．そこで周期を 500 nm に固定して，周期層の厚さを変化させて回折スペクトルを計算した．
図 7 は周期層（図 4 中の d）を 20 nm から 800 nm まで変えて光学応答を計算して得られた回折スペクトルである．いずれも全回折強度をグラフとして示している．
図中のT±1は回折光の±1 次の透過成分，R±1は±1 次の反射成分，SPP(1Q) は周期層と石英基板の界面に生じる（還元）表面プラズモンの 1 次の共鳴エネルギー位置，
SPP(1A) は周期層と空気の界面に生じる表面プラズモンの 1 次の共鳴エネルギー位置を示している．SPP(2Q) は 2 次である．
銀と空気の周期層の厚さが 100 nm 以上になると，大きな反射型の回折が生じる．
しかし，そのピークは厚さの変化にともなって連続的に変わり，ブラッグ回折や表面プラズモンのような特定の共鳴に相関があるわけではない．金属フォトニック結晶の構造によって決まる励起モードが回折の輻射強度に寄与していると考えられるが，その励起モードを直接解析することに成功した報告例はいまのところない．散乱行列のカーネルが励起モードを表すと考えられ，これを数値的に効率よく求めることは今後の課題である．
今回の研究では，定量的な数値計算によって，金属フォトニック結晶内の励起モードと相関の強い応答が回折スペクトルに現れることを示し，今後の励起モードの解析に資するデータを得た．

図 6 回折スペクトル．(a)全回折強度スペクトル．(b)回折の反射成分（点線）と透過成分（破線）への分解．

図 7 構造変化（図 4 での厚さ d を変えること）にともなう回折スペクトルの変化．記号は本文中の説明を参照のこと．
なお，この節のより詳細な結果については文献[4]として公開している．

2.2. 溝の深い回折格子型金属フォトニック結晶スラブ
もう 1 つ具体例をあげる．前節では周期層が単層で表現できる場合であったが，今回は金属フォトニック結晶の表面が滑らかである場合について述べる．このとき，計算を実行するために周期層を図 8 のように多層に分割する必要がある．

図 8 溝の深い金属フォトニック結晶の 1 周期の構造断面図．周期層を多層に分割して滑らかな表面を近似している．周期 660 nm にたいして，金属部の深さがのべ 800 nm ある．この例では金属は金(Au)である．空気側からの入射光の波数ベクトルを k としている．

0.040.0300.020.010.00
550600650700750Wavelength (nm)800

Incident -k10.050.040.0300.020.010.00
550600650700750800

Wavelength (nm)

図 9 入射光が逆の向きでの光学スペクトル．入射光の向きの定義に関しては図 8 参照のこと．入射角は 4.5º で，赤線が反射，青線が透過，灰色線が吸収を示している．

Transmittance

0.05

Reflectance

Incident k1

Transmittance

Reflectance

図 8 は実際に加工して得られた試料を近似した構造で，以下に示す計算結果と実験による測定結果は定量的，定性的に一致をみている．
図 9 は空気側から波数ベクトル k の入射光がきたときの光学スペクトルと基板側から反対向きの入射光（記号として，–k の向きと表す）がきたときの光学スペクトルを示している．赤線が反射スペクトル，青線が透過スペクトル，灰色線が吸収スペクトルを表している．

とくに興味深いこととして，反射や吸収が大きく異なっているにもかかわらず，透過スペクトルのみがいずれの入射光の条件下でも等しい．
この透過率不変性は，相反定理を応用することで回折や吸収が生じる一般的な条件下でも成り立つ性質であることを証明し，2006 年 3 月の日本物理学会で口頭発表した[5]．したがって，透過率不変性は透過率に関する相反性ということができる．相反定理の前提となる電磁応答の局所線形性は一般的に広く成り立つ条件であることから，透過率の相反性はかなり普遍性のある性質であるといえる．
数値計算で精度よく，透過率の不変性を確認できたことから，不変性を証明する動機付けが明確になり，実際に証明できた．この結果は，数値実験から支持される物理的な性質を実際に証明したというよい例になった．

3 . 数値計算の効率
最後に散乱行列法を用いた計算を実行する際の CPU 時間について述べる．
はじめにも書いたように散乱行列のアルゴリズムは基本的に行列演算（掛け算，足し算）の繰り返しがほとんどであるために SX-7 上ではベクトル化率が高く（99%以上），自動並列化できる．
しかし，各周期層での固有モードをもとめるために固有値，固有ベクトルをもとめる手順が 1 回あり，その部分が ASL 関数では並列化できないために並列化率が下がってしまう．現状ではこれは回避できないが，1000×1000 程度の大きな倍精度複素行列を使うと，固有値問題を解くために並列化率が大きく下がる．ちなみに 300×300 程度の行列にたいしては，ほぼ理想的な並列化率を得ている．これは情報シナジーセン
ターの職員の方々のご協力の賜物である．
今回は 1 次元系の結果について報告したが，今後 2 次元系に展開していく計画であり，扱う行列サイズが 2000×2000 程度にまで大きくなっていくと予想される．近い将
来，ASL 関数で固有値問題を解く関数の並列化版が開発されることを期待したい．

謝辞
情報シナジーセンター・大泉氏にはとくにお世話になり，並列化にともなうコードの改良など多くの面で支援いただきました．ここに謝意を表します．

参考文献[1] 野田, 応用物理 74, 147 (2005).
[2] L. Li, J. Opt. Soc. Am. A 13, 1024 (1996).
[3] S. G. Tikhodeev, A. L. Yablonskii, E. A. Muljarov, N. A. Gippius, and T. Ishihara, Phys.
Rev. B 66, 045102 (2002).
[4] M. Iwanaga and T. Ishihara, cond-mat/0510560.
[5] 岩長, Vengurlekar, 石原, 日本物理学会講演概要集 61 巻(1) 第 4 分冊, p. 767(2006).