incorrect neff in long wavelength

[elizaleung830]

Hi,

I was trying to recreate figure 3b(W = 6 um, H = 0.8 um) from this paper. Strangely, GVD is doest not match when wavelength is greater than 3100nm. This happen when i am trying to recreate other curves in figure 3b as well. I wonder if i missed any parameters for the solver. Any help would be appreciate!

Here is my code:

from femwell.maxwell.waveguide import compute_modes
from skfem import Basis, ElementTriP0
from tqdm import tqdm
from femwell.visualization import plot_domains
from skfem.io import from_meshio
import shapely
import matplotlib.pyplot as plt
from refractive_index import n_MgF2, n_Si3N4, n_Air
from collections import OrderedDict
from numpy.polynomial import Polynomial
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import UnivariateSpline
import pandas as pd
from scipy.constants import speed_of_light
from femwell.mesh import mesh_from_OrderedDict

# waveguide parameters
width = 6  # um
height = 0.8  # um

n2 = 2.5e-19  # m^2/W n2 is the nonlinear refractive index at the center
Alpha = 0.7  # loss (dB/cm)

wavelength_range = [310, 5100]
wavelegnth_step = 70  #  steps

n_core = n_Si3N4
n_lower_cladding = lambda w: n_MgF2(w, ray="o") 
n_air = n_Air

# Construct waveguide geometry
core = shapely.geometry.box(-width / 2, 0, +width / 2, height)
lower_cladding = shapely.geometry.box(-6, -6, 6, 0)
air = shapely.geometry.box(-6, 0, 6, 6)
polygons = OrderedDict(
    core=core,
    lower_cladding=lower_cladding,
    air= air
)

# Define material property and resolution of waveguide
resolutions = dict(core={"resolution": 0.04, "distance": 0.1},
                   lower_cladding={"resolution": 0.15, "distance": 0.2},
                   air={"resolution": 0.2, "distance": 0.2})

n_dict = {"core": n_core, "lower_cladding": n_lower_cladding, "air": n_air}

# Calculate dispersion and gamma
mesh = from_meshio(mesh_from_OrderedDict(polygons, resolutions))
mesh.draw().show()
plot_domains(mesh)
plt.show()

print("start")
# Calculate dispersion and gamma
basis0 = Basis(mesh, ElementTriP0())
epsilon = basis0.zeros()
wavelength_list = np.linspace(wavelength_range[0], wavelength_range[1], wavelegnth_step)
neff_list = []
aeff_list = []
for wavelength in tqdm(wavelength_list):
    wavelength = wavelength * 1e-3
    for subdomain, n in n_dict.items():
        epsilon[basis0.get_dofs(elements=subdomain)] = n(wavelength) ** 2
    modes = compute_modes(basis0, epsilon, wavelength=wavelength, num_modes=3, order=1)
    modes_sorted = modes.sorted(key=lambda mode: -np.real(mode.n_eff))
    mode = modes_sorted[0]
    #mode.show(mode.E.real, direction ="x")
    neff_list.append(np.real(mode.n_eff))
    aeff_list.append(mode.calculate_effective_area())

neff_list = np.array(neff_list)
aeff_list = np.array(aeff_list)
wls = np.array(wavelength_list)

#--------calculate GVD-----------------
fig, axs = plt.subplots(3, 1, figsize=(9, 20))

y_spl = UnivariateSpline(wls, neff_list, s=0, k=3)
x_range = np.linspace(wls[0], wls[-1], 1000)
print(y_spl(1000))

axs[0].set_xlabel("Wavelength / µm")
axs[0].set_ylabel("neff")
axs[0].set_title(" neff vs wavelength fit")
axs[0].semilogy(x_range, y_spl(x_range))
axs[0].semilogy(wls, neff_list, 'ro', label='data')
axs[0].legend()

y_spl_2d = y_spl.derivative(n=2)
axs[1].set_xlabel("Wavelength / µm")
axs[1].set_ylabel("neff''")
axs[1].set_title("second derivative fit")
axs[1].plot(x_range, y_spl_2d(x_range))

ref_gvd = pd.read_csv("../reference_data/GVD.csv", dtype=np.float64)
ref_gvd_x, ref_gvd_y = np.split(ref_gvd.values, 2, axis=1)
axs[2].plot(ref_gvd_x, ref_gvd_y, c="green", label="paper")

GVD = (-wls / (2.99792e-7) * y_spl_2d(wls))
axs[2].plot(wls, GVD, label="calculated", c="red")

axs[2].set_ylabel("GVD")
axs[2].set_ylim(-1500, 250)
axs[2].set_xlim(500, 6000)
axs[2].set_title("GVD paramter")
axs[2].legend()

plt.tight_layout()
plt.show()

[HelgeGehring]

I'd guess you just need to make the cladding/air larger - at >=3µm I cannot really imagine that the field is already decayed at the boundaries as those are just a wavelength away. (The wavelengths in the plots seem to be off by 1000 :D )

