単一開口スカラー回折

はじめに	3D領域における完全整合層
電磁気学におけるヘルムホルツ(Helmholtz)方程式	用語集
2Dのスカラー回折モデル	参照
2Dにおける完全整合層

はじめに

回折現象は，波動方程式で完全に記述される．ホイヘンス・フレネル(Huygens-Fresnel)の原理によると，波が障害物や開口部を通過するとき，その障害物を囲む，または開口部内部のすべての点が球面波の点源としての役割を果たす．これらの波の重合せにより，特徴的形態を持つ波面が生成される[Born & Wolf, 1999]．波面の強度プロファイル ℐ [W/m²]は，電場と磁場の外積の大きさ|E×H|²に比例し，回折パターンと呼ばれる．

回折はホログラフィー顕微鏡法，結晶学，天文台等の感知や画像化に使われる最新システムを支える基本現象である．波動方程式を解くことは，特にフーリエ変換等の積分法が正確とは言えない，電磁波長と同程度の大きさの系において，重要な工学的問題である[Born & Wolf, 1999]．

このノートブックで示したシミュレーション結果のアニメーションの多くはRasterizeの呼出しを使って生成されている．これは，このノートブックに必要なディスクスペースを削減するためである．欠点は，アニメーションの見た目が呼出しを使わないときほどくっきりしないという点である．高品質のグラフィックスを得るためにはRasterizeの呼出しを削除するかコメントアウトするかするとよい．

高忠実度の可視化を得るためには，ラスタライズ処理をコメントアウトする．

このチュートリアルで使われる記号および対応する単位は，用語集セクションに要約してある．

有限要素法とOpenCascadeLinkパッケージをロードする．

電磁気学におけるヘルムホルツ(Helmholtz)方程式

等方性の非分散媒体では，回折はベクトル値のsource-freeヘルムホルツ偏微分方程式でモデル化することができる[Jackson, 1999]：

ここでは媒体の中の波数であり，で波長に関連している．材料における波数は，（は相対浸透率，は材料の比誘電率）で自由空間波数に関連している．とは材料の磁化および偏光の大きさである[Jackson, 1999]．一般にとは二階テンソルである．等方性材料では，とはスカラー定数である．

方程式(1)を数値的に解く最初のアプローチとして，電場の一つの成分だけを考える．直交座標において，ベクトル場のラプラシアンは各成分のラプラシアン∇²A=(∇²A_x,∇²A_y,∇²A_z)なので，方程式(2)は3つの独立スカラー微分方程式の集合である．これは対角のおよびテンソルを持つ材料の場合にのみ成り立つ．非対角成分を持つ異方性材料の場合，PDEは混合された従属変数を持つ可能性がある．

場の振動方向を記述する電場ベクトルは，常に電磁波の伝搬方向に垂直である[Born & Wolf, 1999]．つまり，方向に伝搬する電磁波は，および成分にのみエネルギーを含むということである．

場が，線形偏光の電磁波等，1つの成分だけで記述できる場合，その場はスカラー量と言われる．

2Dのスカラー回折モデル

この応用例では，開口幅で不透明なスクリーンに光を当てる，自由空間波長のヘリウムネオンレーザーを考える．開口部は切込みで方向には無限に大きいと想定されるため，方向の回折波へのかく乱はない．したがってこのシステムを2Dでモデル化することができる．レーザーと開口部自体は，放射境界条件でモデル化される．不透明なスクリーンは，光が伝わったり反射されたりすることのない完全吸収材料でできている．境界はオープンシステムを模倣する，入射光を吸収する人為境界である．オープンシステムでは，光は反射することなく無限に伝搬する．したがって，他の境界の場はなので唯一の光源は開口部である．シミュレーション領域の寸法は，である．レーザーは方向に偏波を発する．偏光電磁場では，成分にエネルギーはない．領域内部の屈折率は1である．

上の配置図では，見やすいように90度回転させて伝搬軸が垂直になるようにしている．また，配置は方向について対称であり，方向はスクリーンの外側を向いている．

パラメータの設定

回折スカラー波について，方程式(3)の入射スカラーの電磁波パラメータとヘルムホルツ微分演算子を設定する．

モデルパラメータを設定する．

