Python Scriptで実行するtidy3Dシミュレーションと電場分布解析の一例

本記事では、FDTD法に基づく電磁場シミュレーションツール「Tidy3D (https://www.flexcompute.com/tidy3d/solver/ )」を用いて、簡単な構造の光学応答を再現・解析した一連の手順を紹介する。
Tidy3Dはブラウザベースで動作するFDTDソルバーであり、Pythonによる構造定義・実行・演算・可視化をすべて同一環境内で完結できる点が特徴である。一定の計算量までは無料で使用でき、ローカル環境へのインストールも不要という利点がある。
今回は、実際にTidy3DのNotebook上で
・構造とシミュレーション条件の定義
・シミュレーションの実行
・得られた電場データの解析と可視化
を行った一例を紹介する。
なお、対象とした物理現象は、すでに別のFDTDソフトを用いて確認されている、金属ナノ構造における円偏光励起下での電場モードの空間的なねじれを対象とした。
https://journals.aps.org/prb/abstract/10.1103/PhysRevB.109.035428

Tidy3Dの特徴には、実行したデータの解析からplotまでを環境内で一貫して行えることが挙げられる。このため、従来であればシミュレーションデータを一度エクスポートして他の解析ツールに渡していた処理（例えば電場分布の左右非対称性の評価）も、Tidy3D環境内だけで完結させることが可能である。今回はその実用性を確認する目的も兼ねている。
基本的には各コードブロックは、順番通りにNotebookへ入力・実行していくことで同じシミュレーション結果が再現されるはずである。

• 実行環境の起動
  まずTidy3Dにログインし、画面左側のメニューから「Python Notebook」を選択する。
  
  

「Blank File」を選ぶと以下のようなNotebook環境が起動する。

このNotebook上で、Pythonコードを順に入力していくことでシミュレーションを実行できる。以下にその各ブロックを簡単に説明しながら紹介する。

以下、blockごとに簡単な説明を加えた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

• 変数の決定
  使用する変数などを準備。波長はデフォルトではμmスケール。monitorなどで波長400 nm ~ 1.4 μmを50 nm間隔で設定したい場合などにはあらかじめここでnp.arange(0.4, 1.41, 0.05)などを使用しておくと便利。
  今回は以下の値を設定。
  ・FDTD simulation全域の大きさ
  ・plane waveの位置
  ・計算時間
  ・構造の大きさ
  ・使用する金属の誘電関数
  今回は非常に単純な構造を対象とするためここでまとめてしまったが構造が複雑な場合はobjectの作成箇所で定義したほうが良いかもしれない。
  使用可能な誘電関数は[https://docs.flexcompute.com/projects/tidy3d/en/v2.0.2/material_library.html]を参照。

lmbda0 = 0.633
freq0 = td.C_0 / lmbda0
wavelengths = np.arange(0.4, 1.41, 0.05)
freqs = [td.C_0 / (wl) for wl in wavelengths]

#geometry for simulation
FDTD_x = 1
FDTD_y = 1
FDTD_z = 4
posi_source = -1.5
t_run = 3e-13

#geometry for object
#rectangle
rect_x = 0.3
rect_y = rect_x
rect_z = 0.1

gold = material_library['Au']['RakicLorentzDrude1998']

• 構造を作製
  今回は正方形構造を対象とするためgeometryはtd.Boxを使用。球や円柱を使用する場合はsphereやcylinderなどを使用する[http://docs.flexcompute.com/projects/tidy3d/en/latest/api/geometry.html]。

#Make object
rect = td.Structure(
    geometry = td.Box(center=(0, 0, 0), size=(rect_x, rect_y, rect_z)),  # μm 単位
    medium   = gold
)

• 励起波源の定義
  今回は円偏光を定義したかったため位相をずらした直線偏光を二つ導入する。円偏光を直接定義する方法がある可能性もあるがまだ見つけていない。中心波長や半値などで必要なパルスを定義[https://docs.flexcompute.com/projects/tidy3d/en/latest/api/_autosummary/tidy3d.GaussianPulse.html]。

phase_shift = np.pi/2
pulse_y = td.GaussianPulse(freq0=freq0, fwidth=freq0 / 5.0)
pulse_x = td.GaussianPulse(phase=phase_shift, freq0=freq0, fwidth=freq0 / 5.0)

source_x = td.PlaneWave(
    center=(0, 0, posi_source),
    size=(td.inf, td.inf, 0),
    source_time=pulse_x,
    direction='+',
    pol_angle=np.pi/2,  # x偏光
    name="source_x"
    )
source_y = td.PlaneWave(
    center=(0, 0, posi_source),
    size=(td.inf, td.inf, 0),
    source_time=pulse_y,
    direction='+',
    pol_angle=0,  # y偏光
    name="source_y"
)

• 電場分布や透過/反射成分を観測するためのmonitorの設置
  td.FieldMonitorで電磁場の空間分布を取得、flux monitorでエネルギーの流れを取得する。

# create monitor
monitor_t = 1.4

monitor_xy = td.FieldMonitor(
    center=(0, 0, 0),  # center of the monitor
    size=(td.inf, td.inf, 0),  # size of the monitor
    freqs=freqs,  # frequency points to record the fields at
    name="fields_xy",
)
monitor_yz = td.FieldMonitor(
    center=(0, 0, 0),  # center of the monitor
    size=(0,td.inf, td.inf),  # size of the monitor
    freqs=freqs,  # frequency points to record the fields at
    name="fields_yz",
)

flux_t = td.FluxMonitor(
    center=(0, 0, monitor_t),         
    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"
)

• simulationの作成
  structures=[rect],
  sources=[source_x,source_y],
  monitors=[monitor_xy,monitor_yz,flux_t,flux_r],
  の箇所で設定しておいたobjectやsource, monitorを読み込む。ここに記入しなければ上で定義してもsimulationには反映されないため注意。
  境界条件やmeshサイズなどもここで指定する。必要に応じて変数にして一番上で定義してしまってもよい。

sim = td.Simulation(
    size=(FDTD_x, FDTD_y, FDTD_z),  # 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=td.Boundary.periodic(),
    #y=td.Boundary.periodic(),
    x=td.Boundary.pml(),
    y=td.Boundary.pml(),
    z=td.Boundary.pml(),
    ),
    structures=[rect],
    sources=[source_x,source_y],
    monitors=[monitor_xy,monitor_yz,flux_t,flux_r],
    run_time=t_run,  # physical simulation time in second
) 

• 3D plotで指定した構造が作成できているかやmeshの数を確認する。

sim.plot_3d()
print(f"simulation grid is shaped {sim.grid.num_cells} for {int(np.prod(sim.grid.num_cells) / 1e6)} million cells.")

[sim.plot_3d()]で上に示すような画像が表示される。
今回の構造では
・黄色: monitor
・緑色: plane wave source(緑矢印: 伝搬方向, 桃色: 電場の方向)
・灰色: object
となっている。
また、二行目に対応してMeshサイズが表示されるため、必要に応じてsimulation作成部分におけるgrid_spec=td.GridSpec.auto(min_steps_per_wvl=25)の値を修正する。

• simulationの実行、logの表示

# run simulation
data = td.web.run(sim, task_name="square_Au", path="data/data.hdf5", verbose=True)
# see the log
print(data.log)

• 得られたデータを用いた演算、グラフの表示
  今回は
  ・透過成分のスペクトル
  ・z=0 monitorにおける電場の平均値
  ・z=0 monitorにおける電場分布の左右非対称性
  を波長の変数としてプロットした。
  flux monitorのスペクトルは厳密には透過スぺクトルとは呼んではいけないと思うがおおむね対応すると考えられるため今回はそのまま用いる。

flux_trans = data["flux_trans"].flux
freqs = flux_trans.f.values
wavelengths_nm = 1e9 * td.C_0 / freqs

fd = data["fields_xy"]  
E2 = fd.Ex.abs**2 + fd.Ey.abs**2 + fd.Ez.abs**2  # |E|² を計算
E2_np = E2.values.copy()
E2_flipped_np = E2_np[::-1, :, :, :]
E2_flipped = xr.DataArray(
    E2_flipped_np,
    dims=E2.dims,
    coords=E2.coords,
    attrs=E2.attrs
)

E2_sum = E2 + E2_flipped
E2_diff = abs(E2 - E2_flipped)
eps = 1e-20
E2_sum_safe = E2_sum + eps
E2_asym = E2_diff / E2_sum_safe

E2_avg_vs_wl = [E2.isel(f=i).mean().item() for i in range(E2.sizes["f"])]
E2_asym_vs_wl = [E2_asym.isel(f=i).mean().item() for i in range(E2_asym.sizes["f"])]

fig2, axes2 = plt.subplots(1, 3, figsize=(10, 3))

# 透過スペクトル
axes2[0].plot(wavelengths_nm, flux_trans.values, label="Transmittance", color="C0")
axes2[0].set_xlabel("Wavelength (nm)")
axes2[0].set_ylabel("Flux (linear)")
axes2[0].set_title("Transmission Spectrum")
axes2[0].legend()
axes2[0].grid(True)

# 平均E2
axes2[1].plot(wavelengths_nm, E2_avg_vs_wl, label="⟨|E|²⟩", color="C1")
axes2[1].set_xlabel("Wavelength (nm)")
axes2[1].set_ylabel("Average |E|²")
axes2[1].set_title("Avg. Field Intensity vs Wavelength")
axes2[1].legend()
axes2[1].grid(True)

# 非対称度
axes2[2].plot(wavelengths_nm, E2_asym_vs_wl, label="Asymmetry Factor", color="C2")
axes2[2].set_xlabel("Wavelength (nm)")
axes2[2].set_ylabel("Normalized Difference")
axes2[2].set_title("Field Asymmetry vs Wavelength")
axes2[2].legend()
axes2[2].grid(True)

plt.tight_layout()
plt.show()

上記のcodeの実行で以下の様なグラフが表示できる。
透過波と電場の増強度のグラフより、波長1.4 μm付近に低次mode、500 nm付近に高次のmodeが立っていることが推定できる。
非対称性グラフからは低次-高次モードの間である600 nm付近において電場分布の非対称性がピークを持つことが確認でき、simulation結果のexportなどを行うことなく非対称性が大きくなる波長の探索が可能であることがわかる。

• ピーク波長における電場強度などの空間分布表示
  i = 5で配列中の何番目の周波数に該当するデータを表示するか選択している。

i = 5 # λを選択
wl = wavelengths[i]*1e3
Hz = fd.Hz.real 

vmin = E2.min().item()
vmax = E2.max().item()

# 描画
fig, axes = plt.subplots(1, 3, figsize=(9, 2.5))

E2.isel(f=i).plot(ax=axes[0], cmap="inferno", vmin=0, vmax=vmax/5)
axes[0].set_title(f"Original |E|² at λ = {wl:.0f} nm")
axes[0].set_xlim(-0.3, 0.3)
axes[0].set_ylim(-0.3, 0.3)
E2_asym.isel(f=i).plot(ax=axes[1], cmap="inferno", vmin=0, vmax=1)
axes[1].set_title("difference")
axes[1].set_xlim(-0.3, 0.3)
axes[1].set_ylim(-0.3, 0.3)
# Magnetic field
Hz.isel(f=i).plot(ax=axes[2], cmap="RdBu_r")
axes[2].set_title("Hz (real)")
axes[2].set_xlim(-0.3, 0.3)
axes[2].set_ylim(-0.3, 0.3)

plt.tight_layout()
plt.show()

左に示す電場の平均強度より、波長650 nmにおいて大きな非対称性が生じていることが確認できる(中央の画像は非対称性の計算に用いたもとの電場の左右を反転させてから元の画像との差分をとって絶対値をとったもの)。磁場分布なども同時に表示することで生じているモードの確認が行える。

• 複数波長の電場強度分布の表示

fig, axes = plt.subplots(1, 5, figsize=(20, 3))  # サイズおよび個数の指定

step =5  # 表示間隔
indices = list(range(0, len(wavelengths), step))
n_plot = len(indices)

for plot_i, data_i in enumerate(indices):
    wl_nm = wavelengths[data_i] * 1e3
    ax = axes[plot_i]

    img = E2.isel(f=data_i, z=0).plot(
        ax=ax,
        cmap="inferno",
        # vmin=vmin,
        # vmax=vmax/5,
        add_colorbar=False
    )
    ax.set_xlim(-0.3, 0.3)
    ax.set_ylim(-0.3, 0.3)
    ax.set_title(f"λ = {wl_nm:.0f} nm")

    if plot_i != 0:
        ax.set_ylabel('')
        ax.set_yticklabels([])
        
cbar = fig.colorbar(img, ax=axes.ravel().tolist(), shrink=0.9)
cbar.set_label("|E|² [a.u.]")
plt.show()

波長650 nmにおいて非対称分布が最大化している様子が確認できる。