[nannos]オープンソフトで試すRigorous Coupled Wave Analysis(RCWA) ～Blazed Gratingの光シミュレーション～

はじめに

前回の記事ではRCWA(Fourier modal methodとも呼ばれる)のオープンソースのシミュレータnannosをご紹介しました。

nannosのサイトはこちら

今回はnannosを用いてブレーズド回折格子の透過・反射スペクトルをシミュレーションしてみたので、ご紹介します。

今回の計算対象

こちらが今回計算する構造です。
基板上にとげとげのグレーティング（ブレーズド回折格子）が形成された構造の透過率、反射率(及び各回折次数のスペクトル)のシミュレーションを実施します。

以下が1周期分を切り出した誘電率分布です。

RCWAでは斜目の傾斜構造をそのまま取り込んで計算することはできず、上の図のように階段近似して計算することになります。
今回は1次元周期のブレーズド回折格子で、
- 周期サイズは1 um
- 三角構造の高さ0.8 um, 底辺の長さ1 umとしました。
- 基板の屈折率は1. 5
- 三角構造の屈折率は2.0
- それ以外の媒質は空気に設定

計算するにあたり、
- グレーティングの奥行き方向に振動する電場で入射する場合をS偏光
- 周期方向に振動する電場で入社する場合をp偏光
としました。

シミュレーション結果

今回は計算に含めるハーモニクスの数(構造内のモード解析で含める逆格子ベクトルの数)を51として計算してみました。
こちらが回折光のスペクトルになります。

続いてこちらが透過率・反射率のシミュレーション結果です。

計算に用いたスクリプト

import numpy as np
import matplotlib.pyplot as plt
import nannos as nn

# ============================================================
# 1. 基本パラメータ
# ============================================================
period = 1.0
height = 0.8
nslices = 40

n_sup = 1.0
n_sub = 1.45
n_gr  = 2.0

eps_sup = n_sup**2
eps_sub = n_sub**2
eps_gr  = n_gr**2

# ---- ブレーズ設定
blaze_mode = "both"     # "right", "left", "both"
peak_pos  = 0.5
etched_groove = False

# 入射条件（上から）
theta_deg = 0.0
phi_deg   = 0.0

polarizations = {
    "P": 0.0,
    "S": 90.0,
}

wavelengths = np.linspace(1.3, 1.7, 61)

nh = 51
orders_to_plot = [-2, -1, 0, 1, 2]

Nx, Ny = 1024, 16
Ly = period

# ============================================================
# 2. Lattice
# ============================================================
lattice = nn.Lattice(
    [[period, 0.0], [0.0, Ly]],
    discretization=(Nx, Ny)
)
X, Y = lattice.grid
frac = (X / period) % 1.0

# ============================================================
# 3. ブレーズ構造生成（★ z方向反転 ★）
# ============================================================
def make_blazed_layers():
    layers = [lattice.Layer("Superstrate", epsilon=eps_sup)]
    dz = height / nslices

    for k in range(nslices):
        zc = (k + 0.5) * dz
        zh = zc / height
        zh_inv = 1.0 - zh   # ★ z方向反転

        if blaze_mode == "right":
            in_ramp = frac < zh_inv

        elif blaze_mode == "left":
            in_ramp = frac > (1.0 - zh_inv)

        elif blaze_mode == "both":
            dist = np.abs(frac - peak_pos)
            max_dist = max(peak_pos, 1.0 - peak_pos)
            dist_norm = dist / max_dist
            in_ramp = dist_norm < zh_inv

        else:
            raise ValueError("Invalid blaze_mode")

        if not etched_groove:
            eps = lattice.ones() * eps_sup
            eps[in_ramp] = eps_gr
        else:
            eps = lattice.ones() * eps_gr
            eps[in_ramp] = eps_sup

        layers.append(
            lattice.Layer(
                f"Blaze_{k:02d}",
                thickness=dz,
                epsilon=eps
            )
        )

    layers.append(lattice.Layer("Substrate", epsilon=eps_sub))
    return layers

