Scattering parameters#
What are scattering parameters?#
For a given photonic device, while knowing the actual shape of the device is useful, what we really care about is how light propagates through it. Scattering parameters (S-parameters) are complex numbers that represent the magnitude and phase multiplier acting on the light between every port in a device. S-parameters are collected into an s-matrix which then represents the complete 1st order (linear) input-output response of the device.
For a device with N ports, there will be N2 parameters. The notation for S-Matrices is S(output port)(input port), so S13 will be for the light entering port 3 and exiting port 1. Parameters with the input and output port represent reflections from the device back into the same port, while different output and input ports represent the transmission from the input port into the output ports.
It is easiest to represent s-parameters in polar coordinate form.
Example#
A device has the s-parameter \(S_{21} = 0.98 e^{j\frac{\pi}{2}}\). This means the light accumulates \(90^\circ\) of phase, and its intensity/power will be 96% (\(0.98^2 = 0.9604\)) of the input when we measure the output on port 2 and put input light in port 1.
Why are S-Parameters useful?#
S-parameters allow us to represent a potentially complex photonic component as a matrix, which is much simpler to store and use in computation than other options (e.g. nonlinear functions, FDTD simulations, etc.).
It allows us to connect the device arbitrarily to other components and simulate its behavior in any photonic circuit.
Wavelength dependence (dispersion) can be represented simply by adding an extra dimension to the s-matrix.
S-Parameters are used in circuit-level simulation software, e.g. simphony, sax, Lumerical Interconnect, etc.
How to find S-Parameters for an arbitrary device using a meep simulation#
In this section we’ll demonstrate how to find the s-params for any device using a meep simulation. We’ll show both 2D and 3D simulations, but keep in mind that 2D s-params are not very accurate, but the simulations run much faster. It is often recommended to set up the simulation in 2D to ensure everything works, then extend the simulation into 3D to run the final simulations. This can save us a lot of debugging time since the 2D simulations run faster.
We’ll run through the steps using gdsfactory’s default 2x2 mmi, and explain on the way how the principles can be applied to any component. There are other ways to find s-parameters, one such method is shown on the Gds To Meep page.
Step 1 - Import Geometry#
This tutorial is adapted from the GDSII Import tutorial.
If we have a gds of the device whose parameters we want to simulate, we can add all of our sources and detectors into our gds before we import it into meep. This will create the sources and monitors in meep that we need for the simulation.
Here is a picture of the gds we will be using.
We add a box around the entire component (purple) which will become the simulation region, and 8 lines, 2 on each end that will become the sources and detectors. We also add straight waveguides on the ports of the gdsfactory mmi so that the ports will be exactly where the light would enter and exit the component in an actual circuit.
So, now that we have a gds, lets import it into meep!
# Imports
import meep as mp
import numpy as np
import matplotlib.pyplot as plt
import os
from pathlib import Path
from gplugins.gmeep.get_meep_geometry import get_meep_geometry_from_component
from gdsfactory.read import import_gds
import gdsfactory as gf
mp.verbosity(0)
2024-03-01 12:46:09.088 | INFO | gplugins.gmeep:<module>:39 - Meep '1.28.0' installed at ['/home/andeloth/miniconda3/envs/photonics/lib/python3.11/site-packages/meep']
0
Next we’ll define a whole bunch of variables so our code is more readable down the line.
res = 50 # the resolution of the simulation in pixels/um
sim_is_3D = False # Turn this to false for a 2D simulation
pwd = Path(os.path.abspath(''))
gds_file = pwd.parent / "files/mmi2x2.gds" # The name of our gds file
# The thickness of each material in our simulation (only used in 3D simulations)
t_oxide = 1.0
t_Si = 0.22
t_air = 0.78
dpml = 1 # Diameter of perfectly matched layers
cell_thickness = dpml + t_oxide + t_Si + t_air + dpml # Cell thickness
# Materials used in the simulation
Si = mp.Medium(index=3.45)
SiO2 = mp.Medium(index=1.444)
# Sets the min and max values for the cell and the silicon. Our simulation will be centered at y=0
cell_zmax = 0.5*cell_thickness if sim_is_3D else 0
cell_zmin = -0.5 * cell_thickness if sim_is_3D else 0
# Create a 2D array to hold the S-Parameters for the device
n_ports = 4 # The number of ports, also the size of our array
s_params = np.zeros((n_ports, n_ports), dtype=np.complex128)
mode_parity = mp.NO_PARITY if sim_is_3D else mp.EVEN_Y + mp.ODD_Z
Now, we’ll import our geometry from the gds file.
from gdsfactory.technology import LayerLevel, LayerStack
layers = dict(core=LayerLevel(
layer=(1,0),
thickness=t_Si,
zmin=-t_Si/2,
material="si",
mesh_order=2,
sidewall_angle=0,
width_to_z=0.5,
orientation="100",)
)
layer_stack = LayerStack(layers=layers)
lcen = 1.55
mmi_comp = import_gds(gds_file)
geometry = get_meep_geometry_from_component(mmi_comp, is_3d=sim_is_3D, wavelength=lcen, layer_stack=layer_stack)
# Use this to modify the material of the loaded geometry if needed.
# geometry = [mp.Prism(geom.vertices, geom.height, geom.axis, geom.center, material=mp.Medium(index=3.45)) for geom in geometry]
mmi_comp.plot()
plt.show()
2024-03-01 12:46:09.387 | INFO | gdsfactory.pdk:activate:309 - 'generic' PDK is now active
2024-03-01 12:46:09.559 | WARNING | gdsfactory.config:showwarning:281 - Unnamed cells, 1 in 'Unnamed_0da7224c'
/home/andeloth/miniconda3/envs/photonics/lib/python3.11/site-packages/gdsfactory/component.py:1604: UserWarning: Unnamed cells, 1 in 'Unnamed_0da7224c'
gdspath = component.write_gds(logging=False)
# ###################################################
# Now we actually import the geometry
cell_x = 32
cell_y = 6
cell_z = 3
port_xsize = 0
port_ysize = 1
port_zsize = 0.5
port_xdisp = 13
port_ydisp = (0.75+0.5)/2
port_size = mp.Vector3(port_xsize, port_ysize, port_zsize) if sim_is_3D else mp.Vector3(port_xsize, port_ysize, 0)
cell = mp.Vector3(cell_x, cell_y, cell_z) if sim_is_3D else mp.Vector3(cell_x, cell_y, 0)
port1 = mp.Volume(center=mp.Vector3(-port_xdisp,-port_ydisp,0), size=port_size)
port2 = mp.Volume(center=mp.Vector3(-port_xdisp,port_ydisp,0), size=port_size)
port3 = mp.Volume(center=mp.Vector3(port_xdisp,-port_ydisp,0), size=port_size)
port4 = mp.Volume(center=mp.Vector3(port_xdisp,port_ydisp,0), size=port_size)
source1 = mp.Volume(center=port1.center-mp.Vector3(x=0.5),size=port_size)
source2 = mp.Volume(center=port4.center+mp.Vector3(x=0.5), size=port_size)
Step 2: Run simulation for a single source#
In order to find all the S-Parameters for a device, we have to run a number of simulations where we set the source to be in one of the four ports and then record the output field in all ports. Normally we would have to do this as many times as number of ports, setting a different one to be the input in each time, and we will demonstrate how we can do this in this notebook.
Note, however, that the MMI we are using has mirror symmetry about both the x and y axes, which means we could just do the simulation once and then assign the s-parameters using symmetry arguments. This is a useful trick when working with devices with symmetries about the ports that can save a lot of time.
# Set up frequency points for simulation
npoints = 100
lcen = 1.55
dlam = 0.02
wl = np.linspace(lcen - dlam/2, lcen + dlam/2, npoints)
fcen = 1 / lcen
df_source = 3 * dlam / lcen**2 # ONLY USE FOR SOURCES TO ENSURE LARGE ENOUGH FREQUENCY BANDWIDTH FOR SIMULATION
fpoints = 1 / wl # Center frequency
# Set up the first source for the simulation in port 1 (lower left).
sources = [
mp.EigenModeSource(
src = mp.GaussianSource(fcen, fwidth=df_source),
volume=source1,
eig_band=1,
eig_parity = mode_parity,
eig_match_freq = True,
)
]
# Create Simulation
sim = mp.Simulation(
resolution=res, # The resolution, defined further up
cell_size=cell, # The simulation size, taken from the gds
boundary_layers=[mp.PML(dpml)], # the boundary layers to absorb fields that leave the simulation
sources = sources, # The sources
geometry = geometry, # The geometry
default_material=SiO2
)
After creating the simulation, we can add mode monitors at each of the ports. These will collect the fields on each port and accumulate the fourier transform of the time-domain fields, giving us frequency-domain field information at the end.
After our simulation is completely set up for the first source, we can plot it to make sure everything looks right before running the simulation. We do that using the plot2D() method.
# Adds mode monitors at each of the ports to track the energy that goes in or out
mode_monitor_1 = sim.add_mode_monitor(fpoints, mp.ModeRegion(volume=port1))
mode_monitor_2 = sim.add_mode_monitor(fpoints, mp.ModeRegion(volume=port2))
mode_monitor_3 = sim.add_mode_monitor(fpoints, mp.ModeRegion(volume=port3))
mode_monitor_4 = sim.add_mode_monitor(fpoints, mp.ModeRegion(volume=port4))
# Plot the simulation
plot_plane = mp.Volume(center=mp.Vector3(z=0), size=mp.Vector3(cell.x, cell.y, 0))
sim.plot2D(output_plane=plot_plane if sim_is_3D else None) # No parameters are needed for a 2D simulation.
plt.show()
Show code cell output
As we see from the output of sim.plot2D, our simulation is set up correctly. The red line is our source, the blue are our 4 monitors. We are ready to run the simulation! Actually running the simulation is the most computationally intense part of this, so it may take some time if you have a high resolution and are running a 3D simulation. The until_after_sources=mp.stop_when_dft_decayed parameter for sim.run() means meep will start checking how much the cumulative Fourier Transform of all the fields in the simulation changes only after the sources have turned off. Then, it will wait until the DFT’s change over one time step is below a certain threshold (1e-11 relative change by default), at which point it will stop the simulation
# Runs the simulation
sim.run(until_after_sources=mp.stop_when_dft_decayed(tol=1e-4))
Now that our simulation has been run, the monitors contain the s-parameters which we can retrieve. The function sim.get_eigenmode_coefficients() returns an array that holds the parameters for every mode, every frequency, and in the forward and backward direction. In our simulation we only simulated the fundamental TE mode and at only our center frequency. For most applications, only the fundamental TE mode matters, but often we do want to see how the S-parameters change with frequency. To do that, simply change the df and nfreq parameters for the mode monitors. However, in our simulation we do care about the direction. In our simulation, light will go in port 1 going forward (towards positive X) and exit ports 3 and 4 in the same direction, but any light exiting port 1 and 2 will be going backwards (towards negative X). Therefore, we retrieve the index [0,0,1] for port 1 and 2 which is the index for backward propagating mode.
# Finds the S parameters
norm_mode_coeff = sim.get_eigenmode_coefficients(mode_monitor_1, [1], eig_parity=mode_parity).alpha[0, :, 0]
port1_coeff = sim.get_eigenmode_coefficients(mode_monitor_1, [1], eig_parity=mode_parity).alpha[0, :, 1] / norm_mode_coeff
port2_coeff = sim.get_eigenmode_coefficients(mode_monitor_2, [1], eig_parity=mode_parity).alpha[0, :, 1] / norm_mode_coeff
port3_coeff = sim.get_eigenmode_coefficients(mode_monitor_3, [1], eig_parity=mode_parity).alpha[0, :, 0] / norm_mode_coeff
port4_coeff = sim.get_eigenmode_coefficients(mode_monitor_4, [1], eig_parity=mode_parity).alpha[0, :, 0] / norm_mode_coeff
We’ll also compute and print the transmittance through each port and the total insertion loss. We won’t do this with the other ports as inputs, but it’s a good exercise here.
# Calculates the transmittance based off of the S parameters
port1_trans = abs(port1_coeff) ** 2
port2_trans = abs(port2_coeff) ** 2
port3_trans = abs(port3_coeff) ** 2
port4_trans = abs(port4_coeff) ** 2
# Calculates the Insertion loss as a percent and in dB
insertionLoss = 1-(port2_trans + port3_trans + port4_trans)
insertionLoss_dB = 10*np.log10(insertionLoss)
fig = plt.figure(figsize=(16,10))
ax_trans1 = fig.add_subplot(1,3,1)
ax_trans1.set_title("Transmission per Port")
ax_trans1.plot(wl, port1_trans, label=r's11^2')
ax_trans1.plot(wl, port2_trans, label=r's21^2')
ax_trans1.plot(wl, port3_trans, label=r's31^2')
ax_trans1.plot(wl, port4_trans, label=r's41^2')
ax_trans1.set_xlabel(r"Wavelength (\mu m)")
ax_trans1.set_ylabel(r"Transmission")
ax_trans1.legend()
ax_phase1 = fig.add_subplot(1,3,2)
ax_phase1.set_title("Phase per Port")
ax_phase1.plot(wl, np.unwrap(np.angle(port1_coeff)*180/np.pi), label=r'phase of s11')
ax_phase1.plot(wl, np.unwrap(np.angle(port2_coeff)*180/np.pi), label=r'phase of s21')
ax_phase1.plot(wl, np.unwrap(np.angle(port3_coeff)*180/np.pi), label=r'phase of s31')
ax_phase1.plot(wl, np.unwrap(np.angle(port4_coeff)*180/np.pi), label=r'phase of s41')
ax_phase1.set_xlabel(r"Wavelength (\mu m)")
ax_phase1.set_ylabel("Angle (deg)")
ax_phase1.legend()
ax_loss = fig.add_subplot(1,3,3)
ax_loss.set_title("Loss")
ax_loss.plot(wl, insertionLoss_dB)
ax_loss.set_xlabel(r"Wavelength (\mu m)")
ax_loss.set_ylabel("Loss (dB)")
plt.show()
We find that this component is actually terrible and would never be used in an actual photonic design. Over a third of the light is lost, and the light that isn’t lost is not split in any specific ratio. Fortunately, it was never meant to be used, and just exists as an example of the basic shape of an mmi.
Finally, before we move on to the next step, we’ll run the simulation again and create a plot of the steady state of the mmi.
sim.reset_meep()
whole_dft = sim.add_dft_fields([mp.Ez], fcen, 0, 1, center=mp.Vector3(), size=cell)
sim.run(until_after_sources=mp.stop_when_dft_decayed(tol=1e-4))
eps_data = sim.get_epsilon() # Epsilon Data / The Geometry / An array that holds what materials are where
ez_data = sim.get_dft_array(whole_dft, mp.Ez, 0) # Values for the component of the E-field in the z direction (in/out of screen)
# Creates the plot
plt.figure(dpi=200)
plt.imshow(np.transpose(eps_data), interpolation="spline36", cmap="binary")
plt.imshow(
np.flipud(np.transpose(np.real(ez_data))),
interpolation="spline36",
cmap="RdBu",
alpha=0.9,
)
plt.axis("off")
plt.show()
Step 3: Run the same simulation for each port#
Now, all that’s left to do is to run the simulation for each source. Here, I’ll use a for loop to run the other three simulations. We could have done this for all four, but hopefully this makes the code make more sense.
For the sources going into port 3 and 4, we have to specify that the light will be propagating backwards. Below, I’ve done this by setting the direction to mp.NO_DIRECTION and using the eig_kpoint parameter to specify the direction. You can use this same method to launch sources in any direction.
It should be noted, that since the 2x2 mmi we are using is symmetrical across both the x and y axies, the S-parameters should be the same for all of the other ports. But for the sake of demonstration, we’ll simulate all of them. Then at then end we can check to see the the S-parameters are actually the same.
# Set up the rest of the sources for the simulation.
sources = [
mp.EigenModeSource(
src = mp.GaussianSource(fcen, fwidth=df_source),
size=port1.size, # Here we input the geometry for our first source
center=port1.center + mp.Vector3(-0.5,0,0),
eig_band=1,
eig_parity = mode_parity,
eig_match_freq = True,
),
mp.EigenModeSource(
src = mp.GaussianSource(fcen, fwidth=df_source),
size=port2.size, # Here we input the geometry for our first source
center=port2.center + mp.Vector3(-0.5,0,0),
eig_band=1,
eig_parity = mode_parity,
eig_match_freq = True,
),
mp.EigenModeSource(
src = mp.GaussianSource(fcen, fwidth=df_source),
size=port3.size, # Here we input the geometry for our first source
center=port3.center + mp.Vector3(0.5,0,0),
eig_band=1,
eig_parity = mode_parity,
eig_match_freq = True,
eig_kpoint = mp.Vector3(-1,0,0),
direction = mp.NO_DIRECTION
),
mp.EigenModeSource(
src = mp.GaussianSource(fcen, fwidth=df_source),
size=port4.size, # Here we input the geometry for our first source
center=port4.center + mp.Vector3(0.5,0,0),
eig_band=1,
eig_parity = mode_parity,
eig_match_freq = True,
eig_kpoint = mp.Vector3(-1,0,0),
direction = mp.NO_DIRECTION
)
]
s_params = []
for i in range(n_ports) :
source = sources[i]
# Create Simulation
sim = mp.Simulation(
resolution=res, # The resolution, defined further up
cell_size=cell, # The cell size, taken from the gds
boundary_layers=[mp.PML(dpml)], # the perfectly matched layers, with a diameter as defined above
sources = [source], # The source(s) we just defined
geometry = geometry # The geometry, from above
)
# Adds mode monitors at each of the ports to track the energy that goes in or out
mode_monitor_1 = sim.add_mode_monitor(fpoints, mp.ModeRegion(volume=port1))
mode_monitor_2 = sim.add_mode_monitor(fpoints, mp.ModeRegion(volume=port2))
mode_monitor_3 = sim.add_mode_monitor(fpoints, mp.ModeRegion(volume=port3))
mode_monitor_4 = sim.add_mode_monitor(fpoints, mp.ModeRegion(volume=port4))
# Runs the simulation
sim.run(until_after_sources=mp.stop_when_dft_decayed(tol=1e-4))
#############################################################
# Finds the S parameters
mode_monitors = [mode_monitor_1, mode_monitor_2, mode_monitor_3, mode_monitor_4]
norm_mode = sim.get_eigenmode_coefficients(mode_monitors[i], [1], eig_parity=mode_parity)
norm_mode_coeff = norm_mode.alpha[0,:,0] if i < 2 else norm_mode.alpha[0,:,1]
port1_coeff = sim.get_eigenmode_coefficients(mode_monitor_1, [1], eig_parity=mode_parity).alpha[0, :, 1] / norm_mode_coeff
port2_coeff = sim.get_eigenmode_coefficients(mode_monitor_2, [1], eig_parity=mode_parity).alpha[0, :, 1] / norm_mode_coeff
port3_coeff = sim.get_eigenmode_coefficients(mode_monitor_3, [1], eig_parity=mode_parity).alpha[0, :, 0] / norm_mode_coeff
port4_coeff = sim.get_eigenmode_coefficients(mode_monitor_4, [1], eig_parity=mode_parity).alpha[0, :, 0] / norm_mode_coeff
# Store the S parameters in s_params
index = i
s_params.append([port1_coeff, port2_coeff, port3_coeff, port4_coeff])
s_params = np.array(s_params)
We can now show that due to the symmetry, some of the s-parameters are the same, e.g. s31 and s42.
fig = plt.figure(figsize=(10,10))
ax_trans1 = fig.add_subplot(2,2,1)
ax_trans1.set_title("Transmission for s31")
ax_trans1.plot(wl, np.abs(s_params[0,2])**2, label=r's31^2')
ax_trans1.set_xlabel(r"Wavelength (\mu m)")
ax_trans1.set_ylabel(r"Transmission")
ax_trans1.legend()
ax_phase1 = fig.add_subplot(2,2,2)
ax_phase1.set_title("Phase for s31")
ax_phase1.plot(wl, np.unwrap(np.angle(s_params[0,2])*180/np.pi), label=r'phase of s31')
ax_phase1.set_xlabel(r"Wavelength (\mu m)")
ax_phase1.set_ylabel("Angle (deg)")
ax_phase1.legend()
ax_trans2 = fig.add_subplot(2,2,3)
ax_trans2.set_title("Transmission for s42")
ax_trans2.plot(wl, np.abs(s_params[1,3])**2, label=r's42^2')
ax_trans2.set_xlabel(r"Wavelength (\mu m)")
ax_trans2.set_ylabel(r"Transmission")
ax_trans2.legend()
ax_phase2 = fig.add_subplot(2,2,4)
ax_phase2.set_title("Phase for s42")
ax_phase2.plot(wl, np.unwrap(np.angle(s_params[1,3])*180/np.pi), label=r'phase of s42')
ax_phase2.set_xlabel(r"Wavelength (\mu m)")
ax_phase2.set_ylabel("Angle (deg)")
ax_phase2.legend()
<matplotlib.legend.Legend at 0x7efb5fa43490>
Above, we can see all of the S-parameters from the simulations, and we can verify that the S-parameters are in fact the same for every port, which is what we expected. We also calculated the transmission which is also the same no matter what port we use.
S-Parameters across different Frequencies#
For most applications, you’ll want to know the S-parameters across a range of frequencies. This can be done using the simulations created here. When we get the S-parameters using the sim.get_eigenmode_coefficients().alpha[] we get a matrix. The first entry in the matrix specifies the wavelength. For brevity’s sake, we won’t show that here, but a good exercise would be to find and plot the s-paramters across a range of frequencies for the component shown here.
Multiple frequencies can be simulated at once by setting the second and third arguments in the mode monitors to the desired frequency width and number of points, e.g.
sim.add_mode_monitor(fcen, df, n_freq_points, mp.ModeRegion(volume=port))
Using gplugins.gmeep.write_sparameters_meep#
An alternative method to the one outlined above is using the gplugins package. The write_sparameters_meep function can automate the entire process that we outlined above, but it is less customizable in exchange for its ease of use. This is how we could use it to find the s-parameters for the 2x2 mmi.
from gplugins.gmeep import write_sparameters_meep
pwd = Path(os.path.abspath(''))
gds_file = pwd.parent / "files/mmi2x2.gds"
mmi_comp = import_gds(gds_file)
res = 50
sim_is_3D = False
port_xpos = 17.75
port_ypos = (0.75 + 0.5) / 2
port_ysize = 1
# Enabling this symmetry works for specifically a device with mirror symmetries about the x and y axes and will make the simulation faster
port_symmetries_crossing = {
"o1@0,o1@0": ["o2@0,o2@0", "o3@0,o3@0", "o4@0,o4@0"],
"o2@0,o1@0": ["o1@0,o2@0", "o3@0,o4@0", "o4@0,o3@0"],
"o3@0,o1@0": ["o1@0,o3@0", "o2@0,o4@0", "o4@0,o2@0"],
"o4@0,o1@0": ["o1@0,o4@0", "o2@0,o3@0", "o3@0,o2@0"],
}
mmi_comp.add_port(name="o1", center=[-port_xpos, -port_ypos], width=port_ysize, orientation=0, cross_section=gf.cross_section.cross_section(width=port_ysize, layer='WG'))
mmi_comp.add_port(name="o2", center=[-port_xpos, port_ypos], width=port_ysize, orientation=0, cross_section=gf.cross_section.cross_section(width=port_ysize, layer='WG'))
mmi_comp.add_port(name="o3", center=[port_xpos, -port_ypos], width=port_ysize, orientation=180, cross_section=gf.cross_section.cross_section(width=port_ysize, layer='WG'))
mmi_comp.add_port(name="o4", center=[port_xpos, port_ypos], width=port_ysize, orientation=180, cross_section=gf.cross_section.cross_section(width=port_ysize, layer='WG'))
run_sim = True
sparams = write_sparameters_meep(
component=mmi_comp,
resolution=res,
extend_ports_length=0,
port_symmetries=port_symmetries_crossing,
port_margin=0,
port_source_offset=4.7,
port_monitor_offset=4.5,
run=run_sim,
xmargin=0,
ymargin=1,
zmargin=2 if sim_is_3D else 0,
is_3d=sim_is_3D,
tpml=dpml,
clad_material="SiO2",
wavelength_start=1.54,
wavelength_stop=1.56,
wavelength_points=100,
)
if run_sim:
plt.figure()
plt.plot(sparams['wavelengths'], np.abs(sparams['o3@0,o1@0']))
plt.show()
