Semi-vectorのBeam Propagation MethodをPythonで自作して光導波路の伝搬シミュレーションを行う①~TEモード編

1.はじめに

光導波路の伝搬解析では、Beam Propagation Method（BPM）がよく用いられます。BPMは導波路内の光の進行を数値的に計算する手法として広く知られており、フォトニックデバイス設計や光集積回路の解析などで頻繁に利用されています。実際の設計現場では、有償シミュレーションソフトウェアに実装されている BPM を利用するケースが多く、高機能で安定した解析環境を提供してくれます。

一方で、BPMのアルゴリズムそのものはそれほど複雑ではなく、基本的な数値解析の知識があれば比較的シンプルに実装することができます。特に近年は生成AIを活用することで、数値計算コードの雛形を短時間で作成できるようになり、自作シミュレーターのハードルは以前よりもかなり下がってきています。アルゴリズムの理解を深めたり、特定の用途に合わせて自由にカスタマイズしたりする目的で、自分でBPMコードを実装してみるのも面白い試みだと思います。

本記事では、Pythonを用いて Semi-vector Beam Propagation Method（semi-vector BPM） を実装し、光導波路中の伝搬シミュレーションを行う方法を紹介します。BPMにはスカラー近似（Scalar BPM）やフルベクトル法（Full-vector BPM）などいくつかの種類がありますが、今回は計算コストと精度のバランスが良い Semi-vector BPM を対象とします。これは主偏波成分のみを解きつつ、屈折率境界での電磁場の振る舞いをある程度正しく扱える手法であり、導波路解析では実用的な近似として広く使われています。

記事ではまず、semi-vector BPM の基本式を簡単に整理し、その後 Python による実装方法を紹介します。さらに、自作したコードを用いて簡単な光導波路の伝搬シミュレーションを行い、モードの励振と伝搬の様子を可視化してみます。

なお、本記事で紹介するコードは CPU 上で動作するシンプルな実装になっていますが、BPMの計算は並列化と非常に相性が良いため、GPUを用いた高速化も比較的容易に行うことができます。次回の記事では、今回作成したコードをベースに GPU（CuPy など）を利用した高速化 にも挑戦してみたいと思います。

2. semi-vector BPM

フォトニック導波路の解析では Beam Propagation Method (BPM) がよく使われます。
BPMにはいくつかの近似レベルがあります。

• Scalar BPM
• Semi-vector BPM
• Full-vector BPM

この記事では semi-vector BPM を式から導出しながら説明します。

対象は主に

• リブ導波路
• チャネル導波路
• 光集積回路

などです。

---

BPMについての岡本先生の書籍にしっかり書かれているので、説明は省くとして、

semi-vector BPMでは

主偏波のみを解きます。

例えば TE-like モードでは

E_x

を主成分とします。

このとき横方向の演算子は

\frac{\partial}{\partial x}
\left(
\frac{1}{\varepsilon}
\frac{\partial E_x}{\partial x}
\right)
+
\frac{\partial^2 E_x}{\partial y^2}

になります。

つまり

\partial_x
\left(
\frac{1}{\varepsilon}
\partial_x
\right)

という weighted operator が現れます。

---

3. semi-vector BPM の最終式

BPM形式にすると

\frac{\partial E_x}{\partial z}
=
\frac{i}{2\beta_{ref}}
\left[
\partial_x
\left(
\frac{1}{\varepsilon}
\partial_x E_x
\right)
+
\partial_y^2 E_x
+
k_0^2(n^2-n_{ref}^2)E_x
\right]

になります。

ここで

\varepsilon=n^2

です。

---

4. 数値計算（ADI法）

実際のBPMでは

\partial_x
\left(
\frac{1}{\varepsilon}
\partial_x
\right)
+
\partial_y^2

をそのまま解くと計算が重くなります。

そのため ADI法 (Alternating Direction Implicit) を使います。

1ステップ

z \rightarrow z+\Delta z

は

(I-\alpha L_x)(I-\alpha L_y)E^{n+1}
=
(I+\alpha L_x)(I+\alpha L_y)E^n

のように分解します。

ここで

\alpha = \frac{\Delta z}{4i\beta_{ref}}

です。

この方法を使うと

• 計算量が O(N)
• 数値安定性が高い

