Reflectance of a Waveguide Taper

This example involves computing the reflectance of the fundamental mode of a linear waveguide taper. The structure and the simulation parameters are shown in the schematic below. We will verify that computing the reflectance, the fraction of the incident power which is reflected, using two different methods produces nearly identical results: (1) mode decomposition and (2) Poynting flux. Also, we will demonstrate that the scaling of the reflectance with the taper length is quadratic, consistent with analytical results from Optics Express, Vol. 16, pp. 11376-92, 2008.

The structure, which can be viewed as a two-port network, consists of a single-mode waveguide of width 1 μm (w1) at a wavelength of 6.67 μm and coupled to a second waveguide of width 2 μm (w2) via a linearly-sloped taper of variable length Lt. The material is silicon with ε=12. The taper geometry is defined using a single Prism object with eight vertices. PML absorbing boundaries surround the entire cell. An eigenmode current source with Ez polarization is used to launch the fundamental mode. The dispersion relation (or "band diagram") of the single-mode waveguide is shown in Tutorial/Eigenmode Source/Index-Guided Modes in a Ridge Waveguide. There is an eigenmode-expansion monitor placed at the midpoint of the first waveguide. This is a line monitor which extends beyond the waveguide in order to span the entire mode profile including its evanescent tails. The Fourier-transformed fields along this line monitor are used to compute the basis coefficients of the harmonic modes. These are computed separately via the eigenmode solver MPB. This is described in Mode Decomposition where it is also shown that the squared magnitude of the mode coefficient is equivalent to the power (Poynting flux) in the given eigenmode. The ratio of the complex mode coefficients can be used to compute the S parameters. In this example, we are computing |S11|2 which is the reflectance (shown in the line prefixed by "refl:,"). Another line monitor could have been placed in the second waveguide to compute the transmittance or |S21|2 into the various guided modes (since the second waveguide is multi mode). The scattered power into the radiative modes can then be computed as 1-|S11|2-|S21|2. Following usual practice, a normalization run is required involving a straight waveguide to compute the power in the source.

The structure has mirror symmetry in the $y$ direction which can be exploited to reduce the computation size by a factor of two. This requires using add_flux rather than add_mode_monitor (which is not optimized for symmetry) and specifying the keyword argument eig_parity=mp.ODD_Z+mp.EVEN_Y in the call to get_eigenmode_coefficients. Alternatively, the waveguide could have been oriented along an arbitrary oblique direction which would require specifying direction=mp.NO_DIRECTION and kpoint_func as the waveguide axis. For an example, see Tutorials/Eigenmode Source/Index-Guided Modes in a Ridge Waveguide.

In [1]:
import meep as mp
import matplotlib.pyplot as plt

resolution = 25   # pixels/μm

w1 = 1.0          # width of waveguide 1
w2 = 2.0          # width of waveguide 2
Lw = 10.0         # length of waveguides 1 and 2

# lengths of waveguide taper
Lts = [2**m for m in range(4)]

dair = 3.0        # length of air region
dpml_x = 6.0      # length of PML in x direction
dpml_y = 2.0      # length of PML in y direction

sy = dpml_y+dair+w2+dair+dpml_y

Si = mp.Medium(epsilon=12.0)

boundary_layers = [mp.PML(dpml_x,direction=mp.X),
                   mp.PML(dpml_y,direction=mp.Y)]

lcen = 6.67       # mode wavelength
fcen = 1/lcen     # mode frequency

symmetries = [mp.Mirror(mp.Y)]

R_coeffs = []
R_flux = []

