Photo-thermal optical control in silicon waveguides#
This example demonstrates the usage of the Tidy3D to calculate absorbed optical power, its resultant heating, and the corresponding optical heat perturbation on the silicon-on-insulator (SOI) platform. Our setup is based on All-optic control using a photo-thermal heater in Si photonics by L. Li, T. Tamanuki, and T. Baba. We build two parallel silicon rib waveguides separated by a gap of varying length. One waveguide (the βcontrolβ waveguide) has a triangular-shaped heavily doped region. This waveguide then absorbs an optical βcontrolβ signal, and the absorbed optical power heats both waveguides. The second waveguide (the βsignalβ waveguide) then has its refractive index perturbed by the change in temperature, and we analyze the resulting phase change of a signal sent through this waveguide. This analysis is done for three gap lengths: 0.5 ΞΌm, 1.0 ΞΌm, and 1.5 ΞΌm. We first create the perturbation mediums with Tidy3D heat and charge perturbations in order to capture the properties of doped silicon. We then record the absorbed power and permittivity of the control waveguide, plot the resulting temperature change along the propagation (z) direction and along the y direction, and finally plot the heat-perturbed refractive index along z as well as the signal waveguideβs phase change. Thus our workflow will be: Control optical simulation \(\rightarrow\) Heat simulation \(\rightarrow\) Signal optical simulation
[1]:
import matplotlib.pyplot as plt
import numpy as np
import tidy3d as td
import tidy3d.web as web
First we describe some parameters that will remain constant throughout our simulations.
[2]:
wvl0 = 1.55
freq0 = td.C_0 / wvl0
length = 30
wg_width = 0.45
wg_height = 0.11
slab_height = 0.11
Perturbation Medium Setup#
Because perturbing loss can result in numerical instabilities in our simulations, Tidy3D adds a warning about this when we define perturbation mediums. Since our perturbations will be small enough for stability, we will turn off Tidy3D warnings with the following command.
[3]:
td.config.logging_level = "ERROR"
Here we define the `MultiPhysicsMediums <https://docs.flexcompute.com/projects/tidy3d/en/v2.9.0/api/_autosummary/tidy3d.components.material.multi_physics.MultiPhysicsMedium.html#tidy3d.components.material.multi_physics.MultiPhysicsMedium>`__ of our simulation. To capture the absorption of doped silicon, we will perturb our silicon according to Tidy3Dβs charge perturbation. Here we will use the
`Nedeljkovic-Soref-Mashanovich <https://docs.flexcompute.com/projects/tidy3d/en/v2.9.0/api/_autosummary/tidy3d.NedeljkovicSorefMashanovich.html>`__ model with Li et. al.βs doping concentration.
[4]:
ref_temp = 300
doping_concentration = 2.273e20
Si_heat_coeff = 0.00018
NSM = td.NedeljkovicSorefMashanovich(ref_freq=freq0)
# for maximum generality, we include charge in our index perturbations
Si = td.MultiPhysicsMedium(
name="Si",
optical=td.PerturbationMedium.from_nk(
n=3.5,
k=0,
freq=freq0,
perturbation_spec=td.IndexPerturbation(
delta_n=td.ParameterPerturbation(
heat=td.LinearHeatPerturbation(temperature_ref=ref_temp, coeff=Si_heat_coeff),
charge=NSM.delta_n().charge,
),
delta_k=td.ParameterPerturbation(charge=NSM.delta_k().charge),
freq=freq0,
),
),
heat=td.SolidSpec(
conductivity=148e-6, # W / (um * K)
capacity=0.711, # J / (kg * K)
),
)
Si_doped = td.MultiPhysicsMedium(
name="Si doped",
optical=td.PerturbationMedium.from_nk(
n=3.5 + NSM.delta_n().charge.sample(doping_concentration, 0),
k=0 + NSM.delta_k().charge.sample(doping_concentration, 0),
freq=freq0,
perturbation_spec=Si.optical.perturbation_spec,
),
heat=Si.heat,
)
SiO2 = td.MultiPhysicsMedium(
name="SiO2",
optical=td.PerturbationMedium(
permittivity=1.44**2,
perturbation_spec=td.IndexPerturbation(
delta_n=td.ParameterPerturbation(
heat=td.LinearHeatPerturbation(temperature_ref=ref_temp, coeff=0.000012)
),
freq=freq0,
),
),
heat=td.SolidSpec(
conductivity=1.38e-6, # W / (um * K)
capacity=709, # J / (kg * K)
),
)
First Optical Simulation#
Here we create a function that takes in a gap length and returns a Tidy3D simulation with the corresponding two-waveguide structure. The simulations returned will be our first optical simulations - that is, a simulation where we launch a TE mode and record the field and permittivity along the control waveguide, as this will be used to compute the absorbed power. We also add a field monitor along the propagation direction to visualize the optical power being absorbed along the waveguide.
[5]:
def make_sim(gap):
z_min, z_max = -length / 2.0, length / 2.0
sim_buffer = 4
# Waveguides (above slab) in y; separated in x
wg_signal_x_min, wg_signal_x_max = gap / 2.0, gap / 2.0 + wg_width
wg_control_x_min, wg_control_x_max = -gap / 2.0 - wg_width, -gap / 2.0
wg_signal_box = td.Box.from_bounds(
rmin=(wg_signal_x_min, 0.0, -td.inf), rmax=(wg_signal_x_max, wg_height, td.inf)
)
wg_control_box = td.Box.from_bounds(
rmin=(wg_control_x_min, 0.0, -td.inf),
rmax=(wg_control_x_max, wg_height, td.inf),
)
# Slab (below guides)
slab_y_min, slab_y_max = -slab_height, 0.0
slab_x_min = wg_control_x_min
slab_x_max = wg_signal_x_max + 0.5
slab_box = td.Box.from_bounds(
rmin=(slab_x_min, slab_y_min, -td.inf), rmax=(slab_x_max, slab_y_max, td.inf)
)
# Triangular prism: hypotenuse aligned with z (length) direction
tri_vertices = [
(wg_control_x_min, z_max),
(wg_control_x_min, z_min),
(wg_control_x_max, z_min),
]
tri_geom = td.PolySlab(vertices=tri_vertices, axis=1, slab_bounds=(0.0, wg_height))
slab = td.Structure(geometry=slab_box, medium=Si, name="slab")
wg_signal = td.Structure(geometry=wg_signal_box, medium=Si, name="wg_signal")
wg_control = td.Structure(geometry=wg_control_box, medium=Si, name="wg_control")
tri = td.Structure(geometry=tri_geom, medium=Si_doped, name="doped triangle")
slab_width = slab_x_max - slab_x_min
margin_x = max(0.5 * wvl0 + 0.1, 0.5)
sim_center = (0.5 * (slab_x_max + slab_x_min), 0.5 * (wg_height + slab_y_min), 0.0)
sim_size = (
slab_width + 2 * margin_x,
wvl0 + (wg_height - slab_y_min),
length + sim_buffer,
)
tri_min, tri_max = tri.geometry.bounds
tri_center = tuple(0.5 * (a + b) for a, b in zip(tri_min, tri_max))
tri_size = tuple(a - b for a, b in zip(tri_max, tri_min))
tri_monitor = td.FieldMonitor(
center=tri_center,
size=tri_size,
fields=["Ex", "Ey", "Ez", "Hx", "Hy", "Hz"],
freqs=[freq0],
name="field",
)
perm_monitor = td.PermittivityMonitor(
center=tri_center, size=tri_size, freqs=[freq0], name="tri_permittivity"
)
plane_monitor = td.FieldMonitor(
center=(sim_center[0], wg_height / 2, sim_center[2]),
size=(sim_size[0], 0, sim_size[2]),
fields=["Ex", "Ey", "Ez", "Hx", "Hy", "Hz"],
freqs=[freq0],
name="midplane_field",
)
clearance = min(0.05, 0.25 * gap) if gap > 0 else 0.01
clearance = max(clearance, 1e-3)
buffer_x = max(0.2, 0.5 * wg_width)
buffer_y = max(0.1, 0.2 * wg_height)
extra_x = max(0.3, 0.5 * wg_width)
extra_y = max(0.5, 1.5 * wg_height)
source_x_min = wg_control_x_min - (buffer_x + extra_x)
source_x_max = wg_signal_x_min - clearance
if source_x_max <= source_x_min:
mid = 0.5 * (wg_control_x_max + wg_signal_x_min)
source_x_min = min(wg_control_x_min - (0.1 * wg_width + extra_x), mid - 0.05)
source_x_max = max(mid, wg_signal_x_min - 1e-3)
source_y_min = -buffer_y - extra_y
source_y_max = wg_height + buffer_y + extra_y
mode_center = (
0.5 * (source_x_min + source_x_max),
0.5 * (source_y_min + source_y_max),
z_max + 0.1,
)
mode_size = (source_x_max - source_x_min, source_y_max - source_y_min, 0.0)
mode_source = td.ModeSource(
center=mode_center,
size=mode_size,
direction="-",
source_time=td.GaussianPulse(freq0=freq0, fwidth=freq0 / 10),
mode_spec=td.ModeSpec(num_modes=1),
)
simulation = td.Simulation(
center=sim_center,
size=sim_size,
structures=[slab, wg_signal, wg_control, tri],
sources=[mode_source],
monitors=[tri_monitor, perm_monitor, plane_monitor],
grid_spec=td.GridSpec.auto(min_steps_per_wvl=20, wavelength=wvl0),
medium=SiO2.optical,
run_time=5e-13,
)
return simulation
[6]:
sim_0p5 = make_sim(0.5)
sim_1p0 = make_sim(1.0)
sim_1p5 = make_sim(1.5)
We visualize a constructed simulation to ensure the setup is correct.
[7]:
fig, ax = plt.subplots(1, 2, figsize=(12, 5))
sim_0p5.plot(y=wg_height / 2, ax=ax[0], monitor_alpha=0.1)
sim_0p5.plot(z=sim_0p5.center[0], ax=ax[1], monitor_alpha=0.1)
ax[0].set_aspect(0.1)
plt.tight_layout()
We run the three constructed simulations.
[8]:
sims = {"0.5": sim_0p5, "1.0": sim_1p0, "1.5": sim_1p5}
run_results = {}
for label, sim in sims.items():
result = web.run(sim, task_name=f"sim_gap_{label}")
run_results[label] = result
β simulation_data.hdf5.gz ββββββββββββ 100.0% β’ 25.1/25.1 β’ 3.4 MB/s β’ 0:00:00 MB
11:15:44 UTC Loading results from simulation_data.hdf5
We will visualize the mode getting absorbed as it propagates through the control waveguide for each simulation.
[9]:
fig, axes = plt.subplots(1, 3, figsize=(15, 10))
for ax, label in zip(axes, ["0.5", "1.0", "1.5"]):
run_results[label].plot_field("midplane_field", field_name="S", val="abs", f=freq0, ax=ax)
ax.set_title(f"Midplane |S| (gap {label})")
ax.set_aspect(0.4)
plt.tight_layout()
Heat Simulation#
We will now use the recorded field and permittivity along the control waveguide to compute the absorbed power. The absorbed power is given by the formula
\(P = \frac{1}{2}\omega|E|^2\) \(im(\epsilon)\)
We construct a function that creates a heat source as a `SpatialDataArray <https://docs.flexcompute.com/projects/tidy3d/en/latest/api/_autosummary/tidy3d.SpatialDataArray.html>`__, using the coordinates of the field monitor and interpolating the permittivity monitor data to these points.
[20]:
def absorbed_power(field_data, permittivity_data, freq):
"""Return 0.5 * 2Ο * freq * |E|^2 * imag(permittivity) on the field grid."""
coords = field_data.Ex.coords
coords_dict = {dim: coords.get(dim) for dim in "xyz"}
perm_interp_x = permittivity_data.eps_xx.interp(**coords_dict)
perm_interp_y = permittivity_data.eps_yy.interp(**coords_dict)
perm_interp_z = permittivity_data.eps_zz.interp(**coords_dict)
perm_imag_x = np.imag(perm_interp_x).sel(f=freq, method="nearest") * td.EPSILON_0
perm_imag_y = np.imag(perm_interp_y).sel(f=freq, method="nearest") * td.EPSILON_0
perm_imag_z = np.imag(perm_interp_z).sel(f=freq, method="nearest") * td.EPSILON_0
absorption_x = (
0.5
* 2
* np.pi
* freq
* np.abs(field_data.Ex).sel(f=freq, method="nearest") ** 2
* perm_imag_x
)
absorption_y = (
0.5
* 2
* np.pi
* freq
* np.abs(field_data.Ey).sel(f=freq, method="nearest") ** 2
* perm_imag_y
)
absorption_z = (
0.5
* 2
* np.pi
* freq
* np.abs(field_data.Ez).sel(f=freq, method="nearest") ** 2
* perm_imag_z
)
absorption = absorption_x + absorption_y + absorption_z
return td.SpatialDataArray(absorption, coords=coords_dict)
[11]:
heat_rate_0p5 = absorbed_power(
run_results["0.5"]["field"], run_results["0.5"]["tri_permittivity"], freq0
)
heat_rate_1p0 = absorbed_power(
run_results["1.0"]["field"], run_results["1.0"]["tri_permittivity"], freq0
)
heat_rate_1p5 = absorbed_power(
run_results["1.5"]["field"], run_results["1.5"]["tri_permittivity"], freq0
)
Using the absorbed power data we calculated, we create the heat source by assigning it to the optical waveguide structure.
[12]:
heat_source_0p5 = td.HeatSource(structures=["wg_control"], rate=heat_rate_0p5)
heat_source_1p0 = td.HeatSource(structures=["wg_control"], rate=heat_rate_1p0)
heat_source_1p5 = td.HeatSource(structures=["wg_control"], rate=heat_rate_1p5)
Next we will construct a temperature monitor around the signal waveguide so we can heat-perturb the signal structure.
[13]:
def make_temp_monitor(sim):
rmin, rmax = sim.structures[-1].geometry.bounds
size_z = rmax[2] - rmin[2]
temp_monitor = td.TemperatureMonitor(
center=sim.structures[1].geometry.center,
size=(sim.sources[0].size[0], sim.sources[0].size[1], size_z),
name="temp",
unstructured=True,
)
return temp_monitor
temp_monitor_0p5 = make_temp_monitor(sims["0.5"])
temp_monitor_1p0 = make_temp_monitor(sims["1.0"])
temp_monitor_1p5 = make_temp_monitor(sims["1.5"])
We also define the discretization for our heat simulations.
[14]:
# parameters for meshing
dl_min = slab_height / 3
dl_max = 3 * dl_min
grid_spec = td.DistanceUnstructuredGrid(
dl_interface=dl_min,
dl_bulk=dl_max,
distance_interface=3 * dl_min,
distance_bulk=slab_height,
)
We define boundary specs for our heat simulations.
[15]:
boundaries = td.HeatBoundarySpec(
placement=td.SimulationBoundary(surfaces=["y-", "y+", "x-", "x+"]),
condition=td.TemperatureBC(temperature=ref_temp),
)
We now combine the structures, computed heat source, temperature monitor, discretization, and boundary conditions into Tidy3D Heat simulations.
[16]:
heat_sim_0p5 = td.HeatChargeSimulation(
center=sims["0.5"].center,
size=[sims["0.5"].size[0] * 1.5, sims["0.5"].size[1] * 1.5, sims["0.5"].size[2]],
medium=SiO2,
structures=sims["0.5"].structures,
sources=[heat_source_0p5],
monitors=[temp_monitor_0p5],
grid_spec=grid_spec,
boundary_spec=[boundaries],
)
heat_sim_1p0 = td.HeatChargeSimulation(
center=sims["1.0"].center,
size=[sims["1.0"].size[0] * 1.5, sims["1.0"].size[1] * 1.5, sims["1.0"].size[2]],
medium=SiO2,
structures=sims["1.0"].structures,
sources=[heat_source_1p0],
monitors=[temp_monitor_1p0],
grid_spec=grid_spec,
boundary_spec=[boundaries],
)
heat_sim_1p5 = td.HeatChargeSimulation(
center=sims["1.5"].center,
size=[sims["1.5"].size[0] * 1.5, sims["1.5"].size[1] * 1.5, sims["1.5"].size[2]],
medium=SiO2,
structures=sims["1.5"].structures,
sources=[heat_source_1p5],
monitors=[temp_monitor_1p5],
grid_spec=grid_spec,
boundary_spec=[boundaries],
)
We run the heat simulations.
[17]:
heat_data_0p5 = web.run(heat_sim_0p5, task_name="heat_sim_0p5")
heat_data_1p0 = web.run(heat_sim_1p0, task_name="heat_sim_1p0")
heat_data_1p5 = web.run(heat_sim_1p5, task_name="heat_sim_1p5")
β simulation_data.hdf5.gz βββββββββββ 100.0% β’ 16.5/16.5 β’ 19.9 MB/s β’ 0:00:00 MB
11:16:43 UTC Loading results from simulation_data.hdf5
To visualize our temperature distribution from the heat source, we can add an interactive plot to show a cross-section of the signal waveguide sliding along the propagation (z) axis.
[36]:
# adjust for different sims
temp_data = heat_data_0p5
import ipywidgets as widgets
zs = np.linspace(
temp_data["temp"].monitor.center[2] - temp_data["temp"].monitor.size[2] / 2,
temp_data["temp"].monitor.center[2] + temp_data["temp"].monitor.size[2] / 2,
61,
)
hmin = temp_data["temp"].monitor.center[0] - temp_data["temp"].monitor.size[0] / 2
hmax = temp_data["temp"].monitor.center[0] + temp_data["temp"].monitor.size[0] / 2
vmin = temp_data["temp"].monitor.center[1] - temp_data["temp"].monitor.size[1] / 2
vmax = temp_data["temp"].monitor.center[1] + temp_data["temp"].monitor.size[1] / 2
field_slider = widgets.FloatSlider(
value=4,
min=zs[0],
max=zs[-1],
step=None,
disabled=False,
continuous_update=True,
orientation="horizontal",
readout=True,
readout_format=".3f",
)
def update_wg_plot(x):
if x not in zs:
x = min(zs, key=lambda zz: abs(zz - x))
fig, ax = plt.subplots(figsize=(6, 5))
slice_temp = temp_data["temp"].temperature.plane_slice(axis=2, pos=x)
slice_temp.plot(
cmap="jet", grid=False, vmax=np.max(temp_data["temp"].temperature.values), ax=ax
)
temp_data.simulation.plot_structures(
z=x, fill=False, ax=ax, hlim=[hmin, hmax], vlim=[vmin, vmax]
)
ax.set_title(f"Temperature at z = {x:.3f} Β΅m")
plt.tight_layout()
plt.show()
widgets.interactive(update_wg_plot, x=field_slider)
[36]:
We will also plot the temperature profile through the middle of the structure to show the distribution along the waveguide all at once.
[19]:
fig, ax = plt.subplots(1, 3, figsize=(6, 5))
heat_data_0p5["temp"].temperature.plane_slice(axis=1, pos=0).plot(grid=False, cmap="jet", ax=ax[0])
heat_data_1p0["temp"].temperature.plane_slice(axis=1, pos=0).plot(grid=False, cmap="jet", ax=ax[1])
heat_data_1p5["temp"].temperature.plane_slice(axis=1, pos=0).plot(grid=False, cmap="jet", ax=ax[2])
ax[0].set_title("0.5 Β΅m gap")
ax[1].set_title("1.0 Β΅m gap")
ax[2].set_title("1.5 Β΅m gap")
fig.suptitle("Temperature profiles at y=0")
plt.tight_layout()
plt.show()
From this temperature distribution, we see that the middle section should have more of a perturbation. To verify this, we will visualize the change in refractive index along the propagation direction of the signal waveguide.
[20]:
heat_line_0p5 = heat_data_0p5["temp"].temperature.line_slice(
axis=2, pos=sims["0.5"].structures[1].geometry.center
)
heat_line_1p0 = heat_data_1p0["temp"].temperature.line_slice(
axis=2, pos=sims["1.0"].structures[1].geometry.center
)
heat_line_1p5 = heat_data_1p5["temp"].temperature.line_slice(
axis=2, pos=sims["1.5"].structures[1].geometry.center
)
plt.plot(
heat_line_0p5.z,
np.squeeze(heat_line_0p5.data) * Si_heat_coeff,
color="blue",
label="0.5 Β΅m gap",
)
plt.plot(
heat_line_1p0.z,
np.squeeze(heat_line_1p0.data) * Si_heat_coeff,
color="mediumblue",
label="1.0 Β΅m gap",
)
plt.plot(
heat_line_1p5.z,
np.squeeze(heat_line_1p5.data) * Si_heat_coeff,
color="darkblue",
label="1.5 Β΅m gap",
)
plt.xlabel("z (Β΅m)")
plt.ylabel("Ξn")
plt.title("Refractive Index Change in Signal Waveguide vs. z")
plt.legend()
plt.show()
Final Optical Simulation#
We are now in a position to see the effect on the signal waveguide from the control signal through the control waveguide. We create a function that modifies our first optical simulations into ones that launch and measure the TE mode through the signal waveguide. Then, using the temperature data of the previous heat simulations, we create heat-perturbed copies to compute the phase changes.
[21]:
def update_sim_signal_source(sim):
new_source = sim.sources[0].updated_copy(
center=(
sim.structures[1].geometry.center[0],
sim.sources[0].center[1],
sim.sources[0].center[2],
)
)
new_monitor = td.ModeMonitor(
center=(
sim.structures[1].geometry.center[0],
sim.structures[1].geometry.center[1],
sim.center[2] - sim.size[2] / 2 + 1,
),
size=new_source.size,
name="mode",
freqs=[freq0],
mode_spec=new_source.mode_spec,
)
new_sim = sim.updated_copy(sources=[new_source], monitors=[new_monitor], medium=SiO2.optical)
return new_sim
sims_opt, sims_opt_perturbed = {}, {}
sims_opt["0.5"] = update_sim_signal_source(sims["0.5"])
sims_opt["1.0"] = update_sim_signal_source(sims["1.0"])
sims_opt["1.5"] = update_sim_signal_source(sims["1.5"])
sims_opt_perturbed["0.5"] = sims_opt["0.5"].perturbed_mediums_copy(
temperature=heat_data_0p5["temp"].temperature
)
sims_opt_perturbed["1.0"] = sims_opt["1.0"].perturbed_mediums_copy(
temperature=heat_data_1p0["temp"].temperature
)
sims_opt_perturbed["1.5"] = sims_opt["1.5"].perturbed_mediums_copy(
temperature=heat_data_1p5["temp"].temperature
)
We create batches for the perturbed and unperturbed simulations and run them.
[22]:
opt_unperturbed_batch = td.web.Batch(simulations=sims_opt)
opt_perturbed_batch = td.web.Batch(simulations=sims_opt_perturbed)
# Run all the simulations and get the results.
opt_unperturbed = opt_unperturbed_batch.run()
opt_perturbed = opt_perturbed_batch.run()
0.5 β success βββββββββββββββββββββββββ 100% 0:00:02 1.0 β success βββββββββββββββββββββββββ 100% 0:00:02 1.5 β success βββββββββββββββββββββββββ 100% 0:00:02
11:20:12 UTC Batch complete.
We now compute the phase changes between the heated and unheated simulations for each gap, and plot the phase change vs gap length.
[23]:
phase_0p5_unp = np.angle(opt_unperturbed["0.5"]["mode"].amps.sel(direction="-").values)
phase_0p5_p = np.angle(opt_perturbed["0.5"]["mode"].amps.sel(direction="-").values)
phase_1p0_unp = np.angle(opt_unperturbed["1.0"]["mode"].amps.sel(direction="-").values)
phase_1p0_p = np.angle(opt_perturbed["1.0"]["mode"].amps.sel(direction="-").values)
phase_1p5_unp = np.angle(opt_unperturbed["1.5"]["mode"].amps.sel(direction="-").values)
phase_1p5_p = np.angle(opt_perturbed["1.5"]["mode"].amps.sel(direction="-").values)
phase_shift_0p5 = np.squeeze(phase_0p5_p - phase_0p5_unp)
phase_shift_1p0 = np.squeeze(phase_1p0_p - phase_1p0_unp)
phase_shift_1p5 = np.squeeze(phase_1p5_p - phase_1p5_unp)
phase_shifts = [phase_shift_0p5, phase_shift_1p0, phase_shift_1p5]
plt.plot(
[0.5, 1.0, 1.5],
phase_shifts,
marker="o",
linestyle="-",
linewidth=2,
color="tab:blue",
)
plt.xticks([0.5, 1.0, 1.5])
plt.xlabel("Gap (Β΅m)")
plt.ylabel("Phase Shift (rad)")
plt.title("Phase Shift vs. Gap")
plt.xlim(0.4, 1.6)
plt.show()
[ ]: