tidy3D: Kretschmann配置によるSPPの励起

内容

• 新奇な内容: for loopでの変数探索
• 2D simulation
• 扱った現象: kretschmann配置におけるSPPの励起

目的

本記事では前回(https://qiita.com/linear_tree/items/e95c2c1ce3dd5c7e2db2 )に引き続き、Tidy3D (https://www.flexcompute.com/tidy3d/solver/ )を用いたsimulationを行う。
基本的には前回と同様にscriptを用いたシミュレーションの一例と位置付けられる。
前回は変数を全て固定したうえでmodelを作成したために構造の条件を変えるためには手動で値を変更する必要があったが、FDTD simulationを実行する際にはしばしば最適値の探索のために自動で変数を少しずつ変えたsimulationを実行したい場合がある。この動作は一般にはparameter sweepと呼ばれ、多くのソフトウェアに機能が実装されており、tidy3dのGUI versionにもparameter sweep機能はついている。

今回はこのsweepをscriptで行うため、for loopでの変数指定と評価に用いる値を格納するためのcodeを用いる。なお、tidy3Dは一月あたりにの使用可能な計算コストが決まっている。使用プランには依存するがFree planおよびacademic annual planで使用可能な一月あたり10 creditは複雑な構造を大量に回すとすぐに枯渇してしまう量であるため、今回のようなテスト用構造ではできる限り単純な構造にしておくことが好ましい。このような理由から、今回のシミュレーションでは単純多層膜を対象とした2D simulationを行った。

対象

対象とした物理現象は、COMSOLやLumericalなどの電磁界シミュレーションソフトでしばしばチュートリアルに用いられている表面プラズモンポラリトン(Surface Plasmon Polariton: SPP)の励起とした。
絶縁体基板上に厚み50 nmのAu層が存在しており、さらに上に100 nmの絶縁体層薄膜がある構造において基板側から平面波を斜め入射する典型的なkretschmann配置を想定する。Kretscmann配置では励起されるSPPの波長は構造の構成及び励起光の入射角度に依存するため、単色波レーザーを用いてSPPを励起するには特定の角度で入射する必要がある。ここでは基板の材質や膜厚などの構造は固定し、入射角度だけを変数とすることで指定した波長においてSPPを励起する入射角度の探索を行う。
以下、使用したcodeと簡単な説明を記載する。

使用するlibralyなどのimport

import h5py

import numpy as np
import matplotlib.pyplot as plt
import xarray as xr

import tidy3d as td
from tidy3d import web
from tidy3d import material_library

変数の決定

使用する変数などを準備。今回は以下の値を事前に設定。

• 対象とする波長(周波数)
• monitorで記録する波長(周波数)幅、間隔
• FDTD領域の大きさ
• 波源の位置
• monitorの位置
• 計算時間
• 金属層および絶縁体層の厚み
• 各構造の誘電率

なお、構造に対して十分に大きいサイズとしてinf_effに適当な値を置いている。
前回は長方形構造を
geometry = td.Box(center=(0, 0, 0), size=(rect_x, rect_y, rect_z)),
のように指定していたが、多層膜構造のような構造は境界値を指定する方が定義しやすいため
geometry=td.Box.from_bounds(rmin=(-inf_eff, -inf_eff, Metal_t),rmax=(inf_eff, inf_eff, Metal_t+diel_t)),
のように定義している。今回もxy方向は周期的な構造としているため、時xy方向の値がFDTD領域より大きければ無限に一様な構造として扱われる。

lmbda0 = 0.63
freq0 = td.C_0 / lmbda0
min_wl = 400  # nm
max_wl = 900  # nm
step_wl = 5  # nm
wavelengths = np.arange(min_wl, max_wl + step_wl, step_wl) / 1e3
freqs = [td.C_0 / (wl) for wl in wavelengths]
#geometry for simulation
FDTD_x = lmbda0*2
FDTD_y = 0
FDTD_z = 4
posi_source = -1.5
posi_monitor = 1.4
t_run = 1e-12

#geometry for object
Metal_t = 0.05
diel_t = 0.1

inf_eff = 10

gold = material_library['Au']['RakicLorentzDrude1998']
subst = td.Medium(permittivity=1.43**2)    #基板
coat = td.Medium(permittivity=1.25**2)     #絶縁体層

関数を作製

Simulationを実施する構造の作成を関数化する。ここで引数として励起光の入射角度を選択することでforloopで探索ができるようになる。
今回の例では
source_x = td.PlaneWave(
center=(0, 0, posi_source),
size=(td.inf, td.inf, 0),
source_time=pulse_x,
direction='+',
pol_angle= 0,#np.pi/2,
angle_theta=np.deg2rad(inc_angle),
angle_phi=0,
name="source_x"
)
におけるangle_theta=np.deg2rad(inc_angle),
が変数部分に相当するため、def mkobjct(inc_angle):のように引数をinc_angleとしておく。
[構造][波源][モニター]を個別に関数化することも可能であることは確認できたが今回はそこまで複雑な構造ではないためまとめてしまった。より複雑な構造を作成したい時などは分離した方が整理しやすい。

今回も透過/反射成分の確認用にflux monitorを二つ、電磁場の空間分布をみるためにfield monitorを設置した。最後に金属平面近傍のみの電磁場の演算を行いたかったため、全景を見るための大きなmonitor(monitor_xy)と絶縁体層領域のみをカバーする小さなmonitor(monitor_d)を二つ設定している。

def mkobjct(inc_angle):
    
    substrate = td.Structure(
        geometry=td.Box.from_bounds(rmin=(-inf_eff, -inf_eff, -inf_eff),rmax=(inf_eff, inf_eff, 0)),
        medium = subst,
    )
    film_m = td.Structure(
        geometry=td.Box.from_bounds(rmin=(-inf_eff, -inf_eff, 0),rmax=(inf_eff, inf_eff, Metal_t)),
        medium   = gold,
    )
    film_d = td.Structure(
        geometry=td.Box.from_bounds(rmin=(-inf_eff, -inf_eff, Metal_t),rmax=(inf_eff, inf_eff, Metal_t+diel_t)),
        medium   = coat,
    )
    
    pulse_x = td.GaussianPulse(phase=0, freq0=freq0, fwidth=freq0 / 5)
    source_x = td.PlaneWave(
        center=(0, 0, posi_source),
        size=(td.inf, td.inf, 0),
        source_time=pulse_x,
        direction='+',
        pol_angle= 0,#np.pi/2,
        angle_theta=np.deg2rad(inc_angle),
        angle_phi=0,
        name="source_x"
    )

    # create monitor
    monitor_xz = td.FieldMonitor(
        center=(0, 0, 0),  # center of the monitor
        size=(td.inf, 0, td.inf),  # size of the monitor
        freqs=freqs,  # frequency points to record the fields at
        name="fields_xz",
    )
    monitor_d = td.FieldMonitor(
        center=(0, 0, Metal_t+diel_t/2),  # center of the monitor
        size=(td.inf, 0, diel_t),  # size of the monitor
        freqs=freqs,  # frequency points to record the fields at
        name="fields_d",
    )

    flux_t = td.FluxMonitor(
        center=(0, 0, posi_monitor),
        size=(td.inf, td.inf, 0),
        freqs=freqs,
        name="flux_trans"
    )
    flux_r = td.FluxMonitor(
        center=(0, 0, posi_source-0.05), 
        size=(td.inf, td.inf, 0),
        freqs=freqs,
        name="flux_refl"
    )

    sim_size = (FDTD_x, FDTD_y, FDTD_z)
    bloch_x = td.Boundary.bloch_from_source(source=source_x, domain_size=sim_size[0], axis=0, medium=subst)
    bloch_y = td.Boundary.bloch_from_source(source=source_x, domain_size=sim_size[1], axis=1, medium=subst)

    sim = td.Simulation(
        size=sim_size,  # simulation domain size
        grid_spec=td.GridSpec.auto(min_steps_per_wvl=25),  # automatic nonuniform FDTD grid with 25 grids per wavelength in the material
        boundary_spec=td.BoundarySpec(
            x=bloch_x,
            y=td.Boundary.periodic(),
            z=td.Boundary.pml()
        ),
        structures=[substrate,film_m,film_d],
        sources=[source_x],
        monitors=[monitor_xz,monitor_d,flux_t,flux_r],
        run_time=t_run,  # physical simulation time in second
    )
    return sim

周期構造で斜め入射を用いるときは通常のPeriodic構造を用いることができないためBloch条件を用いることに注意が必要。

構造およびmeshの確認

2D plotを用いて想定した構造や適切なmeshが作成できているかを確認する。
前回は3D plotで三次元モデルを呼び出してしまっていたがはるかに簡便かつ視認性の良い2D plotが可能なことが分かったため今回からこちらを使用する。
左側のoverviewで全体の構造や入射方向や電場方向を確認し、右に表示した拡大図でmeshが適切に切られているかを確認する。

sim = mkobjct(inc_angle=60)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4), tight_layout=True)