for Lt in Lts:
    sx = dpml_x+Lw+Lt+Lw+dpml_x
    cell_size = mp.Vector3(sx,sy,0)

    src_pt = mp.Vector3(-0.5*sx+dpml_x+0.2*Lw)
    sources = [mp.EigenModeSource(src=mp.GaussianSource(fcen,fwidth=0.2*fcen),
                                  center=src_pt,
                                  size=mp.Vector3(y=sy-2*dpml_y),
                                  eig_match_freq=True,
                                  eig_parity=mp.ODD_Z+mp.EVEN_Y)]

    # straight waveguide
    vertices = [mp.Vector3(-0.5*sx-1,0.5*w1),
                mp.Vector3(0.5*sx+1,0.5*w1),
                mp.Vector3(0.5*sx+1,-0.5*w1),
                mp.Vector3(-0.5*sx-1,-0.5*w1)]

    sim = mp.Simulation(resolution=resolution,
                        cell_size=cell_size,
                        boundary_layers=boundary_layers,
                        geometry=[mp.Prism(vertices,height=mp.inf,material=Si)],
                        sources=sources,
                        symmetries=symmetries)

    mon_pt = mp.Vector3(-0.5*sx+dpml_x+0.7*Lw)
    flux = sim.add_flux(fcen,0,1,mp.FluxRegion(center=mon_pt,size=mp.Vector3(y=sy-2*dpml_y)))

    sim.run(until_after_sources=mp.stop_when_fields_decayed(50,mp.Ez,mon_pt,1e-9))

    res = sim.get_eigenmode_coefficients(flux,[1],eig_parity=mp.ODD_Z+mp.EVEN_Y)
    incident_coeffs = res.alpha
    incident_flux = mp.get_fluxes(flux)
    incident_flux_data = sim.get_flux_data(flux)

    sim.reset_meep()

    # linear taper
    vertices = [mp.Vector3(-0.5*sx-1,0.5*w1),
                mp.Vector3(-0.5*Lt,0.5*w1),
                mp.Vector3(0.5*Lt,0.5*w2),
                mp.Vector3(0.5*sx+1,0.5*w2),
                mp.Vector3(0.5*sx+1,-0.5*w2),
                mp.Vector3(0.5*Lt,-0.5*w2),
                mp.Vector3(-0.5*Lt,-0.5*w1),
                mp.Vector3(-0.5*sx-1,-0.5*w1)]

    sim = mp.Simulation(resolution=resolution,
                        cell_size=cell_size,
                        boundary_layers=boundary_layers,
                        geometry=[mp.Prism(vertices,height=mp.inf,material=Si)],
                        sources=sources,
                        symmetries=symmetries)

    flux = sim.add_flux(fcen,0,1,mp.FluxRegion(center=mon_pt,size=mp.Vector3(y=sy-2*dpml_y)))
    sim.load_minus_flux_data(flux,incident_flux_data)

    sim.run(until_after_sources=mp.stop_when_fields_decayed(50,mp.Ez,mon_pt,1e-9))

    res2 = sim.get_eigenmode_coefficients(flux,[1],eig_parity=mp.ODD_Z+mp.EVEN_Y)
    taper_coeffs = res2.alpha
    taper_flux = mp.get_fluxes(flux)

    R_coeffs.append(abs(taper_coeffs[0,0,1])**2/abs(incident_coeffs[0,0,0])**2)
    R_flux.append(-taper_flux[0]/incident_flux[0])
    print("refl:, {}, {:.8f}, {:.8f}".format(Lt,R_coeffs[-1],R_flux[-1]))
-----------
Initializing structure...
Halving computational cell along direction y
time for choose_chunkdivision = 0.00183296 s
Working in 2D dimensions.
Computational cell is 33 x 12 x 0 with resolution 25
     prism, center = (0,0,5e+19)
          height 1e+20, axis (0,0,1), 4 vertices:
          (-17.5,0.5,0)
          (17.5,0.5,0)
          (17.5,-0.5,0)
          (-17.5,-0.5,0)
          dielectric constant epsilon diagonal = (12,12,12)
