Using the mode solver for optical mode analysis

This tutorial shows how to use the mode solver plugin in Tidy3D.

# standard python imports
from __future__ import annotations

import matplotlib.pylab as plt
import numpy as np

# tidy3D import
import tidy3d as td
import tidy3d.web as web
from tidy3d.constants import C_0
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 = 1e-12

# 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()

We can also summarize useful mode information using to_dataframe(). Note that the group index was not computed; this can be included by setting group_index_step=True in the ModeSpec.

mode_data.to_dataframe()

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.75, k_eff = 0.00e+00

The fields stored in mode_data can be visualized using in-built 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 in-built 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 Ez-dominant mode. Finally, next-order 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)))

# In-plane 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")

15:35:29 UTC Created task 'mode_simulation' with resource_id
             'fdve-eeb45fbc-9afc-47a6-b120-a7cfe459b0b5' and task_type 'FDTD'.

             View task using web UI at
             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-eeb45fbc-9af
             c-47a6-b120-a7cfe459b0b5'.

             Task folder: 'default'.

15:35:30 UTC Estimated FlexCredit cost: 0.025. This assumes the FDTD solver runs
             for the full simulation time; if early shutoff is reached, the
             billed cost can be lower. Use 'web.real_cost(task_id)' to get the
             billed FlexCredit cost after a simulation run.

15:35:31 UTC status = queued

             To cancel the simulation, use 'web.abort(task_id)' or
             'web.delete(task_id)' or abort/delete the task in the web UI.
             Terminating the Python script will not stop the job running on the
             cloud.

15:35:43 UTC status = preprocess

15:35:46 UTC status = postprocess

15:35:50 UTC status = success

15:35:52 UTC starting up solver

15:35:53 UTC running solver

             early shutoff detected at 11%, exiting.

15:35:54 UTC status = success

             View simulation result at
             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-eeb45fbc-9af
             c-47a6-b120-a7cfe459b0b5'.

15:35:55 UTC Loading results from data/simulation_data.hdf5

We can now plot the in-plane 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 server-side computed modes

We can also use a ModeSolverMonitor to store the modes as they are computed server-side. 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(
        sort_spec=td.ModeSortSpec(
            filter_key="TM_fraction",
            filter_reference=0.5,
            filter_order="over",
        )
    )
)

# Update mode source to use the highest-tm-fraction 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")

15:35:56 UTC Created task 'mode_simulation' with resource_id
             'fdve-73a38dc0-dfb0-439f-87bf-d18610a1f9bb' and task_type 'FDTD'.

             View task using web UI at
             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-73a38dc0-dfb
             0-439f-87bf-d18610a1f9bb'.

             Task folder: 'default'.

15:35:57 UTC Estimated FlexCredit cost: 0.025. This assumes the FDTD solver runs
             for the full simulation time; if early shutoff is reached, the
             billed cost can be lower. Use 'web.real_cost(task_id)' to get the
             billed FlexCredit cost after a simulation run.

15:35:58 UTC status = queued

             To cancel the simulation, use 'web.abort(task_id)' or
             'web.delete(task_id)' or abort/delete the task in the web UI.
             Terminating the Python script will not stop the job running on the
             cloud.

15:36:12 UTC starting up solver

             running solver

15:36:15 UTC early shutoff detected at 7%, exiting.

             status = postprocess

15:36:17 UTC status = success

15:36:19 UTC View simulation result at
             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-73a38dc0-dfb
             0-439f-87bf-d18610a1f9bb'.

15:36:21 UTC Loading results from data/simulation_data.hdf5

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 Ez-polarized 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 in-built 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()

Subpixel averaging

The mode solver applies subpixel averaging automatically when the tidy3d-extras package is installed (pip install "tidy3d[extras]"). Subpixel averaging smooths the material interfaces on the Yee grid, which increases the computation accuracy, with the whole calculation running locally. Alternatively, the mode solver can also be run in the cloud with web.run(mode_solver, task_name="mode_solver"), in which case subpixel averaging is applied as well.

To see its effect, we can disable subpixel averaging by setting td.config.simulation.use_local_subpixel = False. We use the previous mode solver as an example:

td.config.simulation.use_local_subpixel = False

# solve on a copy of the mode solver, since results from previous solves are cached
mode_data = mode_solver.updated_copy().solve()

# restore subpixel averaging
td.config.simulation.use_local_subpixel = True

15:36:22 UTC WARNING: Use the remote mode solver with subpixel averaging for    
             better accuracy through 'tidy3d.web.run(...)' or the deprecated    
             'tidy3d.plugins.mode.web.run(...)'. Alternatively, you can install 
             the package 'tidy3d-extras' using 'pip install "tidy3d"' and set   
             'config.simulation.use_local_subpixel=True'.                       

fig, ax = plt.subplots(1)
n_eff = mode_data.n_eff
n_eff.plot.line(x="f")
plt.show()

To demonstrate the accuracy improvement, we can show a convergence plot comparing the results with and without subpixel averaging on a waveguide that doesn’t align with the mesh, in this case a dielectric rod.

resolutions = range(10, 31, 2)
n_eff = ([], [])

cylinder = td.Structure(
    geometry=td.Cylinder(axis=1, center=(0, 0, 0), radius=0.7, length=2 * Lz),
    medium=td.Medium(permittivity=wg_permittivity),
)

