Eigenmode Source#


This tutorial demonstrates using the EigenModeSource to launch a single eigenmode propagating in a single direction. Examples are provided for two kinds of eigenmodes in lossless, dielectric media: (1) localized (i.e., guided) and (2) non-localized (i.e., radiative planewave).

Index-Guided Modes in a Ridge Waveguide#

The first structure, shown in the schematic below, is a 2d ridge waveguide with ε=12, width =1 μm, and out-of-plane electric field Ez. The dispersion relation ω(k) for index-guided modes with even mirror symmetry in the -direction is computed using MPB and shown as blue lines. The light cone which denotes radiative modes is the section in solid green. Using this waveguide configuration, we will investigate two different frequency regimes: (1) single mode (normalized frequency of 0.15) and (2) multi mode (normalized frequency of 0.35), both shown as dotted horizontal lines in the figures. We will use the eigenmode source to excite a specific mode in each case — labeled A and B in the band diagram — and compare the results to using a constant-amplitude source for straight and rotated waveguides. Finally, we will demonstrate that a single monitor plane in the -direction is sufficient for computing the total Poynting flux in a waveguide with any arbitrary orientation.

The simulation script is in examples/oblique-source.py. The notebook is examples/oblique-source.ipynb.

The simulation consists of two separate parts: (1) computing the Poynting flux and (2) plotting the field profile. The field profile is generated by setting the flag compute_flux=False. For the single-mode case, a constant-amplitude current source (eig_src=False) excites both the waveguide mode and radiating fields in both directions (i.e., forwards and backwards). This is shown in the main inset of the first of two figures above. The EigenModeSource excites only the forward-going waveguide mode A as shown in the smaller inset. Launching this mode requires setting eig_src=True, fsrc=0.15, and bnum=1. Note that the length of the EigenModeSource must be large enough that it captures the entire guided mode (i.e., the fields should decay to roughly zero at the edges). The line source does not need to span the entire length of the cell but should be slightly larger than the waveguide width. The parameter rot_angle specifies the rotation angle of the waveguide axis and is initially 0° (i.e., straight or horizontal orientation). This enables eig_parity to include EVEN_Y in addition to ODD_Z and the simulation to include an overall mirror symmetry plane in the -direction.

For the multi-mode case, a constant-amplitude current source excites a superposition of the two waveguide modes in addition to the radiating field. This is shown in the main inset of the second figure above. The EigenModeSource excites only a given mode: A (fsrc=0.35, bnum=2) or B (fsrc=0.35, bnum=1) as shown in the smaller insets.

import meep as mp
import numpy as np
import matplotlib.pyplot as plt

resolution = 50 # pixels/μm

cell_size = mp.Vector3(14,14)

pml_layers = [mp.PML(thickness=2)]

# rotation angle (in degrees) of waveguide, counter clockwise (CCW) around z-axis
rot_angle = np.radians(20)

w = 1.0 # width of waveguide

geometry = [mp.Block(center=mp.Vector3(),
                     size=mp.Vector3(mp.inf,w,mp.inf),
                     e1=mp.Vector3(x=1).rotate(mp.Vector3(z=1), rot_angle),
                     e2=mp.Vector3(y=1).rotate(mp.Vector3(z=1), rot_angle),
                     material=mp.Medium(epsilon=12))]

fsrc = 0.15 # frequency of eigenmode or constant-amplitude source
bnum = 1    # band number of eigenmode

kpoint = mp.Vector3(x=1).rotate(mp.Vector3(z=1), rot_angle)

compute_flux = True # compute flux (True) or plot the field profile (False)

eig_src = True # eigenmode (True) or constant-amplitude (False) source

if eig_src:
    sources = [mp.EigenModeSource(src=mp.GaussianSource(fsrc,fwidth=0.2*fsrc) if compute_flux else mp.ContinuousSource(fsrc),
                                  center=mp.Vector3(),
                                  size=mp.Vector3(y=3*w),
                                  direction=mp.NO_DIRECTION,
                                  eig_kpoint=kpoint,
                                  eig_band=bnum,
                                  eig_parity=mp.EVEN_Y+mp.ODD_Z if rot_angle == 0 else mp.ODD_Z,
                                  eig_match_freq=True)]
else:
    sources = [mp.Source(src=mp.GaussianSource(fsrc,fwidth=0.2*fsrc) if compute_flux else mp.ContinuousSource(fsrc),
                         center=mp.Vector3(),
                         size=mp.Vector3(y=3*w),
                         component=mp.Ez)]

sim = mp.Simulation(cell_size=cell_size,
                    resolution=resolution,
                    boundary_layers=pml_layers,
                    sources=sources,
                    geometry=geometry,
                    symmetries=[mp.Mirror(mp.Y)] if rot_angle == 0 else [])

