In this example we use Tidy3D's heat solver to investigate tuning of a straight rectangular waveguide using a heating strip.

Schematic of the thermally tuned waveguide

[1]:

import numpy as np
from matplotlib import pyplot as plt

import tidy3d as td
import tidy3d.web as web

Note that package vtk is required to use full visualization capabilities of heat solver data. To install a compatible version of vtk one can use pip install tidy3d[vtk] command.

Setup Simulation Scene¶

The simulation setup consists of a silicon waveguide in a silicon dioxide cladding on top of a silicon wafer. The heating strip is placed in the cladding above the waveguide. Let us define basic geometric parameters of simulated structures.

[2]:

# in um

w_sim = 16  # simulation width
h_buffer = 1  # vertical simulation buffer

h_clad = 2.8  # thickness of cladding
h_box = 2  # thickness of buried oxide
h_wafer = 0.5  # wafer thickness

# waveguide parameters
w_core = 0.5
h_core = 0.22

# heater parameters
h_heater = 0.14
w_heater = 2
d_heater = 2  # distance between heater and waveguide

We will place the wafer's top surface at location $z=0$. Let us define vertical locations of each structure. Note that we will model buried oxide and cladding as a single structure since they are made of the same material.

[3]:

z_wafer = -h_wafer / 2
z_core = h_box + h_core / 2
z_heater = h_box + h_core + d_heater + h_heater / 2

box_clad_max = h_box + h_clad
box_clad_min = (
    -h_wafer / 2
)  # make structures overlap to ensure no gaps due to numerical roundoff errors
h_box_clad = box_clad_max - box_clad_min
z_box_clad = (box_clad_max + box_clad_min) / 2

z_sim = (h_box + h_clad - h_wafer) / 2
h_sim = h_box + h_clad + h_wafer + 2 * h_buffer

Define optic and thermal properties of materials.

[4]:

Si_n = 3.4777  # Si refraction index
Si_n_slope = 1.86e-4  # Si thermo-optic coefficient dn/dT, 1/K
Si_k = 148e-6  # Si thermal conductivity, W / (um * K)
Si_s = 0.71 * 2.33e-12  # Si volumetric heat capacity, J / (um^3 * K)

SiO2_n = 1.444  # SiO2 refraction index
SiO2_n_slope = 1e-5  # SiO2 thermo-optic coefficient dn/dT, 1/K
SiO2_k = 1.38e-6  # SiO2 thermal conductivity, W/(um*K)
SiO2_s = 0.709 * 2.203e-12  # SiO2 volumetric heat capacity, J / (um^3 * K)

TiN_k = 28e-6  # TiN thermal conductivity W/(um*K)
TiN_s = 0.598 * 5240e-12  # TiN volumetric heat capacity, J / (um^3 * K)
TiN_sigma = 2.3  # Electric conductivity of TiN, S/um

Convert refraction indices into permittivity values.

[5]:

Si_eps = Si_n**2
Si_eps_slope = 2 * Si_n * Si_n_slope

SiO2_eps = SiO2_n**2
SiO2_eps_slope = 2 * SiO2_n * SiO2_n_slope

Create material classes containing relevant information for each medium. Note that we use class PerturbationMedium to specify sensitivity of optic properties to temperature and field heat_spec to specify thermal properties. Optic characteristics of the heating strip material is modeled as a perfect electric conductor.

[6]:

Si = td.PerturbationMedium(
    permittivity=Si_eps,
    permittivity_perturbation=td.ParameterPerturbation(
        heat=td.LinearHeatPerturbation(coeff=Si_eps_slope, temperature_ref=300)
    ),
    heat_spec=td.SolidSpec(
        conductivity=Si_k,
        capacity=Si_s,
    ),
    name="Si",
)

SiO2 = td.PerturbationMedium(
    permittivity=SiO2_eps,
    permittivity_perturbation=td.ParameterPerturbation(
        heat=td.LinearHeatPerturbation(coeff=SiO2_eps_slope, temperature_ref=300)
    ),
    heat_spec=td.SolidSpec(
        conductivity=SiO2_k,
        capacity=SiO2_s,
    ),
    name="SiO2",
)

TiN = td.PECMedium(
    heat_spec=td.SolidSpec(
        conductivity=TiN_k,
        capacity=TiN_s,
    ),
    name="TiN",
)

air = td.Medium(heat_spec=td.FluidSpec(), name="air")

Create structures representing waveguide, heater, cladding, and wafer.

[7]:

core = td.Structure(
    geometry=td.Box(center=(0, 0, z_core), size=(w_core, td.inf, h_core)),
    medium=Si,
)

heater = td.Structure(
    geometry=td.Box(center=(0, 0, z_heater), size=(w_heater, td.inf, h_heater)),
    medium=TiN,
    name="heater",
)

box_clad = td.Structure(
    geometry=td.Box(center=(0, 0, z_box_clad), size=(td.inf, td.inf, h_box_clad)),
    medium=SiO2,
)

wafer = td.Structure(
    geometry=td.Box(center=(0, 0, z_wafer), size=(td.inf, td.inf, h_wafer)),
    medium=Si,
    name="wafer",
)

Pack structures into a simulation scene.

[8]:

scene = td.Scene(
    medium=air,
    structures=[box_clad, core, heater, wafer],
)

Visual validation of the setup.

[9]:

scene.plot(
    y=0, hlim=[-w_sim / 2, w_sim / 2], vlim=[z_sim - h_sim / 2, z_sim + h_sim / 2]
)
plt.show()

No description has been provided for this image

Setup Heat Simulation¶

We assume that the bottom of the wafer is maintained at temperature $T=300$ K and ignore heat transfer at the cladding/air boundary. In this simulation setup these boundary conditions can be set using MediumMediumInterface placement. Note since no other boundary conditions are specified, the other sides of the simulation box will be treated as insulating boundary conditions.

Convective cooling/heating at the cladding/air boundary can be easily taken into account using ConvectionBC boundary conditions.

[10]:

bc_bottom = td.HeatBoundarySpec(
    placement=td.MediumMediumInterface(mediums=["Si", "air"]),
    condition=td.TemperatureBC(temperature=300),
)

bc_top = td.HeatBoundarySpec(
    placement=td.MediumMediumInterface(mediums=["SiO2", "air"]),
    condition=td.HeatFluxBC(flux=0),
)

Create a uniform heat source inside the heater structure equivalent to electrical current of $7.3 \times 10^{-3}$ A.

[11]:

current = 7.4e-3  # A
heat_rate = (current / h_heater / w_heater) ** 2 / TiN_sigma  # convert into power
heater_source = td.UniformHeatSource(rate=heat_rate, structures=[heater.name])

We will measure temperature distribution in the entire cross section of the simulation setup. Note that in this example we request to return heat solver results on the original unstructured grid. Moreover, we request the monitor plane to be meshed conformally.

[12]:

temp_mnt = td.TemperatureMonitor(size=(td.inf, 0, td.inf), name="temperature", unstructured=True, conformal=True)

For spatial discretization we will use a distance-based meshing specifying that solid-solid material interfaces should be discretized with mesh size equal to the heater height h_heater divided by 3 and away from interfaces the grid should coarsen to four times of that size. Setting distance_interface to three of minimal mesh size will create approximately three layers of cells with the minimal size around interfaces. Cells outside of this region will have sizes according to linear relation dl_interface * (1 - ratio(d)) + dl_interface * ratio(d), where ratio(d) = (d - distance_interface) / (distance_bulk - distance_interface) and d is the distance to the closest material interface. We choose this transition to happen on the order of 2 * d_heater. We will also exclude using the higher mesh resolution around the wafer since that region is of less interest.

[13]:

# parameters for meshing
dl_min = h_heater / 3
dl_max = 4 * dl_min

grid_spec = td.DistanceUnstructuredGrid(
    dl_interface=dl_min,
    dl_bulk=dl_max,
    distance_interface=3 * dl_min,
    distance_bulk=2 * d_heater,
    non_refined_structures=["wafer"],  # do not refine near wafer
)

Pack components into a heat simulation specification. Note that from the heat transfer perspective the problem is symmetric around $x=0$. Thus, we can set symmetry specification to (1, 0, 0) to reduce simulation size and cost.

Given that the problem is two-dimensional we will use a small simulation size along the invariant $y$-direction.

[14]:

heat_sim = td.HeatSimulation.from_scene(
    scene=scene,
    center=(0, 0, z_sim),
    size=(w_sim, dl_min, h_sim),
    boundary_spec=[bc_bottom, bc_top],
    sources=[heater_source],
    monitors=[temp_mnt],
    symmetry=(1, 0, 0),
    grid_spec=grid_spec,
)

[15]:

heat_sim.plot(y=0.01)
plt.show()

Mesh Inspection¶

Before proceeding to submitting a parameter sweep batch we run one simulation and inspect its mesh and solution.

[16]:

job = web.Job(simulation=heat_sim, task_name="heat_sim_check_mesh")
heat_sim_data = job.run()

19:43:55 CST Created task 'heat_sim_check_mesh' with task_id                    
             'he-af58dbdc-e428-402f-86c5-3aaf0ec8b2c7' and task_type 'HEAT'.

             Tidy3D's heat solver is currently in the beta stage. All heat      
             simulations are charged a flat fee of 0.025 FlexCredit.

Output()

19:43:56 CST Heat solver status: success

Output()

19:43:58 CST loading simulation from simulation_data.hdf5

Since while setting up the temperature monitor we requested to return data on unstructured grid, we can plot the associated monitor data to visualize both the solver mesh and the solution.

[17]:

_, ax = plt.subplots(3, 1, figsize=(10, 10))

heat_sim.plot_structures_heat_conductivity(y=0, alpha=0.5, ax=ax[0])
heat_sim_data["temperature"].temperature.plot(ax=ax[0], field=False)
ax[0].set_xlim([-8, 8])

heat_sim_data["temperature"].temperature.plot(ax=ax[1], grid=False)
ax[1].set_xlim([-8, 8])

heat_sim_data["temperature"].temperature.plot(ax=ax[2])
ax[2].set_xlim([-8, 8])

plt.show()

Batch Processing of Heat Simulations¶

Having ensured that the simulation mesh is suitable for this simulation setup, let us perform a sweep through several values of electric current in the heating strip. To do this we create a dictionary of simulations with different source definitions.

[18]:

currents = np.linspace(1e-10, 7.4e-3, 10)
heat_sims = {}

for ind, current in enumerate(currents):
    heat_rate = (current / h_heater / w_heater) ** 2 / TiN_sigma
    heater_source = td.UniformHeatSource(rate=heat_rate, structures=[heater.name])

    heat_sim_new = heat_sim.updated_copy(sources=[heater_source])
    heat_sims[f"heat_wg_current_{ind}"] = heat_sim_new

And submit this dictionary using Batch functionality.

[19]:

batch = web.Batch(simulations=heat_sims)

19:44:01 CST Created task 'heat_wg_current_0' with task_id                      
             'he-01827531-5afe-4c82-9d0c-0ee0230d7665' and task_type 'HEAT'.

             Tidy3D's heat solver is currently in the beta stage. All heat      
             simulations are charged a flat fee of 0.025 FlexCredit.

Output()

19:44:02 CST Created task 'heat_wg_current_1' with task_id                      
             'he-e49e54fc-d4d1-4478-85c0-bc1d076e2e1b' and task_type 'HEAT'.

             Tidy3D's heat solver is currently in the beta stage. All heat      
             simulations are charged a flat fee of 0.025 FlexCredit.

Output()

19:44:03 CST Created task 'heat_wg_current_2' with task_id                      
             'he-9a553c15-eea2-44e8-afb3-1db905cca9b3' and task_type 'HEAT'.

             Tidy3D's heat solver is currently in the beta stage. All heat      
             simulations are charged a flat fee of 0.025 FlexCredit.

Output()

19:44:04 CST Created task 'heat_wg_current_3' with task_id                      
             'he-4476021e-7b87-437e-ba01-a450f2d11c2e' and task_type 'HEAT'.

             Tidy3D's heat solver is currently in the beta stage. All heat      
             simulations are charged a flat fee of 0.025 FlexCredit.

Output()

19:44:05 CST Created task 'heat_wg_current_4' with task_id                      
             'he-928dd351-5c03-4e40-bf99-14a51f5981e5' and task_type 'HEAT'.

             Tidy3D's heat solver is currently in the beta stage. All heat      
             simulations are charged a flat fee of 0.025 FlexCredit.

Output()

             Created task 'heat_wg_current_5' with task_id                      
             'he-3134bb5a-d6bc-4429-891b-f26fb8bcbd53' and task_type 'HEAT'.

             Tidy3D's heat solver is currently in the beta stage. All heat      
             simulations are charged a flat fee of 0.025 FlexCredit.

Output()

19:44:06 CST Created task 'heat_wg_current_6' with task_id                      
             'he-805d47ab-1b85-4968-9ee6-6e6a979ab8d0' and task_type 'HEAT'.

             Tidy3D's heat solver is currently in the beta stage. All heat      
             simulations are charged a flat fee of 0.025 FlexCredit.

Output()

19:44:07 CST Created task 'heat_wg_current_7' with task_id                      
             'he-82254872-d4d0-47a6-88a6-be4e28351bd5' and task_type 'HEAT'.

             Tidy3D's heat solver is currently in the beta stage. All heat      
             simulations are charged a flat fee of 0.025 FlexCredit.

Output()

19:44:08 CST Created task 'heat_wg_current_8' with task_id                      
             'he-43c91b10-0426-4cff-9fb0-0bc11e308eff' and task_type 'HEAT'.

             Tidy3D's heat solver is currently in the beta stage. All heat      
             simulations are charged a flat fee of 0.025 FlexCredit.

Output()

             Created task 'heat_wg_current_9' with task_id                      
             'he-fe15222d-95a3-4186-abaf-565ec732c3c5' and task_type 'HEAT'.

             Tidy3D's heat solver is currently in the beta stage. All heat      
             simulations are charged a flat fee of 0.025 FlexCredit.

Output()

[20]:

batch_data = batch.run()

19:44:12 CST Started working on Batch.

19:44:16 CST Maximum FlexCredit cost: 0.250 for the whole batch.

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

Output()

19:44:19 CST Batch complete.

Visualize some of the simulation results.

[21]:

# max temperature value when plotting to enforce the same coloring
temp_max_plot = 340

fig, ax = plt.subplots(1, 3, figsize=(17, 2))

batch_data["heat_wg_current_1"].plot_field("temperature", ax=ax[0], vmax=temp_max_plot)
ax[0].set_title(f"Current I = {currents[1]:1.3} A")

batch_data["heat_wg_current_5"].plot_field("temperature", ax=ax[1], vmax=temp_max_plot)
ax[1].set_title(f"Current I = {currents[5]:1.3} A")

batch_data["heat_wg_current_9"].plot_field("temperature", ax=ax[2], vmax=temp_max_plot)
ax[2].set_title(f"Current I = {currents[9]:1.3} A")

plt.show()

Output()

19:44:22 CST loading simulation from                                            
             ./he-e49e54fc-d4d1-4478-85c0-bc1d076e2e1b.hdf5

Output()

19:44:23 CST loading simulation from                                            
             ./he-3134bb5a-d6bc-4429-891b-f26fb8bcbd53.hdf5

Output()

19:44:25 CST loading simulation from                                            
             ./he-fe15222d-95a3-4186-abaf-565ec732c3c5.hdf5

Coupling to Optic Simulations¶

Now let us create optic simulations where material properties are modified according to obtained temperature distributions. One can conveniently do that using method .perturbed_mediums_copy() of classes Scene and Simulation. We will use obtained optic simulations only to investigate waveguide modes, thus we are defining neither sources nor monitors, and only provide grid specification.

Recall that in this tutorial we chose to retrieve heat solver results on an unstructured grid. However, currently method .perturbed_mediums_copy() supports only Cartesian arrays SpatialDataArray. Thus, we need to interpolate the unstructured data array onto a suitable Cartesian grid. To this end, we first create an unperturbed optic simulation.

[22]:

wvl_um = 1.55
freq0 = td.C_0 / wvl_um

grid_spec = td.GridSpec.auto(min_steps_per_wvl=50, wavelength=wvl_um)

optic_sim = td.Simulation.from_scene(
    scene=scene,
    center=(0, 0, z_sim),
    size=(w_sim, dl_min, h_sim),
    run_time=1e-15,
    grid_spec=grid_spec,
)

Now we can use its grid as the target grid to interpolate the temperature fields to.

[23]:

target_grid = optic_sim.grid.boundaries
perturb_sims = []

for _, hs_data in batch_data.items():
    temp_interpolated = hs_data["temperature"].temperature.interp(x=target_grid.x, y=0, z=target_grid.z, fill_value=300)
    psim = optic_sim.perturbed_mediums_copy(temperature=temp_interpolated)
    perturb_sims.append(psim)

Output()

19:44:27 CST loading simulation from                                            
             ./he-01827531-5afe-4c82-9d0c-0ee0230d7665.hdf5

19:44:30 CST loading simulation from                                            
             ./he-e49e54fc-d4d1-4478-85c0-bc1d076e2e1b.hdf5

Output()

19:44:33 CST loading simulation from                                            
             ./he-9a553c15-eea2-44e8-afb3-1db905cca9b3.hdf5

Output()

19:44:37 CST loading simulation from                                            
             ./he-4476021e-7b87-437e-ba01-a450f2d11c2e.hdf5

Output()

19:44:41 CST loading simulation from                                            
             ./he-928dd351-5c03-4e40-bf99-14a51f5981e5.hdf5

19:44:44 CST loading simulation from                                            
             ./he-3134bb5a-d6bc-4429-891b-f26fb8bcbd53.hdf5

Output()

19:44:47 CST loading simulation from                                            
             ./he-805d47ab-1b85-4968-9ee6-6e6a979ab8d0.hdf5

Output()

19:44:51 CST loading simulation from                                            
             ./he-82254872-d4d0-47a6-88a6-be4e28351bd5.hdf5

Output()

19:44:54 CST loading simulation from                                            
             ./he-43c91b10-0426-4cff-9fb0-0bc11e308eff.hdf5

19:44:57 CST loading simulation from                                            
             ./he-fe15222d-95a3-4186-abaf-565ec732c3c5.hdf5

Define plane where waveguide modes will be calculated.

[24]:

mode_plane = td.Box(center=(0, 0, z_core), size=(4, 0, 3))

Visual inspection.

[25]:

fig, ax = plt.subplots(1, 1)
perturb_sims[0].plot(y=0, ax=ax)
mode_plane.plot(y=0, ax=ax, alpha=0.5)
plt.show()

We will use Tidy3D's mode solver plugin.

[26]:

from tidy3d.plugins.mode.web import run as run_mode
from tidy3d.plugins.mode import ModeSolver

Find only the first mode at the frequency corresponding to 1.55 um wavelength. Note that simulation objects that are passed to the mode solver contain custom medium data in this case. Sometimes this can lead to large amounts of uploaded/downloaded data. To avoid that, the mode solver web API by default restricts passed simulation object containing custom mediums to mode solver plane only. When that happens, a warning is displayed. Setting reduce_simulation=True forces such simulation reductions for all cases (even if no custom mediums present) and effectively silences the warning.

[27]:

mode_datas = []

for psim in perturb_sims:
    ms = ModeSolver(
        simulation=psim,
        plane=mode_plane,
        mode_spec=td.ModeSpec(num_modes=1, precision="double"),
        freqs=[freq0],
    )

    md = run_mode(ms, reduce_simulation=True)

    mode_datas.append(md)

19:45:00 CST Mode solver created with                                           
             task_id='fdve-ddfc2c93-46e8-41a3-978b-08f423d158a1',               
             solver_id='mo-c26288f7-770b-43e8-a49e-5cbba93742a3'.

Output()

Output()

19:45:03 CST Mode solver status: success

Output()

19:45:05 CST Mode solver created with                                           
             task_id='fdve-81477046-9a13-45fd-80d7-c916e02abd22',               
             solver_id='mo-fb644cef-bda4-46a8-a34f-5c8e20382e1b'.

Output()

Output()

19:45:07 CST Mode solver status: success

Output()

