mode overlap

[joamatab]

How can I do a mode overlap from the tidy3d modesolver?

https://gdsfactory.github.io/gdsfactory/notebooks/plugins/tidy3d/01_tidy3d_modes.html#Bend-loss

CC @momchil-flex @twhughes @flaport

[flaport]

something like this should work, I think:

    cross = mode1.Ex * mode2.Hy.conj() - mode1.Ey * mode2.Hx.conj()
    return 0.25*np.trapz(np.trapz(cross, y), x)

Just be sure to zero out the phase before doing so. Something like this in emepy could do the job. Although it's not perfect, I think. Maybe better to find the location of maximum EM energy density in the mode and then use that phase to zero out the phase.

[momchil-flex]

Yeah what @flaport describes is pretty much what we do on the backend when you have a ModeMonitor. I will think about moving this part (specifically the mode overlap function) to the frontend.

[joamatab]

yes Momchil, i think the mode-overlap calculations would be very useful for the client

Hi Floris, i get a ValueError

https://github.com/gdsfactory/gdsfactory/compare/master...5112 https://github.com/gdsfactory/gdsfactory/commit/bdd15e02d22cb84add40dcd9da146627f65114d9

ValueError: operands could not be broadcast together with shapes (1,100) (100,9)

@simbilod

[joamatab]

@helgegehring got it working for the latest gdsfactory version

[smartalecH]

Just passing by @flaport -- but you shouldn't need the conjugate. If the materials are lossless (and the modes are normalized right) then the conjugation won't change the results anyway. But the general mode overlap integral (derived from the unconjugated reciprocity theorem in Snyder and Love) is an unconjugated inner product.

The conjugated inner product we have in meep is more due to convention (physicists seem to like conjugated inner products...) and it's actually a poor choice because the current formulation isn't even holomorphic with the fields, but rather conjugate holomorphic... so the corresponding adjoint is a mess and needs to be reformulated.

(we also need to fix the overlaps in EMEpy...)

PS very nice package @momchil-flex et al!

[momchil-flex]

Yeah this is quite a debate, so hopefully you're only interested in lossless waveguides cause it's best to avoid it... :)

Tidy3D uses the complex conjugate version too. It is also the default in Lumerical FDTD by the way. The advantage is that the interpretation is still power carried by each mode; the disadvantage is that orthonormality does not hold, which can introduce all sorts of trickiness.

[flaport]

That's very interesting @momchil-flex, @smartalecH , I'm definitely not an expert on this particular subject but I ran into an issue with the mode overlaps in emepy not so long ago, resulting in S-parameters bigger than one. Here's an issue I opened for emepy which I'm pretty sure is related to this whole discussion

https://github.com/BYUCamachoLab/emepy/issues/12

What would be your approach to solving that problem there?

[momchil-flex]

From the last reply in that thread, it seems like the modes are orthogonal, so I don't know if this discussion applies.

Generally, it seems to me that our intuitive understanding of scattering matrix is very much predicated on assuming lossless ports. Inside, there could be loss, and that would just make some singular values smaller than one, but the ports are lossless, which means that they are also orthogonal w.r.t. the inner product that can be interpreted as power carried by each mode. In the lossy case there may be no perfect answer - the two inner products give two different results, each with its own meaning. The conjugated one gives you the power carried in each mode, but the modes are not orthogonal, so e.g. P1 + P2 can be > P0, where P0 is input power, and P1, P2, is output carried by modes 1 and 2 (which should translate to singular values larger than 1). In the non-conjugated case, the scattering matrix singular values are <= 1, but the physical interpretation is more challenging.

There's probably more on this that I am not familiar with... I'll also ask Shanhui when I get a chance.

[momchil-flex]

FieldData.dot and ModeSolverData.dot implemented in v1.8.0

[aisichenko]

I'd like to take the mode overlap of two FieldData instances that have a different grid, using the outer_dot method.

One field is type FieldData from a monitor. The other field is type ModeSolverData. I get a KeyError that I am missing the mode_index for the FieldData. The outer_dot documentation, however, seems to suggest that this can be done between these. Perhaps I'm missing something? Thanks!

[momchil-flex]

I think this should be possible yeah, so probably something unexpected is happening. We may need a minimal example to debug. Or maybe @lucas-flexcompute knows better.

[aisichenko]

For sure, see this example:

field1.outer_dot(field1)

where field1 is type FieldData. This throws a KeyError: 'mode_index'.

for field2 of type ModeSolverData, this works as expected:

field2.outer_dot(field2)

[jan-david-fischbach]

Hey @momchil-flex, from your reply (https://github.com/flexcompute/tidy3d/issues/423#issuecomment-1172947318) I believe to understand, that the unconjugated inner product of two distinct modes that the mode solver spits out should be zero (even for lossy modes, e.g. due to PML boundaries), is that correct?

FYI @flaport

[momchil-flex]

I believe so, but this is based on my theoretical understanding rather than actual testing.

[jan-david-fischbach]

Ah ok, we have been seeing some nonorthogonality with respect to the unconjugated inner product. Does anything speak against orthogonalizing (e.g. by Gram Schmidt or similar) the modes from your perspective?

[jan-david-fischbach]

I just realized, that I would have no clue how to find the neff of the then orthogonalized modes...

[momchil-flex]

Right, because they won't be eigenmodes.

It would be good to understand what exactly is causing the nonorthogonality. Do you have a minimal working example? Disclaimer: it may be hard for me to find time to look into it right way.

[jan-david-fischbach]

@flaport can you maybe provide the example? Or I'll sit down later and provide one.

[jan-david-fischbach]

A very simple example using the waveguide-plugin: https://colab.research.google.com/drive/1povHWP-lPL8fSQIlFlW_NQf6KIspo6xe?usp=sharing Note, that the concerning mode is a non guided mode.

[momchil-flex]

So the mode that has some nonzero overlap with other modes is completely unphysical and localized around the PML.

I am not sure if the usual arguments that ensure that the un-conjugated dot product should be orthogonal apply to such modes. Maybe there's some assumption on the permittivity which is broken by PML, which is a very special beast. However, what I do expect is that physical modes that have some loss because of the PML (e.g. modes below the light-line which are only partially confined to the waveguide, but still mostly concentrated there) should have ideally zero overlap with other physical modes. Well, I mean if the first part of my statement is true that maybe PML break something, then certainly there may be a continuous transition where if too much of the field is in a PML, the orthogonality breaks. However, I don't think physical modes should have any visible field in the PML...

Just as a test I tried your notebook using cSi with the lossy Green2008 variant at wavelengths 0.6um and 1um, but without PML. All three modes are physical and the highest overlap is on the 1e-6 level (interestingly, for the conjugated dot product it still looks pretty good, but highest overlap is on the 1e-5 scale).

I am not sure what you are trying to do, but if you are using PML, I think it may be best to filter such modes out. For example, compute the integral of the field intensity over the PML region over the integral over the whole region, and discard modes for which this is higher than some threshold value, e.g. 10%.

[jan-david-fischbach]

Back again with a better example: https://colab.research.google.com/drive/1WrWb2t0fnnBbLs88ElCc7vFRGvcNbO1b?usp=sharing the concerning modes are still not well guided. However they are also persent when using metal bcs and are also returned when using Lumerical with PML (the modes from lumerical are orthogonal however)

[flaport]

Complementary to the example Jan-David gave... Attached you can find the modes from lumerical. As you can see, even though visually the modes (at least their high intensity field components) seem identical, their version of the modes seem to be orthogonal.

Overlaps

Fields

Fields Array

lumerical_fields.npy.zip

• Note1: dimensions of the attached array: (modes, fields [Ex, Ey, Ez, Hx, Hy, Hz], x, y, z)
• Note2: I had to convert the lumerical fields to single precision and zip the npy to stay under the memory limit of an attachment in a github thread.

[momchil-flex]

Interesting! Yeah I see the issue in this example, seems a bit nontrivial to figure out, but still worth looking into it.

So did you use some lumerical built-in to compute the overlap, or did you do it manually based on the fields? And what are the off-diagonal values? Interestingly, if I use the conjugated dot product, I get diagonal elements of 1, while the off-diagonal elements become large.

[flaport]

I used the lumerical built-in way to calculate the overlaps and also confirmed that the modes are orthogonal by manually calculating the overlaps with meow. I could probably use tidy3d as well to doubly confirm the lumerical modes are orthogonal, but I think the overlap calculation is done pretty similarly anyway (between meow and tidy3d).

[momchil-flex]

After looking into this a little bit, it seems to me that the Lumerical orthogonality for modes 2 and 3 comes not so much from the actual integration and some accurate interference in the summation to zero everything out, but rather simply by the fact that they have opposite symmetry w.r.t. the x = 0 plane. Case in point, I don't actually need to apply the diff area elements of the non-uniform grid (I don't have those in the data), and the overlap is 0 with both conjugated and non-conjugated product.

I think the reason things don't come out to zero in Tidy3D may be some off-by-one issue at the domain boundaries, which would manifest itself exactly for such modes which do not decay there. In other words I think what may be happening is that one pixel less or more is counted on the left than on the right, which results in the overlap not going to zero due to symmetry.

Now, that said, I would believe that the overlap should go to zero even without such an implicit symmetry, but I'm starting to wonder... I'll look into fixing Tidy3D, but meanwhile I wonder what you'll get in Lumerical if you offset the waveguide to one side.

[flaport]

Hi I did offset the waveguide by a single pixel in lumerical, but the results seem still pretty stable and definitely orthogonal.

But your comment made me thing of something. We're working on a Yee-Grid, so implicitely when we cut off the grid we have the current situation (note, x is vertical and y is horizontal to match with a numpy array layout):

current situation:

Ez[m-2,n-2]      Ey/Hx[m-2,n-2]   Ez[m-2,n-1]      Ey/Hx[m-2,n-1] | Ez[m-2,n]      Ey/Hx[m-2,n]
                                                                  |
Ex/Hy[m-2,n-2]   Hz[m-2,n-2]      Ex/Hy[m-2,n-1]   Hz[m-2,n-1]    | Ex/Hy[m-2,n]   Hz[m-2,n]
                                                                  |