というメリットがあります。
これについては以下の論文を参考にしました。
"BPM-Matlab: an open-source optical propagation simulation tool in MATLAB"
M. Veettikazhy et al. Optics Express Vol. 29, Issue 8, pp. 11819-11832 (2021)

---

5. BPMの境界条件

BPMでは

z方向は

E(x,y,0)

を与える 初期値問題 になります。

つまり

• 入力面 → 入力モード
• 出力面 → 計算結果

です。

一方横(導波路の断面)方向では

E=0

などの境界条件を置きます。

実際には

• absorbing mask
• complex index

を使って反射を抑えることが多いです。

---

6. semi-vector BPM の特徴

メリット

• scalar BPMより高精度
• full-vectorより高速
• シリコン導波路でも使える

デメリット

• 偏波結合・変換は扱えない
• 強い反射は扱えない
• backward wave は存在しない

---

7. BPMのまとめ

semi-vector BPMは

• scalar BPMの高速性
• vector solverの精度

の中間に位置する手法です。

特に

• rib waveguide
• silicon photonics
• integrated optics
  の解析では非常に有効です。

---

8. 今回計算する光導波路

今回はsemi vector BPMの光導波路の計算コードを作成しましたので、ご紹介します。
semi-vector-bpm/
├── run_demo_semivector_ex_r2.py 実行用のファイル
├── mode_solver.py　導波路断面の屈折率分布と計算グリッドを作るためのファイル
├── semivector_ex_mode_solver.py mode_solver.pyで指定した断面の固有モード計算用
└── semi_vector_bpm_ex.py BPMの伝搬シミュレーション用

上のファイルツリーでのexは入射モードとしてTEモードを設定した場合、eyはTMモード入射時の伝搬シミュレーションに対応しています。
各偏波モードそれぞれにスクリプトで分けています。
TMモード計算用のコードは次回に紹介します。
ちなみに
導波路断面方向の境界条件はディリクレ境界条件(E=0)としました。

9.TEモードの伝搬シミュレーション用のBPMスクリプト

こちらが今回自作したスクリプトになります。
run_demo_semivector_ex_r2.py

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

from mode_solver import build_rib_waveguide
from semi_vector_bpm_ex import semivector_ex_step
from semivector_ex_mode_solver import solve_ex_semivector_modes


def make_absorbing_mask(Nx, Ny, width=18, strength=3.5):
    mask = np.ones((Nx, Ny), dtype=np.float64)

    for i in range(Nx):
        if i < width:
            s = (width - i) / width
            mask[i, :] *= np.exp(-strength * s * s)
        elif i >= Nx - width:
            s = (i - (Nx - width - 1)) / width
            mask[i, :] *= np.exp(-strength * s * s)

    for j in range(Ny):
        if j < width:
            s = (width - j) / width
            mask[:, j] *= np.exp(-strength * s * s)
        elif j >= Ny - width:
            s = (j - (Ny - width - 1)) / width
            mask[:, j] *= np.exp(-strength * s * s)

    return mask


def plot_input_mode(phi, x, y, neff, zoom_um=None):
    plt.figure(figsize=(6, 5))
    plt.imshow(
        np.abs(phi).T[::-1, :],
        origin="lower",
        extent=[x[0] * 1e6, x[-1] * 1e6, y[0] * 1e6, y[-1] * 1e6],
        aspect="equal",
    )
    plt.xlabel("x (um)")
    plt.ylabel("y (um)")
    plt.title(f"Input Ex semi-vector mode, neff={neff}")
    plt.colorbar(label="|Ex|")

    if zoom_um is not None:
        plt.xlim(-zoom_um, zoom_um)
        plt.ylim(-zoom_um, zoom_um)

    plt.tight_layout()