time for set_epsilon = 0.447863 s
-----------
MPB solved for omega_1(0.519356,0,0) = 0.176186 after 10 iters
MPB solved for omega_1(0.424206,0,0) = 0.149878 after 7 iters
MPB solved for omega_1(0.424377,0,0) = 0.149925 after 5 iters
MPB solved for omega_1(0.424377,0,0) = 0.149925 after 1 iters
on time step 2379 (time=47.58), 0.00168153 s/step
field decay(t = 50.02): 7.54740931333862e-08 / 7.54740931333862e-08 = 1.0
on time step 4877 (time=97.54), 0.00160183 s/step
field decay(t = 100.04): 0.001735614282175414 / 0.001735614282175414 = 1.0
on time step 7353 (time=147.06), 0.00161556 s/step
field decay(t = 150.06): 0.45397875691800166 / 0.45397875691800166 = 1.0
on time step 9819 (time=196.38), 0.00162256 s/step
field decay(t = 200.08): 1.5637205312092566 / 1.5637205312092566 = 1.0
on time step 12317 (time=246.34), 0.0016016 s/step
field decay(t = 250.1): 1.2378300815914483 / 1.5637205312092566 = 0.7915929073555166
on time step 14830 (time=296.6), 0.0015922 s/step
field decay(t = 300.12): 0.031049012038060042 / 1.5637205312092566 = 0.01985585750034836
on time step 17280 (time=345.6), 0.00163305 s/step
field decay(t = 350.14): 8.860450697395457e-06 / 1.5637205312092566 = 5.66626230234599e-06
on time step 19780 (time=395.6), 0.00160038 s/step
field decay(t = 400.16): 2.7848020268498908e-11 / 1.5637205312092566 = 1.7808821789251224e-11
run 0 finished at t = 400.16 (20008 timesteps)
MPB solved for omega_1(0.519356,0,0) = 0.176186 after 10 iters
MPB solved for omega_1(0.424206,0,0) = 0.149878 after 7 iters
MPB solved for omega_1(0.424377,0,0) = 0.149925 after 5 iters
MPB solved for omega_1(0.424377,0,0) = 0.149925 after 1 iters
Dominant planewave for band 1: (0.424377,-0.000000,0.000000)
-----------
Initializing structure...
Halving computational cell along direction y
time for choose_chunkdivision = 0.00132203 s
Working in 2D dimensions.
Computational cell is 33 x 12 x 0 with resolution 25
     prism, center = (0,0,5e+19)
          height 1e+20, axis (0,0,1), 8 vertices:
          (-17.5,0.5,0)
          (-0.5,0.5,0)
          (0.5,1,0)
          (17.5,1,0)
          (17.5,-1,0)
          (0.5,-1,0)
          (-0.5,-0.5,0)
          (-17.5,-0.5,0)
          dielectric constant epsilon diagonal = (12,12,12)
time for set_epsilon = 0.664247 s
-----------
MPB solved for omega_1(0.519356,0,0) = 0.176186 after 10 iters
MPB solved for omega_1(0.424206,0,0) = 0.149878 after 7 iters
MPB solved for omega_1(0.424377,0,0) = 0.149925 after 5 iters
MPB solved for omega_1(0.424377,0,0) = 0.149925 after 1 iters
on time step 2337 (time=46.74), 0.00171202 s/step
field decay(t = 50.02): 7.55116423679346e-08 / 7.55116423679346e-08 = 1.0
on time step 4808 (time=96.16), 0.00161879 s/step
field decay(t = 100.04): 0.001737568234822459 / 0.001737568234822459 = 1.0
on time step 7266 (time=145.32), 0.00162736 s/step
field decay(t = 150.06): 0.4546391533256169 / 0.4546391533256169 = 1.0
on time step 9726 (time=194.52), 0.00162628 s/step
field decay(t = 200.08): 1.5661544172637174 / 1.5661544172637174 = 1.0
on time step 12196 (time=243.92), 0.00161985 s/step
field decay(t = 250.1): 1.2398082542830435 / 1.5661544172637174 = 0.7916258068914784
on time step 14700 (time=294), 0.00159784 s/step
field decay(t = 300.12): 0.0313140484926616 / 1.5661544172637174 = 0.01999422799405148
on time step 17203 (time=344.06), 0.00159832 s/step
field decay(t = 350.14): 9.723724333805377e-06 / 1.5661544172637174 = 6.208662585643395e-06
on time step 19681 (time=393.62), 0.0016146 s/step
field decay(t = 400.16): 5.154329670663798e-11 / 1.5661544172637174 = 3.291073736949327e-11
run 0 finished at t = 400.16 (20008 timesteps)
MPB solved for omega_1(0.519356,0,0) = 0.176186 after 9 iters
MPB solved for omega_1(0.424206,0,0) = 0.149878 after 7 iters
MPB solved for omega_1(0.424377,0,0) = 0.149925 after 5 iters
MPB solved for omega_1(0.424377,0,0) = 0.149925 after 1 iters
Dominant planewave for band 1: (0.424377,-0.000000,0.000000)
refl:, 1, 0.00033425, 0.00041190
-----------
Initializing structure...
Halving computational cell along direction y
time for choose_chunkdivision = 0.00108004 s
Working in 2D dimensions.
Computational cell is 34 x 12 x 0 with resolution 25
     prism, center = (0,0,5e+19)
          height 1e+20, axis (0,0,1), 4 vertices:
          (-18,0.5,0)
          (18,0.5,0)
          (18,-0.5,0)
          (-18,-0.5,0)
          dielectric constant epsilon diagonal = (12,12,12)
