get_eigenmode_coefficients for degenerate eigenmodes in uniform media

[oskooi]

For a 2d cell in x and y with a non-zero k_z, get_eigenmode seems to be computing an incorrect set of modes. This is demonstrated by the following script.

import meep as mp
import math

resolution = 50        # pixels/μm                                                                                                                                    

dpml = 2               # PML thickness                                                                                                                                
dmat = 10              # length of material                                                                                                                           

sx = dpml+dmat+dpml
sy = 10

cell_size = mp.Vector3(sx,sy,0)
pml_layers = [mp.Absorber(thickness=dpml,direction=mp.X)]

wvl = 0.5              # center wavelength                                                                                                                            
fcen = 1/wvl           # center frequency                                                                                                                             

ng = 1.5
glass = mp.Medium(index=ng)

# out-of-plane (XZ) rotation angle of incident planewave; CCW about Y axis, 0 degrees along in-plane direction                                                        
theta_out = math.radians(20.3)

# in-plane (XY) rotation angle of incident planewave; CCW about Z axis, 0 degrees along +X axis                                                                       
theta_in = math.radians(10.7)

# k (in source medium) with correct length (plane of incidence: XY)                                                                                                   
k = mp.Vector3(math.cos(theta_in)*math.cos(theta_out),math.sin(theta_in)*math.cos(theta_out),math.sin(theta_out)).scale(fcen*ng)

sim = mp.Simulation(resolution=resolution,
                    cell_size=cell_size,
                    boundary_layers=pml_layers,
                    k_point=k,
                    default_material=glass)

sim.init_sim()

mon_pt = mp.Vector3(0.5*sx-dpml-0.2*dmat,0,0)

for num_band in range(1,7):
    sim.get_eigenmode(fcen, mp.X, mp.Volume(center=mon_pt, size=mp.Vector3(0,sy,0)), num_band, k, verbose=True)

The output shows that the dominant wavevectors of the six modes are not unique and instead come in pairs.

Dominant planewave for band 1: (2.813578,0.022404,1.040807)
Dominant planewave for band 2: (2.813578,0.022404,1.040807)
Dominant planewave for band 3: (2.812597,-0.077596,1.040807)
Dominant planewave for band 4: (2.812597,-0.077596,1.040807)
Dominant planewave for band 5: (2.811003,0.122404,1.040807)
Dominant planewave for band 6: (2.811003,0.122404,1.040807)

Even though each pair of modes has the same dominant wavevector, they are in fact different modes. This is demonstrated by the following binary-grating example involving get_eigenmode_coefficients.

import meep as mp
import math
import cmath
import numpy as np

resolution = 50        # pixels/μm                                                                                                                                    

dpml = 1.0             # PML thickness                                                                                                                                
dsub = 1.0             # substrate thickness                                                                                                                          
dpad = 1.0             # length of padding between grating and pml                                                                                                    
gp = 10.0              # grating period                                                                                                                               
gh = 0.5               # grating height                                                                                                                               
gdc = 0.5              # grating duty cycle                                                                                                                           

sx = dpml+dsub+gh+dpad+dpml
sy = gp

cell_size = mp.Vector3(sx,sy,0)

# replace anisotropic PML with isotropic Absorber to attenuate parallel-directed fields of oblique source                                                             
abs_layers = [mp.Absorber(thickness=dpml,direction=mp.X)]

wvl = 0.5              # center wavelength                                                                                                                            
fcen = 1/wvl           # center frequency                                                                                                                             
df = 0.05*fcen         # frequency width                                                                                                                              

ng = 1.5
glass = mp.Medium(index=ng)

# out-of-plane (XZ) rotation angle of incident planewave; CCW about Y axis, 0 degrees along in-plane direction                                                        
theta_out = math.radians(20.3)

# in-plane (XY) rotation angle of incident planewave; CCW about Z axis, 0 degrees along +X axis                                                                       
theta_in = math.radians(10.7)

# k (in source medium) with correct length (plane of incidence: XY)                                                                                                   
k = mp.Vector3(math.cos(theta_in)*math.cos(theta_out),math.sin(theta_in)*math.cos(theta_out),math.sin(theta_out)).scale(fcen*ng)

