Effective area of a Si-NCs waveguide

Contents

Effective area of a Si-NCs waveguide#

In this example we calculate the effective area and the effective refractive index of a Si-NCs waveguide from [1]

Hide code cell source
from collections import OrderedDict

import matplotlib.pyplot as plt
import numpy as np
import shapely.affinity
from matplotlib.ticker import MultipleLocator
from skfem import Basis, ElementTriP0, Functional
from skfem.io.meshio import from_meshio

from femwell.maxwell.waveguide import compute_modes
from femwell.mesh import mesh_from_OrderedDict

We describe the geometry using shapely. In this case it’s simple: we use a shapely.box for the waveguide according to the parameter provided in the paper For the surrounding we buffer the core and clip it to the part above the waveguide for the air Width is a sweep parameter from 50nm to 700nm

wavelength = 1.55
capital_w = 1.4
capital_h = 0.3
h_list = [0.5, 0.7]
t = 0.1
n_silicon = 3.48
n_air = 1
n_nc = 1.6
n_silica = 1.45
w_list = [x for x in range(250, 700, 100)]
neff_dict = dict()
aeff_dict = dict()
aeff1_dict = dict()
aeff3_dict = dict()
tm_dict = dict()
p_dict = dict()
for h in h_list:
    neff_list = []
    aeff_list = []
    aeff1_list = []
    aeff3_list = []
    tm_list = []
    p_list = []
    for width in w_list:
        width = width * 1e-3
        nc = shapely.geometry.box(
            -width / 2, capital_h + (h - t) / 2, +width / 2, capital_h + (h - t) / 2 + t
        )
        silica = shapely.geometry.box(-capital_w / 2, 0, +capital_w / 2, capital_h)
        silicon = shapely.geometry.box(-width / 2, capital_h, +width / 2, capital_h + h)

        polygons = OrderedDict(
            core=nc,
            silicon=silicon,
            silica=silica,
            air=nc.buffer(10.0, resolution=4),
        )

        resolutions = dict(
            core={"resolution": 0.02, "distance": 0.3},
            silicon={"resolution": 0.02, "distance": 0.1},
            silica={"resolution": 0.04, "distance": 0.5},
        )

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

        basis0 = Basis(mesh, ElementTriP0())
        epsilon = basis0.zeros()
        for subdomain, n in {
            "core": n_nc,
            "silicon": n_silicon,
            "air": n_air,
            "silica": n_silica,
        }.items():
            epsilon[basis0.get_dofs(elements=subdomain)] = n**2
        modes = compute_modes(basis0, epsilon, wavelength=wavelength, num_modes=3, order=1)

        for mode in modes:
            if mode.tm_fraction > 0.5:
                # mode.show("E", part="real")
                print(f"Effective refractive index: {mode.n_eff:.4f}")
                print(f"Effective mode area: {mode.calculate_effective_area(field='y'):.4f}")
                print(f"Mode transversality: {mode.transversality}")
                neff_list.append(np.real(mode.n_eff))
                aeff_list.append(mode.calculate_effective_area())
                tm_list.append(mode.transversality)
                p_list.append(mode.Sz)

                @Functional
                def I(w):
                    return 1

                @Functional
                def Sz(w):
                    return w["Sz"]

                @Functional
                def Sz2(w):
                    return w["Sz"] ** 2

                Sz_basis, Sz_vec = mode.Sz

                int_Sz = Sz.assemble(Sz_basis, Sz=Sz_basis.interpolate(Sz_vec))
                print("int(Sz)", int_Sz)  # 1 as it's normalized
                int_I = I.assemble(Sz_basis.with_elements("core"))
                print("int_core(1)", int_I)  # area of core

                int_Sz_core = Sz.assemble(
                    Sz_basis.with_elements("core"),
                    Sz=Sz_basis.with_elements("core").interpolate(Sz_vec),
                )
                print(
                    "int_core(Sz)",
                    int_Sz_core,
                )
                int_Sz2 = Sz2.assemble(Sz_basis, Sz=Sz_basis.interpolate(Sz_vec))
                print("int(Sz^2)", int_Sz2)

                aeff1_list.append(int_Sz**2 / int_Sz2)
                aeff3_list.append(int_I * int_Sz / int_Sz_core)
                break
        else:
            print(f"no TM mode found for {width}")
    neff_dict[str(h)] = neff_list
    aeff1_dict[str(h)] = aeff1_list
    aeff_dict[str(h)] = aeff_list
    aeff3_dict[str(h)] = aeff3_list
    tm_dict[str(h)] = tm_list
    p_dict[str(h)] = p_list
Hide code cell output
Effective refractive index: 1.6358+0.0000j
Effective mode area: 0.1122
Mode transversality: 0.8486578847897985
int(Sz) (1.0000000035076144+0j)
int_core(1) 0.024999999999999994
int_core(Sz) (0.3234886208573645+0j)
int(Sz^2) (6.310950125433597+0j)
Effective refractive index: 1.8309+0.0000j
Effective mode area: 0.1062
Mode transversality: 0.8378664877597226
int(Sz) (0.9999999995843065+0j)
int_core(1) 0.03499999999999999
int_core(Sz) (0.39140610950315846+0j)
int(Sz^2) (6.466905439707694+0j)
Effective refractive index: 1.9606+0.0000j
Effective mode area: 0.1028
Mode transversality: 0.8042780606112331
int(Sz) (1.0000000031054088+0j)
int_core(1) 0.04499999999999999
int_core(Sz) (0.41572126808058796+0j)
int(Sz^2) (5.931250610125743+0j)
Effective refractive index: 2.0360+0.0000j
Effective mode area: 0.1123
Mode transversality: 0.8131366744098938
int(Sz) (1.0000000033379481+0j)
int_core(1) 0.05499999999999999
int_core(Sz) (0.4331627474807713+0j)
int(Sz^2) (5.39769872917229+0j)
Effective refractive index: 2.0930+0.0000j
Effective mode area: 0.1227
Mode transversality: 0.8068500151074046
int(Sz) (1.0000000030375438+0j)
int_core(1) 0.06499999999999999
int_core(Sz) (0.4403261083823158+0j)
int(Sz^2) (4.887389479005456+0j)
Effective refractive index: 2.0369+0.0000j
Effective mode area: 0.1339
Mode transversality: 0.825400266821272
int(Sz) (0.9999999975957938+0j)
int_core(1) 0.024999999999999994
int_core(Sz) (0.17897702102915045+0j)
int(Sz^2) (4.959302836190315+0j)
Effective refractive index: 2.2844+0.0000j
Effective mode area: 0.1496
Mode transversality: 0.8021836137528735
int(Sz) (1.0000000044053323+0j)
int_core(1) 0.03499999999999999
int_core(Sz) (0.179024592786723+0j)
int(Sz^2) (4.66607404632152+0j)
Effective refractive index: 2.4332+0.0000j
Effective mode area: 0.1699
Mode transversality: 0.779855375929004
int(Sz) (1.0000000024387858+0j)
int_core(1) 0.044999999999999984
int_core(Sz) (0.17233163217698483+0j)
int(Sz^2) (4.057757543716009+0j)
Effective refractive index: 2.5214+0.0000j
Effective mode area: 0.1943
Mode transversality: 0.7783595144584896
int(Sz) (1.0000000013221673+0j)
int_core(1) 0.05499999999999999
int_core(Sz) (0.17374432368397402+0j)
int(Sz^2) (3.5863428497694687+0j)
Effective refractive index: 2.5781+0.0000j
Effective mode area: 0.2166
Mode transversality: 0.7750886740287622
int(Sz) (1.000000001462982+0j)
int_core(1) 0.06499999999999999
int_core(Sz) (0.17348381212331038+0j)
int(Sz^2) (3.199788975358846+0j)

Plot the result

path = "../reference_data/Rukhlenko"

fig, axs = plt.subplots(4, 1, figsize=(9, 20))

for t1, t2, ax1, ax2 in [("0.5", "500", *axs[0:2]), ("0.7", "700", *axs[2:4])]:
    reference_neff = np.loadtxt(path + f"/fig_1c_neff/h_{t2}nm.csv", delimiter=",")
    reference_aeff1 = np.loadtxt(path + f"/fig_1b_aeff/{t1}_Eq1.csv", delimiter=",")
    reference_aeff2 = np.loadtxt(path + f"/fig_1b_aeff/{t1}_Eq2.csv", delimiter=",")
    reference_aeff3 = np.loadtxt(path + f"/fig_1b_aeff/{t1}_Eq3.csv", delimiter=",")
    reference_tm = np.loadtxt(path + f"/fig_1c_neff/tm_h_{t2}nm.csv", delimiter=",")

    ax1.scatter(w_list, aeff_dict[t1], c="r", label="stimulated, eq2")
    ax1.plot(reference_aeff2[:, 0], reference_aeff2[:, 1], c="b", label="reference, eq2")

    ax1.scatter(w_list, aeff1_dict[t1], c="green", label="stimulated, eq1")
    ax1.plot(reference_aeff1[:, 0], reference_aeff1[:, 1], c="orange", label="reference, eq1")

    ax1.scatter(w_list, aeff3_dict[t1], c="purple", label="stimulated, eq3")
    ax1.plot(reference_aeff3[:, 0], reference_aeff3[:, 1], c="brown", label="reference, eq3")

    ax1.set_title(f"aeff at h = {t2}nm")
    ax1.yaxis.set_major_locator(MultipleLocator(0.05))
    ax1.yaxis.set_minor_locator(MultipleLocator(0.01))
    ax1.set_ylim(0, 0.3)
    ax1.legend(loc="upper right")

    ax2b = ax2.twinx()
    ax2b.set_ylabel("mode transversality")
    ax2b.scatter(
        reference_tm[:, 0],
        reference_tm[:, 1],
        marker="v",
        c="b",
        label="reference transversality",
    )
    ax2b.plot(w_list, tm_dict[t1], "-v", c="r", label="stimulated transversality")
    ax2b.set_ylim(0.775, 1)
    ax2b.legend()

    ax2.plot(w_list, neff_dict[t1], "-o", c="r", label="stimulated neff")
    ax2.scatter(reference_neff[:, 0], reference_neff[:, 1], c="b", label="reference neff")
    ax2.set_title(f"neff and mode transversality at h = {t2}nm")
    ax2.set_ylabel("neff")
    ax2.yaxis.set_major_locator(MultipleLocator(0.4))
    ax2.yaxis.set_minor_locator(MultipleLocator(0.2))
    ax2.set_ylim(0, 2.8)
    ax2.legend()