if compute_flux:
    tran = sim.add_flux(fsrc, 0, 1, mp.FluxRegion(center=mp.Vector3(x=5), size=mp.Vector3(y=14)))
    sim.run(until_after_sources=50)
    res = sim.get_eigenmode_coefficients(tran,
                                         [1],
                                         eig_parity=mp.EVEN_Y+mp.ODD_Z if rot_angle == 0 else mp.ODD_Z,
                                         direction=mp.NO_DIRECTION,
                                         kpoint_func=lambda f,n: kpoint)
    print("flux:, {:.6f}, {:.6f}".format(mp.get_fluxes(tran)[0],abs(res.alpha[0,0,0])**2))
else:
    sim.run(until=100)
    sim.plot2D(output_plane=mp.Volume(center=mp.Vector3(), size=mp.Vector3(10,10)),
               fields=mp.Ez,
               field_parameters={'alpha':0.9})
    plt.show()

Note that in EigenModeSource as well as get_eigenmode_coefficients, the direction property must be set to NO_DIRECTION whenever eig_kpoint specifies the waveguide axis.

What Happens When the Source Time Profile is a Pulse?#

The eigenmode source launches a fixed spatial mode profile specified by either its frequency (eig_match_freq=True) or wavevector (eig_match_freq=False) multiplied by the time profile. When the time profile of the source has a finite bandwidth, e.g. a Gaussian pulse (which is typical for calculations involving Fourier-transformed fields such as mode coefficients or S-parameters), then the frequency-dependence (dispersion) of the true modal pattern means that the eigenmode source does not match the desired mode exactly over the whole bandwidth. This is described in Section 4.2.2 of the review article Electromagnetic Wave Source Conditions. A more accurate mode profile may be obtained by adding multiple narrow-band eigenmode sources at the same position at several frequencies across the bandwidth, but this has the disadvantage that the runtime increases as you add more frequency points due to the narrower source bandwidths. However, a single broadband eigenmode source is often sufficient for most practical applications (excepting cases with extreme modal dispersion, e.g. near a cutoff frequency).

This can be demonstrated by computing the error in a broadband eigenmode source via the backward-propagating and scattered power (i.e., any fields which are not forward-propagating waveguide modes) for the single and multi mode ridge waveguide.

import meep as mp
import numpy as np
import matplotlib.pyplot as plt

resolution = 50 # pixels/μm

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

dpml = 1
pml_layers = [mp.PML(thickness=dpml)]

w = 1.0 # width of waveguide

geometry = [mp.Block(center=mp.Vector3(),
                     size=mp.Vector3(mp.inf,w,mp.inf),
                     material=mp.Medium(epsilon=12))]

fsrc = 0.35 # frequency of eigenmode source
bnum = 1    # band number of eigenmode

df = 0.1    # frequency width
nfreq = 51  # number of frequencies

sources = [mp.EigenModeSource(src=mp.GaussianSource(fsrc,fwidth=df),
                              center=mp.Vector3(0,0),
                              size=mp.Vector3(0,sxy),
                              eig_parity=mp.EVEN_Y+mp.ODD_Z,
                              eig_band=bnum,
                              eig_match_freq=True)]

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

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

flux_mon_tp = sim.add_flux(fsrc,df,nfreq,mp.FluxRegion(center=mp.Vector3(0,+0.5*sxy-dpml),size=mp.Vector3(sxy-2*dpml,0),weight=+1))
flux_mon_bt = sim.add_flux(fsrc,df,nfreq,mp.FluxRegion(center=mp.Vector3(0,-0.5*sxy+dpml),size=mp.Vector3(sxy-2*dpml,0),weight=-1))
flux_mon_rt = sim.add_flux(fsrc,df,nfreq,mp.FluxRegion(center=mp.Vector3(+0.5*sxy-dpml,0),size=mp.Vector3(0,sxy-2*dpml),weight=+1))
flux_mon_lt = sim.add_flux(fsrc,df,nfreq,mp.FluxRegion(center=mp.Vector3(-0.5*sxy+dpml,0),size=mp.Vector3(0,sxy-2*dpml),weight=-1))

sim.run(until_after_sources=50)

freqs = mp.get_flux_freqs(flux_mon_tp)
flux_tp = np.asarray(mp.get_fluxes(flux_mon_tp))
flux_bt = np.asarray(mp.get_fluxes(flux_mon_bt))
flux_rt = np.asarray(mp.get_fluxes(flux_mon_rt))
flux_lt = np.asarray(mp.get_fluxes(flux_mon_lt))

res = sim.get_eigenmode_coefficients(flux_mon_rt,[bnum],eig_parity=mp.EVEN_Y+mp.ODD_Z)
coeffs = res.alpha[0,:,0]
guided = np.power(np.abs(coeffs),2)

