FEM mode solver

FEM mode solver#

To understand the propagation in a waveguide, it’s often helpful to look at how the effective refractive indices of the modes change when adjusting the width of the waveguide. As the modes continuously change their refractive index, we can track them through their evolution. We gray the area which would indicate a effective refractive index below the refractive index of the underlying layer (commonly referred to as box), as such modes would not be guided. The refractive index of the box called “cutoff”, as it defines the effective refractive index level under which modes stop being guided.

Hide code cell source

from collections import OrderedDict

import matplotlib.pyplot as plt
import numpy as np
from femwell.maxwell.waveguide import compute_modes
from femwell.mesh import mesh_from_OrderedDict
from shapely.geometry import Polygon, box
from shapely.ops import clip_by_rect
from skfem import Basis, ElementTriP0
from skfem.io.meshio import from_meshio
from tqdm import tqdm
wavelength = 1.55  # wavelength of the light
ncore = 1.9963  # refractive index of the core (SiN)
nclad = 1.444  # refractive index of the cladding (SiO2)
num_modes = 8
wg_thickness = 0.33
widths = np.linspace(0.5, 3.5, 100)

all_neffs = np.zeros((widths.shape[0], num_modes))
all_te_fracs = np.zeros((widths.shape[0], num_modes))
for i, width in enumerate(tqdm(widths)):
    core = box(0, 0, width, wg_thickness)
    polygons = OrderedDict(
        core=core,
        box=clip_by_rect(core.buffer(1.0, resolution=4), -np.inf, -np.inf, np.inf, 0),
        clad=clip_by_rect(core.buffer(1.0, resolution=4), -np.inf, 0, np.inf, np.inf),
    )

    resolutions = {"core": {"resolution": 0.1, "distance": 1}}

    mesh = from_meshio(
        mesh_from_OrderedDict(polygons, resolutions, default_resolution_max=0.6)
    )

    basis0 = Basis(mesh, ElementTriP0())
    epsilon = basis0.zeros(dtype=complex)
    for subdomain, n in {"core": ncore, "box": nclad, "clad": nclad}.items():
        epsilon[basis0.get_dofs(elements=subdomain)] = n**2

    modes = compute_modes(basis0, epsilon, wavelength=wavelength, num_modes=num_modes)
    all_neffs[i] = np.real([mode.n_eff for mode in modes])
    all_te_fracs[i, :] = [mode.te_fraction for mode in modes]

Hide code cell source

all_neffs = np.real(all_neffs)
plt.xlabel("Width of waveguide / µm")
plt.ylabel("Effective refractive index")
plt.fill_between(widths, nclad, alpha=0.5, color="gray")
plt.ylim(1.36, np.max(all_neffs) + 0.1 * (np.max(all_neffs) - 1.444))
for lams, te_fracs in zip(all_neffs.T, all_te_fracs.T):
    plt.plot(widths, lams)
    plt.scatter(widths, lams, c=te_fracs, cmap="cool")
plt.colorbar().set_label("TE fraction")
plt.title("Effective refractive index vs width of waveguide of SiN")
plt.show()
../_images/002e1964e94a62f4eaa2fafee273b0f5c4b80eab06d0cde5f3f9d93c5a13db1b.png

The graph shows a TE mode emerging (getting a greater effective refractive index than the box) at a width of ~750nm. The second mode emerges at a width of ~1700nm. Thus, the waveguide is a single mode waveguide with a width within the range from ~750nm to 1700nm.

The second mode is at the width at which it emerges a TM-mode, but shows for slightly wider waveguides an anti-crossing with the emerging TE-mode. After this anti-crossing, the mode with the second highest effective refracteive index is TE10-mode, while the TM00-mode keeps existing in the system with an effective refractive index slightly above the cutoff. At a width of ~2600nm, a third TE-mode starts being guided.

