Using ModeSimulation for optical mode analysis

This tutorial shows how to use ModeSimulation in Tidy3D.

For a tutorial on the legacy ModeSolver plugin, please refer to this notebook.

# standard python imports
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

Setup

We first set up the mode simulation with information about our system. We start by setting the 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)

We then define 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 to the Simulation object. However, for mode solving, this 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(y=0)
plt.show()

Initialize Mode Simulation

With our system defined, we can now create our mode simulation. We first need to specify the plane on which we want to solve the modes using a `td.Box() <https://docs.flexcompute.com/projects/tidy3d/en/latest/api/_autosummary/tidy3d.Box.html>`__ object.

plane = td.Box(center=(0, 0, 0), size=(4, 0, 3.5))

The mode simulation can now compute the modes given a `ModeSpec <https://docs.flexcompute.com/projects/tidy3d/en/latest/api/_autosummary/tidy3d.ModeSpec.html>`__ object that specifies everything about the modes we are 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 `ModeSimulation <https://docs.flexcompute.com/projects/tidy3d/en/latest/api/_autosummary/tidy3d.ModeSimulation.html>`__ and call the local run method.

mode_simulation = td.ModeSimulation.from_simulation(
    simulation=sim,
    plane=plane,
    mode_spec=mode_spec,
    freqs=freqs,
)
sim_data = mode_simulation.run_local()
mode_data = sim_data.modes

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 <https://docs.flexcompute.com/projects/tidy3d/en/latest/api/_autosummary/tidy3d.ModeSpec.html>`__.

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 used in the solver. The effective index data and the field profile data are 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 built-in 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 built-in plot_field method of sim_data, which also allows us to overlay the structures in the simulation. The image also looks slightly different because the plot_field method uses the robust=True option by default, which scales the colorbar between the 2nd and 98th percentiles of the data.

f, (ax1, ax2) = plt.subplots(1, 2, tight_layout=True, figsize=(10, 3))
sim_data.plot_field("Ex", "abs", mode_index=0, f=freq0, ax=ax1)
sim_data.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))
sim_data.plot_field("Ex", "abs", mode_index=mode_index, f=freq0, ax=ax1)
sim_data.plot_field("Ez", "abs", mode_index=mode_index, f=freq0, ax=ax2)
plt.show()

This looks like an Ez-dominant mode. Finally, the next-order mode has mixed polarization.

mode_index = 2
f, (ax1, ax2) = plt.subplots(1, 2, tight_layout=True, figsize=(10, 3))
sim_data.plot_field("Ex", "abs", mode_index=mode_index, f=freq0, ax=ax1)
sim_data.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, or they can be exported directly from the mode simulation, 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 = td.ModeSource(
    center=plane.center,
    size=plane.size,
    source_time=source_time,
    mode_spec=mode_spec,
    mode_index=2,
    direction="-",
)

# Makes a mode monitor with geometry of `plane`.
mode_mon = td.ModeMonitor(
    center=plane.center,
    size=plane.size,
    freqs=freqs,
    mode_spec=mode_spec,
    name="mode",
)
# 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")

# 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],
)

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

10:37:59 UTC Created task 'mode_simulation' with resource_id
             'fdve-1d3e2e4d-b94b-427d-9def-c66f48d01082' and task_type 'FDTD'.

             View task using web UI at
             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-1d3e2e4d-b94
             b-427d-9def-c66f48d01082'.

             Task folder: 'default'.

10:38:00 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.

10:38:01 UTC status = success

10:38:03 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 have just 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 <https://docs.flexcompute.com/projects/tidy3d/en/latest/api/_autosummary/tidy3d.ModeSolverMonitor.html>`__ to store the modes as they are computed on the server. This is illustrated below. We will also request in the mode specification that the modes are sorted by their tm polarization. In this particular simulation, TM refers to Ez polarization. The effect of the sorting is that modes with a tm polarization fraction greater than or equal to 0.5 come first in the list of modes, while still being ordered by decreasing effective index. After that, the predominantly te-polarized modes follow.

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

10:38:04 UTC Created task 'mode_simulation' with resource_id
             'fdve-dfc79572-4fa4-4a6c-9385-7080664d7b11' and task_type 'FDTD'.

             View task using web UI at
             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-dfc79572-4fa
             4-4a6c-9385-7080664d7b11'.

             Task folder: 'default'.

10:38:05 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.

10:38:06 UTC status = success