flux_total = flux_tp + flux_bt + flux_rt + flux_lt
back = flux_lt / flux_total
scat = (flux_total - guided) / flux_total

if mp.am_master():
    fig = plt.figure()
    plt.subplot(1,2,1)
    plt.plot(freqs,back,'bo-')
    plt.xlabel('frequency')
    plt.ylabel('backward power (fraction of total power)')
    plt.subplot(1,2,2)
    plt.plot(freqs,scat,'ro-')
    plt.xlabel('frequency')
    plt.ylabel('scattered power (fraction of total power)')
    fig.subplots_adjust(wspace=0.5, hspace=0)
    fig.suptitle("multi mode waveguide with pulsed eigenmode source\n center frequency = {}, band = {} (mode B)".format(fsrc,bnum))
    plt.savefig('multi_mode_eigsource_B.png',dpi=150,bbox_inches='tight')

Results are shown for the single mode waveguide with one eigenmode A (band 1) using a center frequency of 0.15 and multi mode waveguide with two eigenmodes A (higher-order mode, band 2) and B (fundamental mode, band 1), all with a center frequency of 0.35.

These results demonstrate that in all cases the error is nearly 0 at the center frequency and increases roughly quadratically away from the center frequency. The error tends to be smallest for single-mode waveguides because a localized source excitation couples most strongly into guided modes. Note that in this case the maximum error is ~1% for a source bandwidth that is 67% of its center frequency. For the multi-mode waveguide, a much larger scattering loss is obtained for the higher-order mode A at frequencies below the center frequency, but this is simply because that mode ceases to be guided around a frequency ≈ 0.3, and the mode pattern changes dramatically as this cutoff is approached.

Another thing to keep in mind is that, even if the modes are imperfectly launched (some power leaks into radiation or into the backward direction), the correct result involving the Poynting flux can still be obtained by the standard procedure of normalizing against a separate straight-waveguide run (and subtracting off fields before computing reflected fluxes), as explained in Introduction. More accurate mode launching may be required for multi-mode waveguide systems, however, as power coupled into an undesired guided mode will not typically decay away.

Oblique Waveguides#

The eigenmode source can also be used to launch modes in an oblique/rotated waveguide. The results are shown in the two figures below for the single- and multi-mode case. There is one subtlety: for mode A in the multi-mode case, the bnum parameter is set to 3 rather than 2. This is because a non-zero rotation angle breaks the symmetry in the -direction which therefore precludes the use of EVEN_Y in eig_parity. Without any parity specified for the -direction, the second band corresponds to odd modes. This is why we must select the third band which contains even modes.

Note that an oblique waveguide leads to a breakdown in the PML. A simple workaround for mitigating the PML reflection artifacts in this case is to double the thickness from 1 to 2.

There are numerical dispersion artifacts due to the FDTD spatial and temporal discretizations which create negligible backward-propagating waves by the eigenmode current source, carrying approximately 10-5 of the power of the desired forward-propagating mode. These artifacts can be seen as residues in the field profiles.

We can demonstrate that the total power in a waveguide with arbitrary orientation — computed using two equivalent methods via get_fluxes and mode decomposition — can be computed by a single flux plane oriented along the direction: thanks to Poynting's theorem, the flux through any plane crossing a lossless waveguide is the same, regardless of whether the plane is oriented perpendicular to the waveguide. Furthermore, the eigenmode source is normalized in such a way as to produce the same power regardless of the waveguide orientation — in consequence, the flux values for mode A of the single-mode case for rotation angles of 0°, 20°, and 40° are 1111.280794, 1109.565028, and 1108.759159, within 0.2% (discretization error) of one another. Note that the Poynting flux could have been normalized to unity by setting the EigenModeSource/Source object parameter amplitude=1/src.fourier_transform(fsrc) where fsrc=0.15 and src=mp.GaussianSource(fsrc,fwidth=0.2*fsrc).

Finally, we demonstrate that as long as the line source intersects the waveguide and eig_kpoint is not nearly parallel to the direction of the line source, the mode can be properly launched. As shown in the field profiles below for the single-mode waveguide, there does not seem to be any noticeable distortion in the launched mode as the waveguide approaches glancing incidence to the source plane up to 80°, where the total power in the forward-propagating mode is 97%. Note that the line source spans the entire length of the cell extending into the PML region (not shown). In this example where the cell length is 10 μm (or 10X the width of the waveguide), the maximum rotation angle is ~84°, where the power drops to 59% and backward-propagating fields are clearly visible.