sim.plot(y=0, freq=freq0, ax=ax1)
ax1.set_title("Structure overview @ y=0")

sim.plot(y=0, freq=freq0, ax=ax2)
sim.plot_grid(y=0, freq=freq0, ax=ax2)
ax2.set_xlim(-0.2, 0.2)
ax2.set_ylim(-0.1, 0.3)  # z軸：-100 nm ~ +300 nm
ax2.set_title("Zoomed structure with grid")

plt.show()

実行・値の格納

forloopで引数を変数として複数のsimulationをまわす。
50~70°までを5°間隔で確認するため、theta_values = np.arange(50, 71, 5)とした。
全てのデータを保存すると構造次第では天文学的な容量を保存してしまうため、普段は評価に必要な値のみを格納しておく手法をとっている。
今回は一般にSPPが励起されていることの確認に用いられるreflectance, 励起されたSPPの分布を確認するため、指定した波長におけるEz成分の空間分布を保存する。

theta_values = np.arange(50, 71, 5)  # 入射角度（単位: degree）
wl_save = 630  # 保存したい波長 [nm]
wl_cell = int((wl_save - min_wl) / step_wl)

# 保存用リスト
Ezfield_list = []
flux_refl_list = []

for theta in theta_values:
    print(f"Running simulation for angle = {theta} degree")
    sim = mkobjct(inc_angle=theta) 

    task_name = f"angle_{theta}"
    filepath = f"data/{task_name}.hdf5"

    data = web.run(sim, task_name=task_name, path=filepath, verbose=False)

    # Reflectance 保存
    flux_r = data["flux_refl"].flux.data
    flux_refl_list.append(flux_r)

    # Ez フィールド保存
    fd = data["fields_xz"]
    Ez_real = fd.Ez.real.isel(f=wl_cell)
    Ezfield_list.append(Ez_real)

