How to use the mode solver in Tidy3D FDTD¶
This tutorial shows how to use the mode solver plugin in Tidy3D.
# standard python imports
import numpy as np
import matplotlib.pylab as plt
# tidy3D import
import tidy3d as td
from tidy3d.constants import C_0
import tidy3d.web as web
from tidy3d.plugins.mode import ModeSolver
Setup¶
We first set up the mode solver with information about our system. We start by setting parameters
# size of simulation domain
Lx, Ly, Lz = 6, 6, 6
dl = 0.05
# waveguide information
wg_width = 1.5
wg_height = 1.0
wg_permittivity = 4.0
# central frequency
wvl_um = 2.0
freq0 = C_0 / wvl_um
fwidth = freq0 / 3
# run_time in ps
run_time = 1e12
# automatic grid specification
grid_spec = td.GridSpec.auto(min_steps_per_wvl=20, wavelength=wvl_um)
Then we set up a simulation, in this case including a straight waveguide and periodic boundary conditions. Note that Tidy3D warns us that we have not added any sources in our Simulation
object, however for purposes of mode solving it is not necessary.
waveguide = td.Structure(
geometry=td.Box(size=(wg_width, td.inf, wg_height)),
medium=td.Medium(permittivity=wg_permittivity),
)
sim = td.Simulation(
size=(Lx, Ly, Lz),
grid_spec=grid_spec,
structures=[waveguide],
run_time=run_time,
boundary_spec=td.BoundarySpec.all_sides(boundary=td.Periodic()),
)
ax = sim.plot(z=0)
plt.show()
Initialize Mode Solver¶
With our system defined, we can now create our mode solver. We first need to specify on what plane we want to solve the modes using a td.Box()
object.
plane = td.Box(center=(0, 0, 0), size=(4, 0, 3.5))
The mode solver can now compute the modes given a ModeSpec
object that specifies everything about the modes we're looking for, for example:

num_modes
: how many modes to compute. 
target_neff
: float, default=None, initial guess for the effective index of the mode; if not specified, the modes with the largest real part of the effective index are computed.
The full list of specification parameters can be found here.
mode_spec = td.ModeSpec(
num_modes=3,
target_neff=2.0,
)
We can also specify a list of frequencies at which to solve for the modes.
num_freqs = 11
f0_ind = num_freqs // 2
freqs = np.linspace(freq0  fwidth / 2, freq0 + fwidth / 2, num_freqs)
Finally, we can initialize the ModeSolver
, and call the solve method.
mode_solver = ModeSolver(
simulation=sim,
plane=plane,
mode_spec=mode_spec,
freqs=freqs,
)
mode_data = mode_solver.solve()
Visualizing Mode Data¶
The mode_info
object contains information about the effective index of the mode and the field profiles, as well as the mode_spec
that was used in the solver. The effective index data and the field profile data is in the form of xarray DataArrays.
We can for example plot the real part of the effective index for all three modes as follows.
fig, ax = plt.subplots(1)
n_eff = mode_data.n_eff # real part of the effective mode index
n_eff.plot.line(x="f")
plt.show()
The raw data can also be accessed.
n_complex = mode_data.n_complex # complex effective index as a DataArray
n_eff = mode_data.n_eff.values # real part of the effective index as numpy array
k_eff = mode_data.k_eff.values # imag part of the effective index as numpy array
print(
f"first mode effective index at freq0: n_eff = {n_eff[f0_ind, 0]:.2f}, k_eff = {k_eff[f0_ind, 0]:.2e}"
)
first mode effective index at freq0: n_eff = 1.77, k_eff = 0.00e+00
The fields stored in mode_data
can be visualized using inbuilt xarray methods.
f, (ax1, ax2) = plt.subplots(1, 2, tight_layout=True, figsize=(10, 3))
abs(mode_data.Ex.isel(mode_index=0, f=f0_ind)).plot(x="x", y="z", ax=ax1, cmap="magma")
abs(mode_data.Ez.isel(mode_index=0, f=f0_ind)).plot(x="x", y="z", ax=ax2, cmap="magma")
ax1.set_title("Ex(x, y)")
ax1.set_aspect("equal")
ax2.set_title("Ez(x, y)")
ax2.set_aspect("equal")
plt.show()
Alternatively, we can use the inbuilt plot_field
method of mode_data
, which also allows us to overlay the structures in the simulation. The image also looks slightly different because the plot_field
method uses robust=True
option by default, which scales the colorbar to between the 2nd and 98th percentile of the data.
f, (ax1, ax2) = plt.subplots(1, 2, tight_layout=True, figsize=(10, 3))
mode_solver.plot_field("Ex", "abs", mode_index=0, f=freq0, ax=ax1)
mode_solver.plot_field("Ez", "abs", mode_index=0, f=freq0, ax=ax2)
plt.show()
Choosing the mode of interest¶
We can also look at the other modes that were computed.
mode_index = 1
f, (ax1, ax2) = plt.subplots(1, 2, tight_layout=True, figsize=(10, 3))
mode_solver.plot_field("Ex", "abs", mode_index=mode_index, f=freq0, ax=ax1)
mode_solver.plot_field("Ez", "abs", mode_index=mode_index, f=freq0, ax=ax2)
plt.show()
This looks like an Ezdominant mode. Finally, nextorder mode has mixed polarization.
mode_index = 2
f, (ax1, ax2) = plt.subplots(1, 2, tight_layout=True, figsize=(10, 3))
mode_solver.plot_field("Ex", "abs", mode_index=mode_index, f=freq0, ax=ax1)
mode_solver.plot_field("Ez", "abs", mode_index=mode_index, f=freq0, ax=ax2)
plt.show()
Exporting Results¶
This looks promising!
Now we can choose the mode specifications to use in our mode source and mode monitors. These can be created separately, can be exported directly from the mode solver, for example:
# Makes a modal source with geometry of `plane` with modes specified by `mode_spec` and a selected `mode_index`
source_time = td.GaussianPulse(freq0=freq0, fwidth=fwidth)
mode_src = mode_solver.to_source(mode_index=2, source_time=source_time, direction="")
# Makes a mode monitor with geometry of `plane`.
mode_mon = mode_solver.to_monitor(name="mode", freqs=freqs)
# Offset the monitor along the propagation direction
mode_mon = mode_mon.copy(update=dict(center=(0, 2, 0)))
# Inplane field monitor, slightly offset along x
monitor = td.FieldMonitor(
center=(0, 0, 0.1), size=(td.inf, td.inf, 0), freqs=[freq0], name="field"
)
sim = td.Simulation(
size=(Lx, Ly, Lz),
grid_spec=grid_spec,
run_time=run_time,
boundary_spec=td.BoundarySpec.all_sides(boundary=td.PML()),
structures=[waveguide],
sources=[mode_src],
monitors=[monitor, mode_mon],
)
sim.plot(z=0)
plt.show()
job = web.Job(simulation=sim, task_name="mode_simulation", verbose=True)
sim_data = job.run(path="data/simulation_data.hdf5")
Output()
Output()
Output()
Output()
Output()
We can now plot the inplane field and the modal amplitudes. Since we injected mode 2 and we just have a straight waveguide, all the power recorded by the modal monitor is in mode 2, going backwards.
fig, ax = plt.subplots(1, 2, figsize=(10, 4))
sim_data.plot_field("field", "Ez", f=freq0, ax=ax[0])
sim_data["mode"].amps.sel(direction="").abs.plot.line(x="f", ax=ax[1])
plt.show()
Storing serverside computed modes¶
We can also use a ModeSolverMonitor
to store the modes as they are computed serverside. This is illustrated below. We will also request in the mode specification that the modes are filtered by their tm
polarization. In this particular simulation, TM refers to Ez
polarization. The effect of the filtering is that modes with a tm
polarization fraction larger than or equal to 0.5 will come first in the list of modes (while still ordered by decreasing effective index). After that, the set of predominantly te
polarized modes (tm
fraction < 0.5) follows.
mode_spec = mode_spec.copy(update=dict(filter_pol="tm"))
# Update mode source to use the highesttmfraction mode
mode_src = mode_src.copy(update=dict(mode_spec=mode_spec))
mode_src = mode_src.copy(update=dict(mode_index=0))
# Update mode monitor to use the tm_fraction ordered mode_spec
mode_mon = mode_mon.copy(update=dict(mode_spec=mode_spec))
# New monitor to record the modes computed at the mode decomposition monitor location
mode_solver_mon = td.ModeSolverMonitor(
center=mode_mon.center,
size=mode_mon.size,
freqs=mode_mon.freqs,
mode_spec=mode_spec,
name="mode_solver",
)
sim = td.Simulation(
size=(Lx, Ly, Lz),
grid_spec=grid_spec,
run_time=run_time,
boundary_spec=td.BoundarySpec.all_sides(boundary=td.PML()),
structures=[waveguide],
sources=[mode_src],
monitors=[monitor, mode_mon, mode_solver_mon],
)
job = web.Job(simulation=sim, task_name="mode_simulation", verbose=True)
sim_data = job.run(path="data/simulation_data.hdf5")
Output()
Output()
Output()
Output()
Output()
Note the different ordering of the recorded modes compared to what we saw above, with the fundamental TM mode coming first.
fig, ax = plt.subplots(1)
n_eff = sim_data["mode"].n_eff # real part of the effective mode index
n_eff.plot.line(x="f")
plt.show()
Now the fundamental Ezpolarized mode is injected, and as before it is the only one that the mode monitor records any intensity in.
fig, ax = plt.subplots(1, 2, figsize=(10, 4))
sim_data.plot_field("field", "Ez", f=freq0, ax=ax[0])
sim_data["mode"].amps.sel(direction="").abs.plot.line(x="f", ax=ax[1])
plt.show()
We can also have a look at the mode fields stored in the ModeFieldMonitor
either directly using xarray methods as above, or using the Tidy3D SimulationData
inbuilt field plotting.
fig, ax = plt.subplots(1, 2, figsize=(12, 4))
sim_data.plot_field("mode_solver", "Ex", f=freq0, val="abs", mode_index=0, ax=ax[0])
sim_data.plot_field("mode_solver", "Ez", f=freq0, val="abs", mode_index=0, ax=ax[1])
plt.show()
Mode tracking¶
As typical for eigenvalue type solvers the Tidy3D's mode solver returns calculated modes in the order corresponding to the magnitudes of the found eigenvalues (effective index). Since the effective index of a mode generally depends on frequency and can become larger or smaller in magnitude compared to the effective index values of other modes, the same mode_index
value in the returned mode solver data may correspond to physically different modes for different frequencies. To avoid such a mismatch ModeSolver
by default performs sorting of modes based on their overlap values (see discussion on mode decomposition), that is, modes at one frequency are matched to the most similar modes at the other frequency. This behavior is controlled by the parameter track_freq
of ModeSpec
object. It can be set to either None
, "lowest"
, "central"
(default), or "highest"
, where in the first case no sorting is performed while in the other three cases sorting is performed starting from the specified frequency. Additionally, any unsorted mode solver data (objects of type ModeSolverData
) returned by ModeSolver.solve()
can be sorted afterwards by using function ModeSolverData.overlap_sort()
. This function takes in arguments track_freq
(default: "central"
), which has the same meaning as in ModeSpec
, and overlap_thresh
(default: 0.9) that defines the threshold overlap value above which a pair of spatial fields is considered to correspond to physically the same mode. While this functionality is not expected to be used in the most cases, we still demonstrate it here to provide a more clear understanding of ModeSolver
results. Note that for mode sources and mode monitors sorting is always performed internally.
We start with calculating twelve modes at eleven different frequencies for a layered waveguide with the default sorting turned off.
# bottom layer
bottom = td.Structure(
geometry=td.Box(
center=(0, 0, 0.8 * wg_height), size=(wg_width, td.inf, wg_height)
),
medium=td.Medium(permittivity=wg_permittivity * 2),
)
# top layer
top = td.Structure(
geometry=td.Box(
center=(0, 0, 0.7 * wg_height), size=(wg_width, td.inf, wg_height * 0.8)
),
medium=td.Medium(permittivity=wg_permittivity * 2),
)
# new simulation object
sim = td.Simulation(
size=(Lx, Ly, Lz),
grid_spec=grid_spec,
structures=[waveguide, bottom, top],
run_time=run_time,
boundary_spec=td.BoundarySpec.all_sides(boundary=td.Periodic()),
)
# visualize
sim.plot(y=0)
plt.show()
Again, Tidy3D is warning us that the simulation does not contain a source. However, since this simulation is used to construct the mode solver and will not be run directly, we can ignore this warning.
# set mode spec to find 12 modes and turn off automatic mode tracking
mode_spec = td.ModeSpec(num_modes=12, track_freq=None)
# we will use a larger solver plane for this structure
plane = td.Box(center=(0, 0, 0), size=(10, 0, 10))
# 11 frequencies in +/ fwidth range
num_freqs = 11
freqs = np.linspace(freq0  fwidth, freq0 + fwidth, num_freqs)
# make new mode solver object and solve for modes
mode_solver = ModeSolver(
simulation=sim,
plane=plane,
mode_spec=mode_spec,
freqs=freqs,
)
mode_data = mode_solver.solve()
By inspecting the values of effective indices we note that some of them come very close to each other, which is a sign that some mode index values might not correspond to the same modes throughout the considered frequency range.
fig, ax = plt.subplots(1)
n_eff = mode_data.n_eff # real part of the effective mode index
n_eff.plot.line(".", x="f")
plt.show()
Indeed, by plotting field distributions for the calculated modes we see mixing of modes for mode_index
> 4.
fig, ax = plt.subplots(
mode_spec.num_modes,
len(freqs),
tight_layout=True,
figsize=(len(freqs), mode_spec.num_modes),
)
for j in range(mode_spec.num_modes):
for i, freq in enumerate(freqs):
ax[j, i].imshow(
mode_data.Ex.isel(mode_index=j, f=i, drop=True).squeeze().abs.T,
cmap="magma",
origin="lower",
)
ax[j, i].axis("off")
plt.show()