Ez[m-1,n-2]      Ey/Hx[m-1,n-2]   Ez[m-1,n-1]      Ey/Hx[m-1,n-1] | Ez[m-1,n]      Ey/Hx[m-1,n]
                                                                  |
Ex/Hy[m-1,n-2]   Hz[m-1,n-2]      Ex/Hy[m-1,n-1]   Hz[m-1,n-1]    | Ex/Hy[m-1,n]   Hz[m-1,n]
------------------------------------------------------------------|
Ez[m,n-2]        Ey/Hx[m,n-2]     Ez[m,n-1]        Ey/Hx[m,n-1]     Ez[m,n]        Ey/Hx[m,n]

Ex/Hy[m,n-2]     Hz[m,n-2]        Ex/Hy[m,n-1]     Hz[m,n-1]        Ex/Hy[m,n]     Hz[m,n]

which implies for the input permittivity eps_xx.shape == eps_yy.shape == eps_zz.shape == (m, n)

whereas to have a more symmetrical situation and boundary conditions, we probaby want the same fields on the left vs right and top vs bottom. So maybe something like this:

Ez[m-2,n-2]      Ey/Hx[m-2,n-2]   Ez[m-2,n-1]      Ey/Hx[m-2,n-1]   Ez[m-2,n]    | Ey/Hx[m-2,n]
                                                                                 |
Ex/Hy[m-2,n-2]   Hz[m-2,n-2]      Ex/Hy[m-2,n-1]   Hz[m-2,n-1]      Ex/Hy[m-2,n] | Hz[m-2,n]
                                                                                 |
Ez[m-1,n-2]      Ey/Hx[m-1,n-2]   Ez[m-1,n-1]      Ey/Hx[m-1,n-1]   Ez[m-1,n]    | Ey/Hx[m-1,n]
                                                                                 |
Ex/Hy[m-1,n-2]   Hz[m-1,n-2]      Ex/Hy[m-1,n-1]   Hz[m-1,n-1]      Ex/Hy[m-1,n] | Hz[m-1,n]
                                                                                 |
Ez[m,n-2]        Ey/Hx[m,n-2]     Ez[m,n-1]        Ey/Hx[m,n-1]     Ez[m,n]      | Ey/Hx[m,n]
---------------------------------------------------------------------------------|
Ex/Hy[m,n-2]     Hz[m,n-2]        Ex/Hy[m,n-1]     Hz[m,n-1]        Ex/Hy[m,n]     Hz[m,n]

which implies for the input permittivity: eps_x.shape == (m, n+1) eps_y.shape == (m+1, n) eps_z.shape == (m+1, n+1)

After getting those results, you can (probably?) still truncate those E and H fields. to have equal size but at least the boundary conditions are now equivalent for left vs right and top vs bottom.

I'm not a 100% sure about this, so let me know if I'm talking nonsense, haha

Do you think it could be related to something like this?

[momchil-flex]

Yeah, this is definitely not nonsense, and has been making my life hard every now and then, since the beginning of Tidy3D. The approach we have under the hood is that everything is periodic if nothing special is done, so all field components can be arrays of the same shape. Then boundaries are applied on top: for example, in the mode solver, there is no tangential E-field stored in the max-side of the simulation domain, but the fact that it is zero (PEC) is applied through the derivative matrix. In the FDTD solver where various boundary conditions are supported, it's more complicated...

The need to interpolate fields to the same locations is also an issue, but it is also done under the hood in Tidy3D if you e.g. use the in-built dot method.

I discovered something slightly off specifically with the pml coefficients in the mode solver (most of the time, we don't apply pml, so these had not been as well tested as the rest of the mode solver plugin). I made a PR to fix this #939 and it seemed to improve your example quite a bit. For starters, it seems to produce fewer spurious modes if I reduce the plane size in the x-direction. Previously, I would not get the same modes if I set e.g. the x size to 12 or 16, now I seem to consistently get the same 4 modes. Secondly, the orthogonality is much better:

It even works fairly well when I offset the waveguide by 1 micron to the right:

This is not ideal, but not terrible either. I wonder how Lumerical's modes do. Maybe there's still something small that's off...

[momchil-flex]

By the way @aisichenko note that the outer_dot has been extended to work for FieldData too in our pre/2.3 branch. Related to that @flaport @Jan-David-Black note that the method now returns a regular xr.DataArray rather than a Tidy3D specific data structure, so you need to do np.abs(overlap) instead of overlap.abs in your scripts.

[flaport]

Thanks a lot @momchil-flex . The new version works a lot better. For comparison, here is the overlap integral for tidy3d vs lumerical:

note:

The lumerical modes are normalized differently, so the values in the lumerical overlap matrix were very large. I divided the resulting overlap matrix by the diagonal values as a quick normalization. I did the same for the Tidy3d modes (even though the diagonal is already close to 1) for the fairest comparison.

No x-offset

Tidy3d:

Lumerical:

1um x-offset

Tidy3d:

Lumerical:

As you can see, the lumerical modes are less sensitive to an offset, but I'm not too worried anymore as the tidy3d values are pretty close to the lumerical values now.