SFWM Aeff calculation

[Kunyuan-LI]

Hello,

This is the PR for the issue #130.

Best,

[lucasgrjn]

Hey @Kunyuan-LI cool example!

I am not sure it is necessary to remove the frozen state of the dataclass: since you only use E0 in calculate_sfwm_Aeff(), maybe you could use it as a variable. What do you think?

Also, to make things more easy, you should update femwell to the last version :)

[Kunyuan-LI]

@lucasgrjn thank you!
@HelgeGehring said if it is more often used, we can add a cached property.
For the update I've alredy merge the main branche before the PR, is there still a update of femwell?

[lucasgrjn]

Yeah you are right, it definitely make sense to cache it if you may need it on anther part of the code :)

For the update, I checked your fork history and it looks like you had the last version. However you delete some pieces of the code such as:

• poynting()
• transversality()
• calculate_effective_area()

[HelgeGehring]

It seems you're removing the newest stuff in waveguide.py, so you probably missed the latest commits on main 🤔
Best would be if you could try to not modify waveguide.py and keep everything in the example.

If we find out over time we use that normalization more often, we can put something in waveguide.py, but for now I don't think it makes sense to add the differently normalized fields there.
If we use them more often, it would also be enough to just cache the norm of the fields (i.e. one number) as that allows to directly calculate the field normalized with that norm.

When I see calculations based on only E, I'm also wondering if there's some approximation in there. The current normalization is based on the poynting vector which should be the most accurate way to normalize the field.
Calculations based on a normalized E often make the assumption that the field is completely in one component and thus can lead to inaccurancies (I think often these approximations come from fiber optics where we don't have that sharp neff steps and the modes are almost pure TE/TM)
While these approximations can be good for certain geometries, there's often a way do the same calculation without approximations.

Could you update the PR to have everything in the example? That'll make debugging easier :)

[HelgeGehring]

Hmm, now you merged the most recent changes, but it seems you're still deleting stuff in waveguide.py 🤔
Could you try to move your changes to the example file? And not change waveguide.py?
Then we can start looking for what we need to change :)

[Kunyuan-LI]

Sure. I made a standalone Aeff.py example.
What do you prefer for the uv(x,y) without using mode.E?🤔

[HelgeGehring]

Ups, that's not what I meant!

I meant you can use waveguide.py, but you shouldn't modify it.

I.e. instead of defining mode.E0 you just calculate the integral needed for normalization and divide the values you get by it.

So that you have something like:

-Calculate normalization factor for each mode
-Calculate the integral (in which you just divide the modes directly by the normalization)

The idea is to have it built on waveguide.py, but decide afterwards if we want to move something there

[HelgeGehring]

Hey @Kunyuan-LI ,

I think I found it! 🥳
Probably I got confused as you defined E_xy, _ = basis.interpolate..., which would then require one less [0] in the integral 🤔
I changed it to stay consistent with the rest of the project.

Also it seemed to me that they use a different crosssection in the paper (?) I changed it to what the mention in the last paragraph

With the changes we now get for the effective area ~0.4 while the paper says 0.5...
I'd guess that'll probably require a little playing with the parameters?

Do they mention anywhere what exact wavelengths they use for that integrals?
Did they mention the sellmeier equations? Or is that based on the materials? (I wouldn't Aeff expect to be sensitive on the exact number)

We should also think what exacly we want to use for the calculation - it's now a dot-product between the E-fields (whereas it's a bit random chosen which ones we dot product in Aeff)
This should for modes which are closed to being pure TE / TM be fine, but for everything in-between it would result in something different than chosing one field...... But I'd guess for non-(close-to)-pure modes we'd anyway need to look deeper in how the integrals are calculated....

