フォトニック結晶ナノ共振器の設計・解析

結果の要約

FDTDを用いて、中央の空孔を取り除いた共振器（H1）、または中央の空孔と次のリング状の空孔を取り除いた共振器（H2）の共振スペクトルを計算することができます。

H1共振器については、以下に示すように、一定値の誘電率モデルと物理的に現実的なGaAs/AlGaAs材料モデルを用いて比較しています。ピークスペクトルは公表されている結果とよく一致しています。

H2共振器については、論文にある8つの共振をスペクトルに見つけることができます。次に、対称/反対称境界条件を使って縮退モードを分離します。

最後に、H2共振器のすべてのモードを計算し、プロットする方法を示します。例として、A11モードを以下に示します。

シミュレーションセットアップ

構造の概観

この構造は、空気中にある薄い GaAs / AlGaAs 膜上に形成されます。この解析では、170nmのGaAsを45nmのAlxGa1-xAs層（x=0.3）で挟んだ層構造を仮定します。文献ではAlGaAs層が40nmでGaAs層が5nm追加で上乗せされた構造になっていることに注意して下さい。ただし、この5nmの層は結果を大きく変えるものではありませんし、さらに、より細かいメッシュサイズも必要になります。このため、単純にAlGaAs層の厚さを40nmではなく45nmとしています。合計の層厚は260nmになります。この層構造はストラクチャグループに含まれています。このグループでは、全層にn=3.4(文献のシミュレーションで使用された値と同じ)の一定値の誘電体、あるいはPalikによるGaAsの実験データ、またはx=0.3のAlxGa1-xAsの理論モデルが使用できるようになっています。

理論モデルはスクリプトAlGaAs_Adachi.lsfによって生成され、nとkの値のテキストファイルを保存し、それを新しい材料としてマテリアルデータベースにインポートしました。

フォトニック結晶は、膜に空孔を持つ三角格子で構成されています。格子のピッチaは366nmです。この例では、穴の半径が0.37a = 135.4nmの場合を考えます。

共振器はコンポーネント・ライブラリの構造体の1つを使用して作成されます。この構造は完全にパラメータ化されており、ピッチやH_ナンバーなどの様々なパラメータを選択することができます。穴の材質はエッチングに設定されています。レイヤー構造は3つの矩形で構成されています。これらは前述のようにグループ化されパラメータ化されています。最後に、以下に示すように、主要なパラメータをモデルに追加しました。これらのパラメータはモデルによって設定され、子モデルの直接の設定を上書きします。

共振器は、中央の空孔（H1）または中央の空孔と次のリング（H2）のいずれかを取り除くことによって作られます。 これは、モデルの H_ナンバー を 1 または 2 に設定することで実現できます。H1共振器を以下に示します。

メッシュのコントロール

周期構造の場合、メッシュサイズを構造に対して周期的に設定するのが最適です。 メッシュサイズを小さく設定すると、計算時間が長くなり、結果の精度が上がります。ここで使用するアプローチは、初期結果を得て問題の物理を理解するために粗いメッシュから始め、必要に応じて最終的な答えを得るために細かいメッシュに移行するというものです。メッシュオーバーライド領域を追加し、dxとdyの値を直接設定することで、xとy方向のメッシュを制御します。z方向については、適切なメッシュを選択するために、自動のnon-uniformメッシュを選んでも良いでしょう。最初に、シミュレーション領域のMesh accuracyスライダを2に設定し、zメッシュを制御します。xとyのメッシュには、PCの1周期あたり10点の粗いメッシュを選択します（36.6nmグリッド）。三角格子を使用しているので、yメッシュのセルサイズをxグリッドのセルサイズの√3/2倍に設定します。これにより、FDTDメッシュがy方向に完全な周期性を持つようになり、所定のメッシュセルサイズでより優れた精度が得られます。この設定はモデルによって処理されるため、”grid points per pitch”プロパティを選択し、10に設定するだけで十分です。その後、モデルが対応するdxとdyの値を設定します。

シミュレーション領域