10:38:09 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 for which the mode monitor records any intensity.

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 <https://docs.flexcompute.com/projects/tidy3d/en/latest/api/_autosummary/tidy3d.SimulationData.html>`__ built-in 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()

Computing modes on the server

We can also run `ModeSimulation <https://docs.flexcompute.com/projects/tidy3d/en/latest/api/_autosummary/tidy3d.ModeSimulation.html>`__ on the server. We use the previous mode simulation as an example:

sim_data = web.run(mode_simulation, task_name="mode_simulation")
mode_data = sim_data.modes

10:38:10 UTC Created task 'mode_simulation' with resource_id
             'mos-0d39805a-d360-4e54-ab57-36e9faa2fd5d' and task_type 'MODE'.

10:38:11 UTC Estimated FlexCredit cost: 0.005. For this solver type, the
             estimate is the final billed cost.

10:38:12 UTC status = success

10:38:14 UTC Loading results from simulation_data.hdf5

             WARNING: updating ModeSimulation from 2.12.0.dev2 to 2.12.0.dev3   

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 local and remote solvers on a waveguide that does not align with the mesh, in this case a dielectric rod.

For more detailed information about local subpixel, please refer to this example.

td.config.simulation.use_local_subpixel = False

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",
)

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


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

    return td.ModeSimulation(
        size=(Lx, Ly, Lz),
        plane=td.Box(center=(0, 0, 0), size=(4, 0, 4)),
        grid_spec=grid_spec,
        structures=[cylinder],
        boundary_spec=td.BoundarySpec.all_sides(boundary=td.Periodic()),
        mode_spec=mode_spec,
        freqs=[freq0],
    )


for resolution in resolutions:
    print("Solving for resolution =", resolution, flush=True)
    mode_simulation = create_mode_simulation(resolution)
    local_sim_data = mode_simulation.run_local()
    local_data = local_sim_data.modes
    n_eff[0].append(local_data.n_eff.item())

    mode_simulation = create_mode_simulation(resolution)
    server_sim_data = web.run(mode_simulation, task_name="mode_simulation", verbose=False)
    server_data = server_sim_data.modes
    n_eff[1].append(server_data.n_eff.item())

             WARNING: The default value of 'current_amplitude_definition' for   
             'UniformCurrentSource' will change from 'density' to 'total' in a  
             future release. To avoid this warning and ensure consistent        
             behavior, please explicitly set                                    
             current_amplitude_definition='density' or                          
             current_amplitude_definition='total' when creating the source.     

Solving for resolution = 10

10:38:18 UTC WARNING: updating ModeSimulation from 2.12.0.dev2 to 2.12.0.dev3   

Solving for resolution = 12

10:38:21 UTC WARNING: updating ModeSimulation from 2.12.0.dev2 to 2.12.0.dev3   

Solving for resolution = 14

10:38:25 UTC WARNING: updating ModeSimulation from 2.12.0.dev2 to 2.12.0.dev3   

Solving for resolution = 16

10:38:29 UTC WARNING: updating ModeSimulation from 2.12.0.dev2 to 2.12.0.dev3   

Solving for resolution = 18

10:38:33 UTC WARNING: updating ModeSimulation from 2.12.0.dev2 to 2.12.0.dev3   

Solving for resolution = 20

10:39:01 UTC WARNING: updating ModeSimulation from 2.12.0.dev2 to 2.12.0.dev3   

Solving for resolution = 22

10:39:05 UTC WARNING: updating ModeSimulation from 2.12.0.dev2 to 2.12.0.dev3   

Solving for resolution = 24

10:39:09 UTC WARNING: updating ModeSimulation from 2.12.0.dev2 to 2.12.0.dev3   

Solving for resolution = 26

10:39:12 UTC WARNING: updating ModeSimulation from 2.12.0.dev2 to 2.12.0.dev3   

Solving for resolution = 28

10:39:16 UTC WARNING: updating ModeSimulation from 2.12.0.dev2 to 2.12.0.dev3   

Solving for resolution = 30

10:39:20 UTC WARNING: updating ModeSimulation from 2.12.0.dev2 to 2.12.0.dev3   

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

Local Subpixel Comparison

Below, we enable local subpixel averaging and repeat the local mode simulation sweep. This allows a direct comparison between:

• local without subpixel,

• local with subpixel,

• server result.

With local subpixel enabled, the local result matches the server result.

Note that it is necessary to install Tidy3D extras with pip install "tidy3d[extras]".