time for set_epsilon = 0.462233 s
-----------
MPB solved for omega_1(0.519356,0,0) = 0.176186 after 10 iters
MPB solved for omega_1(0.424206,0,0) = 0.149878 after 7 iters
MPB solved for omega_1(0.424377,0,0) = 0.149925 after 5 iters
MPB solved for omega_1(0.424377,0,0) = 0.149925 after 1 iters
on time step 2327 (time=46.54), 0.00171948 s/step
field decay(t = 50.02): 7.571674211765607e-08 / 7.571674211765607e-08 = 1.0
on time step 4775 (time=95.5), 0.0016345 s/step
field decay(t = 100.04): 0.0017406599773709985 / 0.0017406599773709985 = 1.0
on time step 7173 (time=143.46), 0.00166846 s/step
field decay(t = 150.06): 0.4552933938871408 / 0.4552933938871408 = 1.0
on time step 9613 (time=192.26), 0.0016394 s/step
field decay(t = 200.08): 1.56824628067331 / 1.56824628067331 = 1.0
on time step 12052 (time=241.04), 0.00164018 s/step
field decay(t = 250.1): 1.2414144426139568 / 1.56824628067331 = 0.7915940614129616
on time step 14512 (time=290.24), 0.00162657 s/step
field decay(t = 300.12): 0.031139198594396233 / 1.56824628067331 = 0.019856064049472478
on time step 16935 (time=338.7), 0.00165142 s/step
field decay(t = 350.14): 8.886312411745595e-06 / 1.56824628067331 = 5.666401075684586e-06
on time step 19391 (time=387.82), 0.00162891 s/step
field decay(t = 400.16): 2.7609192234474572e-11 / 1.56824628067331 = 1.7605138028843823e-11
run 0 finished at t = 400.16 (20008 timesteps)
MPB solved for omega_1(0.519356,0,0) = 0.176186 after 10 iters
MPB solved for omega_1(0.424206,0,0) = 0.149878 after 7 iters
MPB solved for omega_1(0.424377,0,0) = 0.149925 after 5 iters
MPB solved for omega_1(0.424377,0,0) = 0.149925 after 1 iters
Dominant planewave for band 1: (0.424377,-0.000000,0.000000)
-----------
Initializing structure...
Halving computational cell along direction y
time for choose_chunkdivision = 0.00079608 s
Working in 2D dimensions.
Computational cell is 34 x 12 x 0 with resolution 25
     prism, center = (0,0,5e+19)
          height 1e+20, axis (0,0,1), 8 vertices:
          (-18,0.5,0)
          (-1,0.5,0)
          (1,1,0)
          (18,1,0)
          (18,-1,0)
          (1,-1,0)
          (-1,-0.5,0)
          (-18,-0.5,0)
          dielectric constant epsilon diagonal = (12,12,12)