2つ目の作業は、シミュレーション領域の大きさを設定することです。シミュレーション領域が大きいほど計算に時間がかかるので、正確な結果が得られる最小の領域を選びたいと思います。FDTD領域の横方向の寸法については、バンドギャップを生成するのに十分なPC格子の列を含める必要があります。これは高屈折率比を持つPCなので、共振器を囲む3列か4列の穴が良いスタート地点です。xスパンは3200nm、yスパンは3200*√3/2で、3列から4列のPCホールを含みます。垂直方向には、空気中に減衰するエバネッセント場があるため、スラブの上下に空気を含める必要があります。まず、スラブの上下に500nmの空気を含めることから始めます。初期結果が得られたら、誤差解析を行ってシミュレーションの寸法が十分に大きいことを確認します。このモデルでは、xスパンとyスパンがピッチの整数倍（yスパンの場合は√3/2倍）であることが強制されます。xスパンを3200nmに設定した場合、オブジェクトを再編集すると、ピッチ366nmの9倍である3294nmに調整されていることに気づくかもしれません。これは、メッシュオーバーライド領域の設定とともに、FDTDメッシュが構造の各周期で同一であることを保証します。これは、周期デバイスのQファクタを正しく測定するために重要です。

シミュレーション時間を2000 fsに設定して下さい。シミュレーション時間に関する内容は、ソースプロパティセクションで説明されます。

有限サイズのオブジェクト(例えばPC共振器など)をシミュレートする際に使用される境界条件は、PML（Perfectly Matched Layer）として知られています。この境界条件は入射してくる輻射を吸収します。PML境界条件は、新しいシミュレーション領域が作成される際にデフォルトで適用されます。デフォルトのPML設定を変更するには、FDTDシミュレーション領域の詳細オプションで2つの主要な変更を行う必要があります：

　・「PMLを介して構造を拡張する」チェックボックスをオフにします。フォトニック結晶の場合、PCをPMLに拡張させるとPMLの性能が向上します。これを行わないと、境界の端にある材料の最後の層が、PMLを完全に通過して、自動的に拡張されます。

　・“pml kappa”の値を5に設定します。分散材料と組み合わせて使用するPMLは、PML内の数値的不安定性の増加につながる場合があります。これは、PML kappaの値を増やすことで、合理的なシミュレーション時間にすることができます。

構造自体が対称性を示す場合、対称境界条件を使用することも可能です。対称境界条件は計算時間を短縮できますが、その代償として、結果に現れる特定のモードが禁止されます（境界条件と同じ対称関係を示さないモードが含まれます）。このPC共振器には、スラブの中心を通る対称平面（z=0平面）があります。この平面に対して対称境界条件を使用すると、TEモードのみが許可され、TMモードは結果から排除されます。この構造には、TMバンドギャップがほとんどまたは全くなく、キャビティモードがTEであると予想されるため、この境界には対称境界条件を使用します。

最後に、共振器を調査する際にはFDTDのオートシャットオフ機能を無効にすることをお勧めします。その理由は、非常に小さな帯域幅で光を閉じ込める高Q値な共振器が、広帯域である元の励起エネルギーのわずかな割合しか閉じ込めない可能性があるためです。その結果、オートシャットオフがあまりにも早く発動することとなり、時間アポダイゼーションのモニターに問題を引き起こす可能性があります。オートシャットオフは、FDTDシミュレーション領域の詳細オプションタブにある “use early shutoff”のチェックを外すことで、無効にすることができます。

パルス光源とシミュレーション時間