symmetries = []
eig_parity = mp.NO_PARITY
if theta_out == 0:
  eig_parity = mp.ODD_Z
  if theta_in == 0:
    k = mp.Vector3(0,0,0)
    symmetries = mp.Mirror(mp.Y)
    eig_parity += mp.EVEN_Y

def pw_amp(k,x0):
  def _pw_amp(x):
    return cmath.exp(1j*2*math.pi*k.dot(x+x0))
  return _pw_amp

src_pt = mp.Vector3(-0.5*sx+dpml+0.3*dsub,0,0)
sources = [mp.Source(mp.GaussianSource(fcen,fwidth=df),
                     component=mp.Ez,
                     center=src_pt,
                     size=mp.Vector3(0,sy,0),
                     amp_func=pw_amp(k,src_pt))]

sim = mp.Simulation(resolution=resolution,
                    cell_size=cell_size,
                    boundary_layers=abs_layers,
                    k_point=k,
                    default_material=glass,
                    sources=sources,
                    symmetries=symmetries)

refl_pt = mp.Vector3(-0.5*sx+dpml+0.5*dsub,0,0)
refl_flux = sim.add_flux(fcen, 0, 1, mp.FluxRegion(center=refl_pt, size=mp.Vector3(0,sy,0)))

sim.run(until_after_sources=100)

input_flux = mp.get_fluxes(refl_flux)
input_flux_data = sim.get_flux_data(refl_flux)

sim.reset_meep()

geometry = [mp.Block(material=glass, size=mp.Vector3(dpml+dsub,mp.inf,mp.inf), center=mp.Vector3(-0.5*sx+0.5*(dpml+dsub),0,0)),
            mp.Block(material=glass, size=mp.Vector3(gh,gdc*gp,mp.inf), center=mp.Vector3(-0.5*sx+dpml+dsub+0.5*gh,0,0))]

sim = mp.Simulation(resolution=resolution,
                    cell_size=cell_size,
                    boundary_layers=abs_layers,
                    geometry=geometry,
                    k_point=k,
                    sources=sources,
                    symmetries=symmetries)

refl_flux = sim.add_flux(fcen, 0, 1, mp.FluxRegion(center=refl_pt, size=mp.Vector3(0,sy,0)))
sim.load_minus_flux_data(refl_flux,input_flux_data)

tran_pt = mp.Vector3(0.5*sx-dpml-0.5*dpad,0,0)
tran_flux = sim.add_flux(fcen, 0, 1, mp.FluxRegion(center=tran_pt, size=mp.Vector3(0,sy,0)))

sim.run(until_after_sources=200)

nm_r = np.floor((math.sqrt(math.pow(fcen*ng,2)-math.pow(k.z,2))-k.y)*gp)-np.ceil((-math.sqrt(math.pow(fcen*ng,2)-math.pow(k.z,2))-k.y)*gp) # number of reflected orde\
rs                                                                                                                                                                    
if theta_in == 0 and theta_out == 0:
  nm_r = nm_r*0.5 # since eig_parity removes degeneracy in y-direction                                                                                                
Rsum = 0
for nm in range(int(nm_r)):
  res = sim.get_eigenmode_coefficients(refl_flux, [nm+1], eig_parity=eig_parity)
  r_coeffs = res.alpha
  r_kdom = res.kdom[0]
  Rmode = (abs(r_coeffs[0,0,1])**2)/input_flux[0]
  r_angle = np.sign(r_kdom.y)*math.acos(r_kdom.x/(ng*fcen))
  print("refl:, {}, {:.2f}, {:.8f}".format(nm,math.degrees(r_angle),Rmode))
  Rsum += Rmode

nm_t = np.floor((math.sqrt(math.pow(fcen,2)-math.pow(k.z,2))-k.y)*gp)-np.ceil((-math.sqrt(math.pow(fcen,2)-math.pow(k.z,2))-k.y)*gp) # number of transmitted orders  
if theta_in == 0 and theta_out == 0:
  nm_t = nm_t*0.5 # since eig_parity removes degeneracy in y-direction                                                                                                