Could you have a look if you find the reason for the small derivation?
And could you also move it to docs/photonics/examples? that way we could have it in the docs and go on improving/completing it there :)
(See #120 for how to add something to the docs, @elizaleung830 did a great job there!)

[Kunyuan-LI]

Great! 🥳

In my understanding, E_xy, _ = basis.interpolate... and E_xy = basis.interpolate...[0] are equivalent, so there is no issue, after all our goal is the integration over the xy-plane, right? Or as usual I misubderstood something?😂

Regarding the "different crosssection in the paper (?)", I believe it should be that SFWM should be completed within the same crosssection because all the electric field distributions uv(x,y) are normalized. That is, only when using the same geometry, these normalized, physically meaningful 'electric field distributions' make sense in calculating the overlap integral(?). Otherwise, normalizing and then multiplying different crosssections seems make no sense here.🤔

After reviewing the literature at https://doi.org/10.1070/QEL16511, I believe in equations (3) and (4):

$$\begin{align*} H(t) = \frac{3 \varepsilon_{0}}{4} \int \mathrm{d} \boldsymbol{r} \chi^{(3)}(\boldsymbol{r}) &E_{\mathrm{p}}^{(+)}(\boldsymbol{r}, t) E_{\mathrm{p}}^{(+)}(\boldsymbol{r}, t) E_{\mathrm{s}}^{(-)}(\boldsymbol{r}, t) E_{\mathrm{i}}^{(-)}(\boldsymbol{r}, t) \exp \left(\mathrm{i} \Phi_{\mathrm{NL}}\right) +\text { H.c. } \tag{3} \end{align*}$$

Here the effective Hamiltonian (eq3, order equivalently modified) gives the content of the integration(p+p+s-i-), then

$$\begin{align*} & E_{v}^{(+)}(\boldsymbol{r}, t)=\mathrm{i} \sqrt{\frac{\hbar \omega_{v 0}}{2 \varepsilon_{0} n_{v 0} n_{\mathrm{g} v 0}}} \frac{u_{v}(x, y)}{\sqrt{L}} \times \int \frac{\mathrm{d} \omega}{\sqrt{2 \pi}} a_{v}(\omega) \exp \left(\mathrm{i} \beta_{v} z-\mathrm{i} \omega t\right), \quad E_{v}^{(-)}=\left[E_{v}^{(+)}\right]^{\dagger}\tag{4} \end{align*}$$

describes the relationship between $E_{v}^{(+)}(r,t)$ and $u_{v}(x,y)$. I believe $u_{v}(x,y)$ as a time-invariant variable is abstracted here as a highly abstract scalar, so $E_{p}^{(+)}(r,t)$ becomes $u_{p}(x,y)$ and $E_{s}^{(-)}(r,t)$, $E_{i}^{(-)}(r,t)$ become $u_{s}^*(x,y)$, $u_{i}^*(x,y)$ in the Aeff calculation. But here comes a point:

$$\begin{align*} E_{v}^{(-)}=\left[E_{v}^{(+)}\right]^{\dagger} \Rightarrow u_{v}^*(x,y)=\left[u_{v}(x,y)\right]^{\dagger} \end{align*}$$

only exists when $u_{v}(x,y)$ is a scalar. Because the dagger $\dagger$ symbolizes conjugate transpose in quantum optics, but for scalars, the transpose equals itself, thus it simplifies to conjugation. However, this premise relies on treating $u_{v}(x,y)$ as a scalar. Yet, for operators (matrices) and like our situation where $u_{v}(x,y)$ represents a matrix distribution, I think we need to remain the conjugate transpose. This was the point that had previously confused me.😂

So we need again to reevaluate the overlap function 🤔. However, for me the result already makes sense. I used a different geometry and material, but for the comparison between Lumerical and theoretical results, in our case, we have ~0.4 in theoretical (Femwell?) versus 0.469 in Lumerical, and in the tutorial, there's 0.5 in theoretical calculation versus 0.59 in Lumerical. In fact, I'm also not clear on how this difference (significant enough) has emerged, nor whether it will have a notable impact on subsequent phase matching experiments. It can only be said that it remains to be assessed.🤔

They did not mention the Sellmeier equation in the calculation of Aeff​, but we indeed used Sellmeier to obtain neff for TE and TM modes and plotted the contour of Δk=0 under the assumption of zero non-linear shifting, also realized by femwell😙

for the degenerate pump scenario, we can plot the p-si contour, which allows us to determine the specific pump wavelengths corresponding to pairs of si wavelengths that satisfy wave k conservation (of course, energy conservation is a prerequisite). The wavelengths in the example were also obtained by this method, and unlike the Aeff​ for a single mode, the Aeff for SFWM, considering the overlap of three modes, which all have corresponding wavelengths, is indeed sensitive to numerical values. To give an extreme example, suppose in the p-si contour, the wavelengths of p and s can be accommodated by the waveguide, but i escapes because its wavelength is too long for the waveguide to accommodate. Then, the integral of i's electric field in the waveguide is very small (although the integral over the 'entire' space is still 1, there's almost none in the core part), my conclusion is that this is based on both the nonlinear material and the three modes.

Regarding 'how the integrals are calculated', indeed, I also believe we need to reconsider how this overlap is calculated to accurately express it. And although for waveguides we generally only consider TE and TM modes, considering that in reality, there are TE-TM mode SFWM, anyway, we need to delve deeper into how the integrals are calculated....

Move to docs/photonics/examples will be completed instantaly. However, I feel like the holes we are digging are getting bigger and bigger(?)😂

[HelgeGehring]

Yes! both are equivalent. What confused me was that E_xy was put in assemble. This way we need in the integral less indices. As the variable was called E_p I assumed that we would need a [0] to get E_xy, which led to some confusion in our discussion :D

Yeah, all should be calculated with the same waveguide geometry. I didn't read the complete paper, I just looked in the end where they mentioned dimensions to calculate their effective area and took those 🤔 but I was also a bit wondering if they used other parameters before.

Hmm, wouldn't that mean that we calculate the dot products between pump/signal and pump/idler? I mean signal and idler both have a dagger, which would mean that we would need to multiply them with something without a dagger 🤔
@DorisReiter probably knows :)

Your explanation to (4) makes totally sense to me! And yes, I agree, if we don't go from conjugate transposed to just conjugated, it should be fine! (I think they assume the modes to be either almost TE or TM and thus can use just one component) so we probably just need to exchange in the code the fields to have the right overlap integrals?
It also makes sense to have the scalar product between pump/idler and pump/signal as those scalar products each describe a conversion of a photon 🤔

Uijj, that graph looks amazing! Happy to see that the mode solver in use 😊

And makes sense! If one of the mode is not really guided, the probability to the conversion to that mode should also be pretty small!

👍 thanks for moving it! The next step would be to add a header (just copy from https://github.com/HelgeGehring/femwell/blob/main/docs/photonics/examples/waveguide_modes.py) and divide the example into cells using

# %% [markdown]
# for explanations

# %%
print("for code!")

That way we can generate a nice docs page!
(If you use vscode you can use the jupytext extension to execute that python file like a jupyter notebook. But this way there's no problem with saving it, as it's pure python code and not a lot of json or whatever)
An we should add it in https://github.com/HelgeGehring/femwell/blob/main/docs/_toc.yml to have it show up in the docs

Digging deeper sounds great! I had during my master's also a look on the theory for SFWM and I'd be happy to give it a deeper drive :)
What about just digging deeper to explain SFWM completely? Would be a lot of interesting discussions 🥳 Would also be great if we could eventually also have the derivations for equations like (3) and (4) to understand the problem completely :)

What about keeping the main discussion in the issue #130 and merging pull requests as soon as we have a part of the problem solved?

[Kunyuan-LI]

Sure, I'll delve deeper into this, but we can include this PR as part of the problem.😄
For the documentation work, I'll handle it later after using this example as a solid standalone demonstration for calculating the SFWM. Always stay in contact, don't worry.

Now we have two problems:

1. In the example at 500x390 nm, based on the given wavelengths of p, s, i, I obtained a negative Aeff? Regarding the verification of the three modes at 500x390, I have placed it in /femwell/tests, which can ensure that all three modes have propagation distribution in the waveguide. So something must be wrong during the overlap calculation🤔
2. Based on 1, dot(w["E_p"][0], w["E_p"][0]) * dot(np.conj(w["E_s"][0]), np.conj(w["E_i"][0])) is used for overlap calculation, but when I change the order, like dot(w["E_p"][0], np.conj(w["E_s"][0])) * dot(w["E_p"][0], np.conj(w["E_i"][0])), it gives me the same result? I've checked the mechanism of the dot function in skfem.helper, it's implemented with np.einsum and not very clear to me, especially in this 4-matrix-multiplication scenario. Can we just use the explicit dot method like np.dot for testing and for better understanding of how it works? As you say, 'calculate the dot products between pump/signal and pump/idider', that could be a (s * p) (i * p) or (p s*) (p i*), or even something else?? 🤔

[HelgeGehring]

Yes, sounds great :)
I was more thinking to get things in the docs, that makes it easier to discuss. And also to see what's the new stuff, in too long PRs it's sometimes hard to see what's the changes since the last discussion :)

(1) what exactly do you test there? Tests should have the option to "fail", i.e. they need to include something like assertAlmostEqual or similar.
I'm also wondering how to do this best, at the moment I think the easiest would be to just include that in the docs and make them fail in case they don't show what we expect.
That way we don't have copies of code somewhere else, make sure the docs are correct and we also show based on what we test our code :)

And to see the location of the modes: Just include them in this script!
You can use

# %% tags=["hide-input", "hide-output"]

to not show it in the first place (but allow to open the input/output) and keep the example not too long. But for double-checking and understanding it's nice to have these plots :)

What do you think about this combined approach?

About the negative numbers: Don't worry!
Mode-solving is essentially just an eigen-solve of a matrix. That means that your eigenmode is just an eigenvector. Via the normalization of the mode we fix the norm of the mode, but the sign is not defined (in the case of a real-valued matrix. In case of a complex-valued matrix we get a random phase).
Thus, modes can randomly have a negative sign compared to the others. If one would have another sign than the others that would cause the area to be negative in this example (no imaginary part in the refractive indices -> real-valued eigenvectors)

We maybe could think of something like fixing the integral over E to be positive/have a zero phase? That would make the modes perfectly reproducible?
What do you think?

(2) From eq (4) it seems (p s*) (p i*) should be the correct one. For this example I'd expect it to only negligibly change the result, as the mode is quite polarized. Even if we just multiply the x-components, it shouldn't really change. (But I'd guess it's best to go for the most general, so the dot products seem like the best solution :) )

[Kunyuan-LI]

Hi Helge,

After rethinking the integration mechanism of skfem and the integrand functions in SFWM, I rewrote the calculation of the overlap integral:

# 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

I believe this is the only way to accurately describe the nature of the integrand functions being multiplied internally. Any approach using the dot(x,y) will introduce addition within the product, which leads to inconsistencies with the theoretical model.🤔

Therefore, this time I directly considered the contributions of Ex, Ey, and Ez in the x and y directions separately and then integrated. It looks better, but I still haven't fully explained why the integral directly becomes negative when the wavelength is longer than a specific value (usually for the lower-energy idler photon, and the larger the waveguide geometry, the higher this specific value).🤔

[HelgeGehring]

Hmm, like this I don't think you get something depending on E_z - the question is also what the last index in E_p_xy[0,:,0] does - the first one gives you either the x or the y component. I'd guess we shouldn't touch that one 🤔

I read a bit more into it and I think I remember the solution to this - the taylor expansion of the optical susceptibility is a tensor! (see https://en.wikipedia.org/wiki/Electric_susceptibility)

That means in the generalized case we'd need to multiply all combinations of the fields and then multiply it by the tensor element and sum it up.
But for example for silicon it reduces due to symmetry to way less numbers (see https://opg.optica.org/oe/fulltext.cfm?uri=oe-15-25-16604&id=148218) and for polarized modes it's just one number.

Thus, I think that's the reason why they just have a scalar? And therefore they also just multiply scalar fields in the paper? That would also mean that $A_{eff}$ is only defined for polarized modes as different combinations would differently be weighted? And it would also explain why we could get cross-polarization SFWM?

And as long as the absolute value of the $A_{eff}$ is correct we shouldn't worry - as said, each vector has a random sign (or complex phase in case of complex n)

In doi.org/10.1088/1612-2011/13/11/115204 is another example with an $A_{eff}$, maybe we can match it, too to check our assumptions?

[Kunyuan-LI]

Sure,

TOSPDC is also a very interesting non-linear optical phenomenon to study. The differences between SFWM and TOSPDC are discussed in detail in this literature: https://doi.org/10.48550/arXiv.1307.3266. However, the process of TOSPDC, where one pump photon splits into three photons, is indeed more counterintuitive compared to SFWM, which we can, to some extent, analogize with the collision of small balls (?).😂

For the E_p_xy[0,:,0], I'm certain there is a contribution in the integral. In fact, for tensors, considering our example (2,x,3), each unit actually has a (2,3) structure (in reality, xyz is (3,3)), which means that Ex, Ey, and Ez on each axis need to be considered. when I didn't consider the Ez component (x-Ez and y-Ez), I obtained an integral result that differed significantly and was quite far from the actual result.

Here is an example for 670x390nm. We selected three sets of p wavelengths and obtained the corresponding s and i wavelengths. Below is a comparison of the calculation results from femwell and ANSYS Lumerical using my previous overlap calculation method:

1. Pump - signal - idler 1.5 - 1.3045 - 1.7644 : femwell 0.3650617907002221μm2 and Lunerical 0.3677μm2
2. Pump - signal - idler 1.45 - 1.1903 - 1.8547 : femwell 0.3715068458043674μm2 and Lunerical 0.38021μm2
3. Pump - signal - idler 1.54 - 1.4735 - 1.6128 : femwell 0.35987430028852846μm2 and Lunerical 0.3587μm2

Considering the differences between our Sellmeier parameter fitting and Lumerical's fitting (which are observable), I believe this is completely reasonable. And this time, we have truly achieved it.🥳

As for TOSPDC, considering that the calculation method for overlap is p+s-i-r-, I believe the same implementation approach should also be completely feasible. However, I suggest creating a new PR and listing a new example because, firstly we need to reconsider the calculation of wavelengths. On the other hand, keywords such as SFWM and (TO)SPDC are also quite important.😂

[HelgeGehring]

Suggested change

	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 np.sum(E_p_xy * E_p_xy * np.conj(E_s_xy) * np.conj(E_i_xy), axis=0)

Does exactly the same :)
The indices are (direction - x/y , cell_number , quadrature point).

You can easily see this if you change the integration order of the basis you use for setting epsilon. The last index depends on the integration order.

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

will for example change the last index to 4.
The quadrature points is I think something we should leave to skfem, I don't see any situation where we would want to care about individual values on the quadrature points ourselves. (We just care about there being enough to be able to solve the system)

[HelgeGehring]

Suggested change

	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
	# Apply normalization factors to the electric fields for the overlap calculation
	def interpolate_and_normalize_xy(mode):
	@Functional
	def E2(w):
	return dot(w.E, np.conj(w.E))
	E = mode.basis.interpolate(mode.E)[0]
	return E / np.sqrt(E2.assemble(mode.basis, E=E))
	@Functional(dtype=np.complex64)
	def sfwm_overlap(w):
	return np.sum(w.E_p_xy * w.E_p_xy * np.conj(w.E_s_xy) * np.conj(w.E_i_xy), axis=0)
	return 1 / sfwm_overlap.assemble(
	basis,
	E_p_xy=interpolate_and_normalize_xy(mode_p),
	E_s_xy=interpolate_and_normalize_xy(mode_s),
	E_i_xy=interpolate_and_normalize_xy(mode_i),
	)

This would make it even shorter and a bit more efficient as we don't have to interpolate twice :)

[HelgeGehring]

I'd guess we should use here

Suggested change

	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
	def n_X(wavelength):
	x = wavelength * 1000
	return (
	1
	+ 3.0249 / (1 - (135.3406 / x) ** 2)
	+ 40314 / (1 - (1239842 / x) ** 2)
	) ** 0.5

Then it matches way better the linked example, there they got 0.59, I get 0.6 :)
(And we should probably call it n_silicon_nitride)