dummy_source = td.UniformCurrentSource(
    size=(0, 0, 0),
    source_time=td.GaussianPulse(freq0=freq0, fwidth=fwidth),
    polarization="Ex",
    current_amplitude_definition="total",
)

mode_spec = td.ModeSpec(
    num_modes=1,
    target_neff=2.0,
)


def create_mode_solver(resolution):
    grid_spec = td.GridSpec.auto(min_steps_per_wvl=resolution, wavelength=wvl_um)

    sim = td.Simulation(
        size=(Lx, Ly, Lz),
        grid_spec=grid_spec,
        structures=[cylinder],
        sources=[dummy_source],
        run_time=run_time,
        boundary_spec=td.BoundarySpec.all_sides(boundary=td.Periodic()),
    )

    return ModeSolver(
        simulation=sim,
        plane=td.Box(center=(0, 0, 0), size=(4, 0, 4)),
        mode_spec=mode_spec,
        freqs=[freq0],
    )


for resolution in resolutions:
    print("Solving for resolution =", resolution, flush=True)
    td.config.simulation.use_local_subpixel = True
    mode_solver = create_mode_solver(resolution)
    n_eff[0].append(mode_solver.solve().n_eff.item())

    td.config.simulation.use_local_subpixel = False
    mode_solver = create_mode_solver(resolution)
    n_eff[1].append(mode_solver.solve().n_eff.item())

# restore subpixel averaging
td.config.simulation.use_local_subpixel = True

Solving for resolution = 10
Solving for resolution = 12
Solving for resolution = 14
Solving for resolution = 16
Solving for resolution = 18
Solving for resolution = 20
Solving for resolution = 22
Solving for resolution = 24
Solving for resolution = 26
Solving for resolution = 28
Solving for resolution = 30

fig, ax = plt.subplots(1)
ax.plot(resolutions, n_eff[0], label="Subpixel on")
ax.plot(resolutions, n_eff[1], label="Subpixel off")
ax.set(xlabel="Resolution", ylabel="Propagation index")
ax.legend()
ax.grid()

Mode tracking

As typical for eigenvalue type solvers, the Tidy3D’s mode solver finds the modes at each frequency 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 could end up corresponding to physically different modes at different frequencies near mode crossings. To avoid such a mismatch, ModeSolver automatically tracks the modes across frequencies based on their overlap values (see discussion on mode decomposition), so that each mode_index corresponds to the same physical mode throughout the frequency range. This behavior is controlled by the parameter track_freq of the ModeSortSpec object, passed to ModeSpec through its sort_spec parameter. It can be set to either None, "lowest", "central" (default), or "highest", where in the first case no tracking is performed while in the other three cases tracking starts from the specified frequency. Additionally, any unsorted mode solver data (objects of type ModeSolverData) can be sorted afterwards by using the function ModeSolverData.overlap_sort(). Note that for mode sources and mode monitors sorting is always performed internally.

We calculate twelve modes at eleven different frequencies for a layered waveguide with the default tracking:

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

# find 12 modes; mode tracking across frequencies is enabled by default
mode_spec = td.ModeSpec(num_modes=12)

# 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, indicating mode crossings. Thanks to the automatic tracking, each mode_index still corresponds to the same physical mode 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 can confirm that the modes are consistently tracked across all frequencies.

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

Note that depending on the particular situation, the mode sorting, whether default or by hands, might not successfully resolve the tracking of all modes throughout the chosen frequency range. Possible reasons for that: - A mode is physically disappearing at a certain frequency point. - A mode is not included in the calculated number of modes, num_modes, for certain frequency points. In this case, increasing num_modes might alleviate the issue. - Presence of degenerate modes at certain frequencies.

Group Index

The mode solver can also be used to directly calculate the group index of the calculated mode. It is important to have mode tracking enabled, to avoid differentiating wrong mode pairs around mode crossings.

# Use the previous mode solver with the appropriate settings
mode_spec = td.ModeSpec(
    num_modes=12, sort_spec=td.ModeSortSpec(track_freq="central"), group_index_step=True
)
mode_solver = mode_solver.copy(update={"mode_spec": mode_spec})

mode_data = mode_solver.solve()
n_eff = mode_data.n_eff
n_group = mode_data.n_group

_, ax = plt.subplots(1, 2, figsize=(10, 5), tight_layout=True)

n_eff.plot.line(".-", x="f", ax=ax[0])
n_group.plot.line(".-", x="f", ax=ax[1])
ax[1].set_ylabel("Group Index")
plt.show()

Notes / Considerations

• This mode solver runs locally, which means it does not require credits to run.

• When the tidy3d-extras package is installed, the local mode solver applies subpixel averaging, so its results match the modes stored in a ModeSolverMonitor during a simulation run. Without tidy3d-extras, or with td.config.simulation.use_local_subpixel = False, subpixel averaging is not applied locally, and the results are less accurate.

• Symmetries are applied if they are defined in the simulation and the mode plane center sits on the simulation center.

• In cases when the number of modes/frequency points and/or the dimensions of mode solver plane are large and tracking of modes is not important then it might be beneficial to turn off the default mode sorting for increased computational efficiency.