Would be amazing if you could start another example for the docs :) (Or even better just replace the current dispersion example as we don´t have data to match it to) (and add yourself to the contributors list https://github.com/HelgeGehring/femwell/blob/main/README.md :smiley:)

[elizaleung830]

Thanks and I more than happy to add this to the example! However, seems like the size of cladding and air are not the main problem,I have adjusted the mesh to this:

width = 6  # um
height = 0.8  # um

n2 = 2.5e-19  # m^2/W n2 is the nonlinear refractive index at the center
Alpha = 0.7  # loss (dB/cm)

wavelength_range = [310, 5100]
wavelegnth_step = 70  #  steps

n_core = n_Si3N4
n_lower_cladding = lambda w: n_MgF2(w, ray="o") # TODO: check effect of this line
n_air = n_Air

# Construct waveguide geometry
core = shapely.geometry.box(-width / 2, 0, +width / 2, height)
lower_cladding = shapely.geometry.box(-18, -18, 18, 0)
air = shapely.geometry.box(-18, 0, 18, 18)
polygons = OrderedDict(
    core=core,
    lower_cladding=lower_cladding,
    air= air
)

# Define material property and resolution of waveguide
resolutions = dict(core={"resolution": 0.04, "distance": 0.1},
                   lower_cladding={"resolution": 0.3, "distance": 0.2},
                   air={"resolution": 0.3, "distance": 0.2})

but the resulting GVD curve remains the same as before. I have also tried with a higher resolution mesh, still run into the same problems.

[HelgeGehring]

Yay, sounds great 🥳

Hmm, okay, what else could remain?

Could you post a picture of a mode at ~5um, I'm wondering if it looks still like a nice mode or if it touches some of the boundaries 🤔 (even through I wouldn't expect that as changing the size of the domain didn't change anything)

What are you using for the refractive indices? Are you sure that those are still correct at those wavelengths? Might be that you have a fit which is just valid for a certain range? (But your comment already says Todo check that :D )

[elizaleung830]

Thanks for the quick response! This is the picture of mode in current mesh at 5000nm, but seems like it hasn't touch boundary of the mesh yet? In regard for the refractive indices, i have used the same sellmeier equation as mentioned in the paper. I have also tried using extraordinary and ordinary refractive index for cladding, but ordiary seems to give a better fit for GVD curve.

[HelgeGehring]

Hmm, okay, then the mesh should be fine :) Are you sure your sellmeier equations are correct? (without them I cannot run the example) https://refractiveindex.info/?shelf=main&book=Si3N4&page=Luke https://refractiveindex.info/?shelf=main&book=MgF2&page=Dodge-o these should be the ones, right? they even have an option to download the python code :)

[elizaleung830]

The sellmier equation of silcon nitride is same as what i use, but for the magnesium fluroide i use the following function which is dervied from a textbook :

def n_MgF2(wavelength, ray="o"):
    """
    valid for  (0.2–7)um
    :param wavelength: in um
    :return: linear refractive index of MgF2
    """
    if ray == "o":
        return math.sqrt(0.48755108*wavelength**2/(wavelength**2 - 0.04338408**2) + 0.39875031 * wavelength**2 / (wavelength**2 - 0.09461442**2) + 2.3120353 * wavelength**2 /(wavelength**2 - 23.793604**2) + 1)
    elif ray == "e":
        return math.sqrt(1+ 0.41344023*wavelength**2 /(wavelength**2 - 0.03684262**2) + 0.50497499*wavelength**2/(wavelength**2-0.09076162**2)+2.4904862*wavelength**2/(wavelength**2-23.771995**2))

Many thanks!

[HelgeGehring]

Compare the numbers, it's the same sellmeier equation ;)

I did a very quick check with another tool and got exactly the same results femwell produced.

So maybe we're using wrong sellmeier equations? Or we understood something wrong about their geometry (maybe a limited MgF2 thickness?) (Maybe there's also an error on their side?)

Would you have an idea what else could be different from the paper?

[elizaleung830]

The sellmeier equation should be correct ,as they are the one referenced by the paper. They didnt mention about a limited MgF2 thickness, but from the mode profile (Figure 2b): They are at least using thickness >= 3um for $\lambda$ = 4um, and the mode seems to just touched the boundary at this point, so i guess it shouldn't be the main reason GVD differ before $\lambda$ = 4um?

One thing could be different is that they are using a mesh with 500 000 first order triangular element, but our mesh only has about 40 000 elements. However, It causes memeory error for me when the mesh resolution is too high, so i can't verify this.

Anyway, I have also reproduced GVD curve from another paper using femwell recently, just need to add some documentation before i add it to the example. :)

[HelgeGehring]

Thanks for the reproduced GVD curve from the other paper, that looks really good :)

What do we do with this one? As I reproduced femwell results with another tool, it seems to me that the numerics is not the problem.

I'm wondering if there's a problem on their side? Maybe they had a too small domain?

I don't have other ideas what could be wrong :/ Should we close this issue or go on hunting for possible things we have different from them? @elizaleung830 what do you think?

[elizaleung830]

It might be a problem from their side or maybe some assumptions are not specified in their paper. But since it's not a problem from Femwell, I am happy to close this issue. Thanks a lot for helping! :)