Increasing the size of the cell improves results at the expense of a larger simulation. The field profiles shown below are for a cell where the length has been doubled to 20 μm. The waveguide power at 84° increases from 59% to 80%. However, as the waveguide mode approaches glancing incidence, sensitivity to discretization errors increases because the mode varies rapidly with frequency on a glancing-angle cross-section, and you will eventually need to increase the resolution as well as the cell size. For waveguide angles much beyond 45° you probably want to simply change the orientation of the line source by 90°.

Planewaves in Homogeneous Media#

The eigenmode source can also be used to launch planewaves in homogeneous media. The dispersion relation for a planewave is where is the angular frequency of the planewave and its wavevector; is the refractive index of the homogeneous medium ( in Meep units). This example demonstrates launching planewaves in a uniform medium with at four incident angles: 0°, 40°, 120°, and 210°.

The simulation script is in examples/oblique-planewave.py. The notebook is in examples/oblique-planewave.ipynb.

import matplotlib.pyplot as plt
import meep as mp
import numpy as np


mp.verbosity(2)

resolution_um = 50
pml_um = 2.0
size_um = 10.0
cell_size = mp.Vector3(size_um + 2 * pml_um, size_um, 0)
pml_layers = [mp.PML(thickness=pml_um, direction=mp.X)]

# Incident angle of planewave. 0 is +x with rotation in
# counter clockwise (CCW) direction around z axis.
incident_angle = np.radians(40.0)

wavelength_um = 1.0
frequency = 1 / wavelength_um

n_mat = 1.5  # refractive index of homogeneous material
default_material = mp.Medium(index=n_mat)

k_point = mp.Vector3(n_mat * frequency, 0, 0).rotate(
    mp.Vector3(0, 0, 1),
    incident_angle
)

if incident_angle == 0:
    direction = mp.AUTOMATIC
    eig_parity = mp.EVEN_Y + mp.ODD_Z
    symmetries = [mp.Mirror(mp.Y)]
    eig_vol = None
else:
    direction = mp.NO_DIRECTION
    eig_parity = mp.ODD_Z
    symmetries = []
    eig_vol = mp.Volume(
        center=mp.Vector3(),
        size=mp.Vector3(0, 1 / resolution_um, 0)
    )

sources = [
    mp.EigenModeSource(
        src=mp.ContinuousSource(frequency),
        center=mp.Vector3(),
        size=mp.Vector3(0, size_um, 0),
        direction=direction,
        eig_kpoint=k_point,
        eig_band=1,
        eig_parity=eig_parity,
        eig_vol=eig_vol
    )
]

sim = mp.Simulation(
    cell_size=cell_size,
    resolution=resolution_um,
    boundary_layers=pml_layers,
    sources=sources,
    k_point=k_point,
    default_material=default_material,
    symmetries=symmetries
)

sim.run(until=23.56)

output_plane = mp.Volume(
    center=mp.Vector3(),
    size=mp.Vector3(size_um, size_um, 0)
)

fig, ax = plt.subplots()
sim.plot2D(fields=mp.Ez, output_plane=output_plane, ax=ax)
fig.savefig("planewave_source.png", bbox_inches="tight")

There are several items to note in this script.

  1. The line source spans the entire length of the cell. Planewave sources extending into the PML region must include is_integrated=True in the source definition.
  2. The wavevector of the planewave is defined using the k_point which specifies Bloch-periodic boundaries. PML boundaries are used only along the direction.
  3. The eig_vol parameter of the EigenModeSource defines a volume that is one-pixel thick in the direction of the line monitor. This is necessary to ensure that theeig_kpoint is interpreted by MPB as the wavevector of a planewave — that is, so it is the wavevector as defined for a homogeneous material (continuous translational symmetry) rather than as a Bloch wavevector for a system with discrete translational symmetry.
  4. Based on the convention used to define the incident angle whereby 0° is the direction with rotation counter clockwise around the axis, a planewave with an incident angle of 120° is propagating towards the (-x, +y) direction and 210° is propagating towards the (-x, -y) direction.

This example involves a continuous-wave (CW) time profile. For a pulse profile, the oblique planewave is incident at a given angle for only a single frequency component of the source as described in Tutorial/Basics/Angular Reflectance Spectrum of a Planar Interface. This is a fundamental feature of FDTD simulations and not of Meep per se. To simulate an incident planewave at multiple angles for a given frequency , you will need to perform separate simulations involving different values of (k_point) since each set of specifying the Bloch-periodic boundaries and the frequency of the source will produce a different angle of the planewave. For more details, refer to Section 4.5 ("Efficient Frequency-Angle Coverage") in Chapter 4 ("Electromagnetic Wave Source Conditions") of Advances in FDTD Computational Electrodynamics: Photonics and Nanotechnology.

Shown below are the steady-state field profiles generated by the planewave for the four incident angles. Residues of the backward-propagating waves due to the discretization are slightly visible.