# Repeat the local sweep with local subpixel enabled
td.config.simulation.use_local_subpixel = True

n_eff_local_subpixel = []

for resolution in resolutions:
    print("Solving with local subpixel for resolution =", resolution, flush=True)
    mode_simulation = create_mode_simulation(resolution)
    local_subpixel_sim_data = mode_simulation.run_local()
    local_subpixel_data = local_subpixel_sim_data.modes
    n_eff_local_subpixel.append(local_subpixel_data.n_eff.item())

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

fig, ax = plt.subplots(1)
ax.plot(resolutions, n_eff[0], label="Local")
ax.plot(resolutions, n_eff[1], label="Server")
ax.plot(resolutions, n_eff_local_subpixel, "--", label="Local + subpixel")
ax.set(xlabel="Resolution", ylabel="Propagation index")
ax.legend()
ax.grid()
plt.show()

Mode tracking

As is typical for eigenvalue solvers, Tidy3D’s mode simulation returns the calculated modes in the order corresponding to the magnitudes of the found eigenvalues, that is, the effective indices. Since the effective index of a mode generally depends on frequency and can become larger or smaller in magnitude relative to the effective indices of other modes, the same mode_index value in the returned mode simulation data may correspond to physically different modes at different frequencies.

To avoid such mismatches, `ModeSimulation <https://docs.flexcompute.com/projects/tidy3d/en/latest/api/_autosummary/tidy3d.ModeSimulation.html>`__ performs tracking across frequencies by default based on mode overlap values. See the discussion on mode decomposition. In other words, modes at one frequency are matched to the most similar modes at another frequency. This behavior is controlled by `ModeSpec.sort_spec.track_freq <https://docs.flexcompute.com/projects/tidy3d/en/latest/api/_autosummary/tidy3d.ModeSortSpec.html>`__, typically configured as sort_spec=td.ModeSortSpec(track_freq=...). It can be set to either None, "lowest", "central" (default), or "highest". In the first case, no tracking across frequencies is performed. In the other three cases, tracking is performed starting from the specified frequency.

Additionally, any unsorted modes stored in the ModeSimulationData.modes dataset returned by ModeSimulation.run_local() can be sorted afterwards using the overlap_sort() function. This function takes the argument track_freq (default: "central"), which uses the same values as ModeSortSpec.track_freq, and overlap_thresh (default: 0.9), which 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 most cases, we demonstrate it here to provide a clearer understanding of `ModeSimulation <https://docs.flexcompute.com/projects/tidy3d/en/latest/api/_autosummary/tidy3d.ModeSimulation.html>`__ results. Note that for mode sources and mode monitors, sorting is always performed internally.

We start by calculating twelve modes at eleven different frequencies for a layered waveguide with the default tracking across frequencies 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),
)

# set mode spec to find 12 modes and turn off automatic mode tracking
mode_spec = td.ModeSpec(num_modes=12, sort_spec=td.ModeSortSpec(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 simulation object and solve for modes
mode_simulation = td.ModeSimulation(
    size=(Lx, Ly, Lz),
    grid_spec=grid_spec,
    structures=[waveguide, bottom, top],
    boundary_spec=td.BoundarySpec.all_sides(boundary=td.Periodic()),
    plane=plane,
    mode_spec=mode_spec,
    freqs=freqs,
)

# visualize
mode_simulation.plot(y=0)
plt.show()

mode_sim_data = mode_simulation.run_local()
mode_data = mode_sim_data.modes

By inspecting the values of the 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()

This situation can be remedied by using the overlap_sort() function as follows.

mode_data_sorted = mode_data.overlap_sort(track_freq="central")

Plotting the same information as previously, one can confirm that the modes are now matched as closely as possible.

n_eff_sorted = mode_data_sorted.n_eff  # real part of the effective mode index
n_eff_sorted.plot.line(".-", x="f")
plt.show()

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_sorted.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, mode sorting, whether the default behavior or performed manually, might not successfully resolve the tracking of all modes throughout the chosen frequency range. Possible reasons include: - A mode physically disappears 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. - The presence of degenerate modes at certain frequencies.

Group Index

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

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

mode_sim_data = mode_simulation.run_local()
mode_data = mode_sim_data.modes
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 simulation runs locally with run_local(), which means it does not require credits.

• The same `ModeSimulation <https://docs.flexcompute.com/projects/tidy3d/en/latest/api/_autosummary/tidy3d.ModeSimulation.html>`__ object can also be run on the server with web.run(...).

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

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