SFWM Aeff calculation

docs/_toc.yml

@@ -33,6 +33,7 @@ parts:
  - file: photonics/examples/crosstalk.py
  - file: photonics/examples/propagation_loss.py
  - file: photonics/examples/effective_area.py
+ - file: photonics/examples/SFWM_effective_area.py
  - file: photonics/examples/metal_heater_phase_shifter.py
  - file: photonics/examples/metal_heater_phase_shifter_transient.py
  - file: photonics/examples/si_heater_phase_shifter.py

docs/photonics/examples/SFWM_effective_area.py

@@ -0,0 +1,233 @@
+# ---
+# jupyter:
+# jupytext:
+# formats: py:percent,md:myst
+# text_representation:
+# extension: .py
+# format_name: percent
+# format_version: '1.3'
+# jupytext_version: 1.15.0
+# kernelspec:
+# display_name: Python 3
+# name: python3
+# ---
+
+# %% [markdown]
+# # Optimisation of spontaneous four-wave mixing in a ring microcavity
+
+# Here we reproduce {cite}`Chuprina2017`
+
+# %% tags=["remove-stderr", "hide-input", "thebe-init"]
+from collections import OrderedDict
+
+import numpy as np
+from scipy.constants import c, epsilon_0
+from shapely.geometry import box
+from shapely.ops import clip_by_rect
+from skfem import Basis, ElementTriP0, Functional
+from skfem.helpers import dot
+from skfem.io.meshio import from_meshio
+
+from femwell.maxwell.waveguide import compute_modes
+from femwell.mesh import mesh_from_OrderedDict
+
+# %% [markdown]
+#
+# $$
+# \int\int\left|u(x,y)\right|^{2}\mathrm{d}x\mathrm{d}y=1
+# $$
+#
+# $$
+# A_{\mathrm{eff}}
+# =
+# \frac{1}{
+# \displaystyle
+# \int\int
+# \mathrm{d}x \mathrm{d}y
+# u_{\mathrm{p}}(x,y)
+# u_{\mathrm{p}}(x,y)
+# u_{\mathrm{s}}^{*}(x,y)
+# u_{\mathrm{i}}^{*}(x,y)
+# }
+# $$
+
+
+# %%
+def calculate_sfwm_Aeff(basis: Basis, mode_p, mode_s, mode_i) -> np.complex64:
+ """
+ Calculates the overlap integral for SFWM process by considering the interactions
+ between pump, signal, and idler modes in the xy plane.
+
+ Args:
+ basis (Basis): Common basis used for all modes in the same geometric structure.
+ mode_p, mode_s, mode_i: Mode instances for pump, signal, and idler, respectively.
+
+ Returns:
+ np.complex64: The Aeff result for the SFWM process(1/overlap integral).
+ """
+
+ def normalization_factor_mode(mode):
+ @Functional
+ def E2(w):
+ return dot(w["E"][0], np.conj(w["E"][0]))
+
+ E = mode.basis.interpolate(mode.E) # [0]=xy [1]=z
+ E_squared_integral = E2.assemble(mode.basis, E=E)
+ normalization_factor = 1 / np.sqrt(E_squared_integral)
+ return normalization_factor # Return the normalization factor instead of modifying the mode
+
+ # Apply normalization factors to the electric fields for the overlap calculation
+ @Functional(dtype=np.complex64)
+ def sfwm_overlap(w):
+
+ E_p_xy = w["E_p"][0] # shape: (2, x, 3)
+ E_s_xy = w["E_s"][0] # shape: (2, x, 3)
+ E_i_xy = w["E_i"][0] # shape: (2, x, 3)
+
+ overlap_Ex = E_p_xy[0, :, 0] * E_p_xy[0, :, 0] * np.conj(E_s_xy[0, :, 0]) * np.conj(
+ E_i_xy[0, :, 0]
+ ) + E_p_xy[1, :, 0] * E_p_xy[1, :, 0] * np.conj(E_s_xy[1, :, 0]) * np.conj(E_i_xy[1, :, 0])
+ overlap_Ey = E_p_xy[0, :, 1] * E_p_xy[0, :, 1] * np.conj(E_s_xy[0, :, 1]) * np.conj(
+ E_i_xy[0, :, 1]
+ ) + E_p_xy[1, :, 1] * E_p_xy[1, :, 1] * np.conj(E_s_xy[1, :, 1]) * np.conj(E_i_xy[1, :, 1])
+ overlap_Ez = E_p_xy[0, :, 2] * E_p_xy[0, :, 2] * np.conj(E_s_xy[0, :, 2]) * np.conj(
+ E_i_xy[0, :, 2]
+ ) + E_p_xy[1, :, 2] * E_p_xy[1, :, 2] * np.conj(E_s_xy[1, :, 2]) * np.conj(E_i_xy[1, :, 2])
+
+ return np.array([overlap_Ex, overlap_Ey, overlap_Ez]).T
+
+ # return dot(w["E_p"][0], w["E_p"][0]) * dot(np.conj(w["E_s"][0]), np.conj(w["E_i"][0]))#?
+ # return dot(w["E_p"][0], np.conj(w["E_s"][0])) * dot(w["E_p"][0], np.conj(w["E_i"][0]))#??
+
+ overlap_result = sfwm_overlap.assemble(
+ basis,
+ E_p=mode_p.basis.interpolate(mode_p.E * normalization_factor_mode(mode_p)),
+ E_s=mode_s.basis.interpolate(mode_s.E * normalization_factor_mode(mode_s)),
+ E_i=mode_i.basis.interpolate(mode_i.E * normalization_factor_mode(mode_i)),
+ )
+
+ return 1 / overlap_result
+
+
+# %% [markdown]
+#
+# Dispersion relations of materials
+# Silicon nitride {cite}`Luke2015`
+
+
+# %%
+def n_X(wavelength):
+ x = wavelength
+ return (
+ 1
+ + 2.19244563 / (1 - (0.20746607 / x) ** 2)
+ + 0.13435116 / (1 - (0.3985835 / x) ** 2)
+ + 2.20997784 / (1 - (0.20747044 / x) ** 2)
+ ) ** 0.5
+
+
+# %% [markdown]
+#
+# Box {cite}`Malitson1965`
+
+
+# %%
+def n_silicon_dioxide(wavelength):
+ x = wavelength
+ return (
+ 1
+ + 0.6961663 / (1 - (0.0684043 / x) ** 2)
+ + 0.4079426 / (1 - (0.1162414 / x) ** 2)
+ + 0.8974794 / (1 - (9.896161 / x) ** 2)
+ ) ** 0.5
+
+
+Clad = 1
+
+# %% [markdown]
+#
+# Waveguide dimensions 1050x500nm
+
+# %%
+
+core = box(0, 0, 1.05, 0.5)
+# core = box(0, 0, .5, 0.39) # 500x390nm
+polygons = OrderedDict(
+ core=core,
+ box=clip_by_rect(core.buffer(1.5, resolution=4), -np.inf, -np.inf, np.inf, 0),
+ clad=clip_by_rect(core.buffer(1.5, resolution=4), -np.inf, 0, np.inf, np.inf),
+)
+
+resolutions = {"core": {"resolution": 0.025, "distance": 2.0}}
+
+mesh = from_meshio(mesh_from_OrderedDict(polygons, resolutions, default_resolution_max=0.2))
+
+num_modes = 2
+
+# For SFWM we have energy conservation and momemtum(k) conservation for 2pumps and signal+idler
+# lambda_p0 = 1.4
+# lambda_s0 = 1.097
+# lambda_i0 = 1.686
+
+lambda_p0 = 1.55
+lambda_s0 = 1.55
+lambda_i0 = 1.55
+
+basis0 = Basis(mesh, ElementTriP0())
+
+epsilon_p = basis0.zeros()
+epsilon_s = basis0.zeros()
+epsilon_i = basis0.zeros()
+
+
+for wavelength, epsilon in zip(
+ [lambda_p0, lambda_s0, lambda_i0], [epsilon_p, epsilon_s, epsilon_i]
+):
+ for subdomain, n_func in {
+ "core": n_X,
+ "box": n_silicon_dioxide,
+ "clad": lambda _: Clad,
+ }.items():
+ n = n_func(wavelength)
+ epsilon[basis0.get_dofs(elements=subdomain)] = n**2
+
+
+modes_p = compute_modes(basis0, epsilon_p, wavelength=lambda_p0, num_modes=num_modes, order=1)
+modes_s = compute_modes(basis0, epsilon_s, wavelength=lambda_s0, num_modes=num_modes, order=1)
+modes_i = compute_modes(basis0, epsilon_i, wavelength=lambda_i0, num_modes=num_modes, order=1)
+
+# modes_p[0].show(modes_p[0].E.real)
+
+mode_p = max(modes_p, key=lambda mode: mode.te_fraction)
+mode_s = max(modes_s, key=lambda mode: mode.te_fraction)
+mode_i = max(modes_i, key=lambda mode: mode.te_fraction)
+
+A_eff = np.real(calculate_sfwm_Aeff(basis0, mode_p, mode_s, mode_i))
+print("Aeff in um2:", A_eff)
+
+# Calculation for non-linear coef
+chi_3 = 5e-21 # m^2/V^2 #7e-20?
+lambda_p0_m = lambda_p0 * 1e-6 # to m
+n_p0 = np.real(mode_p.n_eff)
+A_eff_m2 = A_eff * 1e-12 # to m^2
+
+omega_p0 = 2 * np.pi * c / lambda_p0_m
+
+# %% [markdown]
+# $$
+# \gamma=\frac{3\chi^{(3)}\omega_{\mathrm{p}}}{4\varepsilon_{0}c^{2}n_{\mathrm{p0}}^{2}A_{\mathrm{eff}}}
+# $$
+
+# %%
+
+gamma = (3 * chi_3 * omega_p0) / (4 * epsilon_0 * c**2 * n_p0**2 * A_eff_m2)
+
+print("gamma:", gamma)
+
+# %% [markdown]
+# ## Bibliography
+#
+# ```{bibliography}
+# :style: unsrt
+# :filter: docname in docnames
+# ```

docs/references.bib

@@ -342,3 +342,17 @@ @article{Jansen1978
  month = feb,
  pages = {75--82}
 }
+@article{Chuprina2017,
+ title = {Optimisation of spontaneous four-wave mixing in a ring microcavity},
+ volume = {47},
+ issn = {1468-4799},
+ url = {http://dx.doi.org/10.1070/QEL16511},
+ doi = {10.1070/qel16511},
+ number = {10},
+ journal = {Quantum Electronics},
+ publisher = {IOP Publishing},
+ author = {Chuprina, I N and An, P P and Zubkova, E G and Kovalyuk, V V and Kalachev, A A and Gol'tsman, G N},
+ year = {2017},
+ month = oct,
+ pages = {887-891}
+}

femwell/tests/test_Aeff_mode.py

@@ -0,0 +1,99 @@
+from collections import OrderedDict
+import matplotlib.pyplot as plt
+import numpy as np
+from shapely.ops import clip_by_rect
+from skfem import Basis, ElementTriP0
+from skfem.io.meshio import from_meshio
+from shapely.geometry import box
+from femwell.maxwell.waveguide import compute_modes
+from femwell.mesh import mesh_from_OrderedDict
+
+
+def n_X(wavelength):
+ x = wavelength
+ return (1 
+ + 2.19244563 / (1 - (0.20746607 / x) ** 2)
+ + 0.13435116 / (1 - (0.3985835 / x) ** 2)
+ + 2.20997784 / (1 - (0.20747044 / x) ** 2)
+ ) ** 0.5
+
+# Box
+def n_silicon_dioxide(wavelength):
+ x = wavelength
+ return (
+ 1
+ + 0.6961663 / (1 - (0.0684043 / x) ** 2)
+ + 0.4079426 / (1 - (0.1162414 / x) ** 2)
+ + 0.8974794 / (1 - (9.896161 / x) ** 2)
+ ) ** 0.5
+
+Clad=1
+
+core = box(0, 0, 0.5, 0.39)
+polygons = OrderedDict(
+ core=core,
+ box=clip_by_rect(core.buffer(1.5, resolution=4), -np.inf, -np.inf, np.inf, 0),
+ clad=clip_by_rect(core.buffer(1.5, resolution=4), -np.inf, 0, np.inf, np.inf),
+)
+
+resolutions = {"core": {"resolution": 0.025, "distance": 2.}}
+
+mesh = from_meshio(mesh_from_OrderedDict(polygons, resolutions, default_resolution_max=0.6))
+
+num_modes = 2
+
+lambda_p0 = 1.4
+lambda_s0 = 1.097
+lambda_i0 = 1.686
+
+basis0 = Basis(mesh, ElementTriP0())
+
+epsilon_p = basis0.zeros()
+epsilon_s = basis0.zeros()
+epsilon_i = basis0.zeros()
+
+for wavelength, epsilon in zip([lambda_p0, lambda_s0, lambda_i0], [epsilon_p, epsilon_s, epsilon_i]):
+ for subdomain, n_func in {
+ "core": n_X,
+ "box": n_silicon_dioxide,
+ "clad": lambda _: Clad,
+ }.items():
+ n = n_func(wavelength)
+ epsilon[basis0.get_dofs(elements=subdomain)] = n**2
+modes = compute_modes(basis0, epsilon, wavelength=wavelength, num_modes=2, order=2)
+
+
+modes_p = compute_modes(basis0, epsilon_p, wavelength=lambda_p0, num_modes=num_modes, order=1)
+modes_s = compute_modes(basis0, epsilon_s, wavelength=lambda_s0, num_modes=num_modes, order=1)
+modes_i = compute_modes(basis0, epsilon_i, wavelength=lambda_i0, num_modes=num_modes, order=1)
+
+mode_p = max(modes_p, key=lambda mode: mode.te_fraction)
+mode_s = max(modes_s, key=lambda mode: mode.te_fraction)
+mode_i = max(modes_i, key=lambda mode: mode.te_fraction)
+
+print(f"Effective refractive index for p mode: {np.real(mode_p.n_eff):.4f}")
+mode_p.show(mode_p.E.real, colorbar=True, direction="x")
+mode_p.show(mode_p.E.imag, colorbar=True, direction="x")
+fig, ax = plt.subplots()
+mode_p.plot_intensity(ax=ax)
+plt.title("Normalized Intensity for p mode")
+plt.tight_layout()
+plt.show()
+
+print(f"Effective refractive index for s mode: {np.real(mode_s.n_eff):.4f}")
+mode_s.show(mode_s.E.real, colorbar=True, direction="x")
+mode_s.show(mode_s.E.imag, colorbar=True, direction="x")
+fig, ax = plt.subplots()
+mode_s.plot_intensity(ax=ax)
+plt.title("Normalized Intensity for s mode")
+plt.tight_layout()
+plt.show()
+
+print(f"Effective refractive index for i mode: {np.real(mode_i.n_eff):.4f}")
+mode_i.show(mode_i.E.real, colorbar=True, direction="x")
+mode_i.show(mode_i.E.imag, colorbar=True, direction="x")
+fig, ax = plt.subplots()
+mode_i.plot_intensity(ax=ax)
+plt.title("Normalized Intensity for i mode")
+plt.tight_layout()
+plt.show()