Radiative losses examples with ModeSolver

[FilipeFcp]

Hi all,

Here are some examples and benchmarks using the ModeSolver to calculate radiative losses. I tried to address some concerns that I have seen from users, such as convergence issues and precision. I believe the message is that the results converge and the ModeSolver captures the general behavior of the losses; however, if you want to be picky with precision, especially for low losses, high resolution and plane size are needed.

Let me know your thoughts.

[tomflexcompute]

Thanks @FilipeFcp for the great benchmark. It works very well. Some comments:

• "imaginay" -> "imaginary"

• "if the exact value are important" -> " if the exact values are important"

• "Losses in bended optical fiber" -> "Losses in a bent optical fiber"

• "Bend Si waveguide" -> "Bent Si waveguide"

• "known beforehand.To automatically" -> "known beforehand. To automatically"

• The comment above cell[12] seems unfinished.

• The comment below [24] seems unfinished too. Missing period in the end of the sentence.

• Add plot titles or axis labels for [13] and [15]

• Add plt.show() in all cells with plots to get rid of the text print out.

• When defining the first mode solver, we can add a few comments and discussions and use this opportunity to educate users on the details of the setup such as

1. Add comments about num_pml in ModeSpec. Explain the default boundary is PEC and for lossy modes, PML needs to be used.
2. Also comment on the necessity of using double precision since the loss is very small and want to have high accuracy.
3. PML layers are added inside the domain so the mode solving plane needs to be sufficiently large.

The use of isMode function is a great idea! A lot of users are looking to do filtering by something like this so it's a very good reference.

[tomflexcompute]

Oh yeah don't forget to add the notebook to the tidy3d-notebooks/docs/features/mode.rst so it will be built into the docs and learning center.

Also the metadata like other tutorial notebooks.

[e-g-melo]

Hi @FilipeFcp! Very nice and useful tutorial!

1. I recommend including cell [1] before the "Defining the simulation domain" header to make it uniform with our other notebooks.
2. Not critical, but I wonder why naming that function as mode_solver_substrate :)
3. Suppress the text output in cell [3], [10], and [17].
4. Here, should the parenthesis come just after E.y.max() - E.y.min()?

    center = (
        E.y.min() + (E.y.max() - E.y.min() / 2),
        E.z.min() + (E.z.max() - E.z.min()) / 2,
    )

5. In the simulation size convergence test you are adding an extra size to the mode plane, but the simulation size is constant in 6 at y and z directions. So, at some point the mode plane will be larger than the simulation domain and the extra_size will not make difference on results, is that correct?
6. I suggest setting the mode plane size directly (maybe making it squared) instead of setting an extra_size so that users don't need to calculate it by themselves. Unless you are reproducing the paper methodology.
7. Something is missing here "index. however, for the imaginay part, a". Additionally, replace "however" with "However".
8. "tune the simulation size" -> "tune the mode plane size"
9. We need axis labels in plots from cells [13], [15], and [22].
10. In [19] and [21]: ax.set_xlabel("Loss (dB/cm)") -> ax.set_ylabel("Loss (dB/cm)")
11. This sentence seems not so clear "It is possible to note that, the results are converging for values around. For resolutions up to x, a value difference of about x can be seen.". And the sentence after it seems duplicated.
12. "We will benchmark the losses for a 30 mm bend radius" -> "We will benchmark the losses for a 30 um bend radius"

[FilipeFcp]

Hi @tomflexcompute and @e-g-melo,

Thanks for the comments. Sorry, I was working on the last part of the notebook and thought I had already revised the second part, so some text was left incomplete.

Anyway, I believe I’ve addressed all the comments. Just to confirm the last point from @e-g-melo: it is indeed 30 mm (or 30,000 µm).

Thank you all!

[momchil-flex]

Hey, thanks for this!

I'm surprised by this slow convergence of the substrate loss with resolution.

Our backend test for the same thing gets Computed complex n: (2.409505649017421+2.9855138851557985e-08j) at resolution 40. Here's the code. Could you see how this compares to your definitions?

def mode_solver_substrate(res, npml):
    """Mode solver for a waveguide with leakage to the substrate as analyzed in  Bienstman et al.,
    "Modelling leaky photonic wires: A mode solver comparison".
    """

    wg = td.Structure(
        geometry=td.Box(size=(1, 0.5, 0.22)),
        medium=td.Medium(permittivity=3.5**2)
    )
    box = td.Structure(
        geometry=td.Box.from_bounds(rmin=(-td.inf, -td.inf, -100), rmax=(td.inf, td.inf, -0.11)),
        medium=td.Medium(permittivity=1.45**2)
    )
    substrate = td.Structure(
        geometry=td.Box.from_bounds(rmin=(-td.inf, -td.inf, -100), rmax=(td.inf, td.inf, -1.11)),
        medium=td.Medium(permittivity=3.5**2)
    )

    sim = td.Simulation(
        size=(4, 6, 6),
        grid_spec=td.GridSpec.auto(min_steps_per_wvl=res, wavelength=1.55),
        # grid_spec=td.GridSpec.uniform(dl=0.010),
        structures=[wg, box, substrate],
        run_time=1e-12,
        boundary_spec=td.BoundarySpec.all_sides(boundary=td.PML()),
        symmetry=(0, -1, 0)
    )
    bms = BackendModeSolver(
        simulation=sim,
        # plane=td.Box.from_bounds(rmin=(0, -1.5, -2.7), rmax=(0, 1.5, 1.3)),
        plane=td.Box.from_bounds(rmin=(0, -1.5, -2.3), rmax=(0, 1.5, 1)),
        mode_spec=td.ModeSpec(num_modes=3, target_neff=2.41, num_pml=(npml, npml), precision="double"),
        freqs=[td.C_0 / 1.55]
    )
    return bms

def test_mode_pml_substrate():
    """Test mode solver pml for leakage into the substrate."""
    bms = mode_solver_substrate(res=40, npml=10)
    modes = bms.data
    wg_mode_ind = np.argmin(np.abs(modes.n_eff.values.squeeze() - 2.41))
    n_complex = modes.n_complex.isel(mode_index=wg_mode_ind).values.squeeze()
    print("Computed complex n: ", n_complex)
    assert np.allclose(n_complex.real, 2.41, rtol=1e-2)
    # There is no consensus result for the imaginary index in the paper.
    # However, 2.9e-8 is a common result among multiple solvers (see Table 6)
    assert np.allclose(n_complex.imag, 2.9e-8, rtol=1e-2)

For the bend waveguide, I think you should use SiN instead of Si. Silicon is too good at confining the light, and you need way too small of a bend radius to get some losses. I don't even know if the mode solver is expected to produce meaningful results if the size of the mode plane is larger than the bend radius, although it seems to be doing well. I think it may be fine as long as the mode is well-localized to the waveguide region. But still, this study is usually much more relevant for lower-index waveguides?

For the optical fiber, is there some value we can compare to?

[FilipeFcp]

Hi @momchil-flex,

The function I used is the same, I just changed the number of PMLs and plane size for running the benchmark. For a resolution of 40, with the parameters that I used, the result is 3.098014e-08. My goal was to see if we could get closer to 2.91e-8, which is the value that the methods described in the paper seem to converge to; that is why I changed the size and number of PMLs.
With the parameters of the test function, the value seems to increase linearly with the resolution:

In my opinion it is a bit picky to try to get that close to the paper value, but sometimes it seems that users can be excessively concerned with numerical error for these low quantities. My idea was just to show that it does converge if we fine-tune the parameters, so the only issue is discretization, which no numerical method can completely avoid.

I will try with the SiN waveguide. I was not sure how meaningful the results would be for a high bend radius, but since they follow some constancy, it looked ok to me.

For the fiber value, there are some analytical formulations. I will take a look and see if we can at least agree on the order of magnitude.

[e-g-melo]

Thanks @FilipeFcp!

A final comment: In the mode plane size convergence chart the x-label should be "mode plane size (um)".

[momchil-flex]

Hi @momchil-flex,

The function I used is the same, I just changed the number of PMLs and plane size for running the benchmark. For a resolution of 40, with the parameters that I used, the result is 3.098014e-08. My goal was to see if we could get closer to 2.91e-8, which is the value that the methods described in the paper seem to converge to; that is why I changed the size and number of PMLs. With the parameters of the test function, the value seems to increase linearly with the resolution:

In my opinion it is a bit picky to try to get that close to the paper value, but sometimes it seems that users can be excessively concerned with numerical error for these low quantities. My idea was just to show that it does converge if we fine-tune the parameters, so the only issue is discretization, which no numerical method can completely avoid.

I see. Pretty finicky huh. Sounds good then.

I will try with the SiN waveguide. I was not sure how meaningful the results would be for a high bend radius, but since they follow some constancy, it looked ok to me.

Thanks.

For the fiber value, there are some analytical formulations. I will take a look and see if we can at least agree on the order of magnitude.

Yeah, that would be good.

[momchil-flex]

I will try with the SiN waveguide. I was not sure how meaningful the results would be for a high bend radius, but since they follow some constancy, it looked ok to me.

Yeah actually I think we are going to validate this out in future releases (at some point). I just realized there was a recent user error because of this. If the coordinates shifted such that the center of the mode plane is at the radial distance happen to contain a coordinate at 0, there will be a division by zero error in the bend transformation.

[tomflexcompute]

@FilipeFcp all looks good but I think metadata is still missing?

[FilipeFcp]

Hi all,

I made some final adjustments and additional tests. I found a paper with an analytical solution for the fiber losses, and after some fine-tuning, the results seem to match quite well.

@momchil-flex, since the fiber example already deals with low-index contrast, do you think it would be okay to keep the Si waveguide instead of the SiN one? That way, we would have one example with a high-confined and one with a low-confined waveguide.

[momchil-flex]

@momchil-flex, since the fiber example already deals with low-index contrast, do you think it would be okay to keep the Si waveguide instead of the SiN one? That way, we would have one example with a high-confined and one with a low-confined waveguide.

I have actually added a PR in our upcoming release which will prohibit the bend radius to be smaller than half the mode plane size along the radial direction. This will break your example, right?

I did this because this can error in some cases. One user hit such an error.

[FilipeFcp]

Yes, it will. I will change them to SiN.

[FilipeFcp]

@momchil-flex I switched to SiN. The size is 20 µm, and the minimum bend radius is also 20 µm, so I believe it shouldn't be a problem.

[e-g-melo]

@FilipeFcp

[FilipeFcp]

Thank you @e-g-melo

[momchil-flex]

Thanks @FilipeFcp ! I think this notebook is looking really good now. Just a few last comments on the SiN part:

• Maybe add a cross-sectional plot of the mode solver. It will just look like a waveguide (without the bend explicitly visible), but at least it's clear. Also these cross-sections are plotted for the other two examples.
• You have changed the index to 2 but still calling the medium si: si = td.Medium(permittivity=2**2)
• Is such a large size (20) really needed? And even if so, I think you should definitely be able to keep the z size fixed to some much smaller value, i.e. it should only be needed in the radial direction.

[e-g-melo]

@FilipeFcp, I will use this image for the GUI version. Maybe you could also include it.

[tomflexcompute]

@FilipeFcp, I will use this image for the GUI version. Maybe you could also include it.

Very nice! @FilipeFcp you can add it to the notebook as well as metadata just like the regular example notebook.

[FilipeFcp]

@momchil-flex That is true, the z-size can be smaller.
This is a plot of the losses (db/cm) varying the sizes, for a bend radius of 130 um:

This is a extreme case, the k_eff is at the order of $10^{-13}$