The co-planar waveguide (CPW) is a transmission line commonly used in RF photonics due to its compact planar nature, ease of fabrication, and resistance to external EM interference. In this notebook series, we will examine the CPW in the context of a Mach-Zehnder modulator (MZM).
In part one (this notebook), we will start with a simple CPW on a dielectric substrate. The goal is to demonstrate the typical 2D and 3D workflow in Tidy3D's RF solver. We will also benchmark key results against other commercial RF simulation software.
In part two, we will simulate the CPW with the MZM optical waveguide. A 2D mode analysis will first be performed on the conventional layout, followed by a 3D simulation of segmented electrode design.
import matplotlib.pyplot as plt
import numpy as np
import tidy3d as td
import tidy3d.rf as rf
from tidy3d import web
td.config.logging.level = "ERROR"
16:43:38 CST WARNING: Using canonical configuration directory at '/Users/yuanshenli/.config/tidy3d'. Found legacy directory at '~/.tidy3d', which will be ignored. Remove it manually or run 'tidy3d config migrate --delete-legacy' to clean up.
Building the simulation¶
Key parameters¶
We start by defining some key parameters.
# Frequency range
f_min, f_max = (1e9, 65e9)
f0 = (f_max + f_min) / 2
freqs = np.linspace(f_min, f_max, 201)
# Key dimensions (default units: um)
# -- CPW --
G = 5 # CPW gap
WS = 30 # Signal trace width
WG = 350 # Ground trace width
W_cpw = WS + 2 * (G + WG) # Total CPW width
T = 1 # Conductor thickness
L = 2000 # Line length (distance between wave ports)
# -- Substrate --
L_sub, W_sub, H_sub = (L + 200, 2400, 1000) # Substrate length, width, thickness
Material and Geometry¶
The substrate is assumed to be lossless. The metal is assumed to have constant conductivity over the frequency range.
# Material properties
cond = 41 # Conductivity in S/um
eps = 4.5 # Relative permittivity, substrate
# Define EM mediums
med_sub = td.Medium(permittivity=eps)
med_metal = rf.LossyMetalMedium(conductivity=cond, frequency_range=(f_min, f_max))
The CPW geometry is created below. The substrate is assumed to occupy the lower half-space of the simulation domain.
# Create substrate
str_sub_layer = td.Structure(
geometry=td.Box(center=(0, -H_sub / 2, 0), size=(W_sub, H_sub, L_sub)),
medium=med_sub,
)
# Create CPW
geom_sig = td.Box(center=(0, T / 2, 0), size=(WS, T, L))
geom_gnd1 = td.Box.from_bounds(rmin=(-WS / 2 - G - WG, 0, -L / 2), rmax=(-WS / 2 - G, T, L / 2))
geom_gnd2 = geom_gnd1.reflected((1, 0, 0))
geom_cpw = td.GeometryGroup(geometries=[geom_sig, geom_gnd1, geom_gnd2])
str_cpw = td.Structure(geometry=geom_cpw, medium=med_metal)
# Full structure list
structure_list = [str_sub_layer, str_cpw]
Grid and Boundary¶
The simulation domain is surrounded by Perfectly Matched Layers (PMLs) by default. The simulation domain size is defined below.
# Define simulation size
sim_LX = W_sub
sim_LY = 2 * H_sub
sim_LZ = L_sub
The simulation grid is automatically generated based on the minimum wavelength in the respective dielectric medium. In RF problems, the field is especially enhanced near metallic structures. For more accurate results, we use the LayerRefinementSpec feature to automatically refine the grid near metallic edges and corners.
To learn more about the LayerRefinementSpec feature, please refer to the API documentation page or the tutorial page.
# Define Layer refinement spec
lr_spec = rf.LayerRefinementSpec.from_structures(
structures=[str_cpw],
min_steps_along_axis=5,
corner_refinement=td.GridRefinement(dl=T / 5, num_cells=2),
)
# Overall grid specification
grid_spec = td.GridSpec.auto(
min_steps_per_wvl=15,
wavelength=td.C_0 / f_max,
layer_refinement_specs=[lr_spec],
)
Monitors¶
We define a field monitor in the CPW plane for visualization purposes below.
mon_1 = td.FieldMonitor(
center=(0, 0, 0), size=(td.inf, 0, td.inf), freqs=[f0], name="field cpw plane"
)
monitor_list = [mon_1]
Wave Ports¶
Wave ports are used to excite the CPW. It is important to ensure that the width of the ports are slightly smaller than the overall width spanned by the CPW side ground traces.
# Define wave ports
WP1 = rf.WavePort(
center=(0, 0, -L / 2),
size=(W_cpw - 50, W_cpw - 50, 0),
direction="+",
name="WP1",
extrude_structures=True,
)
WP2 = WP1.updated_copy(
center=(0, 0, L / 2),
direction="-",
name="WP2",
)
# List of all ports
port_list = [WP1, WP2]
Define Simulation and TerminalComponentModeler¶
The base Simulation object contains information about the overall simulation environment. The TerminalComponentModeler object is used in 3D simulations to conduct a port and frequency sweep in order to obtain the full S-parameter matrix.
# Define base simulation
sim = td.Simulation(
size=(sim_LX, sim_LY, sim_LZ),
grid_spec=grid_spec,
structures=structure_list,
monitors=monitor_list,
run_time=1e-9,
shutoff=1e-8,
)
# Define TerminalComponentModeler
tcm = rf.TerminalComponentModeler(
simulation=sim,
ports=port_list,
freqs=freqs,
run_only=["WP1@0"],
)
The run_only parameter allows the user to sweep only a subset of the available port excitations. In this case, we only want to run wave port WP1 and mode index 0 (the only mode). This is because the structure is symmetric, so there is no need to solve for port WP2.
Visualization¶
We can inspect the grid to ensure a suitable level of refinement. The plots below show the full CPW (left) and the CPW gap (right). The green region is the wave port and the grid lines are in gray. The purple lines indicate regions of grid refinement.
# Inspect transverse grid
fig, ax = plt.subplots(1, 2, figsize=(10, 6))
tcm.plot_sim(z=-L / 2, ax=ax[0], monitor_alpha=0)
sim.plot_grid(z=-L / 2, ax=ax[0], hlim=(-sim_LX / 2, sim_LX / 2), vlim=(-sim_LY / 2, sim_LY / 2))
sim.plot(z=-L / 2, ax=ax[1], monitor_alpha=0)
sim.plot_grid(z=-L / 2, ax=ax[1], hlim=(10, 30), vlim=(-10, 10))
plt.show()
# Plot top view and wave ports (green/blue)
fig, ax = plt.subplots(figsize=(12, 5))
sim.plot_grid(y=0, ax=ax)
tcm.plot_sim(
y=0.0, ax=ax, monitor_alpha=0, vlim=(-sim_LZ / 2, sim_LZ / 2), hlim=(-W_cpw * 0.6, 0.6 * W_cpw)
)
ax.set_aspect(0.2)
plt.show()
2D Analysis¶
In this section, we will perform a 2D mode analysis on the CPW and extract key transmission line parameters.
Run Mode Solver¶
The 2D mode analysis is performed with a ModeSolver object. We can easily create one from the wave port using the to_mode_solver() method.
# Obtain mode solver from wave port
mode_solver = WP1.to_mode_solver(simulation=sim, freqs=np.linspace(f_min, f_max, 51))
We run the mode solver below.
# Run mode solver
mode_data = web.run(mode_solver, task_name="CPW mode solver", path="data/WP1_mode_data.hdf5")
16:43:40 CST Loading simulation from local cache. View cached task using web UI at 'https://tidy3d.simulation.cloud/workbench?taskId=mo-aa473da2-3fef- 4bec-9efa-ce1455250c9b'.
Mode Profile¶
First, let us examine the mode profile. The left plot depicts $\text{Re}(E_x)$ while the right plot shows the field magnitude in dB. (The divide-by-zero warnings are due to the field magnitude being zero in the metallic regions.)
# Preview mode field
fig, ax = plt.subplots(1, 2, figsize=(10, 4), tight_layout=True)
mode_solver.plot_field(field_name="Ex", val="real", f=f_max, mode_index=0, ax=ax[0])
mode_solver.plot_field(
field_name="E", val="abs^2", scale="dB", f=f_max, mode_index=0, ax=ax[1], vmin=-40, vmax=0
)
for axis in ax:
axis.set_xlim(-200, 200)
axis.set_ylim(-200, 200)
plt.show()
Effective Index, Loss, and Line Impedance¶
The key properties that determine a guided mode are effective index $n_\text{eff}$, attenuation $\alpha$, and characteristic impedance $Z_0$. We demonstrate how to obtain those values below.
We extract the first two quantities, $n_\text{eff}$ and $\alpha$, directly from the mode solver data. Then, we use them to compute the complex propagation constant $\gamma = \alpha + j\beta$.
# Gather alpha, neff from mode solver results
neff_mode = mode_data.modes_info["n eff"].squeeze()
alphadB_mode = mode_data.modes_info["loss (dB/cm)"].squeeze()
# Get alpha, beta
alpha_mode = mode_data.alpha.isel(mode_index=0)
beta_mode = mode_data.beta.isel(mode_index=0)
# Gamma in 1/mm
gamma_mode = (alpha_mode + 1j * beta_mode) * 1e-3
The characteristic line impedance $Z_0$ is retrieved from the mode solver data as well. Note the use of np.conjugate() to convert from the physics phase convention (used in Tidy3D) to the standard electrical engineering phase convention.
# Retrieve the characteristic impedance from the mode impedance.
Z0_mode = np.conjugate(mode_data.transmission_line_data.Z0.isel(mode_index=0)).squeeze()
For comparison purposes, we import benchmark data obtained from other commercial RF simulation software.
# Import benchmark data
freqs_ben1, neff_ben1, alphadB_ben1, Zre_ben1 = np.genfromtxt(
"misc/cpw_fem_mode.csv", delimiter=",", unpack=True
)
freqs_ben2, neff_ben2, alphadB_ben2, Zre_ben2 = np.genfromtxt(
"misc/cpw_fit_mode.csv", delimiter=",", unpack=True
)
The three quantities are plotted below. The results compare very well with the benchmark data.
def comparison_plot(data, data_ben1, data_ben2, err_band, ax):
"""Generate comparison plot of three datasets with fractional error band based on dataset"""
ax.plot(data.f / 1e9, data, label="RF Solver")
ax.plot(freqs_ben1, data_ben1, ls="--", color="#555555", label="Commercial FEM")
ax.plot(freqs_ben2, data_ben2, ls=":", color="#555555", label="Commercial FIT")
ax.fill_between(
data.f / 1e9,
(1 + err_band) * data,
(1 - err_band) * data,
alpha=0.2,
color="gray",
label=f"$\\pm${err_band * 100:.2f}%",
)
# Comparison plots
fig, ax = plt.subplots(3, 1, figsize=(10, 10), tight_layout=True)
ax[0].set_title("Effective index")
comparison_plot(neff_mode, neff_ben1, neff_ben2, 0.005, ax=ax[0])
ax[1].set_title("Attenuation (dB/cm)")
comparison_plot(alphadB_mode, alphadB_ben1, alphadB_ben2, 0.05, ax=ax[1])
ax[2].set_title("Real impedance (Ohm)")
comparison_plot(np.real(Z0_mode), Zre_ben1, Zre_ben2, 0.02, ax=ax[2])
for axis in ax:
axis.legend()
axis.grid()
axis.set_xlabel("f (GHz)")
plt.show()
Distributed RLCG transmission line parameters¶
Using the complex propagation constant $\gamma$ and the line impedance $Z_0$, we can calculate the RLCG parameters using the following relationships: $$ R = \text{Re}(\gamma Z_0), \quad L = \text{Im}(\gamma Z_0)/\omega, \quad G = \text{Re}(\gamma/Z_0), \quad C = \text{Im}(\gamma/Z_0)/\omega $$ where $\omega = 2\pi f$.
R_mode = np.real(gamma_mode * Z0_mode) # Ohm/mm
G_mode = np.real(gamma_mode / Z0_mode) # S/mm
L_mode = np.imag(gamma_mode * Z0_mode) / (2 * np.pi * gamma_mode.f) # H/mm
C_mode = np.imag(gamma_mode / Z0_mode) / (2 * np.pi * gamma_mode.f) # F/mm
The RLCG parameters are shown below. Note that since the dielectric is lossless, $G$ is actually zero and the plotted values are due to small amounts of numerical noise.
# Plot RLCG
fig, ax = plt.subplots(2, 2, figsize=(10, 8), tight_layout=False)
ax[0, 0].plot(R_mode.f / 1e9, R_mode)
ax[0, 1].plot(L_mode.f / 1e9, L_mode * 1e9)
ax[1, 0].plot(G_mode.f / 1e9, G_mode)
ax[1, 1].plot(C_mode.f / 1e9, C_mode * 1e12)
titles = ["R (Ohm/mm)", "L (nH/mm)", "G (S/mm)", "C (pF/mm)"]
for ii, title in enumerate(titles):
ax[ii // 2, ii % 2].set_title(title)
ax[ii // 2, ii % 2].set_xlabel("f (GHz)")
ax[ii // 2, ii % 2].grid()
plt.show()
3D Analysis¶
Run Simulation¶
Let us simulate a short section of the CPW to demonstrate the 3D workflow. First, the grid resolution can be reduced from the 2D case. This is because S-parameters generally converge much more quickly than modal parameters. Below, we redefine the grid specification for the 3D case.
# Layer refinement spec
lr_spec_3d = rf.LayerRefinementSpec.from_structures(
structures=[str_cpw],
min_steps_along_axis=2,
corner_refinement=td.GridRefinement(dl=G / 3, num_cells=2),
)
# Overall grid specification
grid_spec_3d = td.GridSpec.auto(
min_steps_per_wvl=15,
wavelength=td.C_0 / f_max,
layer_refinement_specs=[lr_spec_3d],
)
We make an updated copy of the base Simulation and TerminalComponentModeler using the reduced grid specification.
sim_3d = sim.updated_copy(grid_spec=grid_spec_3d)
tcm_3d = tcm.updated_copy(simulation=sim_3d)
The simulation is executed below.
tcm_data = web.run(tcm_3d, task_name="CPW 3D", path="data/tcm_data_cpw.hdf5")
16:43:41 CST Loading simulation from local cache. View cached task using web UI at 'https://tidy3d.simulation.cloud/rf?taskId=sid-b2c5a113-caf3-4597-8 6e4-38e476b21828'.
S-parameters¶
The S-matrix can be calculated using the smatrix() method from the TerminalComponentModelerData instance tcm_data returned by web.run().
smat = tcm_data.smatrix()
Below, we define some convenience functions to extract the individual $S_{ij}$ entries. The port_in and port_out coordinates are used to specify the port name/index. Note the use of np.conjugate() to convert the S-parameter from the physics phase convention (Tidy3D default) to the standard electrical engineering convention.
def sparam(i, j):
return np.conjugate(smat.data.isel(port_in=j - 1, port_out=i - 1))
def sparam_abs(i, j):
return np.abs(sparam(i, j))
def sparam_angle(i, j):
return np.unwrap(np.angle(sparam(i, j))) / np.pi * 180
def sparam_dB(i, j):
return 20 * np.log10(sparam_abs(i, j))
We import some benchmark data for comparison below.
# Import benchmark data
freqs_ben1, S11abs_ben1, S11arg_ben1, S21abs_ben1, S21arg_ben1 = np.genfromtxt(
"misc/cpw_fem_sparam.csv", delimiter=",", unpack=True
)
freqs_ben2, S11abs_ben2, S11arg_ben2, S21abs_ben2, S21arg_ben2 = np.genfromtxt(
"misc/cpw_fit_sparam.csv", delimiter=",", unpack=True
)
The insertion loss magnitude and angle are compared below. They agree very well.
def comparison_plot_sparam(data, data_ben1, data_ben2, err_band, ax, err_band_on=True):
"""Comparison plot of three datasets with fractional error band based on data1"""
ax.plot(freqs / 1e9, data, label="RF Solver")
ax.plot(freqs_ben1, data_ben1, ls="--", color="#555555", label="Commercial FEM")
ax.plot(freqs_ben2, data_ben2, ls=":", color="#555555", label="Commercial FIT")
if err_band_on:
ax.fill_between(
freqs / 1e9,
data + err_band,
data - err_band,
alpha=0.2,
color="gray",
label=f"$\\pm${err_band:.2f} delta S",
)
# Compare insertion loss
fig, ax = plt.subplots(2, 1, figsize=(10, 8), tight_layout=True)
ax[0].set_title("Insertion loss magnitude")
comparison_plot_sparam(sparam_abs(2, 1), S21abs_ben1, S21abs_ben2, 0.01, ax=ax[0])
ax[1].set_title("Insertion loss angle (degrees)")
comparison_plot_sparam(
sparam_angle(2, 1),
np.unwrap(S21arg_ben1),
np.unwrap(S21arg_ben2),
0.0,
ax=ax[1],
err_band_on=False,
)
for axis in ax:
axis.legend()
axis.grid()
axis.set_xlabel("f (GHz)")
plt.show()
Since this is a straight transmission line with no longitudinal variation, the return loss should theoretically be zero. Below, we plot the return loss magnitude in dB. The RF solver is able to achieve a very low level of numerical reflection. The propagating mode is very well matched to the modal absorbing boundary present in each wave port.
fig, ax = plt.subplots(figsize=(10, 4), tight_layout=True)
ax.set_title("Return loss (dB)")
ax.plot(freqs / 1e9, sparam_dB(1, 1))
ax.grid()
ax.set_xlabel("f (GHz)")
plt.show()
Field profile¶
To access the field monitor data, we need to first load the simulation dataset from the TerminalComponentModelerData instance. The simulation data is stored in a dictionary indexed by keys in the format <port name>@<mode number>. For example, to obtain data for the first (and only) excited mode in wave port 1, we use WP1@0.
# Load monitor data
sim_data = tcm_data.data["WP1@0"]
Below, we plot the field magnitude in dB just below the CPW metal plane. (The divide-by-zero error is due to zero field magnitude in some metallic regions.)
fig, ax = plt.subplots(figsize=(10, 8), tight_layout=True)
sim_data.plot_field(
"field cpw plane", field_name="E", val="abs", scale="dB", f=f0, ax=ax, vmin=-60, vmax=10
)
ax.set_xlim(-W_cpw / 2, W_cpw / 2)
ax.set_ylim(-L / 2, L / 2)
plt.show()
Conclusion¶
In this notebook, we used the basic CPW as an example to demonstrate the typical Tidy3D RF workflow for modeling transmission lines. We calculated and benchmarked key quantities such as effective index, attenuation, characteristic impedance, and S-parameters.
In the next notebook, we will model the CPW in the context of a Mach-Zehnder modulator.
TERMS & CONDITIONS
How the process works:
Try this real URL to see an example: