はじめに
散乱光分布を計算していて、ふとこれレンズで結像するとどう見えるんだろう?と思ったので波動光学に基づいて計算してみたシリーズの第3回。
本記事では波動光学から少し離れてpyGDM2を用いて、シリコン基板上の200nmガラス球からの散乱光分布を計算します。近接場光の計算から遠方散乱光への変換、散乱光分布の可視化まで一通り実施します。
ここでの計算結果を次回以降でも活用します。
過去記事
問題設定
こんな状況で散乱光の分布を求めてみます。
- シリコンの上に直径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
特徴
DDA(Discrete Dipole approximation)と呼ばれる計算方法があります。これは散乱体を離散的な双極子(Dipole)を並べることで表現して計算します。pyGDMはDDA系の手法であり、Green関数を用いて双極子間相互作用を扱う点が特徴です。
メリット
- Pythonだけで計算できる
通常のDDAはLINUX環境が必要だったり、fortranの知識が必要になったりします。地獄。 - 任意形状を計算可能 球以外も計算できます。
- 比較的高速(数秒~数分)に計算可能
- 電磁場解析なので、位相も計算可能
デメリット
- dipole間の距離の設定に知識が必要
- 1µm超えるような大きい形状はメモリ不足で計算出来ないことが多い
- 散乱光を特定開口で積算しようとすると自前で計算する必要がある
インストール
pipでインストール可能です。
pip install pyGDM2計算コード
散乱体のオブジェクト、照明条件のオブジェクト、周辺環境の屈折率のオブジェクトを作成します。
これら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でも数秒で計算できます。
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散乱と呼ばれ前方への散乱が強くなります。この例では、異物大きさが波長に対して無視できない領域なので、Mie散乱の特徴が現れ始めているものと思います。

注意事項
dipoleの距離の設定が重要です。
波長くらいまで離れるとうまく干渉を計算できませんし、距離を近くしすぎるとdipole数が増えて計算不可になったり、dipole数が問題なくても計算が不安定になったりします。
距離の決め方はDDAの基準と同じです。pyGDMの計算方法ではdipoleをすべて記憶する必要があります。
必要なメモリ数はdipole数に比例、構造体の一辺の長さの3乗に比例して必要なメモリが増えます。
感覚的に1µmくらいの大きさが限界です。dipoleの距離はx,y,z方向は全部同じなのでアスペクト比が大きい構造だとdipole数が増えやすいです。
dipole距離を大きくすると、一辺が短い方向では構造を再現できなくなるためです。
まとめ
pyGDMはこれ以外にも多くの機能があり、散乱光以外でも様々な電磁場解析が可能です。HPに解析例が載っています。
散乱光の求め方はここまでにします。次回は今回の計算結果を使って、偏光の効果を取り込みつつ結像計算を行いたいと思います。
参考文献
より体系的に学びたい場合は、電磁気学や計算電磁場解析の書籍も参考になります。
短いですが、pyGDMでなくDDAの章があります。
次回
今回求めた散乱光を使って結像計算をしてみたいと思います。





コメント