Avoiding unphysical loss in low-index-contrast waveguides

When working with low-index-contrast waveguides, such as SiN embedded in SiO₂, the optical field can become highly delocalized, especially in narrow structures such as the tip of an edge coupler, resulting in a large mode field size. In such cases, it’s important to carefully define the mode plane to prevent extended evanescent fields from reaching the Perfectly Matched Layer (PML), as this can lead to unphysical losses and inaccurate simulation results.

In this notebook, we will demonstrate how an insufficiently large plane size can result in unphysical losses.

We will also demonstrate how to assemble a `ModeSolver <https://docs.flexcompute.com/projects/tidy3d/en/latest/api/_autosummary/tidy3d.plugins.mode.ModeSolver.html#tidy3d.plugins.mode.ModeSolver>`__ simulation to quickly determine whether the simulation plane is large enough. Furthermore, we will illustrate how to create a grid to efficiently run a simulation with a sufficiently large plane size.

import matplotlib.pyplot as plt
import numpy as np
import tidy3d as td
import tidy3d.web as web
from tidy3d.plugins.mode import ModeSolver

Defining the Geometry

# global parameters
wl = 1.5
num_freqs = 1

freq0 = td.C_0 / wl
fwidth = 0.2 * freq0
freqs = [freq0] if num_freqs == 1 else np.linspace(freq0 - fwidth, freq0 + fwidth, num_freqs)

source_time = td.GaussianPulse(freq0=freq0, fwidth=fwidth)

symmetry = [0, 1, -1]

sim_center = (0, 0, 0)

run_time = 20e-12

# function for returning the waveguide as a function of length and plane size
def get_sim(length=10, plane_size=6, min_steps_per_wvl=15):
    sim_size = (length, plane_size, plane_size)

    structure_index = 1.99
    background_index = 1.44

    width = 0.33
    height = 0.2

    structure_medium = td.Medium(permittivity=structure_index**2)
    background_medium = td.Medium(permittivity=background_index**2)

    structure = td.Structure(
        geometry=td.Box(center=(0, 0, 0), size=(td.inf, width, height)), medium=structure_medium
    )

    # source
    mode_spec = td.ModeSpec(num_modes=1, precision="double")

    source = td.ModeSource(
        center=(-sim_size[0] / 2 + 0.5, 0, 0),
        size=(0, plane_size, plane_size),
        mode_spec=mode_spec,
        direction="+",
        mode_index=0,
        source_time=source_time,
    )
    # flux monitor
    flux_mon = td.FluxMonitor(
        name="flux_mon",
        size=(0, plane_size, plane_size),
        center=(sim_size[0] / 2 - 0.5, 0, 0),
        freqs=freqs,
    )
    # field monitor
    field_mon = td.FieldMonitor(
        name="field_mon",
        size=(td.inf, td.inf, 0),
        center=(0, 0, 0),
        freqs=freqs,
    )

    # grid spec
    grid_spec = td.GridSpec.auto(min_steps_per_wvl=min_steps_per_wvl)

    sim = td.Simulation(
        size=sim_size,
        center=sim_center,
        grid_spec=grid_spec,
        structures=[structure],
        sources=[source],
        monitors=[flux_mon, field_mon],
        run_time=run_time,
        symmetry=symmetry,
        medium=background_medium,
    )

    ms = ModeSolver(
        mode_spec=source.mode_spec.updated_copy(num_pml=(12, 12)),
        freqs=freqs,
        simulation=sim,
        plane=td.Box(size=(0, plane_size, plane_size), center=source.center),
    )

    return sim, ms

# plotting the simulation domain

sim, _ = get_sim()

fig, ax = plt.subplots(1, 2, figsize=(8, 8))
sim.plot(ax=ax[0], z=0)
sim.plot(ax=ax[1], x=0)

plt.show()

Mode Solver Analysis

First, we will carry out ModeSolver simulations varying the plane size.

