Mode Decomposition


This tutorial demonstrates Meep's mode-decomposition feature which is used to separate a given field profile into a superposition of harmonic basis modes. The example involves computing the reflection coefficient of the fundamental mode of a linear waveguide taper as shown in the schematic below. We will verify that the scaling of the reflection coefficient with the taper length is quadratic, consistent with analytical results from Optics Express, Vol. 16, pp. 11376-92, 2008.

The 2d structure consists of a single-mode waveguide of width 1 μm (w1) with ε=12 in vacuum coupled to a second waveguide of width 2 μm (w2) via a linearly-sloped taper of length Lt. PML absorbing boundaries surround the computational cell. An eigenmode source is used to launch the fundamental mode. 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 capture the entire mode profile including its evanescent tails. The Fourier-transformed fields along the line monitor are used to compute the basis coefficients of the harmonic modes which are obtained separately via the eigenmode solver MPB. The details of this calculation are provided in Mode Decomposition. An alternative method for computing the reflection coefficient involves calculating the Poynting flux. This approach requires a separate normalization run to compute the incident fields due to the source and does not provide phase information. The mode-decomposition feature uses only a single run to compute the complex reflection coefficient. This enables the computation of S parameters. The simulation script is shown below.

import meep as mp
import math
import argparse
import numpy as np

def main(args):

    resolution = args.res

    w1 = 1            # width of waveguide 1
    w2 = args.w2      # width of waveguide 2
    Lw = 10           # length of waveguide 1 and 2
    Lt = args.Lt      # taper length

    Si = mp.Medium(epsilon=12.0)

    dair = 3.0
    dpml = 3.0

    sx = dpml + Lw + Lt + Lw + dpml
    sy = dpml + dair + w2 + dair + dpml
    cell_size = mp.Vector3(sx, sy, 0)

    geometry = [ mp.Block(material=Si, center=mp.Vector3(0,0,0), size=mp.Vector3(mp.inf,w1,mp.inf)), 
                 mp.Block(material=Si, center=mp.Vector3(0.5*sx-0.5*(Lt+Lw+dpml),0,0), size=mp.Vector3(Lt+Lw+dpml,w2,mp.inf)) ]

    # form linear taper

    hh = w2
    ww = 2*Lt

    # taper angle (CCW, relative to +X axis)
    rot_theta = math.atan(0.5*(w2-w1)/Lt)

    pvec = mp.Vector3(-0.5*sx+dpml+Lw,0.5*w1,0)
    cvec = mp.Vector3(-0.5*sx+dpml+Lw+0.5*ww,0.5*hh+0.5*w1,0)
    rvec = cvec-pvec
    rrvec = rvec.rotate(mp.Vector3(0,0,1), rot_theta)

    geometry.append(mp.Block(material=mp.air, center=pvec+rrvec, size=mp.Vector3(ww,hh,mp.inf),
                             e1=mp.Vector3(1,0,0).rotate(mp.Vector3(0,0,1),rot_theta),
                             e2=mp.Vector3(0,1,0).rotate(mp.Vector3(0,0,1),rot_theta),
                             e3=mp.Vector3(0,0,1)))

    pvec = mp.Vector3(-0.5*sx+dpml+Lw,-0.5*w1,0)
    cvec = mp.Vector3(-0.5*sx+dpml+Lw+0.5*ww,-(0.5*hh+0.5*w1),0)
    rvec = cvec-pvec
    rrvec = rvec.rotate(mp.Vector3(0,0,1),-rot_theta)

    geometry.append(mp.Block(material=mp.air, center=pvec+rrvec, size=mp.Vector3(ww,hh,mp.inf),
                             e1=mp.Vector3(1,0,0).rotate(mp.Vector3(0,0,1),-rot_theta),
                             e2=mp.Vector3(0,1,0).rotate(mp.Vector3(0,0,1),-rot_theta),
                             e3=mp.Vector3(0,0,1)))

    boundary_layers = [ mp.PML(dpml) ]

    # mode frequency
    fcen = 0.15

    sources = [ mp.EigenModeSource(src=mp.GaussianSource(fcen, fwidth=0.5*fcen),
                                   component=mp.Ez,
                                   size=mp.Vector3(0,sy-2*dpml,0),
                                   center=mp.Vector3(-0.5*sx+dpml,0,0),
                                   eig_match_freq=True,
                                   eig_parity=mp.ODD_Z,
                                   eig_kpoint=mp.Vector3(0.4,0,0),
                                   eig_resolution=32) ]

    symmetries = [ mp.Mirror(mp.Y,+1) ]

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

    xm = -0.5*sx + dpml + 0.5*Lw; # x-coordinate of monitor
    mflux = sim.add_flux(fcen, 0, 1, mp.FluxRegion(center=mp.Vector3(xm,0), size=mp.Vector3(0,sy-2*dpml)));

    sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ez, mp.Vector3(xm,0), 1e-10))

    bands = np.array([1],dtype=np.int32) # indices of modes for which to compute expansion coefficients
    num_bands = 1;
    alpha = np.zeros(2*num_bands,dtype=np.complex128) # preallocate array to store coefficients
    vgrp = np.zeros(num_bands,dtype=np.float64)       # also store mode group velocities
    mvol = mp.volume(mp.vec(xm,-0.5*sy+dpml),mp.vec(xm,+0.5*sy-dpml))
    sim.fields.get_eigenmode_coefficients(mflux, mp.X, mvol, bands, alpha, vgrp)

    alpha0Plus  = alpha[2*0 + 0]; # coefficient of forward-traveling fundamental mode
    alpha0Minus = alpha[2*0 + 1]; # coefficient of backward-traveling fundamental mode

    print("refl:, {}, {:.8f}".format(Lt, abs(alpha0Minus)**2))

if __name__ == '__main__':
    parser = argparse.ArgumentParser()
    parser.add_argument('-Lt',  type=float, default=3.0, help='taper length (default: 3.0)')
    parser.add_argument('-w2',  type=float, default=2.0, help='width of outgoing waveguide (default: 2.0)')
    parser.add_argument('-res', type=int, default=20, help='resolution (default: 20)')
    args = parser.parse_args()
    main(args)

To investigate the scaling, we compute the reflection coefficient for a range of taper lengths chosen on a logarithmic scale (i.e., 1, 2, 4, ..., 64). The results are shown below in blue. For reference, a quadratic scaling is shown in black. Consistent with analytical results, the reflection coefficient of the linear waveguide taper decreases quadratically with the taper length.