Eigenmode overlaps from get_eigenmode_coefficients() vary with k-direction

[ianwilliamson]

When using get_eigenmode_coefficients() with direction=mp.NO_DIRECTION and manually specifying the direction via kpoint_func, the returned mode overlap values change when flipping the sign of the k-point. I believe this relates to what was observed in #1593, because the solution there was to change the k-point sign based on the direction of the requested mode overlap.

Consider the following simulation setup:

frequencies = [0.45, 0.5, 0.55]
wg = mp.Block(
    center=mp.Vector3(0, 0, 0),
    material=mp.Medium(index=3.4),
    size=mp.Vector3(10, 0.5, 0),
)
source = mp.EigenModeSource(
    mp.GaussianSource(frequency=onp.mean(frequencies), fwidth=0.2),
    eig_band=1,
    direction=mp.NO_DIRECTION,
    eig_kpoint=mp.Vector3(x=1),
    size=mp.Vector3(0, 1),
    center=mp.Vector3(),
)
simulation = mp.Simulation(
    cell_size=mp.Vector3(3., 3.),
    boundary_layers=[mp.PML(1.)],
    geometry=[wg],
    sources=[source],
    resolution=20,
)
dft_mon = simulation.add_dft_fields(
    [mp.Ez],
    frequencies,
    where=mp.Volume(size=simulation.cell_size),
)
mon1 = simulation.add_mode_monitor(
    frequencies,
    mp.ModeRegion(center=mp.Vector3(+0.4), size=mp.Vector3(0, 1)),
    yee_grid=True,
)
simulation.run(until=500)
simulation.plot2D()

This produces the following DFT fields:

The output of (automatic direction):

data = simulation.get_eigenmode_coefficients(
    mon1,
    [1],
    direction=mp.AUTOMATIC,
)
print(onp.abs(data.alpha))

is:

[[[1.33664502 0.00561725]
  [5.14267297 0.02121229]
  [1.65917472 0.00646342]]]

The output of (no direction with positive k-point):

data = simulation.get_eigenmode_coefficients(
    mon1,
    [1],
    direction=mp.NO_DIRECTION,
    kpoint_func=lambda *not_used: mp.Vector3(1,0,0),
)
print(onp.abs(data.alpha))

is:

[[[1.33664502 0.00561726]
  [5.14267296 0.02121239]
  [1.65917471 0.00646343]]]

The output of (no direction with negative k-point):

