Manual eigenmode-overlap calculation is not matching built-in eigenmode-coefficients feature

[lxvm]

Hi, I've observed an accuracy issue with Meep's built-in get_eigenmode_coefficients feature whose result for the eigenmode coefficient in on a cross section of an x-invariant rectangular waveguide does not match an overlap integral computed from a combined DFT mode monitor and get_eigenmode call on the same cross section.

Specifically, the result of code (1) using the get_eigenmode_coefficients feature does not match the result of code (2) which uses a get_eigenmode call on a DFT fields object

Code (1):

import meep as mp
import numpy as np

verbosity = 0
waveguide_index = 3.476
SiO2_index = 1.428
resolution = 30
wavelen = 1.55
waveguide_height = 0.5
# waveguide_width = 0.5
waveguide_padding = 1.0
relative_fwidth = 0.2
cell_width = 3.0
cell_height = 4.0
pml_padding = 1.0
coupling_interface_gap = 1.0
px_offsets_x = np.linspace(0, 2, 21)
px_offsets = [mp.Vector3(x=px) for px in px_offsets_x]

# indices reported by https://doi.org/10.1364/OSAC.1.000013
SiO2 = mp.Medium(index=SiO2_index)
waveguide_medium = mp.Medium(index=waveguide_index)

fcen = 1 / wavelen
fwidth = relative_fwidth * fcen

Sx = cell_width + 2*pml_padding
Sy = cell_height + 2*pml_padding
Sz = 0.0
cell_size = mp.Vector3(Sx, Sy, Sz)
pml_layers = [mp.PML(pml_padding)]
# needs to be a few wavelengths away from the structure

kpoint = mp.Vector3(-1, 0, 0)
source_center = mp.Vector3(x=coupling_interface_gap/2)
source_size = mp.Vector3(y=waveguide_height + 2 * waveguide_padding)
src = mp.GaussianSource(frequency=fcen, fwidth=fwidth)

source = [
    mp.EigenModeSource(
        src,
        eig_band=1,
        direction=mp.NO_DIRECTION,
        eig_kpoint=kpoint,
        size=source_size,
        center=source_center,
    ),
]
eig_power = source[0].eig_power(fcen)

geometry = [
    mp.Block(
        center=mp.Vector3(), material=waveguide_medium, size=mp.Vector3(Sx, waveguide_height, 0.0)
    ),  # center waveguide
]

sim = mp.Simulation( 
    cell_size=cell_size,
    boundary_layers=pml_layers,
    geometry=geometry,
    sources=source,
    eps_averaging=True,
    resolution=resolution,
    default_material=SiO2,
    # force_complex_fields=True,
)

offset_mode_monitors = [] 
for px_offset in px_offsets:
    guide_fr = mp.FluxRegion(center = -1*source_center + px_offset / resolution, size = source_size)
    # guide_fr = mp.FluxRegion(center = -1 * source_center, size = source_size)
    guide = sim.add_mode_monitor(fcen, 0, 1, guide_fr)
    offset_mode_monitors.append(guide)

sim.run(until_after_sources=mp.stop_when_energy_decayed(dt=1, decay_by=1e-11))
eig_powers = [np.abs(sim.get_eigenmode_coefficients(guide, [1]).alpha.flatten()[1]) ** 2 for guide in offset_mode_monitors]

Code (2):

import meep as mp
import numpy as np

verbosity = 0
waveguide_index = 3.476
SiO2_index = 1.428
resolution = 30
wavelen = 1.55
waveguide_height = 0.5
# waveguide_width = 0.5
waveguide_padding = 1.0
relative_fwidth = 0.2
cell_width = 3.0
cell_height = 4.0
pml_padding = 1.0
coupling_interface_gap = 1.0
px_offsets_x = np.linspace(0, 2, 21)
px_offsets = [mp.Vector3(x=px) for px in px_offsets_x]

# indices reported by https://doi.org/10.1364/OSAC.1.000013
SiO2 = mp.Medium(index=SiO2_index)
waveguide_medium = mp.Medium(index=waveguide_index)

