Use kpoint_func as initial k_point in fields::get_eigenmode

[oskooi]

As described in the user manual, the parameter kpoint_func of get_eigenmode_coefficients "supplies a rough initial guess for the $\vec{k}$ of band number n at frequency $f=\omega/2\pi$." The manual also states that "If direction is set to mp.NO_DIRECTION, then kpoint_func is not only the initial guess and the search direction of the $\vec{k}$ vectors, but is also taken to be the direction of the waveguide."

This PR fixes a bug whereby if direction = mp.NO_DIRECTION then the initial $\vec{k}$ which had been set to kpoint_func via the argument _kpoint of fields::get_eigenmode was actually being overwritten by the components of k_point of the Simulation object in the periodic directions of the cell for which the mode monitor extends to the cell edge.

The fix ensures that the initial $\vec{k}$ is set to kpoint_func when direction = mp.NO_DIRECTION consistent with what is described in the manual.

[oskooi]

I have added a unit test which is currently failing on master. The test is actually based on my tutorial presentation at MeepCon and involves verifying that the sum of all the reflectance and transmittance of all the diffracted orders of a binary grating are equivalent to one (i.e., conservation of energy).

[stevengj]

This makes no sense to me, because the kpoint corresponding to the boundary conditions is equivalent in the mode solver to the kpoint you are passing, because they are in different Brillouin zones.

In particular, you are claiming that ω₁(kx, 0) ≠ ω₁(kx, integer × 2π/period), or that the mode fields are different, which violates Bloch's theorem. (And we've seen this to work many times in MPB.) So, I don't understand why this PR would fix your test.

As a simple test, you could set up the same unit cell in MPB (an a × no-size × no-size) cell, and check that ω₁(kx, 0) = ω₁(kx, integer × 2π/period) for the kx of your diffracted order.

[mawc2019]

Suppose the periodic direction is the $y$ direction and the incident light is propagating towards the $x$ direction. The k_point in Simulation, which corresponds to the boundary condition, is $(0,0,0)$. Suppose the diffracted light has the wave vector $(\sqrt {f^2-(m/\Lambda)^2},m/\Lambda,0)$, where $f$ is the frequency, $m$ is the order of diffraction, and $\Lambda$ is the periodicity in the $y$ direction. This wave vector is kpoint_func in get_eigenmode_coefficients.

the kpoint corresponding to the boundary conditions is equivalent in the mode solver to the kpoint you are passing, because they are in different Brillouin zones.

Do you mean $(0,0,0)$ and $(\sqrt {f^2-(m/\Lambda)^2},m/\Lambda,0)$ are equivalent? If so, can I ask why? They do not seem equivalent because of their $x$ components.

[stevengj]

The confusion here is that that Bloch wavevector (what you pass to MPB) is not the same thing as the planewave wavevector for a finite period (discrete translational symmetry ≠ continuous translational symmetry).

A quick way to make them the same might be to use the eig_vol parameter to pass an eigenmode volume of 0 in the y direction (or maybe just 1e-20 … not sure MPB supports it being actually zero). In the limit where the period is zero (i.e. continuous translational symmetry), the Bloch wavevector becomes the same thing as the familiar planewave wavevector.

I think it should work as-is with an eig_vol that is smaller than the source volume, but I don't recall ever actually trying this case.

[oskooi]

The confusion here is that that Bloch wavevector (what you pass to MPB) is not the same thing as the planewave wavevector for a finite period (discrete translational symmetry ≠ continuous translational symmetry). A quick way to make them the same might be to use the eig_vol parameter to pass an eigenmode volume of 0 in the $y$ direction

This approach works when the eigenmode volume is a single pixel (1 / resolution) in the periodic $y$ direction. A value of 0 or 1e-20 produces the incorrect result. The unit test in this PR has therefore been updated with this change along with a comment to explain why it is necessary.

[codecov-commenter]

Codecov Report

All modified and coverable lines are covered by tests ✅

Comparison is base (5fa12b4) 73.23% compared to head (7fc28ad) 74.08%.
Report is 144 commits behind head on master.

❗ Your organization needs to install the Codecov GitHub app to enable full functionality.

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #2285      +/-   ##
==========================================
+ Coverage   73.23%   74.08%   +0.85%     
==========================================
  Files          17       18       +1     
  Lines        4931     5399     +468     
==========================================
+ Hits         3611     4000     +389     
- Misses       1320     1399      +79     

see 14 files with indirect coverage changes

[stevengj]

This approach works when the eigenmode volume is a single pixel (1 / resolution) in the periodic direction. A value of 0 or 1e-20 produces the incorrect result.

Maybe it doesn't know how to interpolate with an eigenmode volume that is smaller than a single pixel? This can probably be fixed?

[oskooi]

Turns out I just needed to use a value for the size of the eig_vol in the periodic direction which is larger than machine epsilon to make this work as expected. A value of 1e-20 was too small. I therefore chose 1e-7 which is sufficient for single precision.

Note that a value of 0 produced this MPB error:

CHECK failure on line 393 of eigensolver.c: crazy number detected in trace!!