Ez field is not continuous at the dielectric boundaries when computing TM modes of a simple 2D waveguide

[rafael-fuente]

I'm testing the 2D plane wave expansion of legume computing the TM modes of a simple 2D slab waveguide, a problem that has an analytical solution.

This is my code:

import numpy as np
import matplotlib.pyplot as plt

from legume import PlaneWaveExp, Circle, ShapesLayer, Lattice, viz
import legume


eps_waveguide = 3.5**2  
eps_background = 1

d = 1.0  # waveguide width
w = 0.5
# Initialize lattice
lattice = Lattice(np.array([w,0.]) ,np.array( [0.,8.0])  )

# Initialize layer
layer = ShapesLayer(lattice, eps_b = eps_background)
layer.add_shape(legume.Poly(x_edges=[-0.5*w, 0.5*w, 0.5*w, -0.5*w], y_edges= [-0.5*d, -0.5*d, 0.5*d, 0.5*d], eps=eps_waveguide))

kx = np.linspace(-2.0, 2.0, 2)
ky = kx*0
kpoints = np.c_[kx,ky]
pwe = PlaneWaveExp(layer, gmax=7)
viz.eps_ft(pwe, figsize = (3,8))
path = lattice.bz_path(kpoints, [1])

pwe.run(kpoints=path['kpoints'], pol='tm', numeig=2)
freqs_tm = pwe.freqs

print('TM freqs: ', freqs_tm[-1,:])

The code outputs the correct TM frequencies [0.11859563, 0.19938521], which match the ones computed with the analytical solution.

However, when I visualize the Ez and Dz fields, I found the Ez is discontinuous and Dz is continuous at the waveguide boundaries.

viz.field(pwe, field='D', kind=-1, mind=0,
            component='z', val='abs', N1=100, N2=100, cbar=True, eps=True);
viz.field(pwe, field='E', kind=-1, mind=0,
            component='z', val='abs', N1=100, N2=100, cbar=True, eps=True);

Since these components are perpendicular to the XY plane where the waveguide is defined, from Maxwell equations, it's expected Ez is continuous since the electric field must fulfill at the waveguide boundaries $\hat{n} \times\left(\overline{\mathrm{E}}_1-\overline{\mathrm{E}}_2\right)=0$. From this relation, it can also be seen that Dz must be discontinuous.

Why am I obtaining the opposite results? Is a problem related with legume or am I misinterpreting some output?

[momchil-flex]

I think this may be a numerical artifact when using a mode expansion method (here plane waves). If you want to investigate further, maybe try plotting a 1D cross section, with higher spatial resolution, and zoomed in around the boundary. I wonder if what we're seeing here is related to the Gibbs phenomenon (you're basically doing with Fourier transforms here too). I'm not sure though, because you're using quite high gmax so I would have thought that it should not be so pronounced. Maybe there is a bug in the plotting - if so, I am not sure when I'll find time to investigate in more detail though.

By the way note that legume has a dedicated function to compute (multi-layer) slab modes: https://github.com/fancompute/legume/blob/master/legume/gme/slab_modes.py, which are the basis for the guided-mode expansion. Just in case it is useful for what you are doing (these should be the exact analytical solutions).

[rafael-fuente]

The problem is not a numerical artifact. Due to plane wave expansion, there is always Gibbs phenomenon, but it doesn't cause the discontinuity of Ez. In fact, Ez shouldn't have Gibbs phenomena in the first place because it shouldn't have a discontinuity.

I solved exactly the same problem with MPB and it's outputting the correct fields, with Ez being continuous and Dz discontinuous, so I think it's a legume bug.
MPB and legume frequencies are the same, so the problem isn't there, it's only the fields.

[momchil-flex]

I see, thanks for doing that comparison, I agree that it looks like there's a bug in the field getting or plotting in legume. I'll let you know when I find some time to look into it if I discover what it is or manage to fix it.