fcen = 1 / wavelen
fwidth = relative_fwidth * fcen

Sx = cell_width + 2*pml_padding
Sy = cell_height + 2*pml_padding
Sz = 0.0
cell_size = mp.Vector3(Sx, Sy, Sz)
pml_layers = [mp.PML(pml_padding)]
# needs to be a few wavelengths away from the structure

kpoint = mp.Vector3(-1, 0, 0)
source_center = mp.Vector3(x=coupling_interface_gap/2)
source_size = mp.Vector3(y=waveguide_height + 2 * waveguide_padding)
src = mp.GaussianSource(frequency=fcen, fwidth=fwidth)

source = [
    mp.EigenModeSource(
        src,
        eig_band=1,
        direction=mp.NO_DIRECTION,
        eig_kpoint=kpoint,
        size=source_size,
        center=source_center,
    ),
]
eig_power = source[0].eig_power(fcen)

geometry = [
    mp.Block(
        center=mp.Vector3(), material=waveguide_medium, size=mp.Vector3(Sx, waveguide_height, 0.0)
    ),  # center waveguide
]

sim = mp.Simulation( 
    cell_size=cell_size,
    boundary_layers=pml_layers,
    geometry=geometry,
    sources=source,
    eps_averaging=True,
    resolution=resolution,
    default_material=SiO2,
    # force_complex_fields=True,
)

offset_dft_monitors = []
for px_offset in px_offsets:
    guide_dft = sim.add_dft_fields([mp.Ex, mp.Ey, mp.Ez, mp.Hx, mp.Hy, mp.Hz], fcen, 0, 1, where=mp.Volume(center=-1*source_center + px_offset / resolution, size=source_size))
    offset_dft_monitors.append(guide_dft)

sim.run(until_after_sources=mp.stop_when_energy_decayed(dt=1, decay_by=1e-11))

direction = mp.X
band_num = 1
kpoint = mp.Vector3(-1, 0, 0)

overlap_data = []
for dft_monitor in offset_dft_monitors:
    where = dft_monitor.regions[0].where
    meta = sim.get_array_metadata(dft_cell=dft_monitor)
    Eydft = sim.get_dft_array(dft_monitor, mp.Ey, 0)
    Ezdft = sim.get_dft_array(dft_monitor, mp.Ez, 0)
    wgdata = sim.get_eigenmode(
        fcen,
        direction,
        where,
        band_num,
        kpoint,
    )
    Hympb = [wgdata.amplitude(mp.Vector3(x=x, y=y, z=z), mp.Hy) for x in meta[0] for y in meta[1] for z in meta[2]] 
    Hzmpb = [wgdata.amplitude(mp.Vector3(x=x, y=y, z=z), mp.Hz) for x in meta[0] for y in meta[1] for z in meta[2]] 
    Eympb = [wgdata.amplitude(mp.Vector3(x=x, y=y, z=z), mp.Ey) for x in meta[0] for y in meta[1] for z in meta[2]] 
    Ezmpb = [wgdata.amplitude(mp.Vector3(x=x, y=y, z=z), mp.Ez) for x in meta[0] for y in meta[1] for z in meta[2]] 
    mpbpow = np.abs(np.real(np.sum((np.conj(Eympb) * Hzmpb - np.conj(Ezmpb) * Hympb) * meta[-1])))
    
    mode_overlap = np.sum((Ezdft * (Hympb / np.sqrt(mpbpow)) - Eydft * (Hzmpb / np.sqrt(mpbpow))) * meta[-1])
    mode_power = np.abs(mode_overlap) ** 2
    overlap_data.append(mode_power)

Although these scripts are using 2d simulations, the discrepancy also exists in 3d simulations.

Finally, here is a plotting script to show the discrepancy between the two codes

import matplotlib.pyplot as plt

plt.plot(px_offsets_x, eig_powers, label="get_eigenmode_coefficients")
plt.plot(px_offsets_x, overlap_data, label="manual overlap integral get_eigenmode mode ")
plt.xlabel("fractional shift of a pixel")
plt.ylabel("Power")
plt.legend()

cc @stevengj @oskooi