ms_list = {}
sizes = [10, 15, 20, 25, 30, 35, 40, 45, 50, 60, 70, 80, 90, 100]

for size in sizes:
    _, ms = get_sim(5, plane_size=size)
    ms_list[str(size)] = ms

ms_batch = web.Batch(simulations=ms_list)
ms_data = ms_batch.run()

08:03:51 UTC Started working on Batch containing 14 tasks.

08:04:02 UTC Maximum FlexCredit cost: 0.134 for the whole batch.

             Use 'Batch.real_cost()' to get the billed FlexCredit cost after
             completion.

08:12:23 UTC Batch complete.

Now, we can plot the imaginary part of the effective index as a function of the plane size. As we can see, for planes smaller than 60 µm, there are losses induced by the interaction of the evanescent fields with the simulation boundaries. These losses will lead to unphysical results for a long waveguide.

loss = [ms_data[str(i)].to_dataframe()["k eff"] for i in sizes]

fig, ax = plt.subplots()
ax.set_xlabel(r"Plane Size ($\mu$m)")
ax.set_ylabel(r"$k_{eff}$ (1/$\mu$m)")

ax.plot(sizes, loss, "-o")

plt.show()

To illustrate the importance of a correct simulation domain size, we compare the mode profiles for plane sizes of 20 µm and 80 µm.

As we can see, in the 20 µm case, the fields extending into the PML are not negligible. In contrast, with an 80 µm plane size, the fields decay sufficiently before reaching the PML boundaries.

ms_list["20"].plot_field("E", "abs")
ms_list["80"].plot_field("E", "abs")


plt.show()

FDTD Simulation

Now, we will simulate a relatively long waveguide with two different plane sizes: 10 µm and 60 µm. We will add FluxMonitors every 20 µm and observe the losses as a function of propagation.

sim_10, _ = get_sim(length=150, plane_size=10)
sim_60, _ = get_sim(length=150, plane_size=60)

monitors = []
positions = range(20, 150, 20)
for i in positions:
    flux_mon = td.FluxMonitor(
        name=f"flux_mon{i}",
        size=(0, 50, 50),
        center=(-sim_10.size[0] / 2 + i, 0, 0),
        freqs=freqs,
    )

    monitors.append(flux_mon)

sim_10 = sim_10.updated_copy(monitors=monitors + [sim_10.monitors[-1]])
sim_60 = sim_60.updated_copy(monitors=sim_10.monitors)


batch = web.Batch(simulations={"10": sim_10, "60": sim_60})

batch_data = batch.run()

08:12:46 UTC WARNING: Mode source at sources[0] has a large number (7.59e+05) of
             grid points. This can lead to solver slow-down and increased cost. 
             Consider making the size of the component smaller, as long as the  
             modes of interest decay by the plane boundaries.                   

08:12:48 UTC Started working on Batch containing 2 tasks.

08:12:50 UTC Maximum FlexCredit cost: 29.436 for the whole batch.

             Use 'Batch.real_cost()' to get the billed FlexCredit cost after
             completion.

08:19:11 UTC Batch complete.

As we can see, for the 10 µm plane size simulation, there are significant losses, although increasing the simulation size can significantly raise the simulation cost.

fig, ax = plt.subplots()

loss_10 = [batch_data["10"][f"flux_mon{i}"].flux for i in positions]
loss_60 = [batch_data["60"][f"flux_mon{i}"].flux for i in positions]

ax.plot(positions, loss_10, "o-", label="plane size 10 µm")
ax.plot(positions, loss_60, "o-", label="plane size 60 µm")

ax.legend()
plt.show()

cost_60 = web.real_cost(batch_data.task_ids["60"])
cost_10 = web.real_cost(batch_data.task_ids["10"])
print(f"Cost for plane size 10: {cost_10}")
print(f"Cost for plane size 60: {cost_60}")

08:19:19 UTC Billed flex credit cost: 2.789.

             Note: the task cost pro-rated due to early shutoff was below the
             minimum threshold, due to fast shutoff. Decreasing the simulation
             'run_time' should decrease the estimated, and correspondingly the
             billed cost of such tasks.

             Billed flex credit cost: 0.155.

             Note: the task cost pro-rated due to early shutoff was below the
             minimum threshold, due to fast shutoff. Decreasing the simulation
             'run_time' should decrease the estimated, and correspondingly the
             billed cost of such tasks.

Cost for plane size 10: 0.1547898657015464
Cost for plane size 60: 2.7888519800076015

Custom Grid for Reducing Simulation Cost

To reduce simulation cost, we can use the minimum permitted value for the autogrid min_steps_per_wavelength argument, which is 6, and add a MeshOverrideRegion around the waveguide. This allows enough space for the fields to decay before reaching the simulation boundaries, while maintaining good resolution around the waveguide.

mesh_structure = sim_60.structures[0]

Si02_steps = (wl / 1.45) / 15
SiN_steps = (wl / 1.99) / 15

mesh_override_SiN = td.MeshOverrideStructure(
    geometry=mesh_structure.geometry,
    dl=(SiN_steps,) * 3,
)


mesh_override_SiO2 = td.MeshOverrideStructure(
    geometry=td.Box(center=(0, 0, 0), size=(td.inf, 3, 3)),
    dl=(Si02_steps,) * 3,
)


grid_spec = td.GridSpec.auto(
    min_steps_per_wvl=6, override_structures=[mesh_override_SiO2, mesh_override_SiN]
)

sim_60_newMesh = sim_60.updated_copy(grid_spec=grid_spec)

ax = sim_60_newMesh.plot(x=0)
sim_60_newMesh.plot_grid(ax=ax, x=0)

plt.show()

task_id = web.upload(sim_60_newMesh, "sim_60_newMesh")
web.start(task_id)
web.monitor(task_id)
sim_data = web.load(task_id)

             WARNING: Mode source at sources[0] has a large number (1.47e+05) of
             grid points. This can lead to solver slow-down and increased cost. 
             Consider making the size of the component smaller, as long as the  
             modes of interest decay by the plane boundaries.                   

             Created task 'sim_60_newMesh' with resource_id
             'fdve-92f9260f-391b-47db-9830-27f6b46e66a8' and task_type 'FDTD'.

             View task using web UI at
             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-92f9260f-391
             b-47db-9830-27f6b46e66a8'.

             Task folder: 'default'.

08:19:21 UTC Estimated FlexCredit cost: 7.806. 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.

08:19:23 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.

08:19:32 UTC status = preprocess

08:19:39 UTC starting up solver

08:19:40 UTC running solver

08:21:01 UTC early shutoff detected at 4%, exiting.

             status = postprocess

08:21:04 UTC status = success

08:21:06 UTC View simulation result at
             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-92f9260f-391
             b-47db-9830-27f6b46e66a8'.

08:21:07 UTC Loading results from simulation_data.hdf5

# Here we will use the web API get_tasks to get the cost of the last simulation

cost_60_newMesh = web.real_cost(task_id)
print(f"Cost for plane size 60 with new mesh: {cost_60_newMesh}")

08:21:08 UTC Billed flex credit cost: 0.781.

             Note: the task cost pro-rated due to early shutoff was below the
             minimum threshold, due to fast shutoff. Decreasing the simulation
             'run_time' should decrease the estimated, and correspondingly the
             billed cost of such tasks.

Cost for plane size 60 with new mesh: 0.7805700463461971

As we can see, this meshing strategy returns very similar results, while consuming less than one-third of the Flexcredits.

loss_60_newMesh = [sim_data[f"flux_mon{i}"].flux for i in positions]

fig, ax = plt.subplots()
ax.plot(positions, loss_10, label="plane size 10 µm")
ax.plot(positions, loss_60, label="plane size 60 µm")
ax.plot(positions, loss_60_newMesh, "o", label="plane size 60 µm - new mesh")
ax.legend()

plt.show()