Thermally tuned ring resonator#

Note: the cost of running the entire notebook is larger than 10 FlexCredits.

This example demonstrates the usage of the Tidy3D’s heat solver feature to study the influence of heat propagation on optical properties of a simple ring resonator element.

Our simulation setup will consist of two straight rectangular waveguides coupled by a ring waveguide made of Si. These elements are embedded into a SiO2 cladding on a Si wafer. The bottom of of the wafer is maintained at \(T=300\) K and the top of the cladding is exposed to air, also maintained at \(T=300\) K. The temperature field is controlled by a TiN ring heater above the ring waveguide, which is also embedded in the silicon dioxide cladding.

import numpy as np
from matplotlib import pyplot as plt

import tidy3d as td
from tidy3d import web

This study will involve using two different solvers, thermal and optical ones. The simulation setups for these solvers are defined through two different classes HeatSimulation and Simulation, respectively. However, thanks to a general simulation description class Scene, the geometric setup of the problem can be easily transferred from one simulation to another.

Material Specification#

To create simulation setups that can be used in different solvers we need to create material specifications that contain all relevant information for each of the solvers. Specifically, when performing coupled thermal and optic simulations, each material definition will contain up to three different characteristic: 1. optic properties such as permittivity and conductivity 2. thermal properties such as thermal conductivity 3. response of optic properties to changes in temperature

Let us, for example, create a specification for silicon and silicon dioxide materials that contains all three categories of information. First, we point out that optic mediums containing response models to thermal fields (that is, information from categories 1. and 3.) can be created by using classes PerturbationMedium and PerturbationPoleResidue for non-dispersive and dispersive material models, respectively. Second, thermal properties (that is, information from category 2.) can be specified in any material class by using field heat_spec.

Here, we will use a non-dispersive material model with a linear response of permittivity to temperature. Specification SolidSpec is used to specify thermal properties of the material.

Let us first define values of relevant material properties.

Si_n = 3.4777  # Si refraction index
Si_n_slope = 1.86e-4  # Si thermo-optic coefficient dn/dT, 1/K
Si_k = 148e-12  # 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-12  # SiO2 thermal conductivity, W/(um*K)
SiO2_s = 0.709 * 2.203e-12  # SiO2 volumetric heat capacity, J / (um^3 * K)

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

Convert refraction indices and their thermal sensitivities into permittivity values.

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

Pack material properties into Tidy3D medium classes.

Si = td.PerturbationMedium(
        heat=td.LinearHeatPerturbation(coeff=Si_eps_slope, temperature_ref=300)

SiO2 = td.PerturbationMedium(
        heat=td.LinearHeatPerturbation(coeff=SiO2_eps_slope, temperature_ref=300)

TiN = td.PECMedium(

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

Note that we model TiN as a perfect electric conductor and assume its optic properties are independent of temperature. Thus, a simple non-Perturbation medium class is used.


Now we create structures that constitute our simulation setup. We start with defining basic geometric parameters.

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

# waveguide parameters
w_wg = 0.5  # width
h_wg = 0.22  # thickness
couple_gap = 0.1  # between ring and straight waveguides
ring_radius = 5

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

We will place the bottom surfaces of waveguides at \(z=0\) and the center of the ring resonator at \(x=0\), \(y=0\). Provided that, let us compute auxiliary variable specifying vertical positions of different elements.

# z positions of wafer, heater, and waveguide
z_wafer = -h_box - h_wafer / 2
z_heater = h_wg + d_heater + h_heater / 2
z_wg = h_wg / 2

# y displacement of straight waveguides
y_wg = ring_radius + w_wg / 2 + couple_gap + w_wg / 2

# thickness and z position of box and cladding combined
# make structures overlap to ensure no gaps due to numerical roundoff errors
box_clad_zmax = h_clad
box_clad_zmin = -h_box - h_wafer / 2
h_box_clad = box_clad_zmax - box_clad_zmin
z_box_clad = (box_clad_zmax + box_clad_zmin) / 2

First, we create structures representing the wafer and the SiO2 layer.

wafer = td.Structure(
        center=[0, 0, z_wafer],
        size=[td.inf, td.inf, h_wafer],

box_clad = td.Structure(
        center=[0, 0, z_box_clad],
        size=[td.inf, td.inf, h_box_clad],

Straight waveguide are simply created using Box class.

waveguide_top = td.Structure(
        center=[0, y_wg, z_wg],
        size=[td.inf, w_wg, h_wg],

waveguide_bottom = td.Structure(
        center=[0, -y_wg, z_wg],
        size=[td.inf, w_wg, h_wg],

Ring resonator and ring heater are created by, first, creating cylinders made of required materials and, second, creating smaller cylinders made of SiO2 that will define inner radius of these structures. The gap in the heater ring is also created using an auxiliary structure made of SiO2

# outer ring radius
wg_ring_outer = td.Structure(
        center=[0, 0, z_wg],
        radius=ring_radius + w_wg / 2.0,

# inner ring radius
wg_ring_inner = td.Structure(
        center=[0, 0, z_wg],
        radius=ring_radius - w_wg / 2.0,
        length=w_wg * 1.5,

# heater outer radius
heater_outer = td.Structure(
        center=[0, 0, z_heater],
        radius=ring_radius + w_heater / 2.0,

# heater inner radius
heater_inner = td.Structure(
        center=[0, 0, z_heater],
        radius=ring_radius - w_heater / 2.0,
        length=h_heater * 1.5,

heater_gap = td.Structure(
        center=[0, -ring_radius, z_heater],
            size=(0.5 * w_heater, 1.5 * w_heater, 1.5 * h_heater),


Now the structure can be combined together into a Scene along with simulation size information.

scene = td.Scene(

Visualization to confirm problem setup.

fig, ax = plt.subplots(1, 3, figsize=(10, 3))
scene.plot(x=0, ax=ax[0])
scene.plot(z=z_wg, ax=ax[1])
scene.plot(z=z_heater, ax=ax[2])

We can also visualize structures according to their physical optic and thermal parameters.

fig, ax = plt.subplots(1, 3, figsize=(20,5))
scene.plot(x=0, ax=ax[0])
scene.plot_eps(x=0, ax=ax[1])
scene.plot_heat_conductivity(x=0, ax=ax[2])

Heat Simulation#

To define a heat simulation setup (HeatSimulation) based on the created Scene object we need to define heat related specifications.

We start with specifying the simulation box dimensions. To do that we first define values for buffer space around created structures. We also define variable heat_sim_buffer that define how much larger a heat simulation is compared to an optic one. This is done to make sure that when an optic simulation is constructed we have valid temperature values in the PML regions as well. This will be applied only in x and y directions because the heat simulation is restricted to the region from the bottom of the wafer to the top of the cladding. This is because heat equation is not solved in mediums with heat_spec=FluidSpec(), which is true for medium air defined above.

x_buffer = 2  # buffer in x direction (from ring resonator)
y_buffer = 1  # buffer in y direction (from straight waveguides)
z_buffer = 1  # buffer in z direction

# difference between thermal and optic simulations
heat_sim_buffer = 3

Based on that compute simulation box size and its center along z.

z_sim = (h_clad - h_box - h_wafer) / 2
h_sim = h_box + h_clad + h_wafer + 2 * z_buffer
l_heat_sim = 2 * ring_radius + 2 * x_buffer + 2 * heat_sim_buffer
w_heat_sim = 2 * ring_radius + w_heater + 2 * y_buffer + 2 * heat_sim_buffer

Before proceeding we can inspect the chosen simulation box dimensions using Scene plotting functions with custom vertical and horizontal extents.

fig, ax = plt.subplots(1, 3, figsize=(10, 3))
    x=0, ax=ax[0], hlim=[-w_heat_sim / 2, w_heat_sim / 2], vlim=[z_sim - h_sim / 2, z_sim + h_sim / 2]
    z=z_wg, ax=ax[1], hlim=[-l_heat_sim / 2, l_heat_sim / 2], vlim=[-w_heat_sim / 2, w_heat_sim / 2]
    z=z_heater, ax=ax[2], hlim=[-l_heat_sim / 2, l_heat_sim / 2], vlim=[-w_heat_sim / 2, w_heat_sim / 2]

Having inspected that, let us proceed to other heat simulation specifications.

Create boundary conditions (HeatBoundarySpec) corresponding to temperature (TemperatureBC) boundary conditions on the bottom of the wafer, and convective (ConvectionBC) boundary condition at the cladding/air interface. Note that the convection heat transfer coefficient must be provided in units of W / (um^2 K). To specify the placement of these boundary conditions we use MediumMediumInterface specification.

bc_air = td.HeatBoundarySpec(
    condition=td.ConvectionBC(ambient_temperature=300, transfer_coeff=10 * 1e-12),
bc_wafer = td.HeatBoundarySpec(

Create a uniform heat source (UniformHeatSource) inside the ring heater to represent the Joule heating. For a given electric current in the heater \(I\) we can approximately calculate the volumetric Joule heat generation using formula \(\frac{dP}{dV} = \frac{1}{\sigma_{TiN}} \left( \frac{I}{w_{heater} h_{heater}} \right)^2\). For a current of 50 mA we get:

heat_current = 30e-3  # A
volumetric_heat = (heat_current / w_heater / h_heater) ** 2 / TiN_sigma  # W / um^3
print("Volumetric heat rate: ", volumetric_heat, "W / um^3")
heater_source = td.UniformHeatSource(structures=["heater_outer"], rate=volumetric_heat)
Volumetric heat rate:  4.991126885536821e-09 W / um^3

For spatial discretization we will use class DistanceUnstructuredGrid that specifies desired mesh size based on the distance from material interfaces. Given that h_heater is the smallest feature size of structures, we will select the mesh size near interfaces based on that. Mesh size away from interfaces is selected such that the coarse mesh would have approximately 15 points in the z-direction across the BOX + cladding region. Selecting distance_interface as 2 * dl_interface we are enforcing approximately 2 layers of finest cells around interfaces. And distance_bulk=h_box makes the transition from the finest grid size to the coarsest grid size happen on the order of h_box. We also exclude structure "wafer" from enforcing dl_interafce since we are primarily interested in the temperature distribution in the waveguide region.

dl_interface = h_heater / 4
dl_bulk = (h_box + h_clad) / 15

    distance_interface=2 * dl_interface,

We will use the results of heat simulation in a subsequent optical simulation. Thus we will record a three-dimensional temperature distribution throughout the simulation domain with a TemperatureMonitor.

temp_mnt = td.TemperatureMonitor(size=(td.inf, td.inf, td.inf), name="all")

By default for convenience purposes, heat monitors return data on equivalent Cartesian grids. However, if the unprocessed data and/or the original unstructured grid are of interest, one can request that by setting unstructured=True. Additionally, the option conformal=True allows to request the monitor to be meshed conformally on the unstructured grid. This is useful for visualization purposes, because plots of slices of tetrahedral grids maybe be confusing for visual inspection. Here, we will create monitors with and without conformal meshing for demonstration purposes. In total, we create 3 planar monitors: horizontal through the waveguie structure, horizontal through the heater, and a vertical one.

temp_mnt_ugrid_wg = td.TemperatureMonitor(size=(td.inf, td.inf, 0), center=(0, 0, z_wg), name="ugrid_wg", unstructured=True, conformal=True)
temp_mnt_ugrid_heater = td.TemperatureMonitor(size=(td.inf, td.inf, 0), center=(0, 0, z_heater), name="ugrid_heater", unstructured=True)
temp_mnt_ugrid_vert = td.TemperatureMonitor(size=(0, td.inf, td.inf), center=(1e-3, 0, 0), name="ugrid_vert", unstructured=True)

Combining everything into a HeatSimulation instance. Instead of specifying structures and background medium by hand, we use a convenience class method from_scene(). Note that from the point of view of heat distribution this problem is symmetric with respect to \(x = 0\). Thus, we can use field symmetry to reduce the size, time, and cost for the heat simulation.

# create simulation
sim = td.HeatSimulation.from_scene(
    size=[l_heat_sim, w_heat_sim, h_sim],
    center=(0, 0, z_sim),
    boundary_spec=[bc_air, bc_wafer],
    symmetry=(1, 0, 0),
    monitors=[temp_mnt, temp_mnt_ugrid_wg, temp_mnt_ugrid_heater, temp_mnt_ugrid_vert],

Visualization of HeatSimulation object displays: - thermal conductivity of structures (in grayscale colors), - specified boundary conditions (colored thick lines: yellow for TemperatureBC, green for HeatFluxBC, and red for ConvectionBC), - heat sources (in colored dotted hatching), - and monitors (in transparent yellow color).

# plot the two simulations
fig, ax = plt.subplots(1, 2, figsize=(10, 3))
sim.plot_heat_conductivity(x=0.0, ax=ax[0])
sim.plot_heat_conductivity(z=z_wg, ax=ax[1], monitor_alpha=0)

Now we can submit our heat simulation setup for solving on our servers.

heat_sim_data =, "thermal-ring")
17:52:11 CST Created task 'thermal-ring' with task_id
             'he-f317d26f-67bf-4974-89ae-77e7cbf85b19' 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.
17:52:12 CST Heat solver status: queued
17:52:16 CST Heat solver status: preprocess
17:53:34 CST Heat solver status: running
17:53:39 CST Heat solver status: postprocess
17:54:29 CST Heat solver status: success
17:54:42 CST loading simulation from simulation_data.hdf5

The result of this calculation are returned as an HeatSimulationData instance containing the calculated temperature field. Both unstructured and structured monitor can be used for visualizing resulting temperature fields.

fig, ax = plt.subplots(2, 3, figsize=(20, 10))

heat_sim_data.plot_field("all", x=0, ax=ax[0, 0])
heat_sim_data.plot_field("all", z=z_wg, ax=ax[0, 1])
heat_sim_data.plot_field("all", z=z_heater, ax=ax[0, 2])

heat_sim_data.plot_field("ugrid_vert", ax=ax[1, 0])
heat_sim_data.plot_field("ugrid_wg", ax=ax[1, 1])
heat_sim_data.plot_field("ugrid_heater", ax=ax[1, 2])


The unstructured monitors can also be used to visualize the unstructured grid used by heat solver. Note the difference between resulting grids for conformal and non-conformal monitors. In the former case the monitor plane conforms to finite element cells, while in the latter case the resulting dataset is the result of slicing unstructured tetrahedral grid.

fig, ax = plt.subplots(1, 3, figsize=(15, 5))

sim.plot_structures_heat_conductivity(z=z_wg, ax=ax[0], alpha=0.5)
heat_sim_data["ugrid_wg"].temperature.plot(ax=ax[0], field=False)
ax[0].set_xlim([-2, 2])
ax[0].set_ylim([4, 7])

sim.plot_structures_heat_conductivity(z=z_heater, ax=ax[1], alpha=0.5)
heat_sim_data["ugrid_heater"].temperature.plot(ax=ax[1], field=False)
ax[1].set_xlim([-2, 2])
ax[1].set_ylim([3, 7])

heat_sim_data["ugrid_vert"].temperature.plot(ax=ax[2], cmap="coolwarm")
ax[2].set_xlim([3, 7])


Tidy3D provides a basic functionality for visualization of unstructured data. If a more thorough inspection is required we recommend to export the results into a vtu file using .to_vtu()method, for example:


and processing it using dedicated visualization software like Paraview.

Perturbed Scene#

The obtained temperature field information can be applied to mediums present in the scene so that their optical properties are appropriately modified. This can be done by using function scene.perturbed_mediums_copy(temperature: SpatialDataArray).

perturbed_scene = scene.perturbed_mediums_copy(

We can compare the permittivity fields of the original and perturbed scenes. Note that both ring resonator and waveguides are affected by the heat.

fig, ax = plt.subplots(1, 2, figsize=(10, 5))

eps_lim = [12.08, 12.18]

scene.plot_structures_eps(z=z_wg, eps_lim=eps_lim, ax=ax[0])

perturbed_scene.plot_structures_eps(z=z_wg, eps_lim=eps_lim, ax=ax[1])


Optic simulations#

Now one can run optical simulation on both scenes to investigate the influence of heat on electromagnetic response.

As mentioned above, compared to the heat simulation, the optic simulation will have smaller dimensions in \(x\) and \(y\) directions. Also, since heater and wafer structures are not expected to have significant influence on propagation in the waveguides we will set the size of simulation in z direction to 3 \(\mu\)m.

l_optic_sim = l_heat_sim - 2 * heat_sim_buffer
w_optic_sim = w_heat_sim - 2 * heat_sim_buffer
h_optic_sim = 3

Define the wavelength range of interest and corresponding frequency values.

# wavelength range of interest
lambda_beg = 1.5
lambda_end = 1.6

# define pulse parameters
freq_beg = td.C_0 / lambda_end
freq_end = td.C_0 / lambda_beg
freq0 = (freq_beg + freq_end) / 2
fwidth = (freq_end - freq0) / 1.5

We will inject a waveguide mode into the upper left arm of the structure.

# calculate x-distance for mode source and mode monitors
wg_insert_x = ring_radius + 0.9 * x_buffer

mode_plane = td.Box(
    center=[-wg_insert_x, y_wg, z_wg],
    size=[0, 2, 2],

mode_source = td.ModeSource(
    source_time=td.GaussianPulse(freq0=freq0, fwidth=fwidth),
    num_freqs=7,  # using 7 (Chebyshev) points to approximate frequency dependence

We will monitor the field distribution at the central wavelength in the xy plane passing through waveguide structures and mode amplitudes in the through and drop ports.

# monitor steady state fields at central frequency over whole domain
# frequency is chosen based on simulation results for clearer demonstration
field_monitor = td.FieldMonitor(
    center=[0, 0, z_wg], size=[td.inf, td.inf, 0], freqs=[1.93831331e+14], name="field"

# monitor the mode amps on the output waveguide
lambdas_measure = np.linspace(lambda_beg, lambda_end, 301)
freqs_measure = td.C_0 / lambdas_measure[::-1]

through_mode_mnt = td.ModeMonitor(
    center=[wg_insert_x, y_wg, z_wg],
    name="mode through",

drop_mode_mnt = td.ModeMonitor(
    center=[-wg_insert_x, -y_wg, z_wg],
    name="mode drop",

Define the remaining optic simulation parameters such as grid resolution and total simulation time, and create Simulation instances for both the original and perturbed simulation scenes. Note that as we chose to ignore heater and wafer structures in optic simulations we explicitly exclude them (see structures=scene.structures[0:5]) to avoid accidentally placing them in the PML regions.

min_steps_per_wvl = 15
run_time = 8e-11

optic_sim = td.Simulation(
    size=[l_optic_sim, w_optic_sim, h_optic_sim],
    center=(0, 0, z_sim),
        min_steps_per_wvl=min_steps_per_wvl, wavelength=td.C_0 / freq0
    monitors=[field_monitor, through_mode_mnt, drop_mode_mnt],

perturbed_optic_sim = optic_sim.updated_copy(structures=perturbed_scene.structures[0:5])

Visualization to confirm simulation setup.

<Axes: title={'center': 'cross section at x=0.00'}, xlabel='y', ylabel='z'>
fig, ax = plt.subplots(1, 2, figsize=(10, 5))

optic_sim.plot_eps(z=z_wg, ax=ax[0])

perturbed_optic_sim.plot_eps(z=z_wg, ax=ax[1])


Before proceeding to running actual FDTD simulation, we can take a look at how waveguide modes are affected in the ring resonator.

from tidy3d.plugins.mode import ModeSolver
ring_mode_plane = td.Box(size=(2, 0, 1.9), center=(-ring_radius, 0, z_wg))
mode_solver_freqs = np.linspace(freq_beg, freq_end, 11)
mode_solver = ModeSolver(
    mode_spec=td.ModeSpec(num_modes=2, bend_radius=ring_radius, bend_axis=1),
mode_solver_perturbed = mode_solver.updated_copy(simulation=perturbed_optic_sim)
mode_data = mode_solver.solve()
mode_data_perturbed = mode_solver_perturbed.solve()
17:54:57 CST WARNING: Use the remote mode solver with subpixel averaging for    
             better accuracy through ''.        

In this case additional heat increased the propagation index of the waveguide mode of interest.

plt.legend(["no heat", "with heat"])

Now let us submit both optic simulations for calculation to our servers.

batch = web.Batch(
        "ring_resonator_no_heat": optic_sim,
        "ring_resonator_with_heat": perturbed_optic_sim,

batch_data =
17:55:13 CST WARNING: Simulation has 1.55e+06 time steps. The 'run_time' may be 
             unnecessarily large, unless there are very long-lived resonances.  
             Created task 'ring_resonator_no_heat' with task_id
             'fdve-82e157d9-10b7-4478-8f6b-96f056d75e95' and task_type 'FDTD'.
17:55:14 CST WARNING: Simulation has 1.55e+06 time steps. The 'run_time' may be 
             unnecessarily large, unless there are very long-lived resonances.  
             Created task 'ring_resonator_with_heat' with task_id
             'fdve-806f26c9-1159-4ad7-8ecc-b441f388c89e' and task_type 'FDTD'.
17:56:44 CST Started working on Batch.
17:56:58 CST Maximum FlexCredit cost: 15.492 for the whole batch.
             Use 'Batch.real_cost()' to get the billed FlexCredit cost after the
             Batch has completed.
19:25:40 CST Batch complete.

Visualization of the change in steady-state field distribution at the central frequency.

no_heat_data = batch_data["ring_resonator_no_heat"]
with_heat_data = batch_data["ring_resonator_with_heat"]
19:25:44 CST loading simulation from
19:26:06 CST loading simulation from
fig, ax = plt.subplots(1, 2, figsize=(10, 4))

no_heat_data.plot_field("field", "E", "abs", z=z_wg, ax=ax[0])
with_heat_data.plot_field("field", "E", "abs", z=z_wg, ax=ax[1])


Visualization of the change in mode amplitudes in through and drop ports due to thermal effects. As expected, heating the ring resonator results in shifting of the resonance frequencies.

fig, ax = plt.subplots(1, 2, figsize=(10, 4))

no_heat_data["mode through"].amps.sel(mode_index=0, direction="+").abs.plot(x="f", ax=ax[0])
with_heat_data["mode through"].amps.sel(mode_index=0, direction="+").abs.plot(x="f", ax=ax[0])
ax[0].legend(["no heat", "with heat"])
ax[0].set_title("Through port")
ax[0].set_ylim([0, 1])

no_heat_data["mode drop"].amps.sel(mode_index=0, direction="-").abs.plot(x="f", ax=ax[1])
with_heat_data["mode drop"].amps.sel(mode_index=0, direction="-").abs.plot(x="f", ax=ax[1])
ax[1].legend(["no heat", "with heat"])
ax[1].set_title("Drop port")
ax[1].set_ylim([0, 1])

[ ]: