Coupled mode theory

Contents

Coupled mode theory#

Caution

This example is under construction, results are not yet all correct

https://www.fiberoptics4sale.com/blogs/wave-optics/coupled-mode-theory https://www.fiberoptics4sale.com/blogs/wave-optics/two-mode-coupling

Hide code cell source
from collections import OrderedDict
from itertools import chain

import matplotlib.pyplot as plt
import numpy as np
from scipy.constants import epsilon_0, speed_of_light
from scipy.integrate import solve_ivp
from shapely.geometry import Polygon
from skfem import Basis, ElementTriP0, Mesh
from skfem.io import from_meshio

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

Let’s set up the geometry! It’s the cross-section of two parallel waveguides with different widths:

Hide code cell source
w_sim = 4
h_clad = 1
h_box = 1
w_core_1 = 0.45
w_core_2 = 0.46
gap = 0.4
h_core = 0.22
offset_heater = 2.2
h_heater = 0.14
w_heater = 2

wavelength = 1.55
k0 = 2 * np.pi / wavelength

polygons = OrderedDict(
    core_1=Polygon(
        [
            (-w_core_1 - gap / 2, 0),
            (-w_core_1 - gap / 2, h_core),
            (-gap / 2, h_core),
            (-gap / 2, 0),
        ]
    ),
    core_2=Polygon(
        [
            (w_core_2 + gap / 2, 0),
            (w_core_2 + gap / 2, h_core),
            (gap / 2, h_core),
            (gap / 2, 0),
        ]
    ),
    clad=Polygon(
        [
            (-w_sim / 2, 0),
            (-w_sim / 2, h_clad),
            (w_sim / 2, h_clad),
            (w_sim / 2, 0),
        ]
    ),
    box=Polygon(
        [
            (-w_sim / 2, 0),
            (-w_sim / 2, -h_box),
            (w_sim / 2, -h_box),
            (w_sim / 2, 0),
        ]
    ),
)

resolutions = dict(
    core_1={"resolution": 0.03, "distance": 1},
    core_2={"resolution": 0.03, "distance": 1},
)

mesh = from_meshio(
    mesh_from_OrderedDict(polygons, resolutions, filename="mesh.msh", default_resolution_max=0.2)
)
mesh.draw().show()
../../_images/47675c78f96e45ffbfd6326e9e24f98e2102af19009098d50187c7656b9eefc5.png

First we plot the symmetric and the asymmetric modes of the geometry with both waveguides:

basis0 = Basis(mesh, ElementTriP0(), intorder=4)

epsilon = basis0.zeros() + 1.444**2
epsilon[basis0.get_dofs(elements=("core_1", "core_2"))] = 3.4777**2
# basis0.plot(epsilon, colorbar=True).show()
modes_both = compute_modes(basis0, epsilon, wavelength=wavelength, mu_r=1, num_modes=2)
modes_both[0].show(modes_both[0].E.real, direction="x")
modes_both[1].show(modes_both[1].E.real, direction="x")
print(
    "Refractive index of symmetric and assymetric mode:",
    modes_both[0].n_eff,
    ", ",
    modes_both[1].n_eff,
)
# https://www.fiberoptics4sale.com/blogs/wave-optics/directional-couplers
print(
    f"Maximum power transfer after {np.pi / (2 * np.pi / wavelength * np.real(modes_both[0].n_eff - modes_both[1].n_eff))} um prop length"
)
/tmp/ipykernel_8829/3693127342.py:7: DeprecationWarning: The behavior of passing an array directly to `show` is deprecated and will be removed in the future. Use `plot` instead.
  modes_both[0].show(modes_both[0].E.real, direction="x")
../../_images/df4428700feb3ad5f06839a590d3d89f457b1797d54493cfb381f70b59b0278b.png
/tmp/ipykernel_8829/3693127342.py:8: DeprecationWarning: The behavior of passing an array directly to `show` is deprecated and will be removed in the future. Use `plot` instead.
  modes_both[1].show(modes_both[1].E.real, direction="x")
../../_images/0d61c7e11a89b99dabf70db69d3b0807b4034f507ccde90e70fa01e7feb18e07.png
Refractive index of symmetric and assymetric mode: (2.3708828136066713+0j) ,  (2.349150308782909+0j)
Maximum power transfer after 35.66086865203918 um prop length

And then we plot the modes of each waveguide while setting the shape of the other one to oxide

epsilon = basis0.zeros() + 1.444**2
epsilon[basis0.get_dofs(elements="core_1")] = 3.4777**2
# basis0.plot(epsilon, colorbar=True).show()
modes_1 = compute_modes(basis0, epsilon, wavelength=wavelength, mu_r=1, num_modes=1)
print("Effective refractive index of the mode of the first waveguide", modes_1[0].n_eff)
modes_1[0].show(modes_1[0].E.real, direction="x")

epsilon_2 = basis0.zeros() + 1.444**2
epsilon_2[basis0.get_dofs(elements="core_2")] = 3.4777**2
# basis0.plot(epsilon_2, colorbar=True).show()
modes_2 = compute_modes(basis0, epsilon_2, wavelength=wavelength, mu_r=1, num_modes=1)
print("Effective refractive index of the mode of the second waveguide", modes_2[0].n_eff)
modes_2[0].show(modes_2[0].E.real, direction="x")
Effective refractive index of the mode of the first waveguide (2.349519763537976+0j)
/tmp/ipykernel_8829/185950769.py:6: DeprecationWarning: The behavior of passing an array directly to `show` is deprecated and will be removed in the future. Use `plot` instead.
  modes_1[0].show(modes_1[0].E.real, direction="x")
../../_images/cd1a57842ed59aedeee2fc50f49048cf4c1308f633e097767df308f0353a0abc.png
Effective refractive index of the mode of the second waveguide (2.370423954537899+0j)
/tmp/ipykernel_8829/185950769.py:13: DeprecationWarning: The behavior of passing an array directly to `show` is deprecated and will be removed in the future. Use `plot` instead.
  modes_2[0].show(modes_2[0].E.real, direction="x")
../../_images/b207cd4ffdcbf445d43d65626574907990e2d7e1540929a07074de8772ab1133.png
length = 200
ts = np.linspace(0, length, 1000)
epsilons = [epsilon, epsilon_2]

num_modes = len(modes_1) + len(modes_2)
overlap_integrals = np.zeros((num_modes, num_modes), dtype=complex)
for i, mode_i in enumerate(chain(modes_1, modes_2)):
    for j, mode_j in enumerate(chain(modes_1, modes_2)):
        overlap_integrals[i, j] = mode_i.calculate_overlap(mode_j)

print("overlap", overlap_integrals)
# plt.imshow(np.abs(overlap_integrals))
# plt.colorbar()
# plt.show()

