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.
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.
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).