共振器を励起するために、いくつかの双極子光源を追加しました。これらの光源は、与えられた領域に、双極子のランダムなセットを作成するAnalysis Groupによって追加されます。この領域はモデルによって制御され、選択された H_ナンバーに基づいています(より大きな共振器では、双極子群は拡張されます)。また、対称および非対称境界条件でのシミュレーションを想定して、これらの双極子は第一象限のみに追加されています。これらの光源は、幅広い周波数を励起するために、短いパルスの放射を入射します。共振モードを効果的に励起するためには、双極子光源の偏光が、励起したいモードの偏光と類似していることを確認する必要があります。この場合、TEモードを励起したいので、垂直（z）偏光磁気双極子を使用します。また、モード分布がnullを示さない位置に双極子光源を配置する必要があります。これを達成するために、該当する領域に双極子をランダムにセットします。もしTMモードを励起したいのであれば、垂直方向に偏光している電気双極子を使用し、zmin境界条件をAnti-Symmetricに変更する必要があります。
光源設定の最後のステップは、励起パルスの定義です。これは単純に、シミュレーションしたい適切な波長範囲を選択することで可能です。この場合、1000nmから1400nmを選択します。各ダイポールの設定を個別に変更する必要がないように、グローバルソースプロパティを使用しました。各ダイポールは、デフォルト設定ではないグローバルソースプロパティを用いて設定される必要があることに注意して下さい。
シミュレーション領域のプロパティで定義したシミュレーション時間は2000fsであり、パルス長よりもはるかに長いです。これは、パルスがモードを励起した後の共振器の場をモニターしたいからです。最初のパルスの後の共振器の場を解析することで、共振器モードの周波数が明らかになります。

シミュレーションモニター

モニターには”resonance finder”と呼ばれる解析グループがあります。これらの時間信号をフーリエ変換することで、共振モード周波数を明らかにします。これは、グループ内の解析スクリプトによって自動的に処理されます。デフォルトでは、共振器の領域にランダムに拡散した、合計8つの時間モニターを選択します。これらのモニターの範囲はモデルによって設定され、共振器に選択されたH_ナンバーに依存します。これにより、特定のモードのnullを避けることができます。
また、モード分布を得るための、周波数領域分布モニターもあります。これらのモニターは、mode_profilesと呼ばれる解析グループに格納されています。この解析グループは、最大12モードのモード分布を見つけることができます。ユーザーは、モードの共振周波数とアポダイゼーションの設定を指定する必要があります。アポダイゼーションの設定は、正しいモード分布を見つけるために非常に重要です。アポダイゼーション設定の詳細については、User GuideのApodizationを参照してください。

結果－H1 共振器

H1共振器は、下図に示すように、フォトニック結晶から1つの正孔を取り除くことで形成されます。これは、モデルのH_numberプロパティを1に設定することで実現できます。

Tang_cavity.fspを読み込み、実行して下さい。

次に、”resonance finder”オブジェクトを編集し、”Run Analysis”ボタンを押すと、以下の図と結果が得られます。

スペクトルを規格化単位の関数とした(a/lambda)のプロットすることも、以下のスクリプトコマンドによって可能です。スクリプトは、スクリプトプロンプトへ貼り付け、または新しいスクリプトファイルへの作成も可能です。

select("::model");
a = get("a");
mname = "resonance finder";
runanalysis(mname);
f = getdata(mname,"f");
spectrum = getdata(mname,"spectrum");
plot(f*a/c,log10(spectrum),"frequency (a/lambda)","spectrum (a.u., log scale)");
setplot("x min",0.2);
setplot("x max",0.44);

これは以下に示す図を作ります。

この例では、結果にいくつかのピークが見られます。これらのピークは共振モードとバンドギャップエッジに対応しています。これらの周波数のモードは群速度が非常に低く、すぐに減衰しないため、バンドギャップエッジにピークが観測されることが多いです。共振モードに対応する2つのピークは約0.31と0.37に現れ、スラブ導波路において一定の誘電率を想定して計算したY. Tang論文にあるように、1次モードと2次モードの結果とよく一致しています。

なお、材料に GaAs と AlGaAs のモデルを用いてシミュレーションを再実行することも可能であることに留意して下さい。これは、”layer structure”グループの”use constant dielectric”プロパティを0に設定することで可能です。以下の図は結果を比較しています。

低周波ではスペクトルがよく一致していますが、高周波ではスペクトルがずれていることが分かります。このシフトは、GaAsが260nmの層のうち170nmを占め、以下に示すように高周波では3.4より高い屈折率を持つことを考えれば整合性のある結果であると分かります。

残りの例では、定誘電率モデルを使用します。これはTangらのシミュレーション結果と最もよく一致するはずですが、分散モデルの方が実験結果とよく一致するはずです。

次のステップ

