Calculation of effective area for Spontaneous Four-wave Mixing (SFWM)
Opened this issue Β· 3 comments
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:
Where
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
?According to the comparison with Lumerical's results, the Aeff should be on the order of
Looks like a great start @Kunyuan-LI !
Three quick things on the first look:
-
I'd not unfreeze the
mode
(it's normalized to power=1), but include the factors for the normalization directly in the calculation, i.e. just divide the calculated number by the three normalization factors. There's so many ways to normalize a mode, it might be hard to store all of them π€
If we see in the end that we need them more often, we could add it as acached_property
to the mode object π€
what do you think? -
About the integrals:
@Functional
def E2(w):
return (w["E"][0] * np.conj(w["E"][0]))
I'm actually a bit confused that this doesn't fail π E[0] should be a 2d vector field. I'd guess what you want to do here is
@Functional
def E2(w):
return dot(w["E"][0],np.conj(w["E"][0]))
to get the scalar product of the E.
But this also confuses me a bit - the integral you posted in the end doesn't seem to have any vectors inside but uses scalars.
Do we just need there to use the dominant component? (i.e. ['E'][0][0 or 1]
)
Do you know how that integral is derived?
- As long as you use complex epsilons the E-fields' phases will also be complex with a random global phase. While this is nothing to worry about and it can just be ignored, there's no reason for this example to go to complex epsilons.
Just use
epsilon_p = basis0.zeros()
epsilon_s = basis0.zeros()
epsilon_i = basis0.zeros()
and you'll get real fields. While this can still give a random sign to the integral (as the sign of the eigenmode is not defined), we'll get something real :)
Could you make a PR out of this like @elizaleung830 ? That's easiest to discuss/test/improve and in the end we would also get a nice extra example :)
Thank you for your answer!
Indeed, the mistake about the normalization of E is quite obviousπ . However, maybe getting it right could be helpful for the subsequent calculations?π
In the numerical integration of finite elements, for normalization, what we actually consider is the sum of the
calculated separately on each element of the basis. Considering
needs to be modified as well because function.assemble?
π
And for the complex epsilons, yes I think you are right. As we just concern about the effective refractive index
For the PR sure I'll make it ASAP.π
Hello,
I've just made a PR, you can now critique my coding.π
And for the overlap in SFWM, in this case
So do you have a good idea to multiply them then intergrate them in the xy-plane?π€