time for set_epsilon = 0.693091 s
-----------
MPB solved for omega_1(0.519356,0,0) = 0.176186 after 10 iters
MPB solved for omega_1(0.424206,0,0) = 0.149878 after 7 iters
MPB solved for omega_1(0.424377,0,0) = 0.149925 after 5 iters
MPB solved for omega_1(0.424377,0,0) = 0.149925 after 1 iters
on time step 2308 (time=46.16), 0.00173349 s/step
field decay(t = 50.02): 7.571396797724319e-08 / 7.571396797724319e-08 = 1.0
on time step 4692 (time=93.84), 0.00167819 s/step
field decay(t = 100.04): 0.0017395257849605922 / 0.0017395257849605922 = 1.0
on time step 7113 (time=142.26), 0.00165247 s/step
field decay(t = 150.06): 0.45384932485125673 / 0.45384932485125673 = 1.0
on time step 9546 (time=190.92), 0.00164416 s/step
field decay(t = 200.08): 1.5552337783783405 / 1.5552337783783405 = 1.0
on time step 11961 (time=239.22), 0.00165646 s/step
field decay(t = 250.1): 1.2260556872108188 / 1.5552337783783405 = 0.7883417298775754
on time step 14416 (time=288.32), 0.00162956 s/step
field decay(t = 300.12): 0.030083252819861385 / 1.5552337783783405 = 0.019343235234531443
on time step 16863 (time=337.26), 0.00163533 s/step
field decay(t = 350.14): 8.582297003142968e-06 / 1.5552337783783405 = 5.518332434942239e-06
on time step 19312 (time=386.24), 0.00163389 s/step
field decay(t = 400.16): 5.083050063922431e-11 / 1.5552337783783405 = 3.268351121606029e-11
run 0 finished at t = 400.16 (20008 timesteps)
MPB solved for omega_1(0.519356,0,0) = 0.176186 after 10 iters
MPB solved for omega_1(0.424206,0,0) = 0.149878 after 7 iters
MPB solved for omega_1(0.424377,0,0) = 0.149925 after 5 iters
MPB solved for omega_1(0.424377,0,0) = 0.149925 after 1 iters
Dominant planewave for band 1: (0.424377,-0.000000,0.000000)
refl:, 2, 0.00007150, 0.00008424
-----------
Initializing structure...
Halving computational cell along direction y
time for choose_chunkdivision = 0.00330091 s
Working in 2D dimensions.
Computational cell is 36 x 12 x 0 with resolution 25
     prism, center = (0,0,5e+19)
          height 1e+20, axis (0,0,1), 4 vertices:
          (-19,0.5,0)
          (19,0.5,0)
          (19,-0.5,0)
          (-19,-0.5,0)
          dielectric constant epsilon diagonal = (12,12,12)
time for set_epsilon = 0.496261 s
-----------
MPB solved for omega_1(0.519356,0,0) = 0.176186 after 10 iters
MPB solved for omega_1(0.424206,0,0) = 0.149878 after 7 iters
MPB solved for omega_1(0.424377,0,0) = 0.149925 after 5 iters
MPB solved for omega_1(0.424377,0,0) = 0.149925 after 1 iters
on time step 2138 (time=42.76), 0.00187101 s/step
field decay(t = 50.02): 7.571674212930823e-08 / 7.571674212930823e-08 = 1.0
on time step 4454 (time=89.08), 0.00172765 s/step
field decay(t = 100.04): 0.0017406599772782168 / 0.0017406599772782168 = 1.0
on time step 6777 (time=135.54), 0.0017222 s/step
field decay(t = 150.06): 0.45529339380861855 / 0.45529339380861855 = 1.0
on time step 9069 (time=181.38), 0.0017456 s/step
field decay(t = 200.08): 1.568246280713512 / 1.568246280713512 = 1.0
on time step 11341 (time=226.82), 0.00176058 s/step
field decay(t = 250.1): 1.2414144439515642 / 1.568246280713512 = 0.791594062245601
on time step 13493 (time=269.86), 0.00185906 s/step
field decay(t = 300.12): 0.031139203524192197 / 1.568246280713512 = 0.01985606719247225
on time step 15734 (time=314.68), 0.00178549 s/step
field decay(t = 350.14): 8.886463070035761e-06 / 1.568246280713512 = 5.666497143543454e-06
on time step 18069 (time=361.38), 0.00171313 s/step
field decay(t = 400.16): 2.7650862223635545e-11 / 1.568246280713512 = 1.7631709103148714e-11
run 0 finished at t = 400.16 (20008 timesteps)
MPB solved for omega_1(0.519356,0,0) = 0.176186 after 10 iters
MPB solved for omega_1(0.424206,0,0) = 0.149878 after 7 iters
MPB solved for omega_1(0.424377,0,0) = 0.149925 after 5 iters
MPB solved for omega_1(0.424377,0,0) = 0.149925 after 1 iters
Dominant planewave for band 1: (0.424377,-0.000000,0.000000)
-----------
Initializing structure...
Halving computational cell along direction y
time for choose_chunkdivision = 0.00128293 s
Working in 2D dimensions.
Computational cell is 36 x 12 x 0 with resolution 25
     prism, center = (0,0,5e+19)
          height 1e+20, axis (0,0,1), 8 vertices:
          (-19,0.5,0)
          (-2,0.5,0)
          (2,1,0)
          (19,1,0)
          (19,-1,0)
          (2,-1,0)
          (-2,-0.5,0)
          (-19,-0.5,0)
          dielectric constant epsilon diagonal = (12,12,12)
time for set_epsilon = 0.755732 s
-----------
MPB solved for omega_1(0.519356,0,0) = 0.176186 after 10 iters
MPB solved for omega_1(0.424206,0,0) = 0.149878 after 7 iters
MPB solved for omega_1(0.424377,0,0) = 0.149925 after 5 iters
MPB solved for omega_1(0.424377,0,0) = 0.149925 after 1 iters
on time step 1883 (time=37.66), 0.00212453 s/step
field decay(t = 50.02): 7.569671091769381e-08 / 7.569671091769381e-08 = 1.0
on time step 4139 (time=82.78), 0.00177353 s/step
field decay(t = 100.04): 0.0017391263700697627 / 0.0017391263700697627 = 1.0
on time step 6309 (time=126.18), 0.00184352 s/step
field decay(t = 150.06): 0.453841743379235 / 0.453841743379235 = 1.0
on time step 8560 (time=171.2), 0.0017771 s/step
field decay(t = 200.08): 1.5555871736720326 / 1.5555871736720326 = 1.0
on time step 10769 (time=215.38), 0.0018108 s/step
field decay(t = 250.1): 1.2260941054634609 / 1.5555871736720326 = 0.7881873328700771
on time step 13056 (time=261.12), 0.00174983 s/step
field decay(t = 300.12): 0.02968036484241766 / 1.5555871736720326 = 0.01907984672588669
on time step 15299 (time=305.98), 0.00178368 s/step
field decay(t = 350.14): 7.2347918961150245e-06 / 1.5555871736720326 = 4.650843114781525e-06
on time step 17653 (time=353.06), 0.00169945 s/step
on time step 19969 (time=399.38), 0.00172778 s/step
field decay(t = 400.16): 8.407569561060486e-11 / 1.5555871736720326 = 5.4047562896870926e-11
run 0 finished at t = 400.16 (20008 timesteps)
MPB solved for omega_1(0.519356,0,0) = 0.176186 after 9 iters
MPB solved for omega_1(0.424206,0,0) = 0.149878 after 7 iters
MPB solved for omega_1(0.424377,0,0) = 0.149925 after 5 iters
MPB solved for omega_1(0.424377,0,0) = 0.149925 after 1 iters
Dominant planewave for band 1: (0.424377,-0.000000,0.000000)
refl:, 4, 0.00004024, 0.00004337
-----------
Initializing structure...
Halving computational cell along direction y
time for choose_chunkdivision = 0.00212097 s
Working in 2D dimensions.
Computational cell is 40 x 12 x 0 with resolution 25
     prism, center = (0,0,5e+19)
          height 1e+20, axis (0,0,1), 4 vertices:
          (-21,0.5,0)
          (21,0.5,0)
          (21,-0.5,0)
          (-21,-0.5,0)
          dielectric constant epsilon diagonal = (12,12,12)
time for set_epsilon = 0.565803 s
-----------
MPB solved for omega_1(0.519356,0,0) = 0.176186 after 10 iters
MPB solved for omega_1(0.424206,0,0) = 0.149878 after 7 iters
MPB solved for omega_1(0.424377,0,0) = 0.149925 after 5 iters
MPB solved for omega_1(0.424377,0,0) = 0.149925 after 1 iters
on time step 1905 (time=38.1), 0.00210044 s/step
field decay(t = 50.02): 7.571674211765524e-08 / 7.571674211765524e-08 = 1.0
on time step 4047 (time=80.94), 0.00186782 s/step
field decay(t = 100.04): 0.0017406599773711756 / 0.0017406599773711756 = 1.0
on time step 6216 (time=124.32), 0.0018447 s/step
field decay(t = 150.06): 0.4552933938912802 / 0.4552933938912802 = 1.0
on time step 8357 (time=167.14), 0.00186832 s/step
field decay(t = 200.08): 1.5682462811448037 / 1.5682462811448037 = 1.0
on time step 10518 (time=210.36), 0.00185132 s/step
field decay(t = 250.1): 1.2414144444964543 / 1.5682462811448037 = 0.7915940623753525
on time step 12633 (time=252.66), 0.00189134 s/step
on time step 14831 (time=296.62), 0.00182053 s/step
field decay(t = 300.12): 0.03113920465258198 / 1.5682462811448037 = 0.01985606790653486
on time step 17006 (time=340.12), 0.00183974 s/step
field decay(t = 350.14): 8.886516913937694e-06 / 1.5682462811448037 = 5.666531475815538e-06
on time step 19110 (time=382.2), 0.0019017 s/step
field decay(t = 400.16): 2.7625870947574233e-11 / 1.5682462811448037 = 1.7615773287476016e-11
run 0 finished at t = 400.16 (20008 timesteps)
MPB solved for omega_1(0.519356,0,0) = 0.176186 after 10 iters
MPB solved for omega_1(0.424206,0,0) = 0.149878 after 7 iters
MPB solved for omega_1(0.424377,0,0) = 0.149925 after 5 iters
MPB solved for omega_1(0.424377,0,0) = 0.149925 after 1 iters
Dominant planewave for band 1: (0.424377,-0.000000,0.000000)
-----------
Initializing structure...
Halving computational cell along direction y
time for choose_chunkdivision = 0.000883102 s
Working in 2D dimensions.
Computational cell is 40 x 12 x 0 with resolution 25
     prism, center = (0,0,5e+19)
          height 1e+20, axis (0,0,1), 8 vertices:
          (-21,0.5,0)
          (-4,0.5,0)
          (4,1,0)
          (21,1,0)
          (21,-1,0)
          (4,-1,0)
          (-4,-0.5,0)
          (-21,-0.5,0)
          dielectric constant epsilon diagonal = (12,12,12)
