pyGDMで散乱光分布を計算する 〜DDA法による電磁場解析〜

導入

散乱光分布を計算していて、ふとこれレンズで結像するとどう見えるんだろう？と思ったので計算してみたシリーズの第5回。

今回は散乱光の計算方法について書きます。
PythonライブラリのpyGDM2を用いた計算方法についてです。
ここでの計算結果を次回以降でも活用します。

過去記事

散乱光分布の状況

こんな状況で散乱光の分布を求めてみます。

• シリコンの上に直径200nmのガラス球
• 波長532nmのレーザー光を斜め45度方向から照射

計算方法

散乱光分布を求めるには、FDTDや解析計算を使うなどいくつか方法があります。
今回用いるpyGDMの特徴とメリット・デメリット書きます。

• pyGDM公式ページ
  https://homepages.laas.fr/pwiecha/pygdm_doc/
• pyGDMの論文1
  https://www.sciencedirect.com/science/article/pii/S001046551830225X?via%3Dihub
• pyGDMの論文2
  https://www.sciencedirect.com/science/article/pii/S001046552100254X?via%3Dihub

pyGDMの特徴

DDA(Discrete Dipole approximation)と呼ばれる計算方法があります。
これは散乱体を離散的な双極子(Dipole)を並べることで表現して計算します。
pyGDMはDDAの計算をする際にGreen関数を用いて各dipole間の電場の相互作用を計算します。

pyGDMのメリット

• Pythonだけで計算できる
  通常のDDAはLINUX環境が必要だったり、fortranの知識が必要になったりします。地獄。
• 任意形状を計算可能
• 比較的高速(数秒~数分)に計算可能
• 電磁場解析なので、位相も計算可能

pyGDMのデメリット

• dipole間の距離の設定に知識が必要
• 1µm超えるような大きい形状はメモリ不足で計算出来ないことが多い

pyGDMのインストール

pipでインストール可能です。

pip install pyGDM2

pyGDM計算コード

散乱体のオブジェクト、照明条件のオブジェクト、周辺環境の屈折率のオブジェクトを作成します。
これら3つを渡してSimオブジェクトを作成してDDA計算を実行します。

まず球状の散乱体のオブジェクトを作成します。

from pyGDM2  import structures, materials, fields, propagators, core, visu, tools,linear
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.path import Path

#ラジアンと度数の変換用定数
deg=np.pi/180

#最初に使う物質の誘電率(屈折率)を設定、球の形状を作成
#誘電率を指定する
si = materials.silicon()
sio2 = materials.sio2()
air = materials.dummy(1.0)
n1 =si
n2 = air

#散乱体の形状作成
#dipoleの距離
step = 20 #nm
#球の半径
radi = 100 #nm

#球形状を表すdipoleの座標をnumpyアレイを生成
geometry = structures.sphere(step, R= radi / step, mesh ='cube')

#球形状を表すオブジェクトの作成,引数で物質の誘電率を指定
struct = structures.struct(step, geometry, sio2)

この時に有効なdipole数が表示されます。以下だと619個

structure initialization - automatic mesh detection: cube
structure initialization - consistency check: 619/619 dipoles valid

今回であればこれで十分な数です。
PC次第ですが7,000 ~ 10,000個前後でメモリ不足で計算できなくなったりし始めます。

次に照明条件のオブジェクトを作成します。P偏光の平面波を設定します。
入射角は0度が下から上の方向に対応するため、180+45度としています。

#平面波
field_generator = fields.plane_wave
#波長
wavelengths = [532]
#偏光、入射角をdictで指定
kwargs = dict(E_p=1, E_s =0, inc_angle = [180+45],theta=0)

#条件を渡して照明条件のオブジェクトを作成
efield = fields.efield(field_generator, wavelengths = wavelengths, kwargs = kwargs)

周辺環境の屈折率を設定します。
光ファイバーを想定して垂直方向に3層設定できます。
今回は最下層がシリコン、中間層は真空(屈折率1)とします。3層目は設定しなくても計算できます。

dyads = propagators.DyadsQuasistatic123(n1 = n1, n2=n2)

simオブジェクトを作成します。計算する前に作成した形状を図示して確認できます。
計算を実行、近接場光が求まります。

sim = core.simulation(struct, efield, dyads)

im = visu.structure(sim,show=False)
ct = visu.structure_contour(sim,show = True)

黒の四角が双極子を表します。
直径が200nmとしましたが、合っていることが分かります。

計算を実施、電場分布を可視化してチェックできます。

sim.scatter()

#電場の実部をプロットする
qv = visu.vectorfield_by_fieldindex(sim, 0,projection='XY', show=False)
ct = visu.structure_contour(sim,  show=True)

左側から照明光くるため、前方の電場強いようです。
この図ではシリコンウェーハを上から見た場合ですが、projection引数を"YZ","ZX"などに変えて見る方向を変えられます。

上で求めたのは近接場光と呼ばれるものです。
近接場光は光の波長よりも短い範囲における電場分布です。

求めたいのは散乱光です。
近接場光は伝播に伴う減衰が早く、私たちが観察する散乱光に殆ど寄与しません。そのため遠方解を別途求める必要があります。

以下では100µm離れており、散乱角方向90点、アジマス角が等間隔に360点、
90 x 360= 32400点の電場が計算されます。

#散乱角のサンプリング点数
Nteta = 90
#アジマス角のサンプリング点数
Nphi = 360

E = linear.farfield(sim, field_index=0,return_value='efield',tetamin=5*deg,Nteta=Nteta,Nphi=Nphi,r=100_000.0 )

Eが計算結果で、要素数5のリストです。
E[0],E[1],E[2]のみを以下の解析で使います。

• E[0] : theta(散乱角)のメッシュグリッド
• E[1] : phi(アジマス角)のメッシュグリッド
• E[2] : 散乱された電場

試しに1点の値を確認してみます。

print(E[2][0])

[-5.315976e-04-1.7147543e-04j -5.577258e-13-1.0330385e-12j4.650876e-05+1.5002157e-05j]

ある一点における電場のx,y,z成分です。
複素数なので位相まで含めた情報が得られました。

散乱光分布の可視化

散乱光の分布を確認してみます。
E[2]の形状に注意してパワーに変換、
瞳面上での分布をコンター図にします。

plt.figure(figsize=(6,5))
plt.contourf(np.sin(E[0])*np.cos(E[1]),np.sin(E[0])*np.sin(E[1]),np.sum(np.abs(E[2])**2, axis=1).reshape(E[0].shape),levels=100)
plt.xlabel("azimuth angle(deg)")
plt.ylabel("elevation angle(deg)")
plt.gca().set_aspect('equal')
plt.show()

散乱角を5度から90度としたので、真ん中に穴が空いていますが
照明光は左からくるので、前方への散乱が強いと分かります。
照明波長と球の大きさが同じくらいなのでMie散乱の特徴が現れ始めているものと思います。

注意事項

dipoleの距離の設定が重要です。
波長くらいまで離れるとうまく干渉を計算できませんし、距離を近くしすぎるとdipole数が増えて計算不可になったり、dipole数が問題なくても計算が不安定になったりします。
距離の決め方はDDAの基準と同じです。

pyGDMの計算方法ではdipoleをすべて記憶する必要があります。
必要なメモリ数はdipole数に比例、構造体の一辺の長さの3乗に比例して必要なメモリが増えます。
感覚的に1µmくらいの大きさが限界です。

dipoleの距離はx,y,z方向は全部同じなのでアスペクト比が大きい構造だとdipole数が増えやすいです。
dipole距離を大きくすると、一辺が短い方向では構造を再現できなくなるためです。

次回

pyGDMはこれ以外にも多くの機能があり、散乱光以外でも様々な電磁場解析が可能です。
HPに解析例が載っているので図を見るだけでも面白いです。

散乱光の求め方はここまでにします。
次回はFFT結像に偏光を取り込んでいきたいと思います。

• 次回記事(最終回)