Tsum = 0
for nm in range(int(nm_t)):
  res = sim.get_eigenmode_coefficients(tran_flux, [nm+1], eig_parity=eig_parity)
  t_coeffs = res.alpha
  t_kdom = res.kdom[0]
  Tmode = (abs(t_coeffs[0,0,0])**2)/input_flux[0]
  t_angle = np.sign(t_kdom.y)*math.acos(t_kdom.x/fcen)
  print("tran:, {}, {:.2f}, {:.8f}".format(nm,math.degrees(t_angle),Tmode))
  Tsum += Tmode

print("mode-coeff:, {:.6f}, {:.6f}, {:.6f}".format(Rsum,Tsum,Rsum+Tsum))

r_flux = mp.get_fluxes(refl_flux)
t_flux = mp.get_fluxes(tran_flux)
Rflux = -r_flux[0]/input_flux[0]
Tflux =  t_flux[0]/input_flux[0]
print("poynting-flux:, {:.6f}, {:.6f}, {:.6f}".format(Rflux,Tflux,Rflux+Tflux))

The output from this script indicates that each of the duplicated modes produces a different value for the mode coefficients. For example, here are the reflectance values for the first six modes as well as their dominant wavevectors:

refl:, 0, 20.30, 0.00000393
refl:, 1, 20.30, 0.00001048
refl:, 2, -20.36, 0.00001130
refl:, 3, -20.36, 0.00000988
refl:, 4, 20.45, 0.00001030
refl:, 5, 20.45, 0.00001799

Dominant planewave for band 1: (2.813578,0.022404,1.040807)
Dominant planewave for band 2: (2.813578,0.022404,1.040807)
Dominant planewave for band 3: (2.812597,-0.077596,1.040807)
Dominant planewave for band 4: (2.812597,-0.077596,1.040807)
Dominant planewave for band 5: (2.811003,0.122404,1.040807)
Dominant planewave for band 6: (2.811003,0.122404,1.040807)

Note that the k_z value for each mode's dominant wavevector is constant and equivalent to that of the k_point (following the bug fix in #602).

Since each pair of modes has the same dominant wavevector but produces different values for the mode coefficients, this likely suggests that the parity of each of the two modes is different (i.e., EVEN_Z and ODD_Z). This would be a bug since eig_parity=mp.NO_PARITY in this example. Unfortunately, the parity of the mode cannot be displayed as part of the output in get_eigenmode when the verbose flag is set to True.

As a result of the incorrect modes computed by get_eigenmode, the sum of the reflectance and transmittance values for all these modes does not equal the Poynting flux values:

mode-coeff:, 0.022270, 1.305896, 1.328166
poynting-flux:, 0.056975, 0.942627, 0.999602

[stevengj]

Just because you don't specify a parity doesn't mean that modes have no polarization — on the contrary, NO_PARITY means to compute modes with both possible polarizations. The two polarizations can certainly have different mode coefficients.

"z parity", however, is a meaningless way to talk about polarization in this case because kz breaks the z=0 symmetry: the solutions aren't going to be even/odd in z.

(As usual, however, there is the problem that if you have degeneracies between the two polarizations, then MPB will give two random superpositions.)

[oskooi]

eig_parity=mp.NO_PARITY was used in the script because k_z is non-zero. The modes obviously do not have any symmetry in z. Then the question is: what else could explain why 2 modes with the same dominant wavevector produce different values for the mode coefficient?

[stevengj]

what else could explain why 2 modes with the same dominant wavevector produce different values for the mode coefficient?

The grating structure breaks the rotational symmetry, so it is correct that it couples differently into two polarizations for the same k.

the sum of the reflectance and transmittance values for all these modes does not equal the Poynting flux values:

I think what is happening here is that, for degenerate modes, MPB is computing a pair of modes that is not power-orthogonal. (It computes modes that are orthogonal in the sense that ∫H₁⋅H₂ = 0, not orthogonal in the Poynting-vector sense.) So, there is a cross term that we are not including when we sum up the mode coefficients squared.

[stevengj]

(Note that this only affects the mode calculation in a uniform medium, where there are degeneracies. For something like a waveguide, the modes are non-degenerate and there should be no problem.)