clarify S-matrix computation in directional-coupler calculation

[stevengj]

Might be nice to update this tutorial to explictly calculate the S-matrix elements like

S12 = p2_coeff / p1_coeff
p2_trans = abs(S12)**2

etcetera, and to clarify why the eigenmode-coefficient feature makes it unnecessary to do a separate normalization run (as long as port 1 is strictly to the right of the source — you probably want it to be at least a couple of pixels to the right to avoid overlaps due to discretization).

If you are computing the reflection coefficient S11, then you will see a "noise floor" below which you cannot measure reflection in this way, due to the slight discretization error in the eigenmode-coefficient extraction — in that case, you would have to do a separate normalization run and subtract the Fourier-transformed fields in port 1 (as we do in flux calculations) in order to accurately extract a tiny backscattered field.

[oskooi]

The following plot shows the S-matrix elements S₁₁ and S₁₂ vs. resolution (computed in a single run). Both quantities seem to reach a "noise floor" in the limit of infinite resolution with S₁₂ being ~2 orders of magnitude smaller than S₁₁. However, it's somewhat surprising that S₁₂ is only ~10^-8 and not much smaller (i.e., closer to machine precision).

The scattering parameters are computed using examples/coupler.py via:

    eig_bands = [1]
    eig_parity = mp.EVEN_Y+mp.ODD_Z

    # mode coefficients
    p1_coeff_plus = sim.get_eigenmode_coefficients(mode1,eig_bands,eig_parity).alpha[0,0,0]
    p1_coeff_minus = sim.get_eigenmode_coefficients(mode1,eig_bands,eig_parity).alpha[0,0,1]
    p2_coeff_minus = sim.get_eigenmode_coefficients(mode2,eig_bands,eig_parity).alpha[0,0,1]

    # S parameters
    S11 = abs(p1_coeff_minus)**2/abs(p1_coeff_plus)**2
    S12 = abs(p2_coeff_minus)**2/abs(p1_coeff_plus)**2

[oskooi]

The following is a demonstration of computing S₁₁ using just a straight waveguide for two cases: (1) with and (2) without normalization of the incident fields. The plot of S₁₁ vs. resolution below shows that at the largest resolution (200) (1) produces a value that is 10²⁴ times smaller than (2), i.e., 10^-31 vs. 10^-7. While this is a large discrepancy in the relative value of the backscattered fields, an absolute value of 10^-7 (the "noise floor") is likely sufficiently small that a separate normalization is not necessary. This is the point that should be emphasized in Tutorial/GDSII ImportS-Parameters of a Directional Coupler.

import meep as mp
import argparse

from meep import mpb
mpb.verbosity(0)

parser = argparse.ArgumentParser()
parser.add_argument('-res', type=int, default=50, help='resolution (default: 50 pixels/um)')
args = parser.parse_args()

resolution = args.res

dpml = 1.0

silicon = mp.Medium(epsilon=12)

fcen = 1/1.55
df = 0.2*fcen

sxy = 5.0

cell_size = mp.Vector3(sxy,sxy,0)

boundary_layers = [mp.PML(thickness=dpml)]

eig_parity = mp.EVEN_Y+mp.ODD_Z

geometry = [mp.Block(material=silicon,
                     center=mp.Vector3(),
                     size=mp.Vector3(mp.inf,1.0,mp.inf))]

sources = [mp.EigenModeSource(src=mp.GaussianSource(fcen,fwidth=df),
                              center=mp.Vector3(-0.5*sxy+dpml,0),
                              size=mp.Vector3(0,sxy,0),
                              eig_band=1,
                              eig_parity=eig_parity,
                              eig_match_freq=True)]

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

mode = sim.add_mode_monitor(fcen, 0, 1,
                            mp.ModeRegion(center=mp.Vector3(0.5*sxy-dpml),size=mp.Vector3(0,sxy,0)))

sim.run(until_after_sources=20)

input_data = sim.get_flux_data(mode)
input_flux = mp.get_fluxes(mode)[0]

sim.reset_meep()

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

mode = sim.add_mode_monitor(fcen, 0, 1,
                            mp.ModeRegion(center=mp.Vector3(0.5*sxy-dpml),size=mp.Vector3(0,sxy,0)))

sim.load_minus_flux_data(mode, input_data) ## <----- normalization
sim.run(until_after_sources=20)

# mode coefficients
coeff_minus = sim.get_eigenmode_coefficients(mode,[1],eig_parity).alpha[0,0,1]

# S parameters
S11 = abs(coeff_minus)**2/input_flux

print("S11:, {}, {}".format(resolution,S11))

[stevengj]

S12 being smaller than S11 by 2–3 orders of magnitude seems reasonable. Both S11 and S12 arise from the finite taper/bend length (i.e. the breaking of translational symmetry, which allows backsacttering), and S12 is smaller because the light also has to jump to the second waveguide.