19:45:09 CST Mode solver created with                                           
             task_id='fdve-052aafa6-345d-42bc-b58f-44d1613adf4b',               
             solver_id='mo-2a7065dc-d79f-441b-9c11-ab09a3a9184a'.

Output()

Output()

19:45:11 CST Mode solver status: success

Output()

19:45:13 CST Mode solver created with                                           
             task_id='fdve-1f4f261b-5ee2-4fd6-87da-1e848f87a5fe',               
             solver_id='mo-9d1affff-0d23-4c49-957e-7c3d7c5c2858'.

Output()

Output()

19:45:16 CST Mode solver status: success

Output()

19:45:20 CST Mode solver created with                                           
             task_id='fdve-974f22bb-9832-41eb-8adc-9ed8435bc714',               
             solver_id='mo-ce916752-8e40-4569-95b8-1f158f70b5d8'.

Output()

Output()

19:45:23 CST Mode solver status: success

Output()

19:45:25 CST Mode solver created with                                           
             task_id='fdve-d5ca94c7-99c2-4e06-88cd-fa84b047eba6',               
             solver_id='mo-f9fe962f-08c6-4ddd-8371-f1c80f3a0f7f'.

Output()

Output()

19:45:28 CST Mode solver status: success

Output()

19:45:30 CST Mode solver created with                                           
             task_id='fdve-d82aae5c-e449-4409-91c1-85c572922b1c',               
             solver_id='mo-87b61823-907c-4ead-9d5e-d7c41dfc3cce'.

Output()

Output()

19:45:33 CST Mode solver status: success

Output()

19:45:35 CST Mode solver created with                                           
             task_id='fdve-63b5090f-8ab3-409b-b13a-c6e3a0a4a0aa',               
             solver_id='mo-a0fe5e89-e067-4c4d-855e-9891630ed317'.

Output()

Output()

19:45:38 CST Mode solver status: success

Output()

19:45:40 CST Mode solver created with                                           
             task_id='fdve-37b3c93b-39c8-45ea-9f12-290497619986',               
             solver_id='mo-f1c02607-d9ef-4af1-b797-88002b71d56e'.

Output()

Output()

19:45:44 CST Mode solver status: success

Output()

19:45:46 CST Mode solver created with                                           
             task_id='fdve-b765c948-1e35-41f3-98de-26a8383d14c2',               
             solver_id='mo-9ec5c4d3-5b36-4b11-becc-69fa406e7bd0'.

Output()

Output()

19:45:49 CST Mode solver status: success

Output()

Extract found propagation indices into a single array.

[28]:

n_eff = np.array([md.n_eff.isel(f=0, mode_index=0) for md in mode_datas])

We will compare to results obtained using different thermal and mode solvers.

[29]:

currents_ref = np.linspace(0.0, 7.4e-3, 10)
neff_ref = [
    2.4401133396065293,
    2.440149458129956,
    2.440257814666438,
    2.4404384121136413,
    2.4406912552977094,
    2.441016350968344,
    2.441413707791896,
    2.4418833363425283,
    2.4424252490914293,
    2.443039460394046,
]

The absolute values of propagation indices are quite sensitive to the type of solver used and parameters of mode solver solution domain. Thus, we compare changes to the propagation indices rather then their absolute values.

[30]:

print(f"Phase shift over 320um: {2 * np.pi / 1.55 * (n_eff[-1] - n_eff[0]) * 320}")
print(
    f"Phase shift over 320um (ref): {2 * np.pi / 1.55 * (neff_ref[-1] - neff_ref[0]) * 320}"
)
plt.plot(currents * 1e3, n_eff[:] - n_eff[0], "o:")
plt.plot(currents_ref * 1e3, np.array(neff_ref) - neff_ref[0])
plt.xlabel("Current / mA")
plt.ylabel("Change in effective refractive index $n_{eff}$")
plt.legend(["tidy3d", "reference"])
plt.show()

Phase shift over 320um: 3.807916187409721
Phase shift over 320um (ref): 3.795687048084378

[ ]:

Learning Center

Thermally tuned waveguide¶

Setup Simulation Scene¶

Setup Heat Simulation¶

Mesh Inspection¶

Batch Processing of Heat Simulations¶

Coupling to Optic Simulations¶

More Examples