stack = make_blazed_layers()

# ============================================================
# 4. 誘電率分布表示
# ============================================================

# ---- ε(x,y)
pick = [1, 1 + nslices // 2, nslices]
fig, axs = plt.subplots(1, 3, figsize=(10, 3), constrained_layout=True)

for ax, idx in zip(axs, pick):
    plt.sca(ax)
    im = stack[idx].plot(cmap="viridis")[0]
    ax.set_title(stack[idx].name)
    plt.colorbar(im, ax=ax, orientation="horizontal")

plt.suptitle("Permittivity ε(x,y)")
plt.show()

# ---- ε(x,z)
eps_xz = np.zeros((nslices, Nx))
for k in range(nslices):
    eps_xz[k, :] = np.real(stack[1 + k].epsilon[:, 0])

plt.figure(figsize=(8, 4))
plt.imshow(
    eps_xz,
    origin="lower",
    aspect="auto",
    extent=[0, period, 0, height],
    cmap="viridis"
)
plt.xlabel("x [um]")
plt.ylabel("z [um]")
plt.title("Blazed grating cross section (high-index z-inverted)")
plt.colorbar(label="epsilon")
plt.tight_layout()
plt.show()

# ============================================================
# 5. スペクトル計算
# ============================================================
def run_spectrum(psi_deg):
    Rtot, Ttot = [], []
    Tm = {m: [] for m in orders_to_plot}
    Rm = {m: [] for m in orders_to_plot}

    for lam in wavelengths:
        pw = nn.PlaneWave(lam, angles=(theta_deg, phi_deg, psi_deg))
        sim = nn.Simulation(stack, pw, nh=nh)

        R, T = sim.diffraction_efficiencies()
        Rtot.append(float(R))
        Ttot.append(float(T))

        Ri, Ti = sim.diffraction_efficiencies(orders=True)
        for m in orders_to_plot:
            Tm[m].append(float(sim.get_order(Ti, (m, 0))))
            Rm[m].append(float(sim.get_order(Ri, (m, 0))))

    return np.array(Rtot), np.array(Ttot), Tm, Rm

results = {pol: run_spectrum(psi) for pol, psi in polarizations.items()}

# ============================================================
# 6. スペクトル表示
# ============================================================
for pol, (Rtot, Ttot, Tm, Rm) in results.items():
    plt.figure(figsize=(8, 5))
    for m in orders_to_plot:
        plt.plot(wavelengths, Tm[m], label=f"T(m={m})")
    plt.title(f"{pol}-pol Transmitted orders")
    plt.xlabel("Wavelength [um]")
    plt.ylabel("Efficiency")
    plt.legend()
    plt.tight_layout()
    plt.show()

# ---- 全R/T
plt.figure(figsize=(7, 5))
for pol, (Rtot, Ttot, *_rest) in results.items():
    plt.plot(wavelengths, Ttot, label=f"{pol}-pol T")
    plt.plot(wavelengths, Rtot, "--", label=f"{pol}-pol R")
plt.xlabel("Wavelength [um]")
plt.ylabel("Efficiency")
plt.title("Total Transmission / Reflection")
plt.legend()
plt.tight_layout()
plt.show()

まとめ

今回オープンソースのRCWAシミュレータであるnannosを使ってブレーズド回折格子の透過率・反射率を求めてみました。今回は1次元周期の構造でしたが、2次元構造も対応しているので次回はより複雑な構造のシミュレーションを試していきます。

※本記事は筆者個人の見解であり、所属組織の公式見解を示すものではありません。

問い合わせフォームのご連絡

問い合わせ

光学シミュレーションソフトの導入や技術相談、
設計解析委託をお考えの方はサイバネットシステムにお問合せください。

光学ソリューションサイトについては以下の公式サイトを参照：
👉 [光学ソリューションサイト(サイバネット)]

光学分野のエンジニアリングサービスについては以下の公式サイトを参照：
👉 [光学エンジニアリングサービス(サイバネット)]