time for set_epsilon = 0.863971 s
-----------
MPB solved for omega_1(0.519356,0,0) = 0.176186 after 10 iters
MPB solved for omega_1(0.424206,0,0) = 0.149878 after 7 iters
MPB solved for omega_1(0.424377,0,0) = 0.149925 after 5 iters
MPB solved for omega_1(0.424377,0,0) = 0.149925 after 1 iters
on time step 1966 (time=39.32), 0.00203556 s/step
field decay(t = 50.02): 7.568972576123606e-08 / 7.568972576123606e-08 = 1.0
on time step 4042 (time=80.84), 0.00192733 s/step
field decay(t = 100.04): 0.0017392532865220022 / 0.0017392532865220022 = 1.0
on time step 6101 (time=122.02), 0.00194315 s/step
field decay(t = 150.06): 0.4541781140774447 / 0.4541781140774447 = 1.0
on time step 8223 (time=164.46), 0.00188516 s/step
field decay(t = 200.08): 1.5597082712382018 / 1.5597082712382018 = 1.0
on time step 10355 (time=207.1), 0.00187692 s/step
on time step 12479 (time=249.58), 0.00188403 s/step
field decay(t = 250.1): 1.231655150453248 / 1.5597082712382018 = 0.7896702051053925
on time step 14552 (time=291.04), 0.00193036 s/step
field decay(t = 300.12): 0.030221199579030786 / 1.5597082712382018 = 0.019376187288562084
on time step 16622 (time=332.44), 0.00193299 s/step
field decay(t = 350.14): 6.594515383229233e-06 / 1.5597082712382018 = 4.228044118785151e-06
on time step 18688 (time=373.76), 0.00193665 s/step
field decay(t = 400.16): 3.0985772121881132e-09 / 1.5597082712382018 = 1.986638956353199e-09
on time step 20781 (time=415.62), 0.00191162 s/step
field decay(t = 450.18): 2.4580489315863303e-13 / 1.5597082712382018 = 1.5759671067429584e-13
run 0 finished at t = 450.18 (22509 timesteps)
MPB solved for omega_1(0.519356,0,0) = 0.176186 after 10 iters
MPB solved for omega_1(0.424206,0,0) = 0.149878 after 7 iters
MPB solved for omega_1(0.424377,0,0) = 0.149925 after 5 iters
MPB solved for omega_1(0.424377,0,0) = 0.149925 after 1 iters
Dominant planewave for band 1: (0.424377,-0.000000,0.000000)
refl:, 8, 0.00001847, 0.00001902

Note that the reflectance is computed for five different geometrically-scaled taper lengths: 1, 2, 4, 8, and 16 μm. A quadratic scaling of the reflectance with the taper length appears as a straight line on a log-log plot. The results are plotted below.

In [2]:
plt.figure(dpi=200)
plt.loglog(Lts,R_coeffs,'bo-',label='mode decomposition')
plt.loglog(Lts,R_flux,'ro-',label='Poynting flux')
plt.loglog(Lts,[0.005/Lt**2 for Lt in Lts],'k-',label=r'quadratic reference (1/Lt$^2$)')
plt.legend(loc='upper right')
plt.xlabel('taper length Lt (μm)')
plt.ylabel('reflectance')
plt.show()

The reflectance values computed using the two methods are nearly identical. For reference, a line with quadratic scaling is shown in black. The reflectance of the linear waveguide taper decreases quadratically with the taper length which is consistent with the analytic theory.

In the reflected-flux calculation, we apply our usual trick of first performing a reference simulation with just the incident field and then subtracting that from our taper simulation with load_minus_flux_data, so that what is left over is the reflected fields (from which we obtain the reflected flux). In principle, this trick would not be required for the mode-decomposition method, because the reflected mode is orthogonal to the forward mode and so the decomposition will separate the forward and reflected coefficients automatically. However, this is only true in the limit of infinite resolution — for a finite resolution, the reflected mode used for the mode coefficient calculation (calculated via MPB) is not exactly orthogonal to the forward mode propagating in Meep (whose discretization scheme is different from that of MPB). In consequence, if you did not subtract the fields of the reference simulation, the mode-coefficient could only calculate the reflected power down to a "noise floor" set by the discretization error. With the subtraction, in contrast, you can compute much smaller reflections (limited by the floating-point precision).