print("Finished all simulations.")

# 結果を numpy 配列に変換
flux_refl_array = np.array(flux_refl_list)    # shape: (n_theta, n_freqs)

グラフの表示

• 入射角度ごとのreflectance
• y=0 monitorにおけるEz成分(λ = 630 nm)
  を表示する。

flux_trans = data["flux_trans"].flux
freqs = flux_trans.f.values
wavelengths_nm = 1e9 * td.C_0 / freqs
view_min = 550  # nm
view_max = 750  # nm
plot_wid = lmbda0

wavelengths_nm = 1e3*td.C_0 / freqs

#######################################################################################
fig, ax = plt.subplots(figsize=(6, 4))

# 各角度に対して反射スペクトルをプロット
for refl, theta in zip(flux_refl_array, theta_values):
    ax.plot(wavelengths_nm, -refl, label=f"{theta}°")

ax.axvline(lmbda0 * 1e3, ls="--", color="k", lw=1, label="λ₀")

# 軸と凡例
ax.set_xlim(view_min, view_max)
ax.set_ylim(0, 1)
ax.set_xlabel("Wavelength (nm)")
ax.set_ylabel("Reflectance")
ax.set_title("Reflectance spectra vs incident angle")
ax.legend(title="Incident angle")
ax.grid(True)
plt.tight_layout()
plt.show()
######################################################################################
fig, axes = plt.subplots(1, len(theta_values), figsize=(3.2 * len(theta_values), 2.8))
vmax = np.max(np.abs(Ezfield_list))  #スケーリング

y_index = Ezfield_list[0].sizes["y"] // 2

for i, theta in enumerate(theta_values):
    ax = axes[i]
    Ez = Ezfield_list[i] 
    img = Ez.isel(y=y_index).transpose("z", "x").plot(
        ax=ax,
        cmap="RdBu_r",
        vmin=-vmax,
        vmax=+vmax,
        add_colorbar=False,
    )
    ax.set_xlim(-lmbda0, lmbda0)
    ax.set_ylim(-0.5, 1)
    ax.set_title(f"{theta}°", fontsize=10)
    ax.set_xlabel("x [μm]")
    ax.set_ylabel("z [μm]")

    if i != 0:
        ax.set_ylabel('')
        ax.set_yticklabels([])

cbar = fig.colorbar(img, ax=axes.ravel().tolist(), shrink=0.9)
cbar.set_label("Re[Ez] [a.u.]")
plt.show()

上記のcodeの実行で上の様なグラフが表示できる。
上のグラフは各入射角における反射率を示している。SPPが励起された波長では反射率が低くなるため、各グラフにお蹴るdip位置がSPPが励起されている角度に対応する。
今回対象とした630 nmは、65°入射において強くSPPが励起されていることが示唆される。
下に示すEz成分の空間分布をみると、確かに65°の場合に最も強くSPPが励起されている様子が確認できる。

再計算

実際に励起されたSPPの電磁場分布を用いた計算が行いたいため、sweepで得られた角度(65°)を選択して個別にもう一度simulationを実行する。

# run simulation
theta = 65
print(f"Running simulation for angle = {theta} degree")
sim = mkobjct(inc_angle=theta)
task_name = f"angle_{theta}"
filepath = f"data/{task_name}.hdf5"

data = web.run(sim, task_name=task_name, path=filepath, verbose=True)
print(data.log)

演算・プロット

SPPが励起されている場合に確認できる

• 界面垂直方向への電場成分
• 伝搬軸と直行する方向へのスピン角運動量
  が確かに波長630 nmにおいて極値をとっているかを確認する
  あらかじめ絶縁体層部分のみに対して設定しておいたmonitor(fields_d)で取得した値を対象とし、

transmission = data["flux_trans"].flux
reflection = -data["flux_refl"].flux
fd = data["fields_xz"]
fd_d = data["fields_d"]
Ez_real = fd.Ez.real
freqs = transmission.f.values
wavelengths_nm = 1e9 * td.C_0 / freqs
view_min = 580  # nm
view_max = 680  # nm
n_plot = 5

epsilon_0 = 8.854187817e-12  # [F/m]
mu_0 = 4 * np.pi * 1e-7      # [H/m]
Ex = fd_d.Ex.isel(y=y_index)
Ez = fd_d.Ez.isel(y=y_index)
Hx = fd_d.Hx.isel(y=y_index)
Hz = fd_d.Hz.isel(y=y_index)
Ez_abs = (np.abs(Ez) ** 2).mean(dim=["x", "z"]).values
Ey_comp = np.imag(epsilon * (np.conj(Ez) * Ex - np.conj(Ex) * Ez))
Hy_comp = np.imag(mu * (np.conj(Hz) * Hx - np.conj(Hx) * Hz))
Sy_mean = (Ey_comp + Hy_comp).mean(dim=["x", "z"]).values
Sy_norm = Sy_mean / np.max(np.abs(Sy_mean))
Ez_norm = Ez_abs / np.max(np.abs(Ez_abs))
####################################################################################
# 横並びプロット
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4), tight_layout=True)

# 反射・透過スペクトル
ax1.plot(td.C_0 / freqs * 1e3, reflection, label="Reflection")
ax1.plot(td.C_0 / freqs * 1e3, transmission, label="Transmission")
ax1.axvline(lmbda0 * 1e3, ls="--", color="k", lw=1, label="λ₀")
ax1.set_xlim(view_min, view_max)
ax1.set_ylim(0, 1)
ax1.set_xlabel("Wavelength (nm)")
ax1.set_ylabel("Reflectance / Transmittance")
ax1.set_title("Spectra")
ax1.legend()
ax1.grid(True)

# Syの平均強度
ax2_right = ax2.twinx()
ax2.plot(td.C_0 / freqs * 1e3, Ez_norm, label="|Ez|²", color="tab:red")
ax2_right.plot(td.C_0 / freqs * 1e3, Sy_norm, label="⟨S_y⟩", color="tab:blue")
ax2.set_xlim(view_min, view_max)
ax2.set_xlabel("Wavelength (nm)")
ax2_right.set_ylabel("Spin AM density [a.u.]", color="tab:blue")
ax2.set_ylabel("Intensity (|Ez|²)", color="tab:red")
ax2.set_ylim(-1.2, 1.2)
ax2_right.set_ylim(-1.2, 1.2)
ax2.set_title("Average Ez and Sy")
#ax2.legend()
#ax2_right.legend()
ax2.grid(True)

plt.show()

#################################################################################################
# 表示する波長のインデックスを指定
number_min = int((view_min - min_wl) / step_wl)
number_max = int((view_max - min_wl) / step_wl)
indices = np.linspace(number_min, number_max, n_plot, dtype=int)

fig, axes = plt.subplots(1, len(indices), figsize=(2 * len(indices), 2.5), constrained_layout=True)

y_index = Ez_real.sizes["y"] // 2 
vmax = np.max([abs(Ez_real.isel(f=i, y=y_index)).max().item() for i in indices])

for plot_i, i in enumerate(indices):
    wl_nm = wavelengths[i] * 1e3
    ax = axes[plot_i]
    img = Ez_real.isel(f=i, y=y_index).transpose("z", "x").plot(
        ax=ax,
        cmap="RdBu_r",
        vmin=-vmax,
        vmax=+vmax,
        add_colorbar=False
    )
    ax.set_title(f"Re[Ez] (λ = {wl_nm:.0f} nm)", fontsize=9)
    ax.set_xlim(-lmbda0, lmbda0)
    ax.set_ylim(-0.5, 1)
    ax.set_xlabel("x [μm]")
    ax.set_ylabel("z [μm]")

    if plot_i != 0:
        ax.set_ylabel('')
        ax.set_yticklabels([])

cbar = fig.colorbar(img, ax=axes.ravel().tolist(), shrink=0.9)
cbar.set_label("Field [a.u.]")
plt.show()

波長630 nmにおいてSPPが励起され反射率が最低になっていること、絶縁体層ないにおける電場強度・Poynting vector, スピン角運動量が極値をとることが確認できる。