・共振モードのモード分布を決定するために、周波数領域分布モニターの設定（中心周波数とアポダイゼーション）を調整し、これらの共振のモード分布を計算することができます。
・ モードの対称性を知り、縮退モードのモード分布を見るには、x=0とy=0の平面に様々な対称と非対称の境界条件を組み合わせて適用し、シミュレーションを繰り返します。
・ 精度の高い結果を得るためには、より細かい格子サイズでシミュレーションを繰り返すことができます。また、PMLの境界を共振器から遠ざける/近づけることによる影響を調べることも有効です。

次のセクションで、H2共振モードを探す際に、これらのステップのほとんどを実証します。

結果－H2共振器

H2をシミュレーションするためには、モデルのH_numberを単純に2に変えることが必要です。

最初の結果

シミュレーションを実行した後に、以下のスクリプトを使用することによってスペクトルを規格化周波数の単位(a/lambda)でプロットすることができます。

select("::model");
a = get("a");
mname = "resonance finder";
runanalysis(mname);
f = getdata(mname,"f");
spectrum = getdata(mname,"spectrum");
plot(f*a/c,log10(spectrum),"frequency (a/lambda)","spectrum (a.u., log scale)");
setplot("x min",0.27);
setplot("x max",0.36);

これらの結果を見ると、目的の周波数範囲に8つのピークがあることがすぐにわかります（ここでは実験結果が存在するおおよその周波数範囲に限定しています）。ピークの周波数は実験結果の数パーセント以内ですが、シミュレーション結果のピークと実験結果の各ピークとが必ずしも一致するか、すぐにはわかりません。シミュレーション結果には、実験データでは確認できなかった余分なピークがあります。一般的に、最初の試みとしては、一致度は非常に高いといえます。また、2つのピークが分割されているように見えますが、これは理論的に縮退しているピークに対応し、デカルトメッシュによって縮退が解除されているためであると思われます。三角格子では、60度の格子対称性に基づいてモードが縮退することがあります。デカルトメッシュはこの対称性を人工的に崩す可能性があり、この効果は粗いメッシュほど顕著になります。より細かいメッシュでシミュレーションを再実行するには
　・モデルの “grid points per pitch”プロパティを、10ではなく20に設定します。これにより、dxとdyのサイズが半分になります。
　・ FDTDシミュレーション領域のメッシュ精度スライダを2ではなく4にします。

シミュレーションを再実行（かなりの時間を要する）すると、細かいメッシュと粗いメッシュの結果を以下のように比較することができます。2つの結果には妥当な一致が見られ、周波数シフトは一般的に周波数が高いほど（波長が短いほど）大きくなると予想されます。また、メッシュを細かくするとピークの分裂がなくなることもわかります。しかし、特に波長1000nmではλ0/(n*dx)~8であることを考えると、粗いメッシュは驚くほど正確です。これは、Lumericalのコンフォーマルメッシュテクノロジーの使用により、より粗いメッシュサイズでも良好な結果を得ることが可能になった結果でもあります。この例の残りは粗いメッシュで作業を続けます。論文用の質の高い結果を得るためには、さらに何度かシミュレーションを繰り返し、結果が許容範囲内に収束したことを確認する必要があります。

詳細な解析

この問題をさらに理解するために、対称性と非対称性を使ってモードを分類します。また、周波数領域分布モニターを使って、モード分布を計算します。最初のステップでは、周波数領域パワーモニターを追加し、周波数とアポダイゼーションの設定を行います。

　・アポダイゼーション　シミュレーションは2000fsで実行されます。時間中心を1000fs、時間幅を200fs、完全なアポダイゼーションを選択する必要があります。これは、アポダイゼーションスペクトルのFWHMが4log(2)/(2pi*200fs) ~ 2 THzになることを意味します。このアポダイゼーションを選択するということは、ピークが2THz以上離れており、正しいピークの約2THz以内にあるモニターの周波数を選択しさえすれば、正しいモード分布を測定できるということを意味します。規格化されていない単位でスペクトルをプロットしてみると、これは問題ないことが分かります。アポダイゼーション設定の詳細については、User GuideのApodizationを参照してください。

　・周波数設定　220～300 THzの周波数を持つモニターを8つ追加したいと思います。実際には、276 THzにある人工的な2つのピークを考慮して、9つのモニターを設置します。

これらのモニターは手作業で追加することもできますが、シミュレーションはmode_profilesという解析グループを含んでおり、x-y平面に最大12個の周波数モニターを追加できるようになっています。 これにより、すべてのモニターのアポダイゼーション設定を簡単に調整することができます。9つの共振周波数の計算を簡単にするために、スクリプトファイルTang_cavity_setup_profile_monitors.lsfが、所望の帯域幅内の9つの共振を計算し、それに応じてmode_profile解析グループの周波数を順次設定します。共振周波数は、メッシュ精度、材料設定、または共振ピークをシフトさせる可能性のあるその他の設定に変更があった場合、いつでもリセットされる必要があります。

モードの対称性

次のステップは, モードの対称性を考慮することです。 x=0とy=0の対称平面上の対称境界条件と反対称境界条件の組み合わせを変えて、シミュレーションを4回繰り返すことで、これを実行します。これらのシミュレーションはそれぞれ、完全なシミュレーションの1/4の時間で実行されます。これは、この例題に含まれている“sweep symmetry”オブジェクトを使用することで達成可能です。このオブジェクトは”Optimization and Sweeps”ウィンドウにあります。

オブジェクトの編集は、4つのパラメータ掃引を実行するように構成されています。各パラメータ掃引では、異なる対称条件をxmin, ymin境界に適用します。

パラメータ掃引オブジェクトを実行後、スクリプトの次の行が全てのスペクトルをプロットします。これらの行はTang_cavity_symmetry.lsfのファイルにも保存されます。

select("::model");
a = get("a");
mname = "resonance finder";
runanalysis(mname);
f = getdata(mname,"f");
spectrum = getdata(mname,"spectrum");
sweepname = "sweep symmetry";
spectrum_sym = pinch(getsweepdata(sweepname,"spectrum"));
plot(f*a/c,log10(spectrum),
           log10(pinch(spectrum_sym,2,1))+8,
           log10(pinch(spectrum_sym,2,2))+16,
           log10(pinch(spectrum_sym,2,3))+24,
           log10(pinch(spectrum_sym,2,4))+32,
           "frequency (a/lambda)","spectrum (a.u., log scale)");
legend("No symmetry",
       "Symmmetric    /Symmetric",
       "Anti-Symmetric/Symmetric",
       "Symmmetric    /Anti-Symmetric",
       "Anti-Symmetric/Anti-Symmetric");
setplot("x min",0.27);
setplot("x max",0.36);

下図は結果のプロットです。凡例は、対称性が使用されなかった場合、またはx最小y最小の境界で対称性が使用された場合を示しています。例えば、Symmetric/Anti-Symmetricは、x minの境界条件がSymmetricで、y minがAnti-Symmetricであったことを意味します。一般的に、あるピークは対称境界（非縮退モード）の1つの組み合わせでのみ存在し、他のピークは対称境界（縮退モード）の複数の組み合わせで現れることを示しています。下図の各ピークには、縮退モードまたは非縮退モードを示す “D”または”N”のマークが付けられています。

これで、Y. Tang の論文で示された実験結果を参照し、シミュレーションで得られたモードと実験で得られたモードとを照らし合わせて特定することができます。次の表は、Y. Tang 論文のモード周波数と、線幅/偏光解析によって決定されたそれらの縮退の一覧を示しています。これで、シミュレーションデータのどのモードが実験データに対応するかを簡単に確認することができます。結果は下の表にまとめられています。

ピーク# (モード)	縮退	実験周波数	シミュレーション周波数
1 (E11) 2 (A11) 3 (E12. E13) 4 (B2)	D N D N	0.285 0.306 0.312 / 0.313 0.323	0.276 0.297 0.307 0.316
5 (A2)	N	0.329	0.328
6 (E14) 7 (A12)	D N	0.341 0.348	0.337/0.339 0.345
8	D	NA	0.354

先の表の結果から、シミュレーションデータの8番目のピークは、実験で量子ドットが励起できる周波数よりもわずかに高い縮退モードであることが分かります。このモードが存在することを知ると、実験データの0.352付近に、このモードに対応すると思われる非常に小さなバンプを見ることができます。