def animate_frames(frames, x, y, dz, step_stride, repeat=True):
    fig, ax = plt.subplots(figsize=(6, 5))

    im = ax.imshow(
        frames[0].T[::-1, :],
        origin="lower",
        extent=[x[0] * 1e6, x[-1] * 1e6, y[0] * 1e6, y[-1] * 1e6],
        aspect="equal",
    )
    cbar = plt.colorbar(im, ax=ax)
    cbar.set_label("Intensity")

    ax.set_xlabel("x (um)")
    ax.set_ylabel("y (um)")
    title = ax.set_title("")

    vmax = max(np.max(frame) for frame in frames)
    if vmax > 0:
        im.set_clim(0, vmax)

    def update(frame_idx):
        im.set_data(frames[frame_idx].T[::-1, :])
        z_now_um = frame_idx * step_stride * dz * 1e6
        title.set_text(f"Ex semi-vector movie, z = {z_now_um:.2f} um")
        return [im, title]

    ani = FuncAnimation(
        fig,
        update,
        frames=len(frames),
        interval=100,
        blit=False,
        repeat=repeat,
    )
    plt.tight_layout()
    plt.show()
    return ani


def compute_mode_overlap(E, E0, dx, dy):
    num = np.abs(np.sum(E * np.conjugate(E0)) * dx * dy) ** 2
    den = (
        np.sum(np.abs(E) ** 2) * dx * dy
        * np.sum(np.abs(E0) ** 2) * dx * dy
    )
    if den == 0:
        return 0.0
    return np.real(num / den)


def run_semivector_ex_mode_launch(
    Nz=1000,
    dz=0.05e-6,
    nmodes=1,
    mode_index=0,
    n_eff_guess=None,
    use_absorbing_mask=False,
    mask_width=18,
    mask_strength=3.5,
    save_movie_frames=True,
    step_stride=20,
):
    wavelength = 1.55e-6

    structure = build_rib_waveguide(
        lam0=wavelength,
        spacer=wavelength * np.array([3.0, 3.0, 2.0, 2.0]),
    )

    x = structure["x"]
    y = structure["y"]
    dx = structure["dx"]
    dy = structure["dy"]
    Nx = structure["Nx"]
    Ny = structure["Ny"]
    n = structure["n_map"]

    if n_eff_guess is None:
        n_eff_guess = np.max(n) - 0.05

    ex_result = solve_ex_semivector_modes(
        n=n,
        wavelength=wavelength,
        dx=dx,
        dy=dy,
        nmodes=nmodes,
        n_eff_guess=n_eff_guess,
    )

    E0 = ex_result["modes"][mode_index].copy()
    E = E0.copy()

    n_ref = np.real(ex_result["neff"][mode_index])

    print("Ex semi-vector self-consistent mode solve:")
    print(f"  n_eff_guess = {n_eff_guess}")
    print(f"  selected neff = {ex_result['neff'][mode_index]}")

    if use_absorbing_mask:
        mask = make_absorbing_mask(Nx, Ny, width=mask_width, strength=mask_strength)
    else:
        mask = np.ones((Nx, Ny), dtype=float)

    z = np.arange(Nz) * dz

    power_hist = np.zeros(Nz, dtype=float)
    center_hist = np.zeros(Nz, dtype=float)
    overlap_hist = np.zeros(Nz, dtype=float)
    xz_map = np.zeros((Nz, Nx), dtype=float)

    frames = []

    p0 = np.sum(np.abs(E0) ** 2) * dx * dy

    for iz in range(Nz):
        E = semivector_ex_step(E, n, wavelength, n_ref, dx, dy, dz)
        E *= mask

        power_hist[iz] = np.sum(np.abs(E) ** 2) * dx * dy
        center_hist[iz] = np.abs(E[Nx // 2, Ny // 2]) ** 2
        overlap_hist[iz] = compute_mode_overlap(E, E0, dx, dy)
        xz_map[iz, :] = np.abs(E[:, Ny // 2]) ** 2

        if iz % 50 == 0:
            y_min_edge = np.sum(np.abs(E[:, :10]) ** 2)
            y_max_edge = np.sum(np.abs(E[:, -10:]) ** 2)
            print(
                f"z={iz*dz*1e6:7.2f} um  "
                f"power={power_hist[iz]:12.5e}  "
                f"overlap={overlap_hist[iz]:10.6f}  "
                f"ymin={y_min_edge:12.5e}  "
                f"ymax={y_max_edge:12.5e}  "
                f"center={center_hist[iz]:12.5e}"
            )

        if save_movie_frames and (iz % step_stride == 0):
            frames.append(np.abs(E) ** 2)

    # input mode plot: zoom to -0.5 ~ 0.5 um
    plot_input_mode(E0, x, y, ex_result["neff"][mode_index], zoom_um=1.0)

    # output mode plot: zoom to -0.5 ~ 0.5 um
    plt.figure(figsize=(6, 5))
    plt.imshow(
        (np.abs(E).T ** 2)[::-1, :],
        origin="lower",
        extent=[x[0] * 1e6, x[-1] * 1e6, y[0] * 1e6, y[-1] * 1e6],
        aspect="equal",
    )
    plt.xlabel("x (um)")
    plt.ylabel("y (um)")
    plt.title("Output intensity (Ex semi-vector)")
    plt.colorbar(label="Intensity")
    plt.xlim(-1, 1)
    plt.ylim(-1, 1)
    plt.tight_layout()

    plt.figure(figsize=(7, 5))
    plt.imshow(
        xz_map.T,
        origin="lower",
        aspect="auto",
        extent=[z[0] * 1e6, z[-1] * 1e6, x[0] * 1e6, x[-1] * 1e6],
    )
    plt.xlabel("z (um)")
    plt.ylabel("x (um)")
    plt.title("xz propagation map (Ex semi-vector)")
    plt.colorbar(label="Intensity")
    plt.tight_layout()

    # normalized power
    normalized_power = power_hist / p0
    normalized_power = np.clip(normalized_power, 0.0, 1.0)

    plt.figure(figsize=(6, 4))
    plt.plot(z * 1e6, normalized_power)
    plt.xlabel("z (um)")
    plt.ylabel("Normalized power")
    plt.title("Power history")
    plt.ylim(0, 1.1)
    plt.grid(True)
    plt.tight_layout()

    # overlap
    overlap_plot = np.clip(overlap_hist, 0.0, 1.0)

    plt.figure(figsize=(6, 4))
    plt.plot(z * 1e6, overlap_plot)
    plt.xlabel("z (um)")
    plt.ylabel("Mode overlap with input")
    plt.title("Mode overlap history")
    plt.ylim(0, 1.1)
    plt.grid(True)
    plt.tight_layout()

    plt.figure(figsize=(6, 4))
    plt.plot(z * 1e6, center_hist)
    plt.xlabel("z (um)")
    plt.ylabel("|E(center)|^2")
    plt.title("Center intensity history")
    plt.grid(True)
    plt.tight_layout()

    plt.show()

    ani = None
    if save_movie_frames and len(frames) > 0:
        ani = animate_frames(frames, x, y, dz, step_stride, repeat=True)

    return {
        "E_in": E0,
        "E_out": E,
        "mode_result": ex_result,
        "x": x,
        "y": y,
        "z": z,
        "power_hist": power_hist,
        "overlap_hist": overlap_hist,
        "center_hist": center_hist,
        "xz_map": xz_map,
        "frames": frames,
        "animation": ani,
    }


if __name__ == "__main__":
    result = run_semivector_ex_mode_launch(
        Nz=1000,
        dz=0.05e-6,
        nmodes=1,
        mode_index=0,
        n_eff_guess=None,
        use_absorbing_mask=False,
        save_movie_frames=True,
        step_stride=20,
    )

semi_vector_bpm_ex.py

import numpy as np


def thomas_solve_complex(a, b, c, d):

    N = len(b)

    ac = a.astype(np.complex128).copy()
    bc = b.astype(np.complex128).copy()
    cc = c.astype(np.complex128).copy()
    dc = d.astype(np.complex128).copy()

    for i in range(1, N):
        w = ac[i-1] / bc[i-1]
        bc[i] -= w * cc[i-1]
        dc[i] -= w * dc[i-1]

    x = np.zeros(N, dtype=np.complex128)

    x[-1] = dc[-1] / bc[-1]

    for i in range(N-2, -1, -1):
        x[i] = (dc[i] - cc[i] * x[i+1]) / bc[i]

    return x


# -----------------------------
# weighted x operator
# -----------------------------
def apply_weighted_x(E, eps, dx):

    Nx, Ny = E.shape
    out = np.zeros_like(E, dtype=np.complex128)

    for i in range(1, Nx-1):

        eps_p = 2/(eps[i,:] + eps[i+1,:])
        eps_m = 2/(eps[i,:] + eps[i-1,:])

        out[i,:] = (
            eps_p * (E[i+1,:] - E[i,:])
            - eps_m * (E[i,:] - E[i-1,:])
        ) / dx**2

    return out


# -----------------------------
# scalar y laplacian
# -----------------------------
def apply_laplacian_y(E, dy):

    out = np.zeros_like(E, dtype=np.complex128)

    out[:,1:-1] = (E[:,2:] - 2*E[:,1:-1] + E[:,:-2]) / dy**2

    return out


# -----------------------------
# implicit x solve (weighted)
# -----------------------------
def implicit_solve_x(E_rhs, eps, dx, alpha):

    Nx, Ny = E_rhs.shape

    E_out = np.zeros_like(E_rhs, dtype=np.complex128)

    for j in range(Ny):

        a = np.zeros(Nx-1, dtype=np.complex128)
        b = np.ones(Nx, dtype=np.complex128)
        c = np.zeros(Nx-1, dtype=np.complex128)

        d = E_rhs[:,j].astype(np.complex128).copy()

        for i in range(1, Nx-1):

            eps_p = 2/(eps[i,j] + eps[i+1,j])
            eps_m = 2/(eps[i,j] + eps[i-1,j])

            lower = eps_m/dx**2
            diag  = -(eps_p + eps_m)/dx**2
            upper = eps_p/dx**2

            a[i-1] = -alpha * lower
            b[i]   = 1 - alpha * diag
            c[i]   = -alpha * upper

        b[0] = 1
        b[-1] = 1

        d[0] = 0
        d[-1] = 0

        if Nx>1:
            c[0] = 0
            a[-1] = 0

        E_out[:,j] = thomas_solve_complex(a,b,c,d)

    return E_out


# -----------------------------
# scalar implicit y
# -----------------------------
def implicit_solve_y(E_rhs, dy, alpha):

    Nx, Ny = E_rhs.shape

    E_out = np.zeros_like(E_rhs, dtype=np.complex128)

    lower = 1/dy**2
    diag  = -2/dy**2
    upper = 1/dy**2

    for i in range(Nx):

        a = np.zeros(Ny-1, dtype=np.complex128)
        b = np.ones(Ny, dtype=np.complex128)
        c = np.zeros(Ny-1, dtype=np.complex128)

        d = E_rhs[i,:].astype(np.complex128).copy()

        for j in range(1, Ny-1):

            a[j-1] = -alpha*lower
            b[j]   = 1-alpha*diag
            c[j]   = -alpha*upper

        b[0] = 1
        b[-1] = 1

        d[0] = 0
        d[-1] = 0

        if Ny>1:
            c[0] = 0
            a[-1] = 0

        E_out[i,:] = thomas_solve_complex(a,b,c,d)

    return E_out


# -----------------------------
# semi-vector Ex BPM step
# -----------------------------
def semivector_ex_step(E, n, wavelength, n_ref, dx, dy, dz):

    k0 = 2*np.pi / wavelength
    beta_ref = k0*n_ref

    eps = n**2
    eps_ref = n_ref**2

    alpha = dz/(4j*beta_ref)

    V = k0**2*(eps-eps_ref)

    E = E.copy()

    E[0,:]=0
    E[-1,:]=0
    E[:,0]=0
    E[:,-1]=0

    phase_half = np.exp(dz*V/(4j*beta_ref))

    E1 = E*phase_half

    Ly = apply_laplacian_y(E1,dy)

    rhs1 = E1 + alpha*Ly

    E2 = implicit_solve_x(rhs1,eps,dx,alpha)

    Lx = apply_weighted_x(E2,eps,dx)

    rhs2 = E2 + alpha*Lx

    E3 = implicit_solve_y(rhs2,dy,alpha)

    E4 = E3*phase_half

    E4[0,:]=0
    E4[-1,:]=0
    E4[:,0]=0
    E4[:,-1]=0

    return E4

semivector_ex_mode_solver.py

import numpy as np
import scipy.sparse as sp
import scipy.sparse.linalg as spla
import matplotlib.pyplot as plt

from mode_solver import build_rib_waveguide


def _idx(i, j, Ny):
    return i * Ny + j


def normalize_mode(phi, dx, dy):
    p = np.sum(np.abs(phi) ** 2) * dx * dy
    if p > 0:
        phi = phi / np.sqrt(p)
    return phi


def build_weighted_x_operator(n, dx):
    """
    Ex 用 semi-vector x 演算子:
        d/dx ( s_x d/dx )
    ここで
        s_x(i+1/2,j) = 2 / (eps(i,j) + eps(i+1,j))
        eps = n^2

    Dirichlet 境界。
    """
    Nx, Ny = n.shape
    M = Nx * Ny
    A = sp.lil_matrix((M, M), dtype=np.complex128)

    eps = n**2
    inv_dx2 = 1.0 / dx**2

    for i in range(1, Nx - 1):
        for j in range(Ny):
            p = _idx(i, j, Ny)

            s_p = 2.0 / (eps[i, j] + eps[i + 1, j])
            s_m = 2.0 / (eps[i, j] + eps[i - 1, j])

            A[p, _idx(i - 1, j, Ny)] += s_m * inv_dx2
            A[p, p] += -(s_p + s_m) * inv_dx2
            A[p, _idx(i + 1, j, Ny)] += s_p * inv_dx2

    return A.tocsc()


def build_scalar_y_operator(Nx, Ny, dy):
    """
    通常の d2/dy2
    Dirichlet 境界。
    """
    M = Nx * Ny
    A = sp.lil_matrix((M, M), dtype=np.complex128)

    inv_dy2 = 1.0 / dy**2

    for i in range(Nx):
        for j in range(1, Ny - 1):
            p = _idx(i, j, Ny)

            A[p, _idx(i, j - 1, Ny)] += inv_dy2
            A[p, p] += -2.0 * inv_dy2
            A[p, _idx(i, j + 1, Ny)] += inv_dy2

    return A.tocsc()


def build_index_operator(n, wavelength):
    """
    対角項 k0^2 n^2
    """
    k0 = 2.0 * np.pi / wavelength
    D = (k0**2) * (n**2)
    return sp.diags(D.reshape(-1), format="csc")


def build_ex_mode_operator(n, wavelength, dx, dy):
    """
    Ex quasi-TE semi-vector 断面演算子:
        A = weighted_x + scalar_y + k0^2 n^2
    固有値問題:
        A phi = beta^2 phi
    """
    Nx, Ny = n.shape
    Ax = build_weighted_x_operator(n, dx)
    Ay = build_scalar_y_operator(Nx, Ny, dy)
    An = build_index_operator(n, wavelength)
    return (Ax + Ay + An).tocsc()


def solve_ex_semivector_modes(
    n,
    wavelength,
    dx,
    dy,
    nmodes=1,
    n_eff_guess=None,
    which="LM",
):
    """
    Ex 用 quasi-TE semi-vector モードを解く。

    Parameters
    ----------
    n : ndarray (Nx, Ny)
        屈折率分布
    wavelength : float
        真空波長 [m]
    dx, dy : float
        格子間隔 [m]
    nmodes : int
        求めるモード数
    n_eff_guess : float or None
        シフトインバート用の neff 推定値
        None なら max(n) の少し下を使う
    which : str
        scipy eigs の which。通常は LM でよい

    Returns
    -------
    result : dict
        {
            "beta2": ...,
            "beta": ...,
            "neff": ...,
            "modes": [phi0, phi1, ...],
            "operator": A
        }
    """
    Nx, Ny = n.shape
    k0 = 2.0 * np.pi / wavelength

    A = build_ex_mode_operator(n, wavelength, dx, dy)

    if n_eff_guess is None:
        n_eff_guess = np.max(n) - 1e-3

    sigma = (k0 * n_eff_guess) ** 2

    beta2, vecs = spla.eigs(A, k=nmodes, sigma=sigma, which=which)

    # beta = sqrt(beta^2)
    beta = np.sqrt(beta2)
    neff = beta / k0

    # guided mode っぽい順に real(neff) 降順でソート
    order = np.argsort(-np.real(neff))
    beta2 = beta2[order]
    beta = beta[order]
    neff = neff[order]
    vecs = vecs[:, order]

    modes = []
    for m in range(nmodes):
        phi = vecs[:, m].reshape((Nx, Ny))

        # グローバル位相を中心で揃える
        cx, cy = Nx // 2, Ny // 2
        ref = phi[cx, cy]
        if np.abs(ref) > 0:
            phi = phi * np.exp(-1j * np.angle(ref))

        phi = normalize_mode(phi, dx, dy)
        modes.append(phi)

    return {
        "beta2": beta2,
        "beta": beta,
        "neff": neff,
        "modes": modes,
        "operator": A,
    }


def plot_ex_modes(result, x, y, title_prefix="Ex semi-vector mode"):
    nmodes = len(result["modes"])

    plt.figure(figsize=(5 * nmodes, 4))

    for m in range(nmodes):
        phi = result["modes"][m]

        plt.subplot(1, nmodes, m + 1)
        plt.imshow(
            np.abs(phi).T[::-1, :],
            origin="lower",
            extent=[x[0] * 1e6, x[-1] * 1e6, y[0] * 1e6, y[-1] * 1e6],
            aspect="equal",
        )
        plt.xlabel("x (um)")
        plt.ylabel("y (um)")
        plt.title(f"{title_prefix} {m}\nneff={result['neff'][m]}")
        plt.colorbar(label="|Ex|")

    plt.tight_layout()
    plt.show()
    


def run_ex_mode_solver_demo(
    wavelength=1.55e-6,
    nmodes=3,
    n_eff_guess=None,
    zoom_um=1.0,
):
    structure = build_rib_waveguide(
        lam0=wavelength,
        spacer=wavelength * np.array([3.0, 3.0, 2.0, 2.0]),
    )

    x = structure["x"]
    y = structure["y"]
    dx = structure["dx"]
    dy = structure["dy"]
    n = structure["n_map"]

    if n_eff_guess is None:
        n_eff_guess = np.max(n) - 0.05

    result = solve_ex_semivector_modes(
        n=n,
        wavelength=wavelength,
        dx=dx,
        dy=dy,
        nmodes=nmodes,
        n_eff_guess=n_eff_guess,
    )

    print("Ex semi-vector mode solve")
    print(f"wavelength = {wavelength*1e6:.3f} um")
    print(f"n_eff_guess = {n_eff_guess}")
    for i, neff in enumerate(result["neff"]):
        print(f"mode {i}: neff = {neff}")

    # 屈折率分布
    plt.figure(figsize=(6, 5))
    plt.imshow(
        n.T[::-1, :],
        origin="lower",
        extent=[x[0] * 1e6, x[-1] * 1e6, y[0] * 1e6, y[-1] * 1e6],
        aspect="equal",
    )
    plt.xlabel("x (um)")
    plt.ylabel("y (um)")
    plt.title("Refractive index distribution")
    plt.colorbar(label="n")
    if zoom_um is not None:
        plt.xlim(-zoom_um, zoom_um)
        plt.ylim(-zoom_um, zoom_um)
    plt.tight_layout()

    # 各モード表示
    nm = len(result["modes"])
    plt.figure(figsize=(5 * nm, 4))

    for m in range(nm):
        phi = result["modes"][m]
        plt.subplot(1, nm, m + 1)
        plt.imshow(
            np.abs(phi).T[::-1, :],
            origin="lower",
            extent=[x[0] * 1e6, x[-1] * 1e6, y[0] * 1e6, y[-1] * 1e6],
            aspect="equal",
        )
        plt.xlabel("x (um)")
        plt.ylabel("y (um)")
        plt.title(f"Mode {m}\nneff={result['neff'][m]}")
        plt.colorbar(label="|Ex|")
        if zoom_um is not None:
            plt.xlim(-zoom_um, zoom_um)
            plt.ylim(-zoom_um, zoom_um)

    plt.tight_layout()
    plt.show()

    return {
        "structure": structure,
        "mode_result": result,
    }


if __name__ == "__main__":
    demo = run_ex_mode_solver_demo(
        wavelength=1.55e-6,
        nmodes=3,
        n_eff_guess=None,
        zoom_um=1.0,
    )    

mode_solver.py

import numpy as np
import scipy.sparse as sp
import scipy.sparse.linalg as spla

#from yeeder2d import yeeder2d


def build_rib_waveguide(
    lam0=1.55e-6,
    rib_n1=1.0,
    rib_n2=3.5,
    rib_n3=1.5,
    rib_h=0.10e-6,
    rib_t=0.12e-6,
    rib_w=0.50e-6,
    nres=10,
    spacer=None,
):
    """
    Build a rib waveguide on a Yee-like doubled grid, then extract Nx x Ny maps.
    """
    if spacer is None:
        spacer = lam0 * np.array([1.5, 1.5, 0.5, 0.5])

    nmax = max(rib_n1, rib_n2, rib_n3)

    dx = lam0 / nmax / nres
    dy = lam0 / nmax / nres

    nx = int(np.ceil(rib_w / dx))
    dx = rib_w / nx
    ny = int(np.ceil(rib_t / dy))
    dy = rib_t / ny

    sx = spacer[0] + rib_w + spacer[1]
    Nx = int(np.ceil(sx / dx))
    sx = Nx * dx

    sy = spacer[2] + rib_h + rib_t + spacer[3]
    Ny = int(np.ceil(sy / dy))
    sy = Ny * dy

    Nx2, Ny2 = 2 * Nx, 2 * Ny
    dx2, dy2 = dx / 2.0, dy / 2.0

    x = (np.arange(1, Nx + 1) * dx)
    x = x - np.mean(x)
    y = (np.arange(1, Ny + 1) * dy)
    y = y - np.mean(y)

    ER2 = (rib_n1**2) * np.ones((Nx2, Ny2), dtype=float)
    UR2 = np.ones((Nx2, Ny2), dtype=float)

    nx1 = 1 + int(round(spacer[0] / dx2))
    nx2 = nx1 + int(round(rib_w / dx2)) - 1

    ny1 = 1 + int(round(spacer[2] / dy2))
    ny2 = ny1 + int(round(rib_h / dy2)) - 1
    ny3 = ny2 + 1
    ny4 = ny3 + int(round(rib_t / dy2)) - 1
    ny5 = ny4 + 1

    nx1p, nx2p = nx1 - 1, nx2 - 1
    ny1p, ny2p = ny1 - 1, ny2 - 1
    ny3p, ny4p = ny3 - 1, ny4 - 1
    ny5p = ny5 - 1

    # Rib region
    ER2[nx1p:nx2p + 1, ny1p:ny2p + 1] = rib_n2**2
    # Slab region
    ER2[:, ny3p:ny4p + 1] = rib_n2**2
    # Lower substrate
    ER2[:, ny5p:Ny2] = rib_n3**2

    ERxx = ER2[1::2, 0::2]
    ERyy = ER2[0::2, 1::2]
    ERzz = ER2[0::2, 0::2]

    URxx = UR2[0::2, 1::2]
    URyy = UR2[1::2, 0::2]
    URzz = UR2[1::2, 1::2]

    return {
        "x": x,
        "y": y,
        "dx": dx,
        "dy": dy,
        "Nx": Nx,
        "Ny": Ny,
        "ERxx": ERxx,
        "ERyy": ERyy,
        "ERzz": ERzz,
        "URxx": URxx,
        "URyy": URyy,
        "URzz": URzz,
        "n_map": np.sqrt(ERzz),
        "core_index": rib_n2,
    }

semivector_ex_mode_solver.py単体でも動くので、リブ型導波路のサイズや、屈折率をかえて試してみるのも良いでしょう。
断面方向の計算グリッドはmode_solver.pyのnresで決めてます。

10.シミュレーション結果

こちらがTEモード入射時の伝搬シミュレーションの結果になります。
今回は50 um伝搬させてみました。
入射モード分布

出射モード分布

光の伝搬の様子(powerをプロット)

光伝搬パワーのグラフ表示(長手方向の導波路パワーを各位置でプロット)

光伝搬重なり積分のグラフ表示(長手方向の分布と入射光との重なり積分を各位置でプロット)

どれも計算できていそうです。これぐらいのプログラムであれば、生成AIを活用すれば、数時間でできてしまいました。

次回以降はTMモードの実施例やGPUを使った自作コード、曲げ導波路のシミュレーション事例も紹介していきます

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

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

問い合わせ

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

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

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

参考

• Okamoto – Fundamentals of Optical Waveguides
• Hadley – Transparent Boundary Condition for BPM
• Feit & Fleck – Beam Propagation Method