Mode sources and mode monitors
Contents
Mode sources and mode monitors#
Here, we look at a simple demonstration of how to excite a specific waveguide mode, and how to decompose the fields recorded in a monitor on the basis of waveguide modes, i.e. how to compute the power carried in each mode.
[1]:
# standard python imports
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
# tidy3D import
import tidy3d as td
from tidy3d import web
[22:14:25] WARNING This version of Tidy3D was pip installed from the 'tidy3d-beta' repository on __init__.py:103 PyPI. Future releases will be uploaded to the 'tidy3d' repository. From now on, please use 'pip install tidy3d' instead.
INFO Using client version: 1.9.0rc1 __init__.py:121
Straight waveguide simulation#
First, we will do a simulation of a straight waveguide, using a silicon ridge waveguide on a silicon oxide substrate. We begin by defining some general parameters.
[2]:
# Unit length is micron.
wg_height = 0.22
wg_width = 0.45
# Permittivity of waveguide and substrate
si_eps = 3.48**2
sio_eps = 1.45**2
# Free-space wavelength (in um) and frequency (in Hz)
lambda0 = 1.55
freq0 = td.C_0 / lambda0
fwidth = freq0 / 10
# Simulation size inside the PML along propagation direction
sim_length = 5
# Simulation domain size and total run time
sim_size = [sim_length, 4, 2]
run_time = 20 / fwidth
# Grid specification
grid_spec = td.GridSpec.auto(min_steps_per_wvl=20, wavelength=lambda0)
Initialize structures, mode source, and mode monitor#
When initializing ModeSource and ModeMonitor objects, one of the three values of the size parameter must be zero. This implicitly defines the propagation direction for the mode decomposition. In this example, the waveguide is oriented along the x-axis, and the mode is injected along the positive-x direction (“forward”). Below, we add a mode monitor that will show us the waveguide transmission at a range of frequencies, as well as a simple frequency monitor to examine the fields in
the xy-plane at the central frequency. Note that we use a rather wide frequency range in the flux and modal monitors. This will be used later to demonstrate the effect of mode mismatch away from the central frequency and discuss ways to remedy that.
[3]:
# Waveguide and substrate materials
mat_wg = td.Medium(permittivity=si_eps)
mat_sub = td.Medium(permittivity=sio_eps)
# Substrate
substrate = td.Structure(
geometry=td.Box(
center=[0, 0, -sim_size[2]],
size=[td.inf, td.inf, 2 * sim_size[2]],
),
medium=mat_sub,
)
# Waveguide
waveguide = td.Structure(
geometry=td.Box(
center=[0, 0, wg_height / 2],
size=[100, wg_width, wg_height],
),
medium=mat_wg,
)
# Modal source parameters
src_pos = -sim_size[0] / 2 + 0.5
src_plane = td.Box(center=[src_pos, 0, 0], size=[0, 3, 2])
# xy-plane frequency-domain field monitor at central frequency
freq_mnt = td.FieldMonitor(
center=[0, 0, wg_height / 2], size=[100, 100, 0], freqs=[freq0], name="field"
)
# frequencies
mon_plane = td.Box(center=[-src_pos, 0, 0], size=[0, 3, 2])
Nfreqs = 17
freqs = np.linspace(freq0 - 2 * fwidth, freq0 + 2 * fwidth, Nfreqs)
fcent_ind = Nfreqs // 2 # index of the central frequency
# flux monitor
flux_mnt = td.FluxMonitor(
center=mon_plane.center, size=mon_plane.size, freqs=list(freqs), name="flux"
)
# Modal monitor at a range of frequencies
mode_mnt = td.ModeMonitor(
center=mon_plane.center,
size=mon_plane.size,
freqs=list(freqs),
mode_spec=td.ModeSpec(num_modes=3),
name="mode",
)
Initialize simulation and visualize two cross-sections to make sure we have set up the device correctly.
[4]:
# Simulation
sim = td.Simulation(
size=sim_size,
grid_spec=grid_spec,
structures=[substrate, waveguide],
sources=[],
monitors=[freq_mnt, mode_mnt, flux_mnt],
run_time=run_time,
boundary_spec=td.BoundarySpec.all_sides(boundary=td.PML()),
)
fig, (ax1, ax2) = plt.subplots(1, 2, tight_layout=True, figsize=(11, 4))
sim.plot(z=0, ax=ax1)
sim.plot(y=0, ax=ax2)
WARNING No sources in simulation. simulation.py:494
[4]:
<AxesSubplot: title={'center': 'cross section at y=0.00'}, xlabel='x', ylabel='z'>
Mode source#
Before we can run a simulation with a mode source, we have to select which of the eigenmodes we would like to inject. To do that, we can first visualize the modes using the in-built eigenmode solver and plotting functions. The modes are computed at the central frequency of the source, and in order of decreasing effective index n, such that the modes that are fully below light-line (if any) should appear first. The solver assumes periodic boundary conditions at the boundaries of the 2D plane.
Thus, for accurate results, the plane should be large enough for the fields to decay at the edges.
[5]:
# Visualize the modes
from tidy3d.plugins import ModeSolver
mode_spec = td.ModeSpec(num_modes=3)
ms = ModeSolver(simulation=sim, plane=src_plane, mode_spec=mode_spec, freqs=[freq0])
modes = ms.solve()
print("Effective index of computed modes: ", np.array(modes.n_eff))
fig, axs = plt.subplots(3, 2, figsize=(12, 12))
for mode_ind in range(3):
ms.plot_field("Ey", "abs", f=freq0, mode_index=mode_ind, ax=axs[mode_ind, 0])
ms.plot_field("Ez", "abs", f=freq0, mode_index=mode_ind, ax=axs[mode_ind, 1])
/Users/twhughes/.pyenv/versions/3.10.9/lib/python3.10/site-packages/numpy/linalg/linalg.py:2139: RuntimeWarning: divide by zero encountered in det
r = _umath_linalg.det(a, signature=signature)
/Users/twhughes/.pyenv/versions/3.10.9/lib/python3.10/site-packages/numpy/linalg/linalg.py:2139: RuntimeWarning: invalid value encountered in det
r = _umath_linalg.det(a, signature=signature)
[22:14:26] WARNING Mode field at frequency index 0, mode index 2 does not decay at the plane mode_solver.py:354 boundaries.
Effective index of computed modes: [[2.333904 1.5509037 1.364464 ]]
WARNING No sources in simulation. simulation.py:494
The waveguide has a single guided TE mode, as well as a TM mode which is very close to the light line (effective index close to substrate index). Finally, the last mode shown here, is below the light-line of the substrate and is mostly localized in that region. However, modes like these should always be considered unphysical, because of the assumption is that they decay by the edges of the mode plane. Thus, for meaningful results, only Mode 0 and Mode 1 should be used by the mode source.
Run simulation#
We set the mode source to the fundamental TE mode. Then, we run the simulation as usual through the web API, wait for it to finish, and download and load the results.
[6]:
source_time = td.GaussianPulse(freq0=freq0, fwidth=fwidth)
mode_source = ms.to_source(mode_index=0, direction="+", source_time=source_time)
sim = sim.copy(update=dict(sources=[mode_source]))
[7]:
job = web.Job(simulation=sim, task_name="mode_tutorial")
sim_data = job.run(path="data/simulation.hdf5")
Let’s first examine the in-plane fields recorded by the frequency monitor. We can already see how the source emits all of its power in the desired direction and waveguide mode.
[8]:
sim_data.plot_field("field", "Ey", z=wg_height / 2, f=freq0, val="real")
plt.show()
Mode monitors#
Mode monitors allow us to decompose the frequency-domain fields recorded in a simulation into the propagating modes of a waveguide. Specifically, we can write the modes of the waveguide at angular frequency \(\omega\) that are propagating in the forward direction, i.e. in the positive x direction in the example above, as
where \(p\) is a discrete mode index, \(k_p = n_p \omega/c\) is the wave-vector, \(n_p\) is the effective index of the \(p\)-th mode, and superscript \(f\) specifies forward propagation. The fields in the backward direction can be obtained via reflection symmetry. In the axes of the simulation we have as \(k_p \rightarrow -k_p\) and \(\mathbf{E}_{p}^b(\omega) = (-E_{p,x}^{f}, E_{p,y}^{f}, E_{p,z}^{f})\), \(\mathbf{H}_{p}^b(\omega) = (H_{p,x}^{f}, -H_{p,y}^{f}, -H_{p,z}^{f})\).
The fields stored in a monitor can then be decomposed on the basis of these waveguide modes. Following 1, 2, we can define an inner product between fields in the 2D plane as
where \(\mathbf{u} = (\mathbf{E}, \mathbf{H})\) combines both electromagnetic fields, the integration is over the plane area \(A\), and \(\mathrm{d}\mathbf{A}\) is the surface normal. If a waveguide mode is normalized such that \((\mathbf{u}_p, \mathbf{u}_p) = 1\), and we denote the fields stored in the mode monitor by \(\mathbf{u}_s\), then the power amplitude carried by mode \(p\) is given by the complex coefficient
while the power is given by \(|c_p|^2\).
We can have a look at the monitor modes, but as expected they should be identical to the source modes.
[9]:
# Visualize the monitor modes
ms = td.plugins.ModeSolver(
simulation=sim, plane=mode_mnt.geometry, mode_spec=mode_mnt.mode_spec, freqs=[freq0]
)
modes = ms.solve()
fig, axs = plt.subplots(3, 2, figsize=(12, 12))
for mode_ind in range(3):
ms.plot_field("Ey", "abs", f=freq0, mode_index=mode_ind, ax=axs[mode_ind, 0])
ms.plot_field("Ez", "abs", f=freq0, mode_index=mode_ind, ax=axs[mode_ind, 1])
/Users/twhughes/.pyenv/versions/3.10.9/lib/python3.10/site-packages/numpy/linalg/linalg.py:2139: RuntimeWarning: divide by zero encountered in det
r = _umath_linalg.det(a, signature=signature)
/Users/twhughes/.pyenv/versions/3.10.9/lib/python3.10/site-packages/numpy/linalg/linalg.py:2139: RuntimeWarning: invalid value encountered in det
r = _umath_linalg.det(a, signature=signature)
[22:14:57] WARNING Mode field at frequency index 0, mode index 2 does not decay at the plane mode_solver.py:354 boundaries.
We note that in Tidy3D, the fields recorded by frequency monitors (and thus also mode monitors) are automatically normalized by the power amplitude spectrum of the source (for multiple sources, the user can select which source to use for the normalization). Furthermore, mode sources are normalized to inject exactly 1W of power at the central frequency.
[10]:
# Flux in the mode monitor (total power through the cross-section)
flux_wg = sim_data["flux"].flux
print("Flux at central frequency: ", flux_wg.isel(f=fcent_ind).values)
Flux at central frequency: 0.99999183
We can also use the mode amplitudes recorded in the mode monitor to reveal the decomposition of the radiated power into forward- and backward-propagating modes, respectively. As we would expect, all of the power is injected into the fundamental waveguide mode, in the forward direction. More precisely, this is true up to some numerical precision that decreases with increasing simulation resolution.
[11]:
# Forward and backward power amplitude coefficients
mode_amps = sim_data["mode"]
coeffs_f = mode_amps.amps.sel(direction="+")
coeffs_b = mode_amps.amps.sel(direction="-")
print(
"Power distribution at central frequency in first three modes",
)
print(" positive dir. ", np.abs(coeffs_f.isel(f=fcent_ind) ** 2).values)
print(" negative dir. ", np.abs(coeffs_b.isel(f=fcent_ind) ** 2).values)
# Free-space wavelength corresponding to the monitor frequencies
lambdas = td.C_0 / freqs
fig, ax = plt.subplots(1, figsize=(6, 4))
ax.plot(lambdas, np.abs(coeffs_f.values) ** 2)
ax.set_xlim([lambdas[-1], lambdas[0]])
ax.set_xlabel("Wavelength (um)")
ax.set_ylabel("Power in mode (W)")
ax.set_title("Mode decomposition (forward-propagating)")
ax.legend(["Mode 0", "Mode 1", "Mode 2"])
plt.show()
Power distribution at central frequency in first three modes
positive dir. [9.99991721e-01 5.50487488e-12 9.86647957e-14]
negative dir. [1.41365978e-09 1.88499613e-14 1.09869501e-14]
We can examine the frequency dependence of the results a bit more closely, and compare them to the total power flux, which can be computed for any frequency monitor. The flux is the area-integrated time-averaged Poynting vector and gives the (signed) total power flowing through the surface. The flux computation and the modal decomposition are done in separate monitors and in a completely different way, but because all the power is in the fundamental mode here, the flux matches really well the zero-mode power from the power decomposition.
[12]:
fig, ax = plt.subplots(1, figsize=(6, 4))
ax.plot(lambdas, flux_wg)
ax.plot(lambdas, np.abs(coeffs_f.sel(mode_index=0)) ** 2)
ax.set_xlim([lambdas[-1], lambdas[0]])
ax.set_xlabel("Wavelength (um)")
ax.set_ylabel("Power (W)")
ax.set_title("Power in mode monitor")
ax.legend(["Total", "Mode 0"])
plt.show()
As we already saw, at the central frequency, the source power is extremely well directed in the waveguide mode. Since the source mode is computed at the central frequency only, away from that it is not perfectly matched, leading to a small decrease of the total radiated power. In certain situations, it is even possible to observe injected power larger than one away from the central frequency. That said, we see that all the radiated power is still emitted into the desired waveguide mode, within the wavelength range of interest. For best accuracy when computing scattering parameters away from the central frequency, we can choose one of the following two options. The first option is to do a “normalization” run with a straight waveguide, like the one we just did, and normalize w.r.t. the computed flux to account for the small frequency dependence of the total radiated power. The other one is to use the Tidy3D’s broadband source functionality, which takes into account the frequency dependence of injected fields. Both approaches are illustrated below.
Waveguide junction (using normalization run)#
We repeat the simulation, but this time introduce a much bigger waveguide in the second half of the domain.
[13]:
# Output waveguide
wgout_width = 1.4
waveguide_out = td.Structure(
geometry=td.Box(
center=[2, 0, wg_height / 2],
size=[4, wgout_width, wg_height],
),
medium=mat_wg,
)
[14]:
sim_jct = td.Simulation(
size=sim_size,
grid_spec=grid_spec,
structures=[substrate, waveguide, waveguide_out],
sources=[mode_source],
monitors=[freq_mnt, mode_mnt, flux_mnt],
run_time=run_time,
boundary_spec=td.BoundarySpec.all_sides(boundary=td.PML()),
)
[15]:
fig = plt.figure(figsize=(11, 4))
gs = mpl.gridspec.GridSpec(1, 2, figure=fig, width_ratios=[1, 1.4])
ax1 = fig.add_subplot(gs[0, 0])
ax2 = fig.add_subplot(gs[0, 1])
sim_jct.plot(z=0.1, ax=ax1)
sim_jct.plot(y=0.1, ax=ax2)
[15]:
<AxesSubplot: title={'center': 'cross section at y=0.10'}, xlabel='x', ylabel='z'>
[16]:
job = web.Job(simulation=sim_jct, task_name="mode_tutorial")
sim_data_jct = job.run(path="data/mode_converter.hdf5")
[17]:
sim_data_jct.plot_field("field", "Ey", z=wg_height / 2, f=freq0, val="real")
plt.show()
This time, the output waveguide is multi-mode, and there is obviously some mode-mixing happening. We can use the mode monitor to exactly quantify this.
[18]:
# Forward and backward power amplitude coefficients normalized to the straight waveguide flux
mode_amps_jct = sim_data_jct["mode"]
coeffs_f_jct = mode_amps_jct.amps.sel(direction="+") / np.sqrt(flux_wg.data[:, None])
coeffs_b_jct = mode_amps_jct.amps.sel(direction="-") / np.sqrt(flux_wg.data[:, None])
fig, ax = plt.subplots(1, figsize=(6, 4))
ax.plot(lambdas, np.sum(np.abs(coeffs_f_jct.values) ** 2, axis=1))
ax.plot(lambdas, np.abs(coeffs_f_jct.values) ** 2)
ax.set_xlabel("Wavelength (um)")
ax.set_xlim([lambdas[-1], lambdas[0]])
ax.set_ylabel("Power in mode (W)")
ax.set_title("Mode decomposition (+ propagating)")
ax.legend(["Mode 0 + 1 + 2", "Mode 0", "Mode 1", "Mode 2"])
plt.show()
Because of the symmetry with respect to the \(y=0\) plane, the power in Mode 0 cannot be converted to Mode 1, but a significant amount of power is converted to Mode 3. Note also that the combined power in the computed modes is smaller than 1W. The missing part is most likely lost in scattering at the sharp waveguide interface.
Waveguide junction (using broadband source)#
To demonstrate the broadband source capabilities more clearly we will re-run both the straight waveguide and waveguide junction cases. To create a broadband source one needs to specify the number of points, num_freqs, that will be used for approximating the frequency dependence. Note that these points will be arranged into a Chebyshev grid, thus a very fast convergence of polynomial approximation is expected and only few points are enough to reach converged results. Another important aspect
is that these points are placed in the interval (freq0 - 1.5 * fwidth, freq0 + 1.5 * fwidth), where freq0 and fwidth are the source’s parameters, and using broadband sources is not expected to significantly improve results outside of this interval. Use freq0 and fwidth source parameters to adjust the frequency window where the best accuracy is desired.
[19]:
broadband_mode_source = mode_source.copy(update={"num_freqs": 7})
Duplicate previously defined simulations and update the source lists.
[20]:
sim_bb = sim.copy(update={"sources": [broadband_mode_source]})
sim_jct_bb = sim_jct.copy(update={"sources": [broadband_mode_source]})
First, let us re-run the straight waveguide case.
[21]:
job = web.Job(simulation=sim_bb, task_name="mode_tutorial")
sim_data_bb = job.run(path="data/simulation_bb.hdf5")
Plotting the total injected power and the power in injected mode one can see that a broadband source injects a chosen mode much more accurately throughout the target frequency range (+/- 1.5 * fwidth). That is, the total injected power is equal to 1 W at all frequencies and all of that power is concentrated in the chosen mode. Note that results degrade outside of the target range as expected.
[22]:
flux_wg_bb = sim_data_bb["flux"].flux
# Forward and backward power amplitude coefficients
mode_amps_bb = sim_data_bb["mode"]
coeffs_f_bb = mode_amps_bb.amps.sel(direction="+")
coeffs_b_bb = mode_amps_bb.amps.sel(direction="-")
fig, ax = plt.subplots(1, figsize=(6, 4))
ax.plot(lambdas, flux_wg)
ax.plot(lambdas, flux_wg_bb)
ax.plot(lambdas, np.abs(coeffs_f.sel(mode_index=0)) ** 2, "x--")
ax.plot(lambdas, np.abs(coeffs_f_bb.sel(mode_index=0)) ** 2, ".--")
ax.set_xlim([lambdas[-1], lambdas[0]])
ax.set_xlabel("Wavelength (um)")
ax.set_ylabel("Power (W)")
ax.set_title("Power in mode monitor")
ax.axvspan(
td.C_0 / (freq0 + 1.5 * fwidth),
td.C_0 / (freq0 - 1.5 * fwidth),
facecolor="k",
alpha=0.1,
)
ax.legend(
[
"Total (single freq)",
"Total (broadband)",
"Mode 0 (single freq)",
"Mode 0 (broadband)",
"+/- 1.5 * fwidth",
]
)
plt.show()
Let us re-run the waveguide junction example now and compare to previous results.
[23]:
job = web.Job(simulation=sim_jct_bb, task_name="mode_tutorial")
sim_data_jct_bb = job.run(path="data/mode_converter_bb.hdf5")
Here we plot the comparison for mode mixing results between three options: (1) the simulation containing a mode source that uses only one frequency point, (2) the same as first one, but normalized by flux derived from straight waveguide simulation, and (3) the simulation containing a broadband mode source that uses mutliple frequency points. As one can see, using the broadband source feature allows to obtain results of the same accuracy without running an additional normalization simulation.
[24]:
# Forward and backward power amplitude coefficients using broadband source
mode_amps_jct_bb = sim_data_jct_bb["mode"]
coeffs_f_jct_bb = mode_amps_jct_bb.amps.sel(direction="+")
coeffs_b_jct_bb = mode_amps_jct_bb.amps.sel(direction="-")
# Forward and backward power amplitude coefficients using single frequency source and without normalization
coeffs_f_jct_nonorm = mode_amps_jct.amps.sel(direction="+")
coeffs_b_jct_nonorm = mode_amps_jct.amps.sel(direction="-")
fig, ax = plt.subplots(1, 3, figsize=(18, 4))
for mode_index in range(3):
ax[mode_index].plot(
lambdas, np.abs(coeffs_f_jct_nonorm.sel(mode_index=mode_index)) ** 2, "x:"
)
ax[mode_index].plot(
lambdas, np.abs(coeffs_f_jct.sel(mode_index=mode_index)) ** 2, ".--"
)
ax[mode_index].plot(
lambdas, np.abs(coeffs_f_jct_bb.sel(mode_index=mode_index)) ** 2
)
ax[mode_index].set_xlabel("Wavelength (um)")
ax[mode_index].set_xlim([lambdas[-1], lambdas[0]])
ax[mode_index].set_ylabel("Power in mode (W)")
ax[mode_index].set_title(f"Mode {mode_index}, direction='+'")
ax[mode_index].axvspan(
td.C_0 / (freq0 + 1.5 * fwidth),
td.C_0 / (freq0 - 1.5 * fwidth),
facecolor="k",
alpha=0.1,
)
ax[mode_index].legend(
["Single freq", "Single freq + norm", "Broadband", "+/- 1.5 * fwidth"]
)
plt.show()
[ ]: