femwell
femwell copied to clipboard
Calculation of effective area for Spontaneous Four-wave Mixing (SFWM)
Hello,
Just as @elizaleung830 I'm trying to calculate the effective area, but for Spontaneous Four-wave Mixing (SFWM), which is a third-order non-linear effect and the Aeff needs to take into account the interaction between the electric field distributions of three modes (we consider a degenerate pump, which means there are two identical pumps).
The goal is to replicate the results found in the Ansys Lumerical tutorial on Spontaneous Four-Wave Mixing (SFWM) in a Microring Resonator Photon Source, as documented at https://optics.ansys.com/hc/en-us/articles/15100783091731, following the methodology outlined in the publication https://doi.org/10.1070/QEL16511:
A_{\mathrm{eff}}=\frac{1}{\iint \mathrm{d} x \mathrm{~d} y u_p(x, y) u_p(x, y) u_s^*(x, y) u_i^*(x, y)}
Where $u_p$, $u_s$, $u_i$ are the mode function describing the transverse spatial distribution of the field and normalised $\int|u(x, y)|^{2} \mathrm{~d} x \mathrm{~d} y=1$.
So, for a mode in femwell, does mode.E correspond to the unnormalized electric field distribution? I added normalization in femwell.maxwell.waveguide and calculated the method for Aeff, but I obtained a rather strange result.
Function added (for testing I unfreeze the class to add a E0 for the normalized electric field):
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 overlap integral result for the SFWM process.
"""
def normalize_mode(mode):
@Functional
def E2(w):
return (w["E"][0] * np.conj(w["E"][0]))
E_xy, _ = mode.basis.interpolate(mode.E) #xy, z
E_squared_integral = E2.assemble(mode.basis, E=E_xy)
normalization_factor = 1 / np.sqrt(E_squared_integral)
mode.E0 = mode.E * normalization_factor
normalize_mode(mode_p)
normalize_mode(mode_s)
normalize_mode(mode_i)
#Normalization needed for uv(x,y)!
E_p, _ = mode_p.basis.interpolate(mode_p.E0) #xy, z
E_s, _ = mode_s.basis.interpolate(mode_s.E0)
E_i, _ = mode_i.basis.interpolate(mode_i.E0)
@Functional(dtype=np.complex64)
def sfwm_overlap(w):
overlap_SFWM = w["E_p"][0] * w["E_p"][0] * np.conj(w["E_s"][0]) * np.conj(w["E_i"][0])
return (overlap_SFWM)
# Assemble the integral over the basis to compute the overlap
overlap_result = sfwm_overlap.assemble(
basis,
E_p=E_p,
E_s=E_s,
E_i=E_i,
)
return 1/overlap_result
And for the main code:
from collections import OrderedDict
import matplotlib.pyplot as plt
import numpy as np
from scipy.constants import speed_of_light
from shapely.geometry import box
from shapely.ops import clip_by_rect
from skfem import Basis, ElementTriP0
from skfem.io.meshio import from_meshio
from tqdm import tqdm
from femwell.maxwell.waveguide import compute_modes
from femwell.maxwell.waveguide import calculate_sfwm_Aeff
from femwell.mesh import mesh_from_OrderedDict
from scipy.constants import c, epsilon_0 ,pi
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(dtype=complex)
epsilon_s = basis0.zeros(dtype=complex)
epsilon_i = basis0.zeros(dtype=complex)
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)
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 = calculate_sfwm_Aeff(basis0, mode_p, mode_s, mode_i)
print(A_eff)
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 * pi * c / lambda_p0_m
gamma = (3 * chi_3 * omega_p0) / (4 * epsilon_0 * c**2 * n_p0**2 * A_eff_m2)
print("gamma:",gamma)
After confirming that the electric field distribution was normalized, I obtained a complex number with a negative real part for
{\iint \mathrm{d} x \mathrm{~d} y u_p(x, y) u_p(x, y) u_s^*(x, y) u_i^*(x, y)}
?According to the comparison with Lumerical's results, the Aeff should be on the order of $0.469 \ \text{um}^2$.