変数とPDE演算子を設定する．

NeumannValue境界条件が使えるように微分演算子はInactive形式である．詳細は古典的な偏微分方程式をご参照いただきたい．

PDEモデルを設定するためにHelmholtzPDEComponentを使うこともできたが，後で完全整合層を使うためHelmholtzPDEComponentはあまり実用的ではないのである．セクション間のワークフローを一貫したものにするために，基本的な偏微分方程式項でPDEをモデル化する．

領域

次にシミュレーション領域を設定する．シミュレーション領域は対称軸をとして対称なので，x方向の半分の領域のシミュレーションで十分である．

縮小されたシミュレーション領域の境界メッシュを設定する．

境界メッシュと境界マーカーを可視化する．

黒い数字は境界要素マーカーである．マーカー1は対称線境界（破線）に割り当てられ，マーカー2は境界（青い線）を示し，マーカー3は（赤い線）であり，マーカー4は開口境界の半分を表し，（黒い線）と呼ぶ．前述のようには開口部であるが，ここでは入射光は境界条件によってモデル化される．

要素メッシュは境界メッシュから生成される．波を解決するのには，通常波長につき5つの要素で十分である[Jin, 1993]．初期の離散制約はが推奨される．システムよっては，正確な結果を提供するためにこれより小さい制約を使わなければならない場合がある．この場合，滑らかな曲線を得るためにを使う．

最大セル長が

の完全メッシュを生成する．

境界条件

後は各境界に対する条件を設定するだけである．理論に従うと，境界は境界内の反射によるすべての要因を排除するために無限でなければならない[Jin, 1993]．無限に大きい領域上の解は計算できないので，反射を避けるために辺に人為吸収境界を置く必要がある．は入射光をすべて吸収する不透明なスクリーンなので，これにも吸収境界条件を使うことができる．一次パデ(Padé)の吸収境界条件を使う[Jin, 1993]．

境界要素マーカー2と3に吸収境界条件を設定する．

吸収境界条件と放射条件については時間領域における音響学で解説してある．

の開口部を照らす入射波は以下の形式の放射境界条件でモデル化することができる：

ここでは境界法線，は入射電磁波の伝搬方向，は電磁波の振幅である．境界条件は光の入射項と吸収境界条件項の2つの項の和である．

2つ目の項は，が電磁波を反射しないで入射光だけが境界から外に出ることを確実にするために必要である．この場合，境界法線と波の伝搬方向の両方がz方向であるため，内積は1になる．は1に設定する．

境界要素マーカー4の開口放射条件を

で設定する．

要素境界マーカー1の対称軸はノイマン(Neumann)0値でモデル化できる．境界上に他に何も指定されていない場合は，これがデフォルトの境界条件になる．

モデルの評価

PDEの数値解をNDSolveで計算する．

方程式を解く．

可視化

数値解の絶対値をプロットする．

アニメーションの質を向上させる方法はこちらでご覧いただきたい．

z軸に沿った電場の強さの伝搬を可視化する．

画像の質を向上させる方法はこちらでご覧いただきたい．

上の強度プロットは，電場が軸周りに集中していることを示しているので，中心領域から遠くで過度に細かいメッシュを使うことは不必要に計算的に高価になる可能性がある．MeshRefinementFunctionで設定した，可変の最大セル寸法を使うと，該当する領域の精度を保ちながら，メッシュ要素の数を削減することができる．

伝搬方向に沿って結果を可視化する．

先に述べた通り，吸収境界条件()は無限領域をモデル化するための近似境界条件である．上のプロットではがに近付くにつれて振動が見られる．これらの振動は吸収境界条件における不要な反射によって起こる．次のセクションでは，オープンシステムの計算領域を打ち切る，より正確な別のアプローチを紹介する．

2Dにおける完全整合層

上で使ったような吸収境界条件は，完全吸収境界ではない．波の入射角が表面法線からあまりに逸れると，いくらかの反射が生じることがある[Jin, 1993]．その結果不要な振動が現れる．計算領域を打ち切るよりよいアプローチは完全整合層(PML)である．

PMLの基本的な考え方は，シミュレーション領域の周囲に非常に吸収力のある追加の層を加えるというものである．波が通常のシミュレーション領域を離れPMLに入ると同時に波は即座に減衰する．PMLの層は非常に吸収力の高い直交異方性不均質材料として理解することができる．直交異方性材料は非等方性材料のサブグループであり，材料の主な方向だけが異なり非対角方向は異ならないのでとなる．ここでは直交異方性材料は位置依存の誘電率テンソルで記述される．吸収媒体において，比誘電率は複素値である．ここでは媒体の偏光度，は材料の吸収特性である．2Dモデルでは，xおよびy方向についてだけを使う．

PMLの詳細については時間領域および周波数領域における音響学のチュートリアルをご参照いただきたい．

領域

PMLを利用するためには，PML領域である追加の層の長さで物理的領域を拡張する必要がある．この領域では，到来波は弱まる．

PML領域を持つ境界メッシュを生成し，境界要素マーカーを設定する．

PML領域と境界要素マーカーを可視化する．

境界1から4は前と同じであるが，境界2の領域をPML領域で切り取っているため境界2は吸収境界条件ではない．マーカー5の境界（緑）はPML領域の端であり，マーカー6の境界（オレンジ）は異なるPMLパラメータで3つの別々の領域を定義することが必要となる．

先で，一定の最大セル寸法を持つ一様な要素メッシュを作成したが，強度プロファイルから軸から遠い領域の強度は無視できることが分かった．MeshRefinementFunctionを使って，特定の座標の領域により細かいメッシュを作成し，PML領域についての領域マーカーを設定することができる．2DではMeshRefinementFunctionは要素面積に使うことができる．メッシュの三角要素面積は正三角形面積で近似することができる．制約を満足する面積はである．

可変の最大セル寸法および6つの領域で要素メッシュを作成する．

領域とPML領域を可視化する．

上では，Rasterizeを使うことによってグラフィックスのサイズを小さくすることができている．見栄えをよくしたい場合はRasterizeを削除するとよい．

パラメータの設定

PML層は対角テンソルなので，まだ(4)を使うことができる．物理的領域およびの内部で，最初の例と同じPDEを得る．PML領域は最大値の増加関数である．

PMLパラメータ

に対応する人為的な電気感受率テンソルでヘルムホルツ演算子を設定する．

すべての領域について吸収係数

を定義する．

要素メッシュの吸収係数を可視化する．

上では，Rasterizeを使うことによってグラフィックスのサイズを小さくすることができている．見栄えをよくしたい場合はRasterizeを削除するとよい．

関数のプロットでは，領域内部の電場を変更することなく，吸収係数が領域外の各方向において壁の働きをしていることが分かる．

境界条件

同じ境界条件を入射光の境界とスクリーンに使う．吸収境界条件は今は必要ない．

の境界要素マーカー4に開口部の放射条件を設定する．

境界要素マーカー3に吸収境界条件を設定する．

境界1，2，5，6はデフォルトのノイマン0値を持つ．

モデルの評価

モデルを設定し，PMLパラメータを置き換える．

PDEをNDSolveで解く．

可視化

数値結果をプロットする．

アニメーションの質を向上させる方法はこちらでご覧いただきたい．

伝搬方向に沿って結果を可視化する．

がに近付いても，PMLアプローチでは振動はない．

吸収境界条件アプローチの場合と比較する．

上から分かるように，2つのアプローチの間には何らかの変動誤差がある．これは吸収境界によって引き起こされる不要な反射の結果である．

モデルの検証

ここまで吸収境界条件アプローチと完全整合層(PML)アプローチの数値解を比較してきた．次に解を解析解に対して検証する．モデルを検証するために，グリーン(Green)関数法で推定された積分回折式を使う[Ishimarou, 1991]．

ここでは関心のある領域を囲む面，は系のグリーン関数，は法線導関数である．積分はに対する3つの異なる種類の境界条件を使って解くことができる．その3つとは，ディリクレ境界条件，ノイマン境界条件（それぞれタイプ1，タイプ2のゾンマーフェルト条件とも呼ばれる），ディリクレ・ノイマン境界条件（キルヒホッフの境界条件と呼ばれる）である[Wolf & Marchand, 1964]．キルヒッホフの境界条件は最も一般的に使われるが，想定される境界条件が解から引き出される訳ではないので，数学的に一貫していない[Wolf & Marchand, 1964]．この3つの境界条件はすべて遠い場では同じ結果を与える[Wolf & Marchand, 1964]．

実際には，2Dのスカラー回折モデルで説明しているように，この積分は表面上全体では行われない．境界だけが積分に関わっている．2Dではグリーン関数は階数0の第2種ハンケル(Hankel)関数である．さらにタイプ1のゾンマーフェルト条件を使う[Ishimarou, 1991]．理論的な解は以下になる：

ここではの任意の点と領域内の点の間のユークリッド距離である．開口部は線なので，法線は方向のベクトルである．

ハンケル関数の導関数を計算する．

解析解を計算するためにヘルパー関数を作成する．

における解析解と数値解の間の誤差を可視化する．

における伝搬方向の解析解と数値解の間の誤差を可視化する．

付近の理論解と数値解の間の差は，解析解が課された境界条件を引き出さないことが原因とみられる．例えば，ゾンマーフェルトのタイプ1積分は平面付近で発散する．また，開口境界条件は光の入射と吸収境界条件の和であることを思い出してほしい．放射境界条件におけるこの近似吸収境界条件項は，から開始していくらかの振動を引き起こす可能性がある．これらの振動は波が領域内深くに進むほど減少する．

解析解およびキルヒホッフ・ゾンマーフェルトアプローチの詳細は[Buchwald & Yeang, 2016][Totzeck, 1991][Heurtley, 1973]を参照されたい．

伝搬方向に沿った誤差プロットにおける高周波数の振動は，メッシュ要素の長さをに減らすことでなくすことができる．

3D領域における完全整合層

上のセクションでは2DのPMLを実装した．ここでは3D領域における2Dの円形開口部による波の回折を計算する．3D領域の計算は高価になる可能性があるので，系の対称性を使って領域の一部しか計算する必要がない場合に限定する．

パラメータの設定

入射スカラー電磁波パラメータをに設定し，3Dヘルムホルツ微分演算子を3DのPML領域についての直交異方性比誘電率で設定する．回折スカラー波はとラベル付けする．

3Dヘルムホルツ演算子とPML変数を定義する．

領域

まずPML領域を含む正の象限で立方体領域を設定し，OpenCascadeLinkを使って内側境界線を保存している領域をつなぐ．

内側面を保存している4つの立方体をつないで，必要な領域を得る．

内側境界線を持つオブジェクトの生成に関しての詳細はOpenCascadeLinkのチュートリアルに記載されている．

OpenCascadeオブジェクトを境界要素メッシュオブジェクトに変換し，開口部の場所を黄色で示した3D領域を可視化する．

緑の境界はPMLの立体を囲んでおり，青い境界は対称平面であり，赤い平面はとの境界のを含んでいる．黄色い円板は開口部の四分の一を表す．場合によっては開口部とスクリーンに別々の境界を生成することが困難なことがあるので，赤い平面に放射境界条件と吸収境界条件を同時に設定する関数を使う．

この場合，MeshRefinementFunctionはメッシュ要素の体積に依存するため，推奨される辺長を使う．TetrahedronElementの体積は，辺長という制約を満たすためにとして計算することができる．

異なるPMLに対する領域マーカーを作成し可視化する．

領域マーカーを使ってElementMeshを生成し，MeshRefinementFunctionを使って

軸を中心とするメッシュ要素の大きさを減少させる．

異なる要素サイズを使い，開口部の位置を黄色で示してElementMeshを可視化する．

画像の質を向上させる方法はこちらでご覧いただきたい．

メッシュは左隅で細かく，右上隅で粗くなっているのが分かる．黄色で描かれた開口部はメッシュの一部ではなく可視化のためだけに使われている．

3Dメッシュの可視化は複雑な場合がある．Manipulateを使うと3D領域と3DのElementMarkerを個別に可視化することが簡単にできる．このノートブックのサイズを減少させるために，2番目の要素メッシュは可視化の目的だけで生成されている．

領域マーカーに基づいて各領域の色を定義する．

可視化のためのメッシュを生成する．

領域の一つが動的に選べるManipulateを作成する．

開口部のElementMarkerを見付けるために，境界要素マーカーを可視化する．