wavelength = 1.55  # wavelength of the light
ncore = 3.47  # refractive index of the core (Si)
nclad = 1.444  # refractive index of the cladding (SiO2)
num_modes = 6
widths = np.array([0.5])
wg_thickness = 0.22
sidewall_angle = 0  # sidewall angle of the waveguide (in degrees)
slab_thickness = 0
default_resolution_max = 0.08

core_resolution = 0.02
core_distance = 1
buffer_size = 1.0


all_neffs = np.zeros((widths.shape[0], num_modes))
all_te_fracs = np.zeros((widths.shape[0], num_modes))
for i, width in enumerate(tqdm(widths)):
    top_width = width
    bottom_width = 2 * np.tan(np.deg2rad(sidewall_angle)) + width
    total_slab_width = 4 + bottom_width

    if slab_thickness > 0:
        core = Polygon(
            [
                (-total_slab_width / 2, 0),
                (+total_slab_width / 2, 0),
                (+total_slab_width / 2, slab_thickness),
                (+bottom_width / 2, slab_thickness),
                (+top_width / 2, wg_thickness),
                (-top_width / 2, wg_thickness),
                (-bottom_width / 2, slab_thickness),
                (-total_slab_width / 2, slab_thickness),
                (-total_slab_width / 2, 0),
            ]
        )
    else:
        core = Polygon(
            [
                (-bottom_width / 2, 0),
                (+bottom_width / 2, 0),
                (+top_width / 2, wg_thickness),
                (-top_width / 2, wg_thickness),
                (-bottom_width / 2, 0),
            ]
        )

    polygons = OrderedDict(
        core=core,
        box=clip_by_rect(
            core.buffer(buffer_size, resolution=1), -np.inf, -np.inf, np.inf, 0
        ),
        clad=clip_by_rect(
            core.buffer(buffer_size, resolution=1), -np.inf, 0, np.inf, np.inf
        ),
    )

    resolutions = {"core": {"resolution": core_resolution, "distance": core_distance}}

    mesh = from_meshio(
        mesh_from_OrderedDict(
            polygons, resolutions, default_resolution_max=default_resolution_max
        )
    )

    resolutions = {"core": {"resolution": core_resolution, "distance": core_distance}}
    mesh = from_meshio(
        mesh_from_OrderedDict(
            polygons, resolutions, default_resolution_max=default_resolution_max
        )
    )

    basis0 = Basis(mesh, ElementTriP0())
    epsilon = basis0.zeros(dtype=complex)
    for subdomain, n in {"core": ncore, "box": nclad, "clad": nclad}.items():
        epsilon[basis0.get_dofs(elements=subdomain)] = n**2

    modes = compute_modes(basis0, epsilon, wavelength=wavelength, num_modes=num_modes)
    all_neffs[i] = np.real([mode.n_eff for mode in modes])
    all_te_fracs[i, :] = [mode.te_fraction for mode in modes]
  0%|          | 0/1 [00:00<?, ?it/s]
100%|██████████| 1/1 [00:02<00:00,  2.57s/it]
100%|██████████| 1/1 [00:02<00:00,  2.57s/it]

modes[0].show("E", part="real")
../_images/3393bf566a31bfcc208d78abd0100fa64de7ad026bd10147edf40c13ed45df0e.png
modes[0].plot_component("E", component="x", part="real", colorbar=True)
<Axes: title={'center': 'Ex (real. part)'}>
../_images/408c97c1ceba69ccdded0536567e2a767ec2cfb07b54e87bd054c0744ecce1e4.png
wavelength = 1.55  # wavelength of the light
ncore = 3.47  # refractive index of the core (Si)
nclad = 1.444  # refractive index of the cladding (SiO2)
num_modes = 6
widths = np.linspace(0.2, 1.0, 100)  # width of the waveguide
wg_thickness = 0.22
sidewall_angle = 0  # sidewall angle of the waveguide (in degrees)
slab_thickness = 0
default_resolution_max = 0.6

