missing unit test for phase of mode coefficients

[oskooi]

The existing unit tests for get_eigenmode_coefficients including test_binary_grating.py, test_mode_coeffs.py, test_mode_decomposition.py, test_special_kz.py, and others all test purely real quantities such as the Poynting flux, e.g. $|S_{11}|^2$, etc., which involves the magnitude of the complex mode coefficient. Currently, there are no tests for just the imaginary part of the mode coefficients (or alternatively the real and imaginary parts tested separately).

A possible test could involve comparing the complex reflection and transmission coefficients of a 1d multilayer stack from Meep (2d and 3d simulations) to that computed analytically using a transfer-matrix approach. (We do not have to write our own TMM solver as part of this test. Instead, we can use hard-coded results obtained using a third-party package such as https://github.com/sbyrnes321/tmm.)

[smartalecH]

all test purely real quantities

Yes because FOMs must be purely real, right?

Doesn't the power of a scattering coefficient implicitly test the imaginary part too?

[smartalecH]

(Ah for some reason I thought this was wrt the gradients of the adjoint solver. But you're speaking generally, right?)

[oskooi]

I am referring to the mode-decomposition feature in general, not specific to its use as part of the objective functions of the adjoint solver.

The tutorial Focusing Properties of a Metasurface Lens is the only example thus far where we actually make use of the complex mode coefficients (i.e., to design a lens using the quadratic formula for the local phase). I suppose we could adapt this example as the unit test but it seems overly complicated (i.e., it uses the near-to-far field transformation feature, etc.).

[stevengj]

If you just use a single interface, you can use the Fresnel formula. If it's just the flat interface, even for a 2d or 3d simulation you can make it 1 pixel wide.

Basically, you would compute:

1. the eigenmode coefficient of (e.g.) the transmitted field at several wavelengths
2. the eigenmode coefficient of the incident field at those wavelengths
3. divide them to get a complex transmission coefficient t(ω)
4. compare t(ω) / exp(i phase(t(ω₀))) to the Fresnel formula, where ω₀ is one of your frequencies. That is, you normalize the phase at ω₀ to be 0 (and do the same for the Fresnel formula). That takes away the arbitrary phase factor that comes e.g. from where you place the eigenmode monitors and MPB's convention.

(Otherwise there is an arbitrary overall phase factor that depends on your conventions; it's not interesting or particularly meaningful to try to match the overallphase convention to that of the Fresnel formulas, particularly since we don't even document the phase choice imposed by MPB.)

[oskooi]

The reference analytic result which we will be validating the simulation results against — the complex Fresnel coefficients for a flat interface of two lossless materials with index $n_1$ and $n_2$ given an incident planewave at angle $\theta$ — are shown in slide 2 of these notes:

Note: in this coordinate convention, TE corresponds to the $E$-field polarized normal to the plane of incidence (i.e., $S$ polarization). TM is the $E$-field parallel to the plane of incidence (i.e., $P$ polarization).

[oskooi]

From the Fresnel equations (shown above), a complex mode coefficient requires $n^2 - \sin^2\theta < 0$. This occurs for total internal reflection (TIR) whereby $n = \frac{n_2}{n_1} < 1$ or equivalently $n_1 > n2$. In this case, $|r{TE,TM}|=1$ which means $r{TE,TM}=|r{TE,TM}|e^{i\phi}=e^{i\phi}$. The phase $\phi$ of the reflected wave spans the $2\pi$ range of $[-\pi,+\pi]$. This topic is reviewed in slides 13, 17-20 of the notes.

For the unit test, we can use the same parameters as the notes: $n_1=1.5$ and $n_2=1.0$ for which the critical angle for TIR is 41.8°. The key result to demonstrate in the unit test are several points on the phase plot: relative phase shift vs. angle of incidence of the planewave for both polarizations. These results are summarized on slide 20:

(On slides 21 and 22, the notes review the operation of the Fresnel rhomb. This could be a useful tutorial demonstration.)

[oskooi]

Here is an initial attempt to validate the complex reflectance coefficient computed using the mode-decomposition feature in Meep with the Fresnel equations. The procedure for computing the reflectance coefficient in Meep is based on steps 1-3 from @stevengj's comment above.

The setup is 2d (cell size is arbitrary) and consists of a flat interface of $n_1=1.5$ and $n_2=1.0$. There are two test cases. A planewave with $S$ polarization ($E$-field pointing out of the $xy$ plane of incidence) is incident at two different angles: (1) 50° and (2) 61°. These correspond to total internal reflected (TIR) modes.

A comparison of the Meep and Fresnel results are shown in the tables below. Note that the magnitude of the reflectance coefficient (second column) is one in all cases. This is expected because of the TIR mode. However, there is an obvious discrepancy in the actual phase. Since the Meep results are converged with resolution, perhaps there is a bug in the calculation of the reflection coefficient itself?

θ = 50°	reflectance coefficient	phase (radians)
Meep	-0.56759-0.82108j	-2.17565
Fresnel	0.48743-0.87316j	-1.06165

θ = 61°	reflectance coefficient	phase (radians)
Meep	-0.65889+0.75084j	2.29106
Fresnel	-0.15385-0.98809j	-1.72526

'''                                                                                                               
Computes the complex reflectance coefficient for a flat interface                                                 
of two lossless materials given a planewave at oblique incidence.                                                 

ref: https://www.patarnott.com/atms749/pdf/FresnelEquations.pdf                                                   
'''
import cmath
import math

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

def mode_coeff(dft_mon: mp.DftFlux) -> np.ndarray:
    """Computes the complex mode coefficients.                                                                    

       Args:                                                                                                      
          dft_mon: the DFT flux monitor object.                                                                   

       Returns:                                                                                                   
          The mode coefficients as a Numpy array of (bands, frequencies,                                          
            forward/backward).                                                                                    
    """
    res = sim.get_eigenmode_coefficients(
        dft_mon,
        bands=[1],
        eig_parity=mp.ODD_Z,
    kpoint_func=lambda *not_used: k,
        direction=mp.NO_DIRECTION,
    )
    return res.alpha

def refl_coeff_te(th: float, n1: float, n2: float) -> float:
    """Computes the reflectance coefficient for TE via the Fresnel equations.                                     

       Args:                                                                                                      
          th: angle of incident planewave (radians).                                                              
          n1: refractive index of source medium.                                                                  
          n2: refractive index of medium adjacent to source medium.                                               

       Returns:                                                                                                   
          The reflectance coefficient (a complex number).                                                         
    """
    return ((math.cos(th) - ((n2/n1)**2 - math.sin(th)**2)**0.5) /
            (math.cos(th) + ((n2/n1)**2 - math.sin(th)**2)**0.5))

resolution = 100.

sx = 10.
sy = 5.
dpml = 1.
cell_size = mp.Vector3(sx + 2*dpml, sy, 0)
pml_layers = [mp.PML(dpml, direction=mp.X)]

n1 = 1.5
n2 = 1.0

# angle of incident planewave at center frequency                                                                 
# 0° is +x; rotated CCW about +z axis                                                                             
theta = np.radians(50.)

fcen = 1. # center frequency                                                                                      
df = 0.1 * fcen
cutoff = 10.0

# k (in source medium) with correct length                                                                        
# plane of incidence is xy                                                                                        
k = mp.Vector3(n1*fcen, 0, 0).rotate(mp.Vector3(0, 0, 1), theta)
print("k:, ({:.6f}, {:.6f}, {:.6f})".format(k.x, k.y, k.z))

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, 0, 0)
sources = [
    mp.Source(
        mp.GaussianSource(fcen, fwidth=df, cutoff=cutoff),
        component=mp.Ez, # S polarization                                                                         
        center=src_pt,
        size=mp.Vector3(0, cell_size.y, 0),
        amp_func=pw_amp(k, src_pt),
    ),
]

sim = mp.Simulation(
    resolution=resolution,
    cell_size=cell_size,
    default_material=mp.Medium(index=n1),
    boundary_layers=pml_layers,
    k_point=k,
    sources=sources,
)

# DFT monitor for reflected fields                                                                                
refl_pt = mp.Vector3(-0.25*sx, 0, 0)
fluxreg = mp.FluxRegion(center=refl_pt, size=mp.Vector3(0, cell_size.y, 0))
refl_mon = sim.add_mode_monitor(fcen, 0, 1, fluxreg)

tol = 1e-8
sim.run(
    until_after_sources=mp.stop_when_fields_decayed(
        50,
        mp.Ez,
        src_pt,
        tol,
    ),
)

input_mode_coeff = mode_coeff(refl_mon)[0, 0, 0]
input_flux_data = sim.get_flux_data(refl_mon)

sim.reset_meep()

geometry = [
    mp.Block(
        material=mp.Medium(index=n1),
        center=mp.Vector3(-0.25*(sx+2*dpml), 0, 0),
        size=mp.Vector3(0.5*(sx+2*dpml), mp.inf, mp.inf),
    ),
    mp.Block(
        material=mp.Medium(index=n2),
        center=mp.Vector3(0.25*(sx+2*dpml), 0, 0),
        size=mp.Vector3(0.5*(sx+2*dpml), mp.inf, mp.inf),
    ),
]

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

refl_mon = sim.add_mode_monitor(fcen, 0, 1, fluxreg)
sim.load_minus_flux_data(refl_mon, input_flux_data)

fig, ax = plt.subplots()
sim.plot2D(ax=ax)
if mp.am_master():
    fig.savefig('refl_coeff_flat_interface.png', dpi=150, bbox_inches='tight')

sim.run(
    until_after_sources=mp.stop_when_fields_decayed(
        50,
        mp.Ez,
        src_pt,
        tol,
    ),
)

# mode coefficient of reflected planewave                                                                         
refl_mode_coeff = mode_coeff(refl_mon)[0, 0, 1]

# reflection coefficient                                                                                          
refl_coeff = refl_mode_coeff / input_mode_coeff

# reflection coefficient (analytic)                                                                               
Fresnel_refl_coeff = refl_coeff_te(theta, n1, n2)

print(f"refl-coeff:, {refl_coeff} (Meep), {Fresnel_refl_coeff} (Fresnel)")

print(f"phase:, {cmath.phase(refl_coeff)} (Meep), "
      f"{cmath.phase(Fresnel_refl_coeff)} (Fresnel)")

[oskooi]

It seems that to demonstrate agreement in the phase of the reflected mode computed by Meep and the Fresnel equations requires placing the mode monitor exactly at the interface. The example above (incorrectly) placed the mode monitor halfway between the source and the interface. The position of the mode monitor is critical in calculations of the phase of the outgoing mode relative to the incident mode.

With the mode monitor positioned at the interface, the phase values are in agreement (see table below). However, the magnitude of the reflectance coefficient computed by Meep is less than one which is incorrect. Could this simply be due to a discretization error?

θ = 50°	reflectance coefficient	phase (radians)
Meep	0.08128-0.13957j	-1.04344
Fresnel	0.48743-0.87316j	-1.06165

θ = 61°	reflectance coefficient	phase (radians)
Meep	-0.00321-0.02276j	-1.710768
Fresnel	-0.15385-0.98809j	-1.72526

[smartalecH]

Wouldn't a simpler test just be free-space propagation of a planewave in a uniform medium (or even a mode propagating down a waveguide)? Place two monitors an arbitrary distance apart and measure their phase difference (arg of the ratio of the complex coefficients). It should approach exp(j*2*pi*n_eff/wavelength * L) where L is the length (the error will accumulate the longer the L due to numerical dispersion etc, but still approaches 0 with infinite resolution).

[oskooi]

Wouldn't a simpler test just be free-space propagation of a planewave in a uniform medium

This could be an alternative test. The main reason for using the example of a flat interface is that its setup involving a line/area DFT monitor can be adapted to the diffraction orders of (metasurface) gratings which are of general interest.

[stevengj]

It seems that to demonstrate agreement in the phase of the reflected mode computed by Meep and the Fresnel equations requires placing the mode monitor exactly at the interface.

If you place it at a distance L from the interface, I forgot to mention that there is a frequency-dependent correction factor exp(i*kx*L) where kx = cos(θ) ωn/c is the surface-normal wavevector (i.e. the phase accumulated as the wave propagates for a distance L in a medium with index n). (Here, L is the total distance from the interface from both monitors.)

[oskooi]

If you place it at a distance L from the interface, I forgot to mention that there is a frequency-dependent correction factor exp(ikxL)

In one of the tests above, the reflectance coefficient measured by Meep is -0.56759-0.82108j for $\theta = 50^{\circ}$. If we apply the phase correction factor exp(i*kx*L) to this result in which the monitor is a distance $L = 2.5$ from the interface, the corrected phase is 1.69709. However this value does not match the Fresnel result of -1.06165.

> r = -0.56759-0.82108j
> fcen = 1.0
> n1 = 1.5
> theta = math.radians(50.0)
> kx = math.cos(theta) * n1 * fcen
> L = 2.5
> corr_factor = cmath.exp(1j * kx * L) # phase correction factor
> cmath.phase(r / corr_factor)
1.6970909295581122

[stevengj]

You forgot a 2π factor to convert f to ω. Also, note that there is also a phase factor from the separation of the source and the interface, so you might need 2L instead of L. With these two fixes I get a phase of -1.050367.

[oskooi]

You forgot a 2π factor to convert f to ω. Also, note that there is also a phase factor from the separation of the source and the interface, so you might need 2L instead of L.

Good catch. I updated the test in #2428 to include this phase-correction factor.

https://github.com/NanoComp/meep/pull/2428/files#diff-f3b088e6eccdd36c8153197ed7245cc1648355d49a443429c011903e75f9d731R838-R839

[stevengj]

Note that it's problematic to measure the reflection coefficient exactly at the interface, because the amplitude of the left-going wave drops discontinuously to zero across the interface.

(Presumably, if you measure the left-going wave exactly at the interface in a discretized system like this, it gives some average of the correct results r and the 0 amplitude on the right of the interface. But that is equivalent to just multiplying r by a real number, which messes up the amplitude but not the phase, as you noted above.)