mismatched frequency between eigenmode and Gaussian pulse launches wrong mode

[oskooi]

A 2d EigenModeSource in a 2d cell does not seem to correctly support y mirror symmetry. This is demonstrated using the parallel waveguide example from TutorialOptical Forces. In this example, eig_parity=mp.ODD_Y is part of the eigenmode source definition.

The E_y field profile is correct when there are no symmetries and when symmetries = [mp.Mirror(mp.X, phase=+1)]. However, whenever the symmetries object includes mp.Mirror(mp.Y, phase=-1), the fields are incorrect.

This does not seem to be related to #1208.

import meep as mp
import numpy as np
import matplotlib
matplotlib.use('agg')
import matplotlib.pyplot as plt

resolution = 30   # pixels/μm                                                                                                                              

Si = mp.Medium(index=3.45)

dpml = 1.0
pml_layers = [mp.PML(dpml)]

sx = 5
sy = 3
cell = mp.Vector3(sx+2*dpml,sy+2*dpml,0)

a = 1.0     # waveguide width/height                                                                                                                       
s = 0.5     # waveguide separation                                                                                                                         

k_point = mp.Vector3(z=0.5)

fcen = 0.2

geometry = [mp.Block(center=mp.Vector3(-0.5*(s+a)),
                     size=mp.Vector3(a,a,mp.inf),
                     material=Si),
            mp.Block(center=mp.Vector3(0.5*(s+a)),
                     size=mp.Vector3(a,a,mp.inf),
                     material=Si)]

symmetries = [mp.Mirror(mp.X, phase=+1),                                                                                                                 
              mp.Mirror(mp.Y, phase=-1)]                                                                                                                 

eig_sources = [mp.EigenModeSource(src=mp.GaussianSource(fcen, fwidth=0.1*fcen),
                                  size=mp.Vector3(sx,sy),
                                  center=mp.Vector3(),
                                  eig_band=1,
                                  eig_kpoint=k_point,
                                  eig_match_freq=False,
                                  eig_parity=mp.ODD_Y)]

sim = mp.Simulation(resolution=resolution,
                    cell_size=cell,
                    geometry=geometry,
                    boundary_layers=pml_layers,
                    sources=eig_sources,
                    symmetries=symmetries,
                    k_point=k_point)

sim.run(until_after_sources=1000)

sim.plot2D(fields=mp.Ey)
plt.savefig("parallel_ey.png")

no symmetries

symmetries = [mp.Mirror(mp.X, phase=+1)]

symmetries = [mp.Mirror(mp.Y, phase=-1)]

symmetries = [mp.Mirror(mp.X, phase=+1),mp.Mirror(mp.Y, phase=-1)]

[stevengj]

Is it finding the same mode? (see the mode frequency output in the simulation output)

[oskooi]

Yes, MPB is finding the same mode in all four cases. Does this mean therefore that the bug is in the way that the y mirror symmetry is implemented in a 2d cell with non-zero kz?

no symmetries

Initializing structure...
Splitting into 2 chunks evenly
time for choose_chunkdivision = 0.00045725 s
Working in 2D dimensions.
Computational cell is 7 x 5 x 0 with resolution 30
     block, center = (-0.75,0,0)
          size (1,1,1e+20)
          axes (1,0,0), (0,1,0), (0,0,1)
          dielectric constant epsilon diagonal = (11.9025,11.9025,11.9025)
     block, center = (0.75,0,0)
          size (1,1,1e+20)
          axes (1,0,0), (0,1,0), (0,0,1)
          dielectric constant epsilon diagonal = (11.9025,11.9025,11.9025)
time for set_epsilon = 0.0224114 s
-----------
Meep: using complex fields.
    iteration   15: trace = 0.05034806646248322 (0.00114643% change)
    iteration   29: trace = 0.05034785978091706 (2.79759e-09% change)
MPB solved for frequency_1(0,0,0.5) = 0.224383 after 31 iters

symmetries = [mp.Mirror(mp.X, phase=+1)]

Initializing structure...
Halving computational cell along direction x
Splitting into 2 chunks evenly
time for choose_chunkdivision = 0.000237776 s
Working in 2D dimensions.
Computational cell is 7 x 5 x 0 with resolution 30
     block, center = (-0.75,0,0)
          size (1,1,1e+20)
          axes (1,0,0), (0,1,0), (0,0,1)
          dielectric constant epsilon diagonal = (11.9025,11.9025,11.9025)
     block, center = (0.75,0,0)
          size (1,1,1e+20)
          axes (1,0,0), (0,1,0), (0,0,1)
          dielectric constant epsilon diagonal = (11.9025,11.9025,11.9025)
time for set_epsilon = 0.010961 s
-----------
Meep: using complex fields.
    iteration   13: trace = 0.05035031580867293 (0.0238154% change)
    iteration   27: trace = 0.05034785978997125 (6.42218e-08% change)
MPB solved for frequency_1(0,0,0.5) = 0.224383 after 31 iters

symmetries = [mp.Mirror(mp.Y, phase=-1)]

Initializing structure...
Halving computational cell along direction y
Splitting into 2 chunks evenly
time for choose_chunkdivision = 0.000313096 s
Working in 2D dimensions.
Computational cell is 7 x 5 x 0 with resolution 30
     block, center = (-0.75,0,0)
          size (1,1,1e+20)
          axes (1,0,0), (0,1,0), (0,0,1)
          dielectric constant epsilon diagonal = (11.9025,11.9025,11.9025)
     block, center = (0.75,0,0)
          size (1,1,1e+20)
          axes (1,0,0), (0,1,0), (0,0,1)
          dielectric constant epsilon diagonal = (11.9025,11.9025,11.9025)
time for set_epsilon = 0.0116178 s
-----------
Meep: using complex fields.
    iteration   15: trace = 0.05034806646248328 (0.00114643% change)
    iteration   30: trace = 0.05034785978082769 (1.77621e-10% change)
MPB solved for frequency_1(0,0,0.5) = 0.224383 after 31 iters

symmetries = [mp.Mirror(mp.X, phase=+1),mp.Mirror(mp.Y, phase=-1)]

Initializing structure...
Halving computational cell along direction x
Halving computational cell along direction y
Splitting into 2 chunks evenly
time for choose_chunkdivision = 0.000155622 s
Working in 2D dimensions.
Computational cell is 7 x 5 x 0 with resolution 30
     block, center = (-0.75,0,0)
          size (1,1,1e+20)
          axes (1,0,0), (0,1,0), (0,0,1)
          dielectric constant epsilon diagonal = (11.9025,11.9025,11.9025)
     block, center = (0.75,0,0)
          size (1,1,1e+20)
          axes (1,0,0), (0,1,0), (0,0,1)
          dielectric constant epsilon diagonal = (11.9025,11.9025,11.9025)
time for set_epsilon = 0.00448635 s
-----------
Meep: using complex fields.
    iteration   14: trace = 0.0503486436687603 (0.00332107% change)
    iteration   27: trace = 0.05034785978997258 (6.42228e-08% change)
MPB solved for frequency_1(0,0,0.5) = 0.224383 after 31 iters

[stevengj]

A CW simulation might be clearer to visualize what mode is being launched?

[oskooi]

I think I've tracked down the cause of the bug: when the center frequency of the GaussianSource does not match the frequency of the eigenmode (computed using eig_match_freq=False), the mode that is launched is incorrect. This does not seem to be related to the y mirror symmetry.

In the example posted above, the frequency of the eigenmode computed by MPB is 0.224383. However, the center frequency of the GaussianSource is 0.2 (fcen) and frequency bandwidth (fwidth) is 0.02 which does not overlap 0.224383.

When fcen is changed to 0.22, the resulting mode is correct:

symmetries = [mp.Mirror(mp.Y, phase=-1)]

symmetries = [mp.Mirror(mp.X, phase=+1),mp.Mirror(mp.Y, phase=-1)]

The issue here is that a user launching an eigenmode via eig_match_freq=False does not know apriori what the frequency of the mode is and if they just arbitrarily specify a center frequency for a Gaussian pulse time profile which does not match the eigenmode frequency then the mode that is launched is incorrect. The correct thing to do would be for the eigenmode frequency to override the center frequency of the pulse (which is simply ignored in this case) and the fwidth parameter should be defined as e.g. 10% of the eigenmode frequency (fixed default).

[stevengj]

It looks like it was supposed to do that already, but there may be a bug. Try:

diff --git a/src/mpb.cpp b/src/mpb.cpp
index cd507e25..66ada5ab 100644
--- a/src/mpb.cpp
+++ b/src/mpb.cpp
@@ -696,7 +696,7 @@ void fields::add_eigenmode_source(component c0, const src_time &src, direction d
   global_eigenmode_data->amp_func = A ? A : default_amp_func;
 
   src_time *src_mpb = src.clone();
-  if (!match_frequency) src_mpb->set_frequency(frequency);
+  if (!match_frequency) src_mpb->set_frequency(global_eigenmode_data->frequency);
 
   /*--------------------------------------------------------------*/
   // step 2: add sources whose radiated field reproduces the      */

[stevengj]

In particular, it looks like it used to work, but was broken by #192.

[stevengj]

Note that this will not only work for a gaussian, but also for a CW source — it will set the frequency to the eigenmode frequency (but won't change the width or any other parameter of the source).

[oskooi]

Tested the patch and it looks like it fixed the bug (i.e., it now does not matter what the center frequency of the GaussianSource is, the value is overwritten by the eigenmode frequency).

[stevengj]

Note that, in get_eigenmode, if match_frequency=false or the k point is ≠ 0 in the given direction then the frequency is ignored (it is only used to for the initial k guess in the Newton search).

[stevengj]

(It won't set the frequency for a Python CustomSource or a Scheme custom-src, however.)