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]
Show 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
Show code cell output
Effective refractive index: 1.6352+0.0000j
Effective mode area: 0.0866
Mode transversality: 0.8364573555269028
int(Sz) (1.0000000027154932+0j)
int_core(1) 0.024999999999999994
int_core(Sz) (0.3461699676246049+0j)
int(Sz^2) (6.849160748546705+0j)
Effective refractive index: 1.8334+0.0000j
Effective mode area: 0.0880
Mode transversality: 0.8313585190086243
int(Sz) (1.0000000015631885+0j)
int_core(1) 0.03499999999999999
int_core(Sz) (0.40701440350805296+0j)
int(Sz^2) (6.733755685964607+0j)
Effective refractive index: 1.9568+0.0000j
Effective mode area: 0.0976
Mode transversality: 0.8266226822935947
int(Sz) (1.0000000034430936+0j)
int_core(1) 0.04499999999999999
int_core(Sz) (0.43015535402434624+0j)
int(Sz^2) (6.073848322289527+0j)
Effective refractive index: 2.0367+0.0000j
Effective mode area: 0.1093
Mode transversality: 0.8228259941564221
int(Sz) (1.0000000020891473+0j)
int_core(1) 0.05499999999999999
int_core(Sz) (0.4395054900178267+0j)
int(Sz^2) (5.409137609393959+0j)
Effective refractive index: 2.0908+0.0000j
Effective mode area: 0.1218
Mode transversality: 0.8197180576035481
int(Sz) (1.000000001167669+0j)
int_core(1) 0.06499999999999999
int_core(Sz) (0.44315250888026103+0j)
int(Sz^2) (4.8325015913504155+0j)
Effective refractive index: 2.0639+0.0000j
Effective mode area: 0.1091
Mode transversality: 0.8165931065816691
int(Sz) (1.0000000039954564+0j)
int_core(1) 0.024999999999999994
int_core(Sz) (0.18288719132675213+0j)
int(Sz^2) (4.966105830433412+0j)
Effective refractive index: 2.3037+0.0000j
Effective mode area: 0.1344
Mode transversality: 0.8009123980406774
int(Sz) (1.0000000025400038+0j)
int_core(1) 0.03499999999999999
int_core(Sz) (0.18210406425445635+0j)
int(Sz^2) (4.583317896938594+0j)
Effective refractive index: 2.4405+0.0000j
Effective mode area: 0.1619
Mode transversality: 0.7942437495221183
int(Sz) (1.0000000020827027+0j)
int_core(1) 0.04499999999999999
int_core(Sz) (0.1791383172238547+0j)
int(Sz^2) (4.063743899371225+0j)
Effective refractive index: 2.5243+0.0000j
Effective mode area: 0.1892
Mode transversality: 0.791011695864681
int(Sz) (1.0000000015770518+0j)
int_core(1) 0.05499999999999999
int_core(Sz) (0.17684374595910685+0j)
int(Sz^2) (3.5904097028603448+0j)
Effective refractive index: 2.5788+0.0000j
Effective mode area: 0.2159
Mode transversality: 0.7891039189018594
int(Sz) (1.000000002643856+0j)
int_core(1) 0.06499999999999999
int_core(Sz) (0.17515582092533805+0j)
int(Sz^2) (3.1943473251511194+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)
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()
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.