[HelgeGehring]

YAY, that's just amazing! 🥳 🎉

I'm wondering if we should generalize the method for calculating the effective area a bit more to consider $\chi^{(3)}$ being a tensor. That would allow also calculating things like cross-polarization mixing 🤔 What do you think? Maybe better in another example? Or for now just explain it in the text that we made that simplification?

Thanks for the link, I'll have a look! And nice analogon with the small balls! :D

The indices are (vector-direction, cell_number, quadrature point). I don't really see an application where we would need to care about the quadrature points in femwell, that's more something deep in FEM. Usually, we just care about having enough quadrature points to have a defined system to solve.
You can see that it's the quadrature points if you increase the intorder in the basis which is defined for epsilon, i.e.

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

Wow, that looks like perfect matching! 🥳 🥳

We should now care about turning the code into a real example by adding some explanations and maybe a graph to look at?

How about showing first the graph for the phase-matching you posted?
And then we show a heatmap for the effective area / gamma for different signal / idler wavelengths? (maybe also a heatmap for phase matching?)
And in the end we could also show $N_{\text{puls}}$ (eq 20) as a heatmap? (maybe in another PR 🤔 )
That would be a somehow very nice complete story?

TOSPDC sounds also super interesting! Let's investigate that deeper in a new PR!

[HelgeGehring]

I've just added some markdown and split the example in multiple cells to have it ready for the docs :)

Have a look on the suggestion to simplify the code, didn't want to just change that :) but I'd guess that would make it easier to read and shorter :)

[Kunyuan-LI]

Hi, I'm back. The material of the figure has a multi-layer structure based on TiO2 and various testing materials, so it's a little complicated. For now, we can treat it as an isotropic material to simplify the analysis. However, you're right that using silicon nitride directly for the example would be better.

As for the TOSPDC, I'll create a 3D contour surface that describes the phase matching condition in x, y, z (signal, idler, rest) space. This will provide a clear visualization of the phase matching.

[HelgeGehring]

Yay! Ah, okay, I thought we were matching the paper :)

Sounds cool!

Have a look on the suggestions above, those will make the code way easier to read :)