data = simulation.get_eigenmode_coefficients(
    mon1,
    [1],
    direction=mp.NO_DIRECTION,
    kpoint_func=lambda *not_used: mp.Vector3(-1,0,0),
)
print(onp.abs(data.alpha)

is:

[[[0.0535472  1.27748126]
  [0.26983257 4.85162955]
  [0.11030892 1.54240248]]]

So, the automatic direction and the positive k-point coefficient values are consistent but the negative k-point coefficient values are slightly different. I would have thought that the negative k-point coefficient values are the same, but with the fwd / bkwd array indexes swapped.

[ianwilliamson]

Is what I am seeing in the results above expected or does it appear to be a bug? If it is not a bug, which mode overlap values are considered to be correct?

[oskooi]

I looked into what could be causing the discrepancy in the mode coefficients in the forward- and backward-propagating in the example above for a single frequency (0.5) and noticed two things.

1. When get_eigenmode_coefficients is passed either kpoint_func=lambda *not_used: mp.Vector3(1,0,0) or kpoint_func=lambda *not_used: mp.Vector3(-1,0,0), the mode computed by MPB has the same dominant wavevector with different sign and group velocity with same sign. The following is the output from MPB using meep.verbosity(2).

(a) kpoint_func=lambda *not_used: mp.Vector3(1,0,0)

Newton step: group velocity v=0.283361, kmatch=1.53045
Dominant planewave for band 1: (1.530449,-0.000000,0.000000)

(b) kpoint_func=lambda *not_used: mp.Vector3(-1,0,0)

Newton step: group velocity v=0.283361, kmatch=1.53045
Dominant planewave for band 1: (-1.530449,0.000000,0.000000)

The problem therefore in (b) is that MPB is computing the wrong mode. The group velocity should be negative rather than positive.

2. When I printed out the actual values for the mode coefficients cplus and cminus from the following lines, the results were different for the two cases. This is expected because the mode profiles are different.

meep/src/mpb.cpp

Lines 945 to 946 in fa23466

	complex<double> cplus = 0.5 * (mode_flux[0] + mode_flux[1]);
	complex<double> cminus = 0.5 * (mode_flux[0] - mode_flux[1]);

(a) kpoint_func=lambda *not_used: mp.Vector3(1,0,0)

cplus=7.66864+6.51408i, cminus=-0.01000-0.00830i

(b) kpoint_func=lambda *not_used: mp.Vector3(-1,0,0)

cplus=0.10096+0.08594i, cminus=7.55769+6.41984i

We need to first figure out why MPB is not computing the correct backward-propagating mode.

[stevengj]

Realize that MPB is computing the group velocity relative to the direction of k, so if k is in the backward direction, and the mode is backward propagating, then it will return a group velocity that is positive.

[oskooi]

I think I have tracked down what could be the cause of the unexpected results reported here (i.e., why the mode coefficients of the backward mode are not identical to those of the forward mode with the order just swapped). It has to do with the phase of the eigenmodes returned by MPB via get_eigenmode which varies with the sign of the user-specified input parameter kpoint:

meep/src/mpb.cpp

Lines 912 to 916 in d12e544

	if (user_kpoint_func) kpoint = user_kpoint_func(flux.freq[nf], band_num, user_kpoint_data);
	am_now_working_on(MPBTime);
	void *mode_data =
	get_eigenmode(flux.freq[nf], d, flux.where, eig_vol, band_num, kpoint, match_frequency,
	parity, eig_resolution, eigensolver_tol, kdom, (void **)&mdata, dp);

For the 2d example posted above, I plotted the eigenmode profile for the H_y and E_z fields computed by MPB at the DFT line monitor and discovered that while the H_y field is identical for the forward (kpoint = (+1,0,0)) and backward (kpoint = (-1,0,0)) modes, E_z for the backward mode is 180° out of phase relative to the forward mode (i.e., there is a sign flip of -1). This phase shift for a single field component is why the resulting overlap integral and thus the mode coefficients are different for the forward and backward modes.

The fix therefore involves ensuring that an identical eigenmode profile is used to compute the mode coefficients for the forward and backward modes.

[stevengj]

That's correct — for the backward mode, you have to flip the relative phase of Hy and Ez, in order to flip the sign of the Poynting vector.

(Note that the transverse fields Hy and Ez can be chosen to be purely real. MPB might be multiplying both by an arbitrary phase, but if one is real then at first glance it seems that the other should be real too.)

[oskooi]

For some reason, when plotting the mode profile of the actual eigenmode from MPB interpolated into the Meep grid and used to compute the overlap integral from this line:

meep/src/dft.cpp

Line 945 in 434f3fb

if (mode1_data) mode1val = eigenmode_amplitude(mode1_data, loc, S.transform(c_conjugate, sn));

The E_z field is complex but H_y is purely real as shown in the following figures for the real and imaginary parts of each field component separately. This should not affect the results but it is unexpected.

Note: the E_z fields of the "forward" (kpoint=(+1,0,0)) and "backward" (kpoint=(-1,0,0)) modes have a sign flip which is expected because the Poynting flux is in opposite directions.

[oskooi]

It seems there is a bug when computing the mode coefficients which involve adding two different cross products conj[E(meep)] x H(mpb) and conj[E(mpb)] x H(meep) on these lines:

meep/src/mpb.cpp

Lines 945 to 946 in 849bf3d

	complex<double> cplus = 0.5 * (mode_flux[0] + mode_flux[1]);
	complex<double> cminus = 0.5 * (mode_flux[0] - mode_flux[1]);

The problem is that E(meep) and H(meep) are located at different Yee grid points and therefore the cross products must each be computed at the same centered grid point before they can be added/subtracted together in the sum/difference.

[stevengj]

Even if there is a bug there (and you would need to check if things are being interpolated elsewhere), that would only be a first-order error and would go away with increasing resolution, which is not the problem you're seeing here.

[ianwilliamson]

The problem is that E(meep) and H(meep) are located at different Yee grid points and therefore the cross products must each be computed at the same centered grid point before they can be added/subtracted together in the sum/difference.

Whether or not this happens is determined by yee_grid=False or yee_grid=True on add_mode_monitor(), no?

[oskooi]

It looks like this issue is related to #1470 and only when the DFT mode monitors used to compute the mode coefficients are defined on the centered grid, the error converges to zero with increasing resolution (see right plot of first figure below). This is not the case when the mode coefficients are computed using a DFT mode monitor on the Yee grid which is the case for the EigenmodeCoefficient class of the adjoint solver as reported in #1793. In this case, the error in the backward-propagating mode seems to be a constant independent of the resolution (see right plot of second figure below).

To demonstrate this, I set up a test case based on a 2d waveguide from the user manual and computed |S₁₁|² and |S₂₁|² using four different cases which are a combination of: (1) yee_grid parameter of add_mode_monitor is either True or False and (2) kpoint_func parameter of get_eigenmode_coefficients is either (+1,0,0) or (-1,0,0). The mode coefficients are computed at resolutions of 25, 50, 100, 200, and 400. The simulation script is here.

The following two figures show these results comparing the centered vs. Yee grid.

It therefore seems that there is nothing fundamentally wrong with the computation of the mode overlap integrals as long as the mode monitor is defined on the centered grid. The fix therefore involves simply changing the EigenmodeCoefficient class to use the centered grid rather than the Yee grid for the mode monitors.

[ianwilliamson]

It looks like this issue is related to #1470

I am confused as to how what we're observing here relates to #1470 and reciprocity, since the location of the sources and monitors are not being changed in this setup. It's the same monitor and the same simulation but the result of the overlap calculation is changing. This seems surprising to me, even at a finite simulation resolution.

[oskooi]

The connection to #1470 relates to the discretization error introduced by the interpolation of the eigenmode from MPB onto the Meep grid when computing the overlap integral for the mode coefficients. #1470 and the results shown above demonstrate that this error converges to zero in the limit of infinite resolution.

Note from the results in the right plot of the first figure above that when using the centered grid, |S₁₁|² is nearly identical (and converging to zero with increasing resolution) when using a kpoint_func of either (+1,0,0) or (-1,0,0). Also note that the example posted at the top of this issue is using yee_grid=True. If you switch this to yee_grid=False, you will find the same behavior that as you increase the resolution |S₁₁|² is the same when kpoint_func is either (+1,0,0) or (-1,0,0).

[ianwilliamson]

when using the centered grid, |S11|2 is nearly identical

The value is nearly identical but, from my calculation, the relative difference between the S11 for the +/-k cases is actually consistent with the difference between S21 for the +/-k cases. It looks like the value varies by more than 10% for a resolution of 20. Here we're considering a straight waveguide so what we're calling S11 will be much smaller than S21, but this won't be the case in general.

I would have thought that changing the sign of the k-point just flips the coordinate system of MPB. It's surprising that you can get significantly different values depending on how you choose to specify your direction. I understand that the difference goes away at infinite resolution, but at a finite resolution it's unclear what the preferred approach / value is from a user POV.

[oskooi]

For accurate computations of the S-parameters, note that you can always do a separate normalization run to compute the incident DFT fields from just the sources and then subtract these from the DFT fields of the "real" run. This approach is described in Section 4.3 "Separating Incident and Scattered Fields" of our book chapter: https://arxiv.org/pdf/1301.5366.pdf.

[oskooi]

The relative difference in |S₂₁|² and |S₁₁|² for +k and -k is actually very similar which is reassuring. This is shown in the following figures. In the figure above, the seemingly large difference in |S₂₁|² relative to |S₁₁|² is actually just an artifact of the different magnitudes: in the limit of infinite resolution, |S₂₁|² is ~1100 whereas |S₁₁|² is 0 (1e-4); a difference of about 7 orders. The choice of the loglog axes obscures the difference in magnitude. For this geometry, the relative difference in +k vs. -k (second figure below) is less than 1% at the smallest resolutions. Nevertheless, the relative difference seems to be geometry dependent and just happens to be small in this case.

[oskooi]

Closing this issue as it is not a bug but rather a feature of computing mode coefficients on the Yee grid vs. centered grid for add_mode_monitor.