plt.show()
/home/runner/miniconda3/lib/python3.12/site-packages/matplotlib/cbook.py:1762: ComplexWarning: Casting complex values to real discards the imaginary part
  return math.isfinite(val)
/home/runner/miniconda3/lib/python3.12/site-packages/matplotlib/collections.py:197: ComplexWarning: Casting complex values to real discards the imaginary part
  offsets = np.asanyarray(offsets, float)
../../_images/ed9c9a68b9154a47a4e8a40da05776701876a51d43a6d6511e34e6e861c8618b.png

We can also plot the modes to compare them with the reference article

# Data figure h
w_fig_h = 500
idx_fig_h = min(range(len(w_list)), key=lambda x: abs(w_list[x] - w_fig_h))
h_fig_h = 0.7
basis_fig_h, Pz_fig_h = p_dict[str(h_fig_h)][idx_fig_h]

# Data figure f
w_fig_f = 600
idx_fig_f = min(range(len(w_list)), key=lambda x: abs(w_list[x] - w_fig_f))
h_fig_f = 0.5
basis_fig_f, Pz_fig_f = p_dict[str(h_fig_f)][idx_fig_f]

fig, ax = plt.subplots(1, 2)

basis_fig_h.plot(np.abs(Pz_fig_h), ax=ax[0], aspect="equal")
ax[0].set_title(
    f"Poynting vector $S_z$\nfor h = {h_fig_h}μm & w = {w_fig_h}nm\n(Reproduction of Fig.1.h)"
)
ax[0].set_xlim(-w_fig_h * 1e-3 / 2 - 0.1, w_fig_h * 1e-3 / 2 + 0.1)
ax[0].set_ylim(capital_h - 0.1, capital_h + h_fig_h + 0.1)
# Turn off the axis wrt to the article figure
ax[0].axis("off")
# Add the contour
for subdomain in basis_fig_h.mesh.subdomains.keys() - {"gmsh:bounding_entities"}:
    basis_fig_h.mesh.restrict(subdomain).draw(
        ax=ax[0], boundaries_only=True, color="k", linewidth=1.0
    )

basis_fig_f.plot(np.abs(Pz_fig_f), ax=ax[1], aspect="equal")
ax[1].set_title(
    f"Poynting vector $S_z$\nfor h = {h_fig_f}μm & w = {w_fig_f}nm\n(Reproduction of Fig.1.f)"
)
ax[1].set_xlim(-w_list[idx_fig_f] * 1e-3 / 2 - 0.1, w_list[idx_fig_f] * 1e-3 / 2 + 0.1)
ax[1].set_ylim(capital_h - 0.1, capital_h + h_fig_f + 0.1)
# Turn off the axis wrt to the article figure
ax[1].axis("off")
# Add the contour
for subdomain in basis_fig_f.mesh.subdomains.keys() - {"gmsh:bounding_entities"}:
    basis_fig_f.mesh.restrict(subdomain).draw(
        ax=ax[1], boundaries_only=True, color="k", linewidth=1.0
    )

fig.tight_layout()
plt.show()
../../_images/3a99a4c4866a456b4976a4f1048db633b1fc08c202d13f051f89772ef1eeed8a.png

Bibliography#

[1]

Ivan D. Rukhlenko, Malin Premaratne, and Govind P. Agrawal. Effective mode area and its optimization in silicon-nanocrystal waveguides. Optics Letters, 37(12):2295, June 2012. URL: http://dx.doi.org/10.1364/OL.37.002295, doi:10.1364/ol.37.002295.