シミュレーションデータは実験と非常によく一致しています。モデルと実験との差は最悪の場合で～3.5%です。これは平面波モデルとの差（実験と比較した）の最悪の場合と一致しています。実験とシミュレーションのモード周波数の平均差は2％未満です。これは、今回使用した単純なモデルとしては非常に良い一致です。この結果は、より細かいメッシュの使用、より現実的な材料モデルの使用、またはシミュレーションの構造寸法が実験で作製されたデバイスと一致していることを確認する(例えば、SEMインポートを使用して構造を定義するなど)ことによって、確実に改善することができます。

結果－H2共振器モード分布

前のステップでは、対称性を変えていくつものシミュレーションを行いました。また、mode_profiles解析グループオブジェクトを設定し、9つの共振周波数のそれぞれの電場を記録するモニターを適切な時間アポダイゼーションで設定しました。
パラメータ掃引オブジェクトは、各対称設定のスペクトルを収集するだけでなく、｜E｜2の行列も収集します。この行列は、各シミュレーション中の mode_profiles 解析グループによって作成されます。パラメータ掃引では、このすべてのデータが1つの行列になるよう並べます。掃引の最後には、E2という4次元の行列(サイズが length(x) * length(y) * 9 * 4)ができます。最初の2次元は、x-y平面におけるx軸とy軸です。3番目の次元は、収集した9つの共振周波数に対応します。最後の次元は、パラメータ掃引した4つの値で、対称/非対称の４つの組み合わせに対応します。
スクリプトファイル(plot_modes.lsf)は、4次元行列E2を取得します。その後、より細かいメッシュに補間し、構造の概観の結果を屈折率モニターからいくつか収集し、結果をよりイメージしやすくするための解像度の高い行列を作成します。そして、特定の共振周波数、特定の対称/非対称条件を選択することにより、さまざまな結果をプロットします。各図は自動的にjpgファイルとしてエクスポートされます。結果は以下に示されています。

縮退モード

非縮退モード

ここで示されたモード分布は、平面波モデルを用いて計算されたY. Tangの論文結果と非常によく一致していることが分かります。

varFDTDソルバーの利用

このセクションの目的は、以下に示すフォトニック結晶ナノ共振器の解析にMODEのvarFDTDソルバーがどのように使用できるか示すことです。これは、Y. Tangによる論文の結果を再現することで達成されます。この論文には、共振器構造に関する詳細な情報と、比較するための実験データおよびシミュレーションデータが含まれています。

セットアップ

シミュレーションファイルはcavities_3Dpcm.lmsとなります。全体として、MODEでのシミュレーションのセットアップ方法はFDTDと非常によく似ています。 構造はFDTDのサンプルファイルから直接コピーできます。 シミュレーション領域、光源、モニターも同様に設定します。セットアップは非常に似ていますが、もちろん、いくつかの違いがあります：

-　Propagator object – Effective index タブ：ここでは、スラブ分散を無視するNarrowbandオプションを含んだデフォルト設定を使用します。 精度を上げるにはBroadbandオプションを使用してください。

結果

シミュレーションは非常に速いので、数秒で終了するはずです。続いて、Qunalysis Objectを実行して共振周波数を求めてください。この処理にはかなりの数のfftが含まれるので、いくらかの時間が必要になります。以下のコマンドを使うことで、結果をスクリプトプロンプトに出力することができます。スペクトル(そして計算に用いたフィルター)のプロットにQanalysis objectを使うことができます。

runanalysis;
Q = getresult("Qanalysis","Q");
?Q.f/1e12;
?Q.Q;
result: 
192.376  
207.302  
197.875  
result: 
202.889  293.017  237.532  

各共振ピークの振幅はそれほど意味がありません。なぜなら、系を励振する光源の位置や、応答を測定するタイムモニターの位置によるからです。

計算された共振周波数から、二つのプロファイルモニターをセットアップすることが出来ます。これで、二つの周波数をモニターし、そのプロファイルを得ます。

参考文献

Tang, Y. Mintairov, A.M. Merz, J.L. Tokranov, V. Oktyabrsky, S., Characterization of 2D-photonic crystal nanocavities by polarization-dependent photoluminescence, 5th IEEE Conference on Nanotechnology, July 2005, vol. 1, pp 35-38