coupling_coefficients = np.zeros((num_modes, num_modes), dtype=complex)
for i, mode_i in enumerate(chain(modes_1, modes_2)):
    for j, mode_j in enumerate(chain(modes_1, modes_2)):
        coupling_coefficients[i, j] = (
            k0
            * speed_of_light
            * epsilon_0
            * mode_i.calculate_coupling_coefficient(
                mode_j, epsilons[(j // len(modes_1) + 1) % 2] - 1.444**2
            )
            * 0.5
        )


print(coupling_coefficients)
# plt.imshow(np.abs(coupling_coefficients))
# plt.colorbar()
# plt.show()

kappas = np.array(
    [
        [
            (
                coupling_coefficients[i, j]
                - overlap_integrals[i, (i + 1) % 2]
                * coupling_coefficients[(i + 1) % 2, j]
                / overlap_integrals[(i + 1) % 2, (i + 1) % 2]
            )
            / (
                1
                - overlap_integrals[0, 1]
                * overlap_integrals[1, 0]
                / (overlap_integrals[0, 0] * overlap_integrals[1, 1])
            )
            for i in range(2)
        ]
        for j in range(2)
    ]
)
print(kappas)

delta = 0.5 * (
    np.real(modes_1[0].n_eff) * k0 + kappas[1, 1] - (np.real(modes_2[0].n_eff) * k0 + kappas[0, 0])
)
print(delta, np.real(modes_1[0].n_eff) * k0, kappas[1, 1])

beta_c = (kappas[0, 1] * kappas[1, 0] + delta**2) ** 0.5

print(np.pi / (2 * beta_c))

eta = np.abs(kappas[1, 0] ** 2 / beta_c**2) * np.sin(beta_c * 1e3)
print("eta", eta, np.abs(kappas[1, 0] ** 2 / beta_c**2))

plt.plot(ts, 1 - np.abs(kappas[1, 0] ** 2 / beta_c**2 * np.sin(beta_c * ts) ** 2))
plt.plot(ts, np.abs(kappas[1, 0] ** 2 / beta_c**2 * np.sin(beta_c * ts) ** 2))
plt.show()
overlap [[ 1.        +0.j -0.01393994+0.j]
 [-0.01393994+0.j  1.        +0.j]]
[[ 0.00072551+0.j -0.01161464+0.j]
 [-0.0127873 +0.j  0.00061889+0.j]]
[[ 0.00054736+0.j -0.01277967+0.j]
 [-0.01160827+0.j  0.00045707+0.j]]
(-0.042414467451351356+0j) 9.524172940122495 (0.0004570730574596199+0j)
(35.595846991119586+0j)
eta (0.01009055164071527+0j) 0.0691980549984048
../../_images/b910148691faf6d7137b66def18f0ea2ade38ff7b5fb619cf72daca28c6803e1.png

see http://home.iitj.ac.in/~k.r.hiremath/research/thesis.pdf , not yet finished

def fun(t, y):
    phase_matrix = [
        [
            np.exp(2j * np.pi / wavelength * (mode_i.n_eff - mode_j.n_eff) * t)
            for mode_j in chain(modes_1, modes_2)
        ]
        for mode_i in chain(modes_1, modes_2)
    ]
    matrix = (
        np.linalg.inv(overlap_integrals * phase_matrix)
        @ (coupling_coefficients * phase_matrix)
        * -1j
    )
    return (matrix @ y).ravel()


result = solve_ivp(fun, [0, length], np.array((1, 0), dtype=complex), t_eval=ts)
ys = result.y.T

plt.plot(ts, np.abs(np.array(ys)[:, 0]) ** 2, "r")
plt.plot(ts, 1 - np.abs(np.array(ys)[:, 0]) ** 2, "r")
# plt.plot(ts, np.array(ys).imag.reshape((-1,)+matrix.shape)@(1,0), 'g')
plt.show()
../../_images/eaee90cae2b31ebcdb68c97f343bcc0f3a431266778ecc1614b7fb7526e80f1c.png

two modes#

R = [np.abs(modes_1[0].calculate_overlap(mode_j) ** 2) for mode_j in modes_both]
print(R)
P = (
    R[0] ** 2
    + R[1] ** 2
    + 2
    * R[0]
    * R[1]
    * np.cos(2 * np.pi / wavelength * (modes_both[0].n_eff - modes_both[1].n_eff) * ts)
)

plt.plot(ts, P)
plt.show()
[0.021346039780330346, 0.9786991534283266]
/home/runner/miniconda3/lib/python3.12/site-packages/matplotlib/cbook.py:1709: ComplexWarning: Casting complex values to real discards the imaginary part
  return math.isfinite(val)
/home/runner/miniconda3/lib/python3.12/site-packages/matplotlib/cbook.py:1345: ComplexWarning: Casting complex values to real discards the imaginary part
  return np.asarray(x, float)
../../_images/aa66cda10973e2108a6cdc032a8793ccb57e7c000a3b2fb06a9c08e61bf26a6c.png
plt.plot(ts, P)


plt.plot(ts, 1 - np.abs(kappas[1, 0] ** 2 / beta_c**2) * np.sin(beta_c * ts) ** 2)

plt.plot(ts, np.abs(np.array(ys)[:, 0]) ** 2, "r")
plt.plot(ts, 1 - np.abs(np.array(ys)[:, 0]) ** 2, "r")
# plt.plot(ts, np.array(ys).imag.reshape((-1,)+matrix.shape)@(1,0), 'g')
plt.show()
../../_images/240a7bebad2c55eba59bf3b9b29c9faa27679f34fba66a25f1b94843af822945.png