Scattering Matrix Plugin
Contents
Scattering Matrix Plugin#
Run this notebook in your browser using Binder.
This notebook will give a demo of the tidy3d ComponentModeler plugin used to compute scattering matrix elements.
[1]:
# make sure notebook plots inline
%matplotlib inline
# standard python imports
import numpy as np
import matplotlib.pyplot as plt
import os
import gdspy
# tidy3D imports
import tidy3d as td
from tidy3d import web
# set tidy3d to only print error information to reduce verbosity
td.logging_level = 'error'
Setup#
We will simulate a directional coupler, similar to the GDS and Parameter scan tutorials.
Let’s start by setting up some basic parameters.
[2]:
# wavelength / frequency
lambda0 = 1.550 # all length scales in microns
freq0 = td.constants.C_0 / lambda0
fwidth = freq0 / 10
# Spatial grid specification
grid_spec = td.GridSpec.auto(min_steps_per_wvl=14, wavelength=lambda0)
# Permittivity of waveguide and substrate
wg_n = 3.48
sub_n = 1.45
mat_wg = td.Medium(permittivity=wg_n**2)
mat_sub = td.Medium(permittivity=sub_n**2)
# Waveguide dimensions
# Waveguide height
wg_height = 0.22
# Waveguide width
wg_width = 1.0
# Waveguide separation in the beginning/end
wg_spacing_in = 8
# length of coupling region (um)
coup_length = 6.0
# spacing between waveguides in coupling region (um)
wg_spacing_coup = 0.05
# Total device length along propagation direction
device_length = 100
# Length of the bend region
bend_length = 16
# Straight waveguide sections on each side
straight_wg_length = 4
# space between waveguide and PML
pml_spacing = 2
Define waveguide bends and coupler#
Here is where we define our directional coupler shape programmatically in terms of the geometric parameters
[3]:
def bend_pts(bend_length, width, npts=10):
""" Set of points describing a tanh bend from (0, 0) to (length, width)"""
x = np.linspace(0, bend_length, npts)
y = width*(1 + np.tanh(6*(x/bend_length - 0.5)))/2
return np.stack((x, y), axis=1)
def arm_pts(length, width, coup_length, bend_length, npts_bend=30):
""" Set of points defining one arm of an integrated coupler """
### Make the right half of the coupler arm first
# Make bend and offset by coup_length/2
bend = bend_pts(bend_length, width, npts_bend)
bend[:, 0] += coup_length / 2
# Add starting point as (0, 0)
right_half = np.concatenate(([[0, 0]], bend))
# Add an extra point to make sure waveguide is straight past the bend
right_half = np.concatenate((right_half, [[right_half[-1, 0] + 0.1, width]]))
# Add end point as (length/2, width)
right_half = np.concatenate((right_half, [[length/2, width]]))
# Make the left half by reflecting and omitting the (0, 0) point
left_half = np.copy(right_half)[1:, :]
left_half[:, 0] = -left_half[::-1, 0]
left_half[:, 1] = left_half[::-1, 1]
return np.concatenate((left_half, right_half), axis=0)
def make_coupler(
length,
wg_spacing_in,
wg_width,
wg_spacing_coup,
coup_length,
bend_length,
npts_bend=30):
""" Make an integrated coupler using the gdspy FlexPath object. """
# Compute one arm of the coupler
arm_width = (wg_spacing_in - wg_width - wg_spacing_coup)/2
arm = arm_pts(length, arm_width, coup_length, bend_length, npts_bend)
# Reflect and offset bottom arm
coup_bot = np.copy(arm)
coup_bot[:, 1] = -coup_bot[::-1, 1] - wg_width/2 - wg_spacing_coup/2
# Offset top arm
coup_top = np.copy(arm)
coup_top[:, 1] += wg_width/2 + wg_spacing_coup/2
# Create waveguides as GDS paths
path_bot = gdspy.FlexPath(coup_bot, wg_width, layer=1, datatype=0)
path_top = gdspy.FlexPath(coup_top, wg_width, layer=1, datatype=1)
return [path_bot, path_top]
Create Base Simulation#
The scattering matrix tool requires the “base” Simulation (without the modal sources or monitors used to compute S-parameters), so we will construct that now.
We generate the structures and add a FieldMonitor so we can inspect the field patterns.
[4]:
gdspy.current_library = gdspy.GdsLibrary()
lib = gdspy.GdsLibrary()
# Geometry must be placed in GDS cells to import into Tidy3D
coup_cell = lib.new_cell('Coupler')
substrate = gdspy.Rectangle(
(-device_length/2, -wg_spacing_in/2-10),
(device_length/2, wg_spacing_in/2+10),
layer=0)
coup_cell.add(substrate)
# Add the coupler to a gdspy cell
gds_coup = make_coupler(
device_length,
wg_spacing_in,
wg_width,
wg_spacing_coup,
coup_length,
bend_length)
coup_cell.add(gds_coup)
# Substrate
[oxide_geo] = td.PolySlab.from_gds(
gds_cell=coup_cell,
gds_layer=0,
gds_dtype=0,
slab_bounds=(-10, 0),
axis=2)
oxide = td.Structure(
geometry=oxide_geo,
medium=mat_sub)
# Waveguides (import all datatypes if gds_dtype not specified)
coupler1_geo, coupler2_geo = td.PolySlab.from_gds(
gds_cell=coup_cell,
gds_layer=1,
slab_bounds=(0, wg_height),
axis=2)
coupler1 = td.Structure(
geometry=coupler1_geo,
medium=mat_wg
)
coupler2 = td.Structure(
geometry=coupler2_geo,
medium=mat_wg
)
# Simulation size along propagation direction
sim_length = 2*straight_wg_length + 2*bend_length + coup_length
# Spacing between waveguides and PML
sim_size = [
sim_length,
wg_spacing_in + wg_width + 2*pml_spacing,
wg_height + 2*pml_spacing]
# source
src_pos = sim_length/2 - straight_wg_length/2
# in-plane field monitor (optional, increases required data storage)
domain_monitor = td.FieldMonitor(
center = [0,0,wg_height/2],
size = [td.inf, td.inf, 0],
freqs = [freq0],
name='field'
)
# initialize the simulation
sim = td.Simulation(
size=sim_size,
grid_spec=grid_spec,
structures=[oxide, coupler1, coupler2],
sources=[],
monitors=[domain_monitor],
run_time=50/fwidth,
boundary_spec=td.BoundarySpec.all_sides(boundary=td.PML())
)
[17:14:19] WARNING No sources in simulation. simulation.py:406
[5]:
f, (ax1, ax2) = plt.subplots(1, 2, tight_layout=True, figsize=(15, 10))
ax1 = sim.plot(z=wg_height/2, ax=ax1)
ax2 = sim.plot(x=src_pos, ax=ax2)
<Figure size 1080x720 with 2 Axes>
Setting up Scattering Matrix Tool#
Now, to use the S matrix tool, we need to defing the spatial extent of the “ports” of our system using Port objects.
These ports will be converted into modal sources and monitors later, so they require both some mode specification and a definition of the direction that points into the system.
We’ll also give them names to refer to later.
[6]:
from tidy3d.plugins.smatrix.smatrix import Port
port_right_top = Port(
center=[src_pos, wg_spacing_in / 2, wg_height / 2],
size=[0, 4, 2],
mode_spec=td.ModeSpec(num_modes=2),
direction='-',
name='right_top')
port_right_bot = Port(
center=[src_pos, -wg_spacing_in / 2, wg_height / 2],
size=[0, 4, 2],
mode_spec=td.ModeSpec(num_modes=2),
direction='-',
name='right_bot')
port_left_top = Port(
center=[-src_pos, wg_spacing_in / 2, wg_height / 2],
size=[0, 4, 2],
mode_spec=td.ModeSpec(num_modes=2),
direction='+',
name='left_top')
port_left_bot = Port(
center=[-src_pos, -wg_spacing_in / 2, wg_height / 2],
size=[0, 4, 2],
mode_spec=td.ModeSpec(num_modes=2),
direction='+',
name='left_bot')
ports = [port_right_top, port_right_bot, port_left_top, port_left_bot]
Next, we will add the base simulation and ports to the ComponentModeler, along with the frequency of interest and a name for saving the batch of simulations that will get created later.
[7]:
from tidy3d.plugins.smatrix.smatrix import ComponentModeler
modeler = ComponentModeler(simulation=sim, ports=ports, freq=freq0)
We can plot the simulation with all of the ports as sources to check things are set up correctly.
[8]:
f, (ax1, ax2) = plt.subplots(1, 2, tight_layout=True, figsize=(15, 10))
ax1 = modeler.plot_sim(z=wg_height/2, ax=ax1)
ax2 = modeler.plot_sim(x=src_pos, ax=ax2)
<Figure size 1080x720 with 2 Axes>
Solving for the S matrix#
With the component modeler defined, we may call it’s .solve() method to run a batch of simulations to compute the S matrix. The tool will loop through each port and create one simulation per mode index (as defined by the mode specifications) where a unique modal source is injected. Each of the ports will also be converted to mode monitors to measure the mode amplitudes and normalization.
[9]:
smatrix = modeler.run(path_dir='data')
[17:14:45] Started working on Batch. container.py:384
[17:19:51] Batch complete. container.py:405
Working with Scattering Matrix#
The scattering matrix returned by the solve is actually a nested dictionarty relating the port names. For example smatrix[name1][name2] gives an array of shape (m, n) relating the coupling between the m modes injected into port name1 with the n modes measured in port name2.
For example, here each waveguide has 2 modes, so we can compute the coupling between the 2 input modes in left_top with the 2 output modes in right_bot as:
[10]:
smatrix['left_top']['right_bot']
array([[0.06405768+0.07115856j, 0.00669322-0.00306965j], [0.00574468-0.00463113j, 0.35245579+0.43931115j]])
Alternatively, we can convert this into a numpy array:
[11]:
blocks_cols = []
for port_name_in, val_in in smatrix.items():
blocks_rows = []
for port_name_out, S_in_out in val_in.items():
blocks_rows.append(S_in_out)
blocks_cols.append(np.block(blocks_rows))
S = np.concatenate(blocks_cols)
print(S.shape)
(8, 8)
We can inspect S and note that the diagonal elements are very small indicating low backscattering.
[12]:
print(abs(S)**2)
[[1.19474563e-04 3.84524744e-05 1.30710199e-06 1.71791191e-05
9.88835458e-01 3.45354262e-04 9.16672365e-03 5.51433494e-05]
[3.84105732e-05 1.81103843e-04 1.75045775e-05 1.33686022e-04
3.46289767e-04 6.66286303e-01 5.53308841e-05 3.17179490e-01]
[1.30942979e-06 1.74195952e-05 1.18416747e-04 3.83326123e-05
9.16690648e-03 5.42348612e-05 9.88832771e-01 3.45066232e-04]
[1.74761613e-05 1.33721545e-04 3.70578836e-05 1.78586736e-04
5.44390178e-05 3.17219231e-01 3.46296752e-04 6.66280941e-01]
[9.88824386e-01 3.47533504e-04 9.16692743e-03 5.42219524e-05
1.19366910e-04 3.83262152e-05 1.32683404e-06 1.74373165e-05]
[3.46757462e-04 6.66369136e-01 5.44486295e-05 3.17219373e-01
3.70781649e-05 1.78675694e-04 1.74779135e-05 1.33136106e-04]
[9.16639762e-03 5.51809273e-05 9.88830888e-01 3.45410542e-04
1.32680003e-06 1.71683153e-05 1.20300400e-04 3.84957392e-05]
[5.53905501e-05 3.17215450e-01 3.46314313e-04 6.66283607e-01
1.75122371e-05 1.33131540e-04 3.84814298e-05 1.79949441e-04]]
Summing each rows of the matrix should give 1.0 if no power was lost.
[13]:
np.sum(abs(S)**2, axis=0)
array([0.9985696 , 0.984358 , 0.99857286, 0.9842704 , 0.99857838, 0.98427243, 0.99857871, 0.98422966])
Finally, we can check whether S is close to unitary as expected.
S times it’s Hermitian conjugate should be the identy matrix.
[14]:
mat = S @ (np.conj(S.T))
[15]:
f, (ax1, ax2, ax3) = plt.subplots(1, 3, tight_layout=True, figsize=(12, 3.5))
imabs = ax1.imshow(abs(mat))
imreal = ax2.imshow(mat.real)
imimag = ax3.imshow(mat.imag)
plt.colorbar(imabs, ax=ax1)
plt.colorbar(imreal, ax=ax2)
plt.colorbar(imimag, ax=ax3)
ax1.set_title('abs{$S^\dagger S$}')
ax2.set_title('real{$S^\dagger S$}')
ax3.set_title('imag{$S^\dagger S$}')
plt.show()
<Figure size 864x252 with 6 Axes>
It looks pretty close, but there seems to indeed be a bit of loss (expected).
Viewing individual Simulation Data#
To verify, we may want to take a look the individual simulation data. For that, we can load up the batch and inspect the SimulationData for each task.
[16]:
batch_data = modeler.batch.load(path_dir='data')
[17]:
f, (ax1, ax2) = plt.subplots(2, 1, tight_layout=True, figsize=(15, 10))
ax1 = batch_data['smatrix_portleft_top_mode0'].plot_field('field', 'int', freq=freq0, z=wg_height/2, ax=ax1)
ax2 = batch_data['smatrix_portleft_top_mode1'].plot_field('field', 'int', freq=freq0, z=wg_height/2, ax=ax2)
<Figure size 1080x720 with 4 Axes>
[ ]: