波動光学シミュレーションにおけるビームエンベロープ法の使い方

2018年 1月 8日

波動光学の分野では, マックスウェルの方程式を厳密に解く方法で大規模な光学システムをシミュレートすることは困難です. これは, システム内に現れる波を十分に細かいメッシュで解像する必要があるためです. この目的のためには, COMSOL Multiphysics® ソフトウェアのビームエンベロープ法が1つの選択肢となります. このブログでは, 電磁波 (ビームエンベロープ)インターフェースの使用方法とその制限の対処方法についてご説明します.

大きな波動光学モデルの求解方法の比較

電磁気シミュレーションでは, マックスウェル方程式の正確な解を見つけるために, 波長は常にメッシュによって解像される必要があります. この要件により, 波長に比べて大きいモデルをシミュレートすることが困難になります. 定常光学問題では, いくつかの手法を使って大きなモデルを扱うことができます. これらの方法には, Fraunhofer, Fresnel-Kirchhoff, Rayleigh-Sommerfeld 回折式などのいわゆる回折式と, 近軸 BPM や角度スペクトル法などのビーム伝搬法 (BPM) が含まれます(参照1).

これらの手法のほとんどは, ヘルムホルツ方程式の近似式を使用します. これらは, ある平面の既知の電場から別の平面の電場を求解する伝播方法に基づいているため, 大きなモデルを処理できます. そのため, ドメイン全体をメッシュ化する必要はなく, 必要な平面の 2D メッシュがあればよいのです.

これらの方法と比較して, COMSOL Multiphysics の電磁波 (ビームエンベロープ) インターフェース(以下, ビームエンベロープインターフェース) は, ドメイン内のヘルムホルツ方程式の厳密な解を求めます. 大規模なモデルにも対応できます. つまり, 特定の制限が満たされている場合, メッシュ要件を大幅に緩和できるということです.

ビームエンベロープ法によるレンズシミュレーション.

1μm の波長のビームに対して, ミリメートル範囲の焦点距離を持つレンズのビームエンベロープシミュレーション.

ビームエンベロープインターフェースについて, 以下で詳しく説明します.

ビームエンベロープインターフェースの背後にある理論

ビームエンベロープインターフェースの内部での計算を見てみましょう. このインターフェースをモデルに追加し, フィジックスインターフェースノードをクリックして, 位相指定のタイプをユーザー定義に変更すると, 数式セクションに以下のように表示されます.

(\nabla-i \nabla \phi_1) \times \mu^{-1}_r (( \nabla-i \nabla \phi_1) \times {\bf E1}) -k_0^2 \left (\epsilon_r -\frac{j \sigma}{\omega \epsilon_0} \right ) {\bf E1}

ここで, \bf E1はインターフェースが解く従属変数で, エンベロープ関数と呼ばれます.

電場のフェーザー表現では, \bf E1が振幅に, \phi_1が位相に相当し, すなわち

{\bf E} = {\bf E1} \exp(-i \phi_1).

となります. 最初の式は, ビームエンベロープインターフェースの支配方程式であり, 電場の第2の定義をヘルムホルツ方程式に代入することで導き出されます. \phi_1がわかっていれば, 未知数は\bf E1だけとなり, それを解くことができます. 問題を求解するには, 位相\phi_1事前に指定しなければなりません.

2番目の式では, 高速振動部分である位相を場から因数分解できるような形式を想定しています. そうすると, エンベロープ\bf E1はゆっくりと変化するので, 波長を解像する必要はありません. 代わりに, エンベロープの緩やかな波を解像するだけで済みます. このプロセスにより, 大規模な波動光学問題をパソコンでシミュレーションすることが可能になったのです.

場自体ではなく, エンベロープの解が欲しいのはどのようなときか?という質問がよくあります. レンズシミュレーションはその一例です. 場合によって, 複素電場ではなく, 強度を解く必要があります. ここで実は, エンベロープのノルムの2乗で強度が得られます. このような場合には, エンベロープ関数を求めれば十分です.

位相関数が正確にわからない場合は?

ビームエンベロープ法の背後にある数学には, さらなる疑問があります.

  • 位相が正確に分からない場合はどうするのか?
  • そのような場合でも, ビームエンベロープインターフェースを使用できるのか?
  • その結果は正しいのか?

これらの疑問に答えるためには, もう少し計算する必要があります.

1次元の例

最も簡単な例として, Ez = \exp(-i k_0 x)という平面波が, k_0 = 2\pi / \lambda_0, 波長\lambda_0= 1 umで, 長さ20um の長方形のドメインを伝搬する場合を考えてみましょう. (ここでは, 説明のために意図的に短いドメインを使用しています.)

面外波は左側の境界から入り, 反射せずに右側の境界を透過します. これは, ビームエンベロープインターフェースで, 左に励振あり, 右に励振なしの適合境界条件を追加し, 上下に完全磁気伝導体境界条件を追加することでシミュレーションできます (y方向は気にしないという意味です).

位相仕様の正しい設定を下図に示します.

COMSOL Multiphysics GUI で波動ベクトルの設定を示したスクリーンショット.

事前に正しい位相関数 k_0 x, または波動ベクトルが (k_0,0) であることがわかっているので, 答えは Ez = \exp(-i k_0 x) となります. 位相関数を2つ目の式に代入すると, 逆に定数関数である E1z = 1 が得られます.

定数関数を求解するために必要なメッシュ要素の数は? たった1つです. (高周波モデリングに関する以前のブログ参照.)

次の結果は, エンベロープ関数 \bf E1 と, \bf E のノルム ewbe.normE(これは |{\bf E1}| とも等しい) を示しています. ここでは, 予想通り, 任意のメッシュ数に対して正確な位相関数, 定数1 を与えれば, 正しいエンベロープ関数が得られることがわかります. 確認のため, \bf E1z の位相である arg(E1z) もプロットしています. これも予想通り, ゼロです.

正しい位相関数のためのエンベロープ関数, 電場ノルム, エンベロープ関数の位相の3つの結果.

さまざまなメッシュサイズに対して計算された, 正しい位相関数 k0x のエンベロープ関数 (赤), 電場ノルム (青), およびエンベロープ関数の位相 (緑).

ここで, 位相関数の推測が少しずれている場合にどうなるかを見てみましょう. たとえば, 正確な (k_0,0)ではなく, (0.95k_0,0) とします. どのような解が得られるでしょうか? 見てみましょう.

誤った位相関数に対するエンベロープ関数, 電場ノルム, エンベロープ関数の位相の3つの結果.

さまざまなメッシュサイズに対して計算された, 誤った位相関数 0.95 k0x のエンベロープ関数 (赤), 電場ノルム (青), およびエンベロープ関数の位相 (緑).

ここでエンベロープ関数に見られるのはいわゆるうなりです. すべてがメッシュサイズに依存していることは明らかです. 何が起こっているのかを理解するには, まず鉛筆と紙と忍耐力が必要です.

答えが Ez = \exp(-i k_0 x) であることはわかっていましたが, COMSOL® ソフトウェアで意図的に誤った位相見積もりを入力してみました. 2番目の式に誤った位相関数を代入すると, \exp(-i k_0 x)={\bf E1z} \exp(-0.95i k_0 x) が得られます. この結果, {\bf E1z} = \exp(-0.05i k_0 x) となり, 定数1ではなくなります. これは \lambda_b= 2\pi/0.05k_0=20um の, いわゆるうなり波長と呼ばれる, 波長を持つ波です.

上のプロットを, 6つのメッシュ要素について見てみましょう. 予想通りの結果 (赤線), つまり {\bf E1z} = \exp(-0.05i k_0 x) が得られました. プロットは自動的に実数部を取り, {\bf E1z} = \cos(-0.05 k_0 x) と表示します. 低解像度のプロットでも, エンベロープ関数の近似解が示されています. 有限要素シミュレーションで予想される通り, メッシュが粗いほど近似が粗くなるという結果が得られるということです.

これによって, 位相関数の推測を誤ると, 誤った (うなりの混じった) エンベロープ関数が得られることがわかります. 誤った推測により, エンベロープ関数はうなりの位相 (緑の線) を加えたものになっており, これは -0.05 k_0 x となります.

\bf E のノルムはどうでしょうか? 上のプロットの青い線を見てください. COMSOL Multiphysics ソフトウェアは ewbe.normE の正しい解 (定数1) を生成しているように見えます. 計算してみましょう. 誤った (解析的な) 位相関数と誤った (うなりが混ざった) エンベロープ関数の両方を代入すると, {\bf Ez} = \exp(-0.05i k_0 x) \times \exp(-0.95i k_0 x) = \exp(-i k_0 x) となりますが, これは正しい高速振動電場です.

\bf E のノルムを取ると, 正しい解である定数1が得られます. これは期待通りの結果です. ドメインが大きい場合には, \bf E 自体を表示することができないこともあることにご注意ください. ただし, \bf E を解析的に求め, 粗いメッシュで \bf E のノルムを表示することは可能です.

これには理由があります. 実際, 位相関数がずれていると, うなりが混ざり, エンベロープ関数もずれてしまうことがわかります. しかし, それでも電場のノルムは依然として正しいのです. したがって, 正しい電場を得るためには, うなりの混ざったエンベロープ関数を正しく計算することが重要です. 上記のプロットは, そのことを明確に示しています. 6要素のメッシュケースでは, うなりの混ざったエンベロープ関数を完全に解像しているため, 完全に正しい電場ノルムが導かれています. 他のメッシュは, メッシュサイズに応じて, うなりの混ざったエンベロープ関数の近似解を出しています. また, 場のノルムについても同様です. これは, 任意のケースに当てはまる一般的な結果です.

COMSOL Multiphysics で使用する位相関数に関係なく, \bf E1 の最初の方程式を正しく解く限り, また位相関数がドメイン全体で連続している限り, 問題ありません. 1つのドメインに複数の材料がある場合, 位相関数の連続性も解の精度に大きく影響します. これについては, 今後のブログでご紹介するかもしれませんが, 高周波モデリングに関するこの前のブログでも解説しています.

2次元の例

これまで, スカラー波数について説明してきました. より一般的には, 位相関数は波数ベクトルによって指定されます. 波数ベクトルが正しく推測されない場合, 結果はベクトル的な影響を受けます. 最初の例と同じ平面波を見てみましょう. ここで, 位相を誤って推測した場合, つまり k_0 x ではなく, k_0(x \cos \theta + y \sin \theta) と推測した場合を考えます. この場合, 波数は正しいですが, 波数ベクトルがずれています. 今回は, うなりは2 次元で起こります.

まず, 1次元の例と同じ計算を行ってみましょう. \exp(-i k_0 x)= {\bf E1z}(x,y) \exp(-i k_0 (x \cos \theta+y \sin \theta) ) となり, エンベロープ関数は {\bf E1z}(x,y) = \exp(-i k_0 (x (1-\cos \theta) -y \sin \theta) ) と計算されます. これは, (1-\cos \theta, -\sin \theta) 方向に伝搬する傾斜波で, うなり波数は k_b = 2 k_0/\sin (\theta/2), うなり波長は \lambda_b=\lambda_0/(2\sin (\theta/2)) となります.

次のプロットは, さまざまな最大メッシュサイズでの3.8637 um x 29.348 um のドメインでのθ = 15° の結果です. 前述の1次元の例のケースと同じ境界条件が与えられています. 唯一の違いは, 左の境界の入射波が{\bf E1z}(0,y) = \exp(i k_0 y \sin \theta) であることです. (位相の推測が誤っているため, 対応する誤った境界条件を与えなければならないことにご注意ください.)

一番細かいメッシュの結果 (右端) では, \bf E1z が先程解析した通りに計算され, \bf Ez のノルムも定数1になるように計算されていることが確認できます. これらの結果は, 1次元の例の場合と一致しています.

誤った位相関数に対する電場ノルムとエンベロープ関数の異なる結果.

異なるメッシュサイズに対して計算された, 誤った位相関数 k_0(x \cos\theta +y \sin\theta ) の電場ノルム (上) とエンベロープ関数 (下). 色の範囲は, -1から1までの値を表します.

ビームエンベロープインターフェースによるレンズのシミュレーション

ここでの最終的な目標は, ビームエンベロープインターフェースを使用して, ミリメートル単位のドメインで光学レンズを通過する電磁ビームをシミュレートすることです. どのようにしてこれを実現するのでしょうか? 正しい解を計算する方法はすでにご説明しました. 以下の例は, 曲率半径500 um, 屈折率 1.5の平凸レンズ (約1mm の焦点距離) に, ハードアパーチャーのフラットトップ入射ビームを入射させた場合のシミュレーションです.

ここでは \phi_1 = k_0 x としていますが, これは全く正確ではありません. レンズの前のドメインでは反射があり, 干渉が発生しています. レンズの中では, 複数の反射があります. レンズを過ぎると, 位相は球形になり, ビームは一点に集束します. つまり, この位相関数は, レンズ周辺で起きていることとは大きく異なるのです. しかし, 手がかりはあります. \bf E1z をプロットしてみると, うなりがかかっているのがわかります.

レンズ内のうなり波長のシミュレーション.

\bf E1zのプロット. 挿入図は, レンズ内の最も細かいうなり波長を示しています.

プロットを見ると分かるように, レンズに顕著なうなりが発生しています (挿入図参照).

注: バージョン5.5以降, レンズ表面に遷移境界条件を導入し, 反射防止膜を模倣することで, 反射によるうなりを解消しています. (以前のブログをご参照ください. また, 反射防止膜を施していても, 法線方向の入射にしか効果がないため, メッシュの微細化が必要ということにご注意ください. レンズの端の周辺の場には, 若干のうなりが見られる場合があります.

実際には, 最も細かいうなり波長は, レンズの前では\lambda_0/2です. これを証明するには, これまでの例と同じ計算をするだけです. 最も細かいうなり波長は, 入射ビームと反射ビームの干渉によるものですが, 順方向の伝搬には寄与しないので, これは無視することができます. また, メッシュがレンズ前のうなりを解像していないことがわかりますが, とりあえずこれは無視しておきましょう.

レンズ内のうなり波長は, n= 1.5の場合, 後方ビームでは 3\lambda_0/2, 前方ビームでは 2\lambda_0 となり, これも前の例と同じ方法で証明することができます. ここでも, 後方のビームは無視します. プロットで見えるのは, 前方ビームの 2\lambda_0 のうなりです. 後方ビームはほんの一部 (n= 1.5の場合, 入射ビームの約4%なので, 見えない) に過ぎません. 次の図は, 10個のメッシュ要素を使用してレンズ内のうなりを解像するメッシュを示しています.

COMSOL Multiphysics で作成されたレンズ内のうなり波長のメッシュ.

レンズ内部のうなり波長. メッシュが10個のメッシュ要素でうなりを解像します.

レンズ内の伝搬ビームのうなり以外では, 後続の空気ドメインでのうなりはかなり大きいので, ここでは粗いメッシュを使います. しかし, より短焦点距離のレンズでは, 2次位相が速くなり, うなり波長が非常に短くなることがあるため, これは当てはまらないかもしれません. この例では, 最速のうなりを解像するために, レンズドメインでのみより細かいメッシュを使用する必要があります.

計算された場のノルムは, このブログの上部に示されています. 結果を検証するために, 周波数領域インターフェースを使用してレンズ出口面での場を計算し, 次にフレネル回折式を使用して焦点での場を計算します. ノルムの結果は非常によく一致しています.

ビームエンベロープインターフェースとフレネル回折式を比較した1D プロット.

ビームエンベロープインターフェースとフレネル回折式の比較. メッシュは10個のメッシュ要素でレンズ内のうなりを解像しています.

次の比較は, メッシュサイズの依存性を示しています. 推奨値 \lambda_b/6 に等しい \lambda_0/3 を使用すると, かなり良い結果が得られます.

場のノルムのメッシュサイズ依存性を示す1次元プロット.

焦点における場のノルムのメッシュサイズ依存性.

COMSOL® ソフトウェアのバージョン5.3a 以降, フレネルレンズのチュートリアルモデルに, ビームエンベロープインターフェースを使った計算が含まれています. フレネルレンズは通常, 非常に薄いです (波長オーダー). レンズ表面の不連続部とその周辺に回折があったとしても, レンズ部分の周りの細かいメッシュはメッシュ要素の総数に大きな影響を与えません.

最後に

このブログでは, ビームエンベロープインターフェースが内部で何をするのか, そして波動光学問題の正確な解をどのように得ることができるのかについて説明しました. たとえうなりが発生したとしても, うなり波長は波長よりもはるかに長くなる可能性があるため, 大規模な光学系をシミュレートすることができます.

うなりを解像するためにメッシュサイズを確認するのは面倒に思えますが, この余分な作業はビームエンベロープインターフェース以外の手法でも必要となります. 有限要素法を使用する場合にも, 正確に計算された解のメッシュサイズ依存性を常に確認する必要があります.

次のステップ

お試しください: 下のボタンをクリックすると, ミリメートル範囲の焦点距離レンズのシミュレーションファイルをダウンロードできます.

参照

  1. J. Goodman, Fourier Optics, Roberts and Company Publishers, 2005.

コメント (0)

コメントを残す
ログイン | 登録
Loading...
COMSOL ブログを探索