core_resolution = 0.08
core_distance = 1
buffer_size = 1.0


all_neffs = np.zeros((widths.shape[0], num_modes))
all_te_fracs = np.zeros((widths.shape[0], num_modes))
for i, width in enumerate(tqdm(widths)):
    top_width = width
    bottom_width = 2 * np.tan(np.deg2rad(sidewall_angle)) + width
    total_slab_width = 4 + bottom_width

    if slab_thickness > 0:
        core = Polygon(
            [
                (-total_slab_width / 2, 0),
                (+total_slab_width / 2, 0),
                (+total_slab_width / 2, slab_thickness),
                (+bottom_width / 2, slab_thickness),
                (+top_width / 2, wg_thickness),
                (-top_width / 2, wg_thickness),
                (-bottom_width / 2, slab_thickness),
                (-total_slab_width / 2, slab_thickness),
                (-total_slab_width / 2, 0),
            ]
        )
    else:
        core = Polygon(
            [
                (-bottom_width / 2, 0),
                (+bottom_width / 2, 0),
                (+top_width / 2, wg_thickness),
                (-top_width / 2, wg_thickness),
                (-bottom_width / 2, 0),
            ]
        )

    polygons = OrderedDict(
        core=core,
        box=clip_by_rect(
            core.buffer(buffer_size, resolution=4), -np.inf, -np.inf, np.inf, 0
        ),
        clad=clip_by_rect(
            core.buffer(buffer_size, resolution=4), -np.inf, 0, np.inf, np.inf
        ),
    )

    resolutions = {"core": {"resolution": core_resolution, "distance": core_distance}}

    mesh = from_meshio(
        mesh_from_OrderedDict(
            polygons, resolutions, default_resolution_max=default_resolution_max
        )
    )

    resolutions = {"core": {"resolution": core_resolution, "distance": core_distance}}
    mesh = from_meshio(
        mesh_from_OrderedDict(
            polygons, resolutions, default_resolution_max=default_resolution_max
        )
    )

    basis0 = Basis(mesh, ElementTriP0())
    epsilon = basis0.zeros(dtype=complex)
    for subdomain, n in {"core": ncore, "box": nclad, "clad": nclad}.items():
        epsilon[basis0.get_dofs(elements=subdomain)] = n**2

    modes = compute_modes(basis0, epsilon, wavelength=wavelength, num_modes=num_modes)
    all_neffs[i] = np.real([mode.n_eff for mode in modes])
    all_te_fracs[i, :] = [mode.te_fraction for mode in modes]
  0%|          | 0/100 [00:00<?, ?it/s]
  1%|          | 1/100 [00:00<00:33,  2.93it/s]
  2%|▏         | 2/100 [00:00<00:32,  3.04it/s]
  3%|▎         | 3/100 [00:00<00:31,  3.08it/s]
  4%|▍         | 4/100 [00:01<00:31,  3.08it/s]
  5%|▌         | 5/100 [00:01<00:30,  3.08it/s]
  6%|▌         | 6/100 [00:01<00:30,  3.07it/s]
  7%|▋         | 7/100 [00:02<00:30,  3.03it/s]
  8%|▊         | 8/100 [00:02<00:30,  3.04it/s]
  9%|▉         | 9/100 [00:02<00:30,  3.03it/s]
 10%|█         | 10/100 [00:03<00:29,  3.02it/s]
 11%|█         | 11/100 [00:03<00:29,  3.03it/s]
 12%|█▏        | 12/100 [00:03<00:29,  3.03it/s]
 13%|█▎        | 13/100 [00:04<00:28,  3.04it/s]
 14%|█▍        | 14/100 [00:04<00:28,  3.05it/s]
 15%|█▌        | 15/100 [00:04<00:27,  3.05it/s]
 16%|█▌        | 16/100 [00:05<00:27,  3.02it/s]
 17%|█▋        | 17/100 [00:05<00:27,  2.97it/s]
 18%|█▊        | 18/100 [00:05<00:28,  2.92it/s]
 19%|█▉        | 19/100 [00:06<00:27,  2.92it/s]
 20%|██        | 20/100 [00:06<00:27,  2.92it/s]
 21%|██        | 21/100 [00:06<00:27,  2.92it/s]
 22%|██▏       | 22/100 [00:07<00:26,  2.92it/s]
 23%|██▎       | 23/100 [00:07<00:26,  2.92it/s]
 24%|██▍       | 24/100 [00:08<00:25,  2.93it/s]
 25%|██▌       | 25/100 [00:08<00:25,  2.92it/s]
 26%|██▌       | 26/100 [00:08<00:25,  2.91it/s]
 27%|██▋       | 27/100 [00:09<00:25,  2.91it/s]
 28%|██▊       | 28/100 [00:09<00:24,  2.90it/s]
 29%|██▉       | 29/100 [00:09<00:24,  2.91it/s]
 30%|███       | 30/100 [00:10<00:24,  2.91it/s]
 31%|███       | 31/100 [00:10<00:25,  2.72it/s]
 32%|███▏      | 32/100 [00:10<00:24,  2.76it/s]
 33%|███▎      | 33/100 [00:11<00:24,  2.79it/s]
 34%|███▍      | 34/100 [00:11<00:23,  2.81it/s]
 35%|███▌      | 35/100 [00:11<00:23,  2.82it/s]
 36%|███▌      | 36/100 [00:12<00:22,  2.81it/s]
 37%|███▋      | 37/100 [00:12<00:22,  2.79it/s]
 38%|███▊      | 38/100 [00:12<00:22,  2.78it/s]
 39%|███▉      | 39/100 [00:13<00:21,  2.79it/s]
 40%|████      | 40/100 [00:13<00:21,  2.79it/s]
 41%|████      | 41/100 [00:14<00:21,  2.78it/s]
 42%|████▏     | 42/100 [00:14<00:20,  2.77it/s]
 43%|████▎     | 43/100 [00:14<00:20,  2.77it/s]
 44%|████▍     | 44/100 [00:15<00:20,  2.78it/s]
 45%|████▌     | 45/100 [00:15<00:19,  2.77it/s]
 46%|████▌     | 46/100 [00:15<00:19,  2.75it/s]
 47%|████▋     | 47/100 [00:16<00:19,  2.73it/s]
 48%|████▊     | 48/100 [00:16<00:19,  2.72it/s]
 49%|████▉     | 49/100 [00:17<00:18,  2.71it/s]
 50%|█████     | 50/100 [00:17<00:18,  2.70it/s]
 51%|█████     | 51/100 [00:17<00:18,  2.69it/s]
 52%|█████▏    | 52/100 [00:18<00:17,  2.68it/s]
 53%|█████▎    | 53/100 [00:18<00:17,  2.67it/s]
 54%|█████▍    | 54/100 [00:18<00:17,  2.66it/s]
 55%|█████▌    | 55/100 [00:19<00:16,  2.67it/s]
 56%|█████▌    | 56/100 [00:19<00:16,  2.64it/s]
 57%|█████▋    | 57/100 [00:20<00:16,  2.64it/s]
 58%|█████▊    | 58/100 [00:20<00:15,  2.63it/s]
 59%|█████▉    | 59/100 [00:20<00:15,  2.62it/s]
 60%|██████    | 60/100 [00:21<00:15,  2.59it/s]
 61%|██████    | 61/100 [00:21<00:15,  2.58it/s]
 62%|██████▏   | 62/100 [00:21<00:14,  2.59it/s]
 63%|██████▎   | 63/100 [00:22<00:14,  2.60it/s]
 64%|██████▍   | 64/100 [00:22<00:13,  2.59it/s]
 65%|██████▌   | 65/100 [00:23<00:13,  2.58it/s]
 66%|██████▌   | 66/100 [00:23<00:13,  2.54it/s]
 67%|██████▋   | 67/100 [00:23<00:13,  2.53it/s]
 68%|██████▊   | 68/100 [00:24<00:12,  2.53it/s]
 69%|██████▉   | 69/100 [00:24<00:12,  2.53it/s]
 70%|███████   | 70/100 [00:25<00:11,  2.53it/s]
 71%|███████   | 71/100 [00:25<00:11,  2.53it/s]
 72%|███████▏  | 72/100 [00:25<00:11,  2.54it/s]
 73%|███████▎  | 73/100 [00:26<00:10,  2.55it/s]
 74%|███████▍  | 74/100 [00:26<00:10,  2.56it/s]
 75%|███████▌  | 75/100 [00:27<00:09,  2.57it/s]
 76%|███████▌  | 76/100 [00:27<00:09,  2.57it/s]
 77%|███████▋  | 77/100 [00:27<00:08,  2.56it/s]
 78%|███████▊  | 78/100 [00:28<00:08,  2.56it/s]
 79%|███████▉  | 79/100 [00:28<00:08,  2.54it/s]
 80%|████████  | 80/100 [00:29<00:07,  2.52it/s]
 81%|████████  | 81/100 [00:29<00:07,  2.51it/s]
 82%|████████▏ | 82/100 [00:29<00:07,  2.51it/s]
 83%|████████▎ | 83/100 [00:30<00:06,  2.50it/s]
 84%|████████▍ | 84/100 [00:30<00:06,  2.51it/s]
 85%|████████▌ | 85/100 [00:31<00:05,  2.52it/s]
 86%|████████▌ | 86/100 [00:31<00:05,  2.50it/s]
 87%|████████▋ | 87/100 [00:31<00:05,  2.49it/s]
 88%|████████▊ | 88/100 [00:32<00:04,  2.47it/s]
 89%|████████▉ | 89/100 [00:32<00:04,  2.45it/s]
 90%|█████████ | 90/100 [00:33<00:04,  2.45it/s]
 91%|█████████ | 91/100 [00:33<00:03,  2.44it/s]
 92%|█████████▏| 92/100 [00:33<00:03,  2.42it/s]
 93%|█████████▎| 93/100 [00:34<00:02,  2.41it/s]
 94%|█████████▍| 94/100 [00:34<00:02,  2.41it/s]
 95%|█████████▌| 95/100 [00:35<00:02,  2.42it/s]
 96%|█████████▌| 96/100 [00:35<00:01,  2.41it/s]
 97%|█████████▋| 97/100 [00:36<00:01,  2.38it/s]
 98%|█████████▊| 98/100 [00:36<00:00,  2.35it/s]
 99%|█████████▉| 99/100 [00:36<00:00,  2.35it/s]
100%|██████████| 100/100 [00:37<00:00,  2.36it/s]
100%|██████████| 100/100 [00:37<00:00,  2.68it/s]

all_neffs = np.real(all_neffs)
plt.title("Effective refractive index vs width of waveguide of Si")
plt.xlabel("Width of waveguide / µm")
plt.ylabel("Effective refractive index")
plt.fill_between(widths, nclad, alpha=0.5, color="gray")
plt.ylim(1.36, np.max(all_neffs) + 0.1 * (np.max(all_neffs) - 1.444))
for lams, te_fracs in zip(all_neffs.T, all_te_fracs.T):
    plt.plot(widths, lams)
    plt.scatter(widths, lams, c=te_fracs, cmap="cool")
plt.colorbar().set_label("TE fraction")
plt.show()
../_images/c0aee6cdd9eb37ed51a5e3decc7297fab48254f40a6f692dc7a5cb0555a1b0d3.png

The second order TE mode (pink) starts to be guided around 0.45 um, which is the single mode condition for these types of waveguides.