境界ElementMarkerのリストを抽出する．

Manipulateを作成して，表面境界を動的に調べる．

境界条件

先で境界と境界を分割したが，より複雑な立体や3D領域を扱うときに適切な境界を設定することは難しい場合がある．このような場合に適したアプローチとして，条件とを区別する関数を使うというものがある．方程式(5)から，をに設定すると，方程式(6)を取り戻すことができるので，開口部で，スクリーンでとなる関数で条件との両方を設定することができる．

次のステップでは直径の円形開口部の関数を作成する．

開口部にHeavisidePi使い，それを可視化する．

これで統一された境界条件が生成できる．

ElementMarker 9の放射条件を設定する．

ElementMarker 1および5の表面境界は，この領域について対称の表面であり，ノイマン0値でモデル化できる．これは，境界に他の指定がない場合のデフォルトの境界条件である．

完全整合層のパラメータの設定

領域マーカーが分かったので，完全整合層の吸収係数が設定できるようになった．

各完全整合層の各方向に対する吸収係数

を設定する．

領域すべてにおける吸収係数

を構築する．

モデルの評価

モデルの完全整合層パラメータと境界条件を置き換える．

NDSolveでPDEを解き，時間と使用メモリを測定する．

可視化

数値解をプロットする．

アニメーションの質を向上させる方法はこちらでご覧いただきたい．

結果を2DのDensityPlotとして可視化する．

アニメーションの質を向上させる方法はこちらでご覧いただきたい．

密度プロットは円形開口部の特徴的なリング状のパターンを示している．伝搬距離が長い場合，このパターンは0次の第一階ベッセル関数に進化する．

モデルの検証

3Dでは，系のグリーン関数は形式の球面波である．方程式(7)のグリーン関数を置換すると，次が得られる：

ここではの任意の点と領域内部の点の間のユークリッド距離である．解析解を利用するために，ヘルパー関数を設定する．

直線

および

における理論解についての関数を作成する．

における解析解と数値解を可視化する．

における解の誤差をプロットする．

まとめると，吸収境界条件は問題を素早く設定するには役立つが，正確さが必須の場合は完全整合層を使った方がよい．

用語集

記号	解説	単位
E	電場	[V/m]
λ_n	波長	[m]
k_n	波数	[1/m]
μ_r	比透磁率	N/A
ϵ_r	比誘電率	N/A
L	伝搬領域の長さ	[μm]
H	伝搬領域の幅	[μm]
d	開口サイズ	[m]
Γ_screen	不透明なスクリーン境界	N/A
Γ_aperture	開口境界	N/A
Γ_ABC	吸収境界	N/A
Ω	計算領域	N/A

参照

1. Born, M. and Wolf, E.(1999). Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light (7th ed.). Cambridge: Cambridge University Press.

2. Jackson, John D. Classical Electrodynamics. New York :Wiley, 1999.

3. Jin, Jian-Ming. The finite element method in electromagnetics. John Wiley & Sons, 2015.

4. A. Ishimarou, Electromagnetic Wave Propagation, Radiation, and Scattering, Prentice Hall, 1991.

5. Wolf, E., & Marchand, E. W. Comparison of the Kirchhoff and the Rayleigh–Sommerfeld Theories of Diffraction at an Aperture. Journal of the Optical Society of America, 1964, 54(5), 587. doi:10.1364/josa.54.000587

6. Buchwald, J.Z., Yeang, CP. Kirchhoff’s theory for optical diffraction, its predecessor and subsequent development: the resilience of an inconsistent theory. Arch. Hist. Exact Sci. 70, 463–511 (2016). https://doi.org/10.1007/s00407-016-0176-1

7. Totzeck, M. (1991). Validity of the scalar Kirchhoff and Rayleigh–Sommerfeld diffraction theories in the near field of small phase objects. Journal of the Optical Society of America A, 1991, 8(1), 27. doi:10.1364/josaa.8.000027

8. Heurtley, J. C. Scalar Rayleigh–Sommerfeld and Kirchhoff diffraction integrals: A comparison of exact evaluations for axial points*. Journal of the Optical Society of America, 1973, 63(8), 1003. doi:10.1364/josa.63.001003

トップへ