TIDY3D
LEARNING CENTER

Lectures
Tutorials
Applications

Example Library / Passive Photonic Integrated Circuit Components / Star couplers and arrayed waveguide gratings

Star couplers and arrayed waveguide gratings

Sign up now to simulate .ipynb

Reference: S. Hulyal et al., “Arrayed waveguide gratings in lithium tantalate integrated photonics,” Optica 12, 978-984 (2025) DOI:10.1364/OPTICA.565570.

A star coupler is a free-propagation slab that connects waveguides arranged along a focal arc. Light from the input waveguide diffracts across the slab, and each array arm samples a different amplitude and phase. The output star coupler performs the reciprocal process, allowing the fields from all arms to expand, interfere, and focus onto the output waveguides.

An arrayed waveguide grating (AWG) adds a fixed path-length increment \(\Delta L_\mathrm{AWG}\) between adjacent arms. At frequency \(f\), this produces the propagation phase

\[\phi(f)=\frac{2\pi f}{c}n_\mathrm{eff}(f)\Delta L_\mathrm{AWG},\]

where \(n_\mathrm{eff}\) is the effective index of waveguides. The array therefore acts as a phased aperture whose focal position depends on frequency. Output waveguides placed along the focal arc sample different spectral channels, enabling wavelength demultiplexing or on-chip spectroscopy.

This notebook combines full-wave simulations of the two star couplers, especially the Rowland type from the paper, with analytical propagation through the array as follows:

  1. Simulation setup and mode analysis: calculate the fundamental mode and its dispersion, and derive the geometry parameters for Rowland couplers.
  2. Input coupler: simulate diffraction into the array arms.
  3. Arrayed waveguide gratings: apply arm-dependent propagation analytically.
  4. Output coupler: simulate wavelength-dependent recombination at the output ports.

All lengths are in μm and all frequencies are in Hz unless stated otherwise. This notebook costs ~ 65 FlexCredits to run the input/output coupler simulations.

schematic

Simulation Setup

The Tidy3D simulation and cloud job-submission APIs, together with the post-processing libraries used throughout the notebook, are imported below.

# Core numerics and plotting
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np

# Tidy3D simulation and cloud job submission
import tidy3d as td
from tidy3d import web

td.config.logging.level = "ERROR"

The design parameters of the device are specified first; the operating band and material models used throughout the notebook are then derived from them.

# --- Simulation parameters (design values for the target device) ---

# Operating band
freq0 = 198.4e12  # center frequency (Hz), 198.4 THz
f_spacing = 100e9  # channel spacing (Hz), 100 GHz

# Waveguide cross-section geometry (um)
t_sub = 4.7  # SiO2 substrate thickness
t_cladding = 0.8  # SiO2 cladding thickness
t_LiTaO3 = 0.3  # LiTaO3 core thickness
w_wg = 1.8  # channel waveguide width
l_wg = 50  # straight waveguide length used in the mode simulation
sidewall_angle_deg = 20  # ridge sidewall angle (deg)
box_factor = (
    1.2,
    3,
)  # (width, height) multiplier on (w_wg, t_LiTaO3) for confinement/monitor boxes

# Materials
sio2_dataset = "Palik_LowLoss"
n_LiTaO3_o = 2.119
n_LiTaO3_e = 2.123


# AWG / Rowland star-coupler design
m_order = 160  # AWG diffraction order
Ra = 75.0  # grating-arc radius (um)
da = 2.3  # array-aperture pitch (um)
n_grating_arms = 17  # number of array arms
n_output_channels = 9  # number of output waveguides

# Taper / access-waveguide geometry (um)
w_taper = 2.0  # taper facet width
L_taper = 5.0  # taper length
wg_extend = 100  # straight waveguide extension length

# Center wavelength and frequency grid spanning the operating band
lda0 = td.C_0 / freq0

fwidth = f_spacing * 11
fmin, fmax = freq0 - fwidth / 2, freq0 + fwidth / 2

freqs = np.linspace(fmin, fmax, 111)
freqs_reduced = np.linspace(fmin, fmax, 23)
source_time = td.GaussianPulse(freq0=freq0, fwidth=fwidth)

lda_max, lda_min = td.C_0 / fmin, td.C_0 / fmax

print(f"Target wavelength = {lda0:.2f} um")
print(f"Target center frequency = {freq0 / 1e12:.2f} THz")
Target wavelength = 1.51 um
Target center frequency = 198.40 THz

SiO2 for cladding is imported from the td.material_library. LiTaO3 for waveguide/slab core is defined using td.AnisotropicMedium with refractive indices in ordinary (\(x, y\)) and extraordinary (\(z\)) directions.

SiO2 = td.material_library["SiO2"][sio2_dataset]

LiTaO3 = td.AnisotropicMedium(
    xx=td.Medium.from_nk(n=n_LiTaO3_o, k=0, freq=freq0),
    yy=td.Medium.from_nk(n=n_LiTaO3_o, k=0, freq=freq0),
    zz=td.Medium.from_nk(n=n_LiTaO3_e, k=0, freq=freq0),
)

Waveguide Geometry

The ridge waveguide consists of a 1.8 μm-wide, 300 nm-thick LiTaO3 core with 20° sidewalls, sandwiched between SiO2 substrate and cladding. A transverse ModeSimulation solves for the fundamental mode using 20 grid points per shortest wavelength and PML boundaries, while a companion bulk-slab ModeSimulation with periodic in-plane boundaries yields the free-propagation-region (FPR) mode used later in the star-coupler design. The resulting cross-section plot serves to verify the waveguide geometry.

# Ridge waveguide cross-section: LiTaO3 slab/core between SiO2 substrate and cladding
sidewall = sidewall_angle_deg * np.pi / 180
mode_size = (w_wg * 4, 0, t_LiTaO3 * 12)

slab = td.Structure(
    geometry=td.Box(center=(0, 0, -t_LiTaO3 / 2), size=(td.inf, td.inf, t_LiTaO3)),
    medium=LiTaO3,
    name="slab",
)
substrate = td.Structure(
    geometry=td.Box(center=(0, 0, -t_LiTaO3 - t_sub / 2), size=(td.inf, td.inf, t_sub)),
    medium=SiO2,
    name="substrate",
)
cladding = td.Structure(
    geometry=td.Box(
        center=(0, 0, (t_cladding + t_LiTaO3) / 2),
        size=(td.inf, td.inf, t_cladding + t_LiTaO3),
    ),
    medium=SiO2,
    name="cladding",
)
base_layers = [slab, substrate, cladding]

# Straight channel waveguide used for the transverse mode solve
wg = td.Structure(
    geometry=td.PolySlab(
        vertices=[
            (w_wg / 2, -l_wg / 2),
            (-w_wg / 2, -l_wg / 2),
            (-w_wg / 2, l_wg / 2),
            (w_wg / 2, l_wg / 2),
        ],
        sidewall_angle=sidewall,
        slab_bounds=(0, t_LiTaO3),
        reference_plane="bottom",
    ),
    medium=LiTaO3,
)

grid_spec = td.GridSpec.auto(min_steps_per_wvl=20, wavelength=lda_min)

# ModeSimulation for the channel waveguide (PML boundaries)
sim_mode_wg = td.ModeSimulation(
    size=mode_size,
    structures=[*base_layers, wg],
    grid_spec=grid_spec,
    boundary_spec=td.BoundarySpec.all_sides(td.PML()),
    freqs=freqs,
    mode_spec=td.ModeSpec(num_modes=1, target_neff=n_LiTaO3_o),
    symmetry=(-1, 0, 0),
)

# Bulk-slab ModeSimulation (periodic in-plane) for the FPR mode
sim_mode_fpr = sim_mode_wg.updated_copy(
    size=(0, 0, t_LiTaO3 * 12),
    boundary_spec=td.BoundarySpec(
        x=td.Boundary(minus=td.Periodic(), plus=td.Periodic()),
        y=td.Boundary(minus=td.Periodic(), plus=td.Periodic()),
        z=td.Boundary(minus=td.PML(), plus=td.PML()),
    ),
    symmetry=(0, 0, 0),
    freqs=[freq0],
)
# Plot the cross-section to verify the geometry
fig, ax = plt.subplots()
sim_mode_wg.plot(y=0, ax=ax)
ax.text(0, 0, "LiTaO$_3$", horizontalalignment="left", verticalalignment="center")
ax.text(
    0, -t_LiTaO3 * 3, "SiO$_2$", horizontalalignment="left", verticalalignment="center"
)
ax.text(
    0,
    t_LiTaO3 * 2 + t_cladding,
    "Air",
    horizontalalignment="left",
    verticalalignment="center",
);

Run the Mode Simulations

Both transverse eigenmode problems are solved, and the resulting fields and propagation constants are stored in mode_data.

# Run both mode simulations (channel waveguide and FPR slab) concurrently in the cloud
mode_data = web.run_async(
    simulations={"wg": sim_mode_wg, "fpr": sim_mode_fpr},
    folder_name="StarCoupler",
    verbose=False,
    path_dir="output/mode",
)

Mode Analysis

\(\mathrm{Re}(H_z)\) and \(|S_y|\) are plotted at the center frequency to verify the mode polarization and confinement, and the power flux is integrated within two cross-sectional boxes to quantify the fraction confined near the core. The frequency-dependent effective index yields the group index

\[n_g = n_\mathrm{eff} + f\,\frac{\mathrm{d} n_\mathrm{eff}}{\mathrm{d}f},\]

which characterizes propagation through the array arms, while the bulk-slab simulation’s \(|H_z|\) field profile (plotted as an inset) and free-propagation-region index \(n_\mathrm{FPR}\) are used in the star-coupler design below.

# Plot mode field (Re Hz) and power flux (|Sy|) at the center frequency
fig, ax = plt.subplots(2, 1, sharex=True, sharey=True, tight_layout=True)
mode_data["wg"].plot_field(
    field_name="Hz",
    val="real",
    f=freq0,
    ax=ax[0],
    cmap=plt.cm.RdBu,
    vmin=-0.17,
    vmax=0.17,
)
mode_data["wg"].plot_field(
    field_name="Sy", val="abs", f=freq0, ax=ax[1], cmap=plt.cm.inferno, vmin=0, vmax=2
)
ax[0].set_title(r"Field, Re($H_z$)")
ax[1].set_title(r"Power flux, $S_y$")


# Time-averaged Sy(x, z, f) for the transverse mode, and the fraction of its integrated
# power at a given frequency contained within |x|<=x_half, |z|<=z_half
def mode_poynting(mode_data_wg):
    Ex, Ez = mode_data_wg.modes.Ex.squeeze(), mode_data_wg.modes.Ez.squeeze()
    Hx, Hz = mode_data_wg.modes.Hx.squeeze(), mode_data_wg.modes.Hz.squeeze()
    return 0.5 * np.real(Ez * np.conj(Hx) - Ex * np.conj(Hz))


def confined_power_fraction(Sy, freq, x_half, z_half):
    x, z = Sy.x, Sy.z
    P_tot = Sy.sel(f=freq).integrate(coord=["x", "z"]).data
    P_box = (
        Sy.sel(f=freq, x=np.abs(x) <= x_half, z=np.abs(z) <= z_half)
        .integrate(coord=["x", "z"])
        .data
    )
    return P_box / P_tot


Sy = mode_poynting(mode_data["wg"])
P_box = confined_power_fraction(
    Sy, freq0, w_wg * box_factor[0] / 2, t_LiTaO3 * box_factor[1]
)
P_box_ref = confined_power_fraction(Sy, freq0, w_wg / 2, t_LiTaO3)

# Overlay the two confinement boxes and their captured-power fractions on the Sy panel
ax[1].plot(
    w_wg / 2 * np.array([1, 1, -1, -1, 1]) * box_factor[0],
    t_LiTaO3 * np.array([1, -1, -1, 1, 1]) * box_factor[1],
    "y--",
    lw=0.75,
)
ax[1].text(
    w_wg / 2 * 1.4, 3 * t_LiTaO3, f"P = {P_box * 100:.1f}%", color="y", fontsize=10
)
ax[1].plot(
    w_wg / 2 * np.array([1, 1, -1, -1, 1]),
    t_LiTaO3 * np.array([1, -1, -1, 1, 1]),
    "w--",
    lw=0.75,
)
ax[1].text(
    w_wg / 2 * 1.4, 1 * t_LiTaO3, f"P = {P_box_ref * 100:.1f}%", color="w", fontsize=10
)

# Inset: FPR (bulk-slab) mode field profile
ax2 = ax[0].inset_axes([1.5, 0, 0.5, 1])
Hz_fpr = mode_data["fpr"].modes.Hz.squeeze()
ax2.plot(np.abs(Hz_fpr), Hz_fpr.z)
ax2.set(
    ylim=(Hz_fpr.z.min(), Hz_fpr.z.max()),
    xlim=(-0.01, 0.15),
    xlabel=r"$|H_z|$",
    yticklabels=[],
)
ax2.fill_between((-0.01, 0.15), t_LiTaO3, -t_LiTaO3, color="gray", alpha=0.2)
ax2.set_title("Field (FPR)");

# Effective and group index at the center frequency, and the FPR effective index
n_eff = mode_data["wg"].modes.n_eff.squeeze()
n_group = n_eff + freqs * np.gradient(n_eff, freqs)
n_eff0 = n_eff.sel(f=freq0, method="nearest").squeeze().data
n_group0 = n_group.sel(f=freq0, method="nearest").squeeze().data
n_FPR = mode_data["fpr"].modes.n_eff.squeeze().data

print(
    f"Effective and group indices at center frequency = {n_eff0:.03f}, {n_group0:.03f}."
)
print(f"Effective free-propagation-region index = {n_FPR:.03f}.")

# Plot the effective index vs. frequency
fig, ax = plt.subplots()
ax.plot(freqs / 1e12, n_eff)
ax.set(
    xlabel="Frequency (THz)",
    ylabel=r"Effective Index, $n_\mathrm{eff}$",
    xlim=(fmin / 1e12, fmax / 1e12),
);
Effective and group indices at center frequency = 1.922, 2.192.
Effective free-propagation-region index = 1.955.

Rowland Geometry

The array apertures lie on a grating arc of radius \(R_a=75~\mu\mathrm{m}\); the access and output waveguides lie on the Rowland focal arc of radius \(R_a/2\).

For points \(\mathbf{P} = R_a(-\sin\alpha, \cos\alpha-1/2)\) and \(\mathbf{Q}=R_a/2 (\sin 2\theta, -\cos2\theta)\) on the grating and Rowland arcs, the distance is \(L(\alpha, \theta) = R_a\left[ 1 + \sin^2\theta + 2\sin\theta \sin(\alpha-\theta)\right]^{1/2}\), which expands near \(\alpha=0\) as \[L(\alpha,\theta) = R_a\left[ \cos\theta + \sin\theta\alpha - \frac{\sin\theta}{6}\alpha^3 + O(\alpha^4) \right] \approx R_a\cos\theta + R_a\alpha\sin\theta.\]

The constant term sets the reference propagation distance, while the linear term provides the phase ramp across the array required for focusing at angle (\(\theta\)) and is included in the array phase-matching condition. Crucially, the quadratic term vanishes identically. The Rowland focal geometry therefore eliminates second-order defocus about the central array aperture.

rowland geometry

Focusing Condition

For angular aperture spacing \(\Delta\alpha\), the arc pitch is \(d_a=R_a\Delta\alpha\). Using \(\Delta L_{\mathrm{FPR}}\approx d_a\sin\theta_f\), the constructive-interference condition \(\delta\phi + n_\mathrm{FPR} k \Delta L_\mathrm{FPR} = 0\) gives

\[\sin\theta_f = -\frac{\delta\phi} {n_{\mathrm{FPR}}k d_a},\]

where \(k = 2\pi f/c\) is the free-space wave number.

AWG Order and Free Spectral Range

Let \(f_0\) be the frequency focused at the center output waveguide (\(\theta_f=0\)), so the phase difference between adjacent array arms is an integer multiple of \(2\pi\):

\[\delta\phi_0 = \frac{2\pi f_0}{c} n_{\mathrm{eff},0}\Delta L_{\mathrm{AWG}} = 2\pi m,\]

giving

\[\boxed{ \Delta L_{\mathrm{AWG}} = \frac{mc}{n_{\mathrm{eff},0}f_0} }.\]

With group index \(n_g = n_{\mathrm{eff}} + f \cdot \mathrm{d}n_{\mathrm{eff}}/\mathrm{d}f\), the phase near \(f_0\) varies as

\[\delta\phi(f)-\delta\phi_0 \approx \frac{2\pi n_g\Delta L_{\mathrm{AWG}}}{c}(f-f_0).\]

A \(2\pi\) phase change defines the free spectral range:

\[\boxed{ f_{\mathrm{FSR}} = \frac{c}{n_g\Delta L_{\mathrm{AWG}}} = \frac{n_{\mathrm{eff},0}}{mn_g}f_0 }.\]

Frequency and Output-Waveguide Spacing

For channel frequencies \(f_l=f_0+l\Delta f\), the residual phase is

\[\delta\phi_l \approx 2\pi m \frac{n_g}{n_{\mathrm{eff},0}} \frac{l\Delta f}{f_0},\]

giving a focal angle

\[\boxed{ \sin\theta_l \approx - \frac{m\lambda_l}{n_{\mathrm{FPR}}d_a} \frac{n_g}{n_{\mathrm{eff},0}} \frac{l\Delta f}{f_0} }.\]

The sign of \(\theta_l\) depends on the array-index and angular-coordinate conventions.

# AWG arm-length increment, adjacent-arm phase increment, and free spectral range
Delta_L = m_order * lda0 / n_eff0
Delta_phi = 2 * np.pi * n_group0 * Delta_L * f_spacing / td.C_0
f_FSR = n_eff0 * freq0 / m_order / n_group0

# Angular positions of the array apertures along the grating arc
Delta_alpha = da / Ra

grating_idx = np.arange(-(n_grating_arms // 2), n_grating_arms // 2 + 1)
alpha_grt = Delta_alpha * grating_idx

# Focal angle of each output waveguide, from the Rowland focusing condition
output_idx = np.arange(-(n_output_channels // 2), n_output_channels // 2 + 1)
sin_theta_out = (
    -(m_order / n_FPR / da)
    * (td.C_0 / (freq0 + f_spacing * output_idx))
    * (n_group0 / n_eff0)
    * (output_idx * f_spacing / freq0)
)
theta_out = np.arcsin(sin_theta_out)

print(f"m order = {m_order}")
print(f"Delta L = {Delta_L:.1f} um")
print(f"Delta phi = {Delta_phi * 180 / np.pi:.2f} deg")
print(f"Free spectral range = {f_FSR / 1e12:.2f} THz")
print(f"AWG spacing, d_a = {da:.2f} um")
print(f"Grating angle, Delta_alpha = {Delta_alpha * 180 / np.pi:.2f} deg")
print(f"Divergence angle, theta_out = {theta_out[5] * 180 / np.pi:.2f} deg")
m order = 160
Delta L = 125.8 um
Delta phi = 33.11 deg
Free spectral range = 1.09 THz
AWG spacing, d_a = 2.30 um
Grating angle, Delta_alpha = 1.76 deg
Divergence angle, theta_out = -1.77 deg

The physical geometry is constructed from these design values: the Rowland slab bounded by the grating and focal arcs, a tapered waveguide template, and rotated and translated copies of it via wg_transform for the input arm, the 17 array arms, and the 9 output waveguides. The resulting input- and output-coupler layouts are plotted in the two panels below for visual inspection.

# Rowland structure: grating-arc + Rowland-focal-arc polygon (the slab region)
theta_grt = np.linspace(60, 120, 121) * np.pi / 180
vert_grt_x = Ra * np.cos(theta_grt)
vert_grt_y = -Ra / 2 + Ra * np.sin(theta_grt)
vert_grt = np.column_stack([vert_grt_x, vert_grt_y])

theta_rowland = np.linspace(180, 360, 181) * np.pi / 180
vert_rowland_x = Ra / 2 * np.cos(theta_rowland)
vert_rowland_y = Ra / 2 * np.sin(theta_rowland)
vert_rowland = np.column_stack([vert_rowland_x, vert_rowland_y])

vert = np.concat([vert_grt, vert_rowland], axis=0)

Rowland = td.Structure(
    geometry=td.PolySlab(
        vertices=vert,
        sidewall_angle=sidewall,
        slab_bounds=(0, t_LiTaO3),
        reference_plane="bottom",
    ),
    medium=LiTaO3,
    name="Rowland",
)

# Tapered waveguide template: narrow taper at the facet widening into the straight channel
y = np.linspace(0, L_taper, 21)
w = (-w_taper + w_wg) / L_taper**2 * y**2 + w_taper
vert_x = np.concat(
    [
        w / 2,
        np.array([w_wg / 2, -w_wg / 2]),
        -w[::-1] / 2,
        np.array([-w_taper / 2, w_taper / 2]),
    ]
)
vert_y = np.concat([y, np.array([wg_extend, wg_extend]), y[::-1], -np.array([5, 5])])
vert = np.column_stack([vert_x, vert_y])

wg_unit = td.Structure(
    geometry=td.PolySlab(
        vertices=vert,
        sidewall_angle=sidewall,
        slab_bounds=(0, t_LiTaO3),
        reference_plane="bottom",
    ),
    medium=LiTaO3,
)


# Rotate/translate a copy of the waveguide template onto the arc at the given angle
def wg_transform(center=(0, -Ra / 2), radius=Ra, angle=0, name="wg"):
    geom = wg_unit.geometry.rotated(axis=2, angle=angle)
    geom = geom.translated(
        x=center[0] - radius * np.sin(angle), y=center[1] + radius * np.cos(angle), z=0
    )
    return wg_unit.updated_copy(geometry=geom, name=name)


# Instantiate the input arm, the 17 array arms, and the 9 output arms
wg_in = wg_transform(center=(0, 0), radius=Ra / 2, angle=np.pi, name="wg_in")
wg_gratings = [
    wg_transform(angle=a, name=f"wg_grating_{n}") for n, a in enumerate(alpha_grt)
]
wg_out = [
    wg_transform(
        center=(0, Ra / 2),
        radius=Ra * np.cos(theta),
        angle=np.pi + theta,
        name=f"wg_out_{n}",
    )
    for n, theta in enumerate(theta_out)
]
# Plot the input- and output-coupler layouts for visual inspection
scene_in = td.Scene(structures=[Rowland, wg_in, *wg_gratings])
scene_out = td.Scene(structures=[Rowland, *wg_out, *wg_gratings])

fig, ax = plt.subplots(
    1, 2, figsize=(8, 8), sharey=True, sharex=True, tight_layout=True
)
scene_in.plot(z=0, ax=ax[0])
scene_out.plot(z=0, ax=ax[1])
ax[0].set(xlim=(-Ra, Ra))
ax[0].set_title("Input coupler")
ax[1].set_title("Output coupler")

for n, alpha in enumerate(alpha_grt[::8]):
    xx, yy = Ra * np.cos(alpha + np.pi / 2), -Ra / 2 + Ra * np.sin(alpha + np.pi / 2)
    ax[0].plot([0, xx], [-Ra / 2, yy], "r--")
    if n == 0:
        ax[0].text(xx + 3, yy + 3, f"grating arm {n * 8}", color="r")
        ax[1].text(xx + 3, yy + 3, "leading pulse", color="r")
    elif n == 2:
        ax[0].text(
            xx - 3,
            yy + 3,
            f"grating arm {n * 8}",
            color="r",
            horizontalalignment="right",
        )
        ax[1].text(
            xx - 3,
            yy + 3,
            "delayed pulse",
            color="r",
            horizontalalignment="right",
        )

for n, theta in enumerate(theta_out[::4]):
    xx, yy = (
        Ra * np.cos(theta) * np.cos(theta + np.pi / 2 * 3),
        +Ra / 2 + Ra * np.cos(theta) * np.sin(theta + np.pi / 2 * 3),
    )
    ax[1].plot([0, xx], [+Ra / 2, yy], "m--")
    if n == 0:
        ax[1].text(xx + 3, yy - 3, "output (low freq.)", color="m")
    elif n == 2:
        ax[1].text(
            xx - 3,
            yy - 3,
            "output (high freq.)",
            color="m",
            horizontalalignment="right",
        )

Input Coupler Simulation

The fundamental mode is launched into the input star coupler. A spatially subsampled planar monitor records the slab field, a small full-resolution monitor covers the array facets for the close-up plots, and transverse monitors capture the complex fields entering each of the 17 array arms.

In-plane absorbers and vertical PML suppress boundary reflections, and mirror symmetry about \(x=0\) is imposed to reduce the computational cost.

# Mode source launching the fundamental mode into the input star coupler
r_mnt = Ra + 5.0 + lda_max * 1.5
source = td.ModeSource(
    center=(0, Ra / 2 - r_mnt, 0),
    size=mode_size,
    source_time=source_time,
    direction="+",
    mode_spec=td.ModeSpec(num_modes=1, target_neff=n_LiTaO3_o),
)

# Planar field monitor spanning the whole slab (subsampled, Ex/Ey/Hz only, to limit size)
mnt_field = td.FieldMonitor(
    center=(0, 0, 0),
    size=(td.inf, td.inf, 0),
    name="field",
    freqs=[freq0],
    interval_space=(4, 4, 1),
    fields=["Hz", "Ex", "Ey"],
)
# Small full-resolution monitor near the array facets, for the zoomed inset plots
mnt_field_facets = td.FieldMonitor(
    center=(-12, Ra / 2 + 2.5, 0),
    size=(24, 12, 0),
    name="field_facets",
    freqs=[freq0],
    fields=["Hz"],
)

# Transverse monitors capturing the complex field entering each array arm
mnts_wg = [
    td.FieldMonitor(
        center=(-r_mnt * np.sin(a), -Ra / 2 + r_mnt * np.cos(a), 0),
        size=(w_wg / np.cos(a) * box_factor[0], 0, 2 * t_LiTaO3 * box_factor[1]),
        freqs=[freq0],
        fields=["Ex", "Ez", "Hx", "Hz"],
        name=f"wg_{n}",
    )
    for n, a in enumerate(alpha_grt)
]

# In-plane absorbers and vertical PML to suppress boundary reflections
boundary_spec = td.BoundarySpec(
    x=td.Boundary.absorber(), y=td.Boundary.absorber(), z=td.Boundary.pml()
)

# Assemble the input-coupler simulation (mirror symmetry about x=0 reduces cost)
run_time_in = 2 / fwidth + 1.5 * Ra / (td.C_0 / n_FPR)
sim_in = td.Simulation(
    structures=[*base_layers, Rowland, wg_in, *wg_gratings],
    size=(Ra, Ra + 10 + lda_max * 5, 5),
    sources=[source],
    monitors=[mnt_field, mnt_field_facets, *mnts_wg],
    symmetry=(-1, 0, 0),
    boundary_spec=boundary_spec,
    grid_spec=grid_spec,
    run_time=run_time_in,
)

fig, ax = plt.subplots(tight_layout=True)
sim_in.plot(z=1e-5, ax=ax);

Run the Input Simulation

The input-coupler simulation is submitted, and the slab-field and array-arm monitor data are collected.

# Submit the input-coupler simulation and retrieve the results
sim_data = web.run(
    simulation=sim_in,
    task_name="Rowland_input",
    folder_name="StarCoupler",
    path="output/sim_data_in.hdf5",
    verbose=True,
)
19:27:35 UTC Created task 'Rowland_input' with resource_id                      
             'fdve-2430ac98-c960-4745-80bf-ca693312cccd' and task_type 'FDTD'.  

             View task using web UI at                                          
             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-2430ac98-c96
             0-4745-80bf-ca693312cccd'.                                         

             Task folder: 'default'.                                            







             Estimated FlexCredit cost: 8.317. This assumes the FDTD solver runs
             for the full simulation time; if early shutoff is reached, the     
             billed cost can be lower. Use 'web.real_cost(task_id)' to get the  
             billed FlexCredit cost after a simulation run.                     







             status = queued                                                    







             To cancel the simulation, use 'web.abort(task_id)' or              
             'web.delete(task_id)' or abort/delete the task in the web UI.      
             Terminating the Python script will not stop the job running on the 
             cloud.                                                             












19:27:47 UTC status = preprocess                                                













19:27:51 UTC starting up solver                                                 







             running solver                                                     












19:30:34 UTC early shutoff detected at 93%, exiting.                            













             status = postprocess                                               












19:30:40 UTC status = success                                                   













19:30:42 UTC View simulation result at                                          
             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-2430ac98-c96
             0-4745-80bf-ca693312cccd'.                                         


















             Loading results from output/sim_data_in.hdf5                       











Helpers


A few small helper functions used throughout the rest of the notebook are defined next: meshgrid_xy builds the 2D coordinate mesh for a monitor’s DataArray, add_inset_colorbar attaches a colorbar on an inset axis, inplane_poynting computes the in-plane time-averaged Poynting magnitude from a planar field monitor, and clip_outside_arc together with draw_outline draw the array-arm and Rowland/grating-arc outlines used to annotate the field plots below.






# 2D (x, y) coordinate mesh for pcolormesh, from a monitor DataArray
def meshgrid_xy(da):
    return np.meshgrid(da.x, da.y, indexing="ij")


# Colorbar on an inset axis attached to `ax`, labeled with `label`
def add_inset_colorbar(ax, mappable, label, extend=None):
    cax = ax.inset_axes([1.04, 0, 0.04, 1])
    cbar = plt.colorbar(mappable, cax=cax, extend=extend)
    cbar.ax.set_ylabel(label)
    return cbar


# In-plane time-averaged Poynting magnitude from the recorded Ex, Ey, Hz
def inplane_poynting(mnt_data):
    Ex, Ey, Hz = mnt_data.Ex.squeeze(), mnt_data.Ey.squeeze(), mnt_data.Hz.squeeze()
    Sy = 0.5 * np.real(-Ex * np.conj(Hz))
    Sx = 0.5 * np.real(Ey * np.conj(Hz))
    return np.sqrt(Sx**2 + Sy**2)


# Keep only the waveguide vertices outside the given center/radius circle, for plotting
def clip_outside_arc(wg, center_y, radius):
    vert_wg = wg.geometry.vertices
    out_of_circle = (
        np.sqrt(vert_wg[:, 0] ** 2 + (vert_wg[:, 1] - center_y) ** 2) >= radius
    )
    return vert_wg[out_of_circle]


# Draw the arm and arc outlines on ax; output=True draws the output arms instead of wg_in
def draw_outline(ax, color="w", output=False):
    for wg in wg_gratings:
        vert_wg = clip_outside_arc(wg, -Ra / 2, Ra)
        ax.plot(*vert_wg.T, color, lw=0.75)
    if output:
        for wg in wg_out:
            vert_wg = clip_outside_arc(wg, 0, Ra / 2)
            ax.plot(*vert_wg.T, color, lw=0.75)
    else:
        vert_wg = clip_outside_arc(wg_in, 0, Ra / 2)
        ax.plot(*vert_wg.T, color, lw=0.75)
    ax.plot(*vert_grt.T, color + "--", lw=0.75)
    ax.plot(*vert_rowland.T, color + "--", lw=0.75)











Field Distribution


The plot_slab_flux and plot_facet_inset helpers defined below (combined in plot_slab_with_facet_inset) draw a slab power-flux map with a zoomed facet inset, and are reused again for the output-coupler plots further down. The power flow through the input coupler is plotted here; the inset overlays \(\mathrm{Re}(H_z)\) with the focal arc and array facets, illustrating the sampled phase fronts.






# Slab power-flux map on `ax`, with the waveguide/arc outlines overlaid; `f` selects the
# frequency to plot when `S` still has a frequency dimension
def plot_slab_flux(ax, S, f=None, vmax_frac=1 / 8, output=False, title=None):
    draw_outline(ax, output=output)
    S_plot = S if f is None else S.sel(f=f)
    x, y = meshgrid_xy(S)
    vmax = S.max()
    pc = ax.pcolormesh(x, y, S_plot, vmin=0, vmax=vmax * vmax_frac, cmap=plt.cm.inferno)
    ax.set(
        aspect="equal",
        xlabel=r"$x$ (um)",
        ylabel=r"$y$ (um)",
        xlim=(-sim_in.size[0] / 2, sim_in.size[0] / 2),
        ylim=(-sim_in.size[1] / 2, sim_in.size[1] / 2),
        title=title,
    )
    add_inset_colorbar(ax, pc, r"$|\mathbf{S}|$", extend="max")
    return pc


# Zoomed Re(Hz) inset above `ax`, near the array facets, linked back to `ax` with a zoom indicator
def plot_facet_inset(ax, Hz_facets, f=None, color="k", output=False):
    axins = ax.inset_axes([0, 1, 1, 0.5])
    draw_outline(axins, color=color, output=output)
    axins.set(
        aspect="equal",
        xticks=[],
        yticks=[],
        xlim=(-24, 0),
        ylim=(Ra / 2 + 2.5 - 6, Ra / 2 + 2.5 + 6),
    )
    indicator = ax.indicate_inset_zoom(axins, edgecolor="gray", lw=3)
    indicator.connectors[0].set_visible(True)
    indicator.connectors[1].set_visible(False)
    indicator.connectors[2].set_visible(True)
    indicator.connectors[3].set_visible(False)

    Hz_plot = Hz_facets if f is None else Hz_facets.sel(f=f)
    x, y = meshgrid_xy(Hz_facets)
    vmax = np.abs(Hz_facets).max()
    pc = axins.pcolormesh(
        x, y, np.real(Hz_plot), vmin=-vmax, vmax=vmax, cmap=plt.cm.RdBu
    )
    add_inset_colorbar(axins, pc, r"Re($H_z$)")
    return axins


# Combined slab flux map + zoomed facet inset, used for both the input- and output-coupler plots
def plot_slab_with_facet_inset(
    ax, S, Hz_facets, f=None, vmax_frac=1 / 8, output=False, title=None, facet_color="k"
):
    pc = plot_slab_flux(ax, S, f=f, vmax_frac=vmax_frac, output=output, title=title)
    plot_facet_inset(ax, Hz_facets, f=f, color=facet_color, output=output)
    return pc











# Plot the power-flux distribution across the slab, with a zoomed Re(Hz) inset near the array facets
fig, ax = plt.subplots(figsize=(8, 8), tight_layout=True)

S = inplane_poynting(sim_data["field"])
Hz_facets = sim_data["field_facets"].Hz.squeeze()
plot_slab_with_facet_inset(ax, S, Hz_facets, title="Rowland input coupler");

























Input-to-Array Coupling


The flux coupled into each array arm is plotted as a function of divergence angle, revealing the angular illumination envelope and its roll-off toward the outer arms.






# Flux coupled into each array arm vs. divergence angle
flux_list = [
    sim_data[f"wg_{n}"].flux.squeeze().data / P_box for n in range(len(alpha_grt))
]

fig, ax = plt.subplots()
ax.plot(alpha_grt * 180 / np.pi, 10 * np.log10(flux_list), "o-")
ax.set(xlabel=r"Divergence angle (deg)", ylabel=r"Transmission (dB)");



























AWG Propagation and Dispersion


The waveguide array is too large to include directly in the FDTD domain. Since its arms have different physical lengths (\(L_p=L_0+p\Delta L_\mathrm{AWG}\)), a pulse propagating through arm \(p\) acquires the relative group delay


\[\Delta t_p=\frac{L_p-L_0}{v_g}
=\frac{p\Delta L}{v_g}
=\frac{p n_g\Delta L}{c},\]


where \(p\) is the arm index and \(v_g=c/n_g\). This propagation delay is applied by shifting the source pulse associated with each arm relative to arm 0 using the offset of its GaussianPulse.






# Extend the run time for the (longer) output simulation to accommodate delayed pulses
dt = sim_in.dt
Delta_T = Delta_L / (td.C_0 / n_group0)
run_time_out = run_time_in + (n_grating_arms - 1) * Delta_T
t = np.arange(0, run_time_out, dt)

# Build a group-delayed Gaussian pulse for each array arm and plot their envelopes
fig, ax = plt.subplots()
delayed_source_time_list = []
cmap = plt.cm.managua
for n, a in enumerate(alpha_grt):
    t_delay = n * Delta_T
    source_time_delayed = td.GaussianPulse(
        freq0=freq0, fwidth=fwidth, offset=5 + t_delay * (2 * np.pi * fwidth)
    )
    delayed_source_time_list.append(source_time_delayed)
    amp = source_time_delayed.amp_time(time=t)
    ax.plot(
        t / 1e-12,
        np.abs(amp),
        color=cmap(n / (n_grating_arms - 1)),
        lw=1.5,
        label=f"arm {n}" if n % 4 == 0 else None,
    )

ax.legend()
ax.set(
    xlabel=r"Time, $t$ (ps)",
    ylabel=r"Envelope, $|a(t)|$",
    ylim=(0, 1.05),
    xlim=(-1, run_time_out / 1e-12 + 2),
);

























Output Coupler Simulation


The output star coupler recombines the array-arm fields. Their relative phases move the constructive-interference focus between output waveguides as the frequency changes.




Backward-Propagating Field Sources


The output simulation reuses the complex fields recorded at the array-arm entrances. For the reciprocal propagation assumed here, time reversal transforms the phasor fields as


\[\mathbf{E}_\mathrm{rev}(\mathbf{r},\omega)=\mathbf{E}^{\ast}(\mathbf{r},\omega),\qquad
\mathbf{H}_\mathrm{rev}(\mathbf{r},\omega)=-\mathbf{H}^{\ast}(\mathbf{r},\omega).\]


Conjugation reverses the phase progression, while the magnetic field changes sign because it is odd under time reversal. The time-averaged Poynting vector therefore reverses:


\[\begin{aligned}
\mathbf{S}_\mathrm{rev}
&=\frac{1}{2}\operatorname{Re}\!\left(\mathbf{E}_\mathrm{rev}\times\mathbf{H}_\mathrm{rev}^{\ast}\right)=\frac{1}{2}\operatorname{Re}\!\left[\mathbf{E}^{\ast}\times(-\mathbf{H})\right]
=-\mathbf{S}.
\end{aligned}\]


This transformation is applied to the tangential field components, and each reversed field is paired with its arm-specific delayed waveform. The plots below show \(\mathrm{Re}(H_z)\) and \(|S_y|\) for two representative arms (indices 0 and 8) as a consistency check, prior to time-reversing every arm and converting it into a source.


Note: The array arms are too closely spaced near the coupler to support separate, nonoverlapping ModeSource. The recorded fields, by contrast, preserve the amplitude and phase delivered to each arm, which is required to reproduce the output interference correctly.






# Plot Re(Hz) and Sy recorded at two representative array-arm monitors
fig, ax = plt.subplots(2, 2, figsize=(10, 6), sharey=True, tight_layout=True)
for ii in range(2):
    sim_data.plot_field(
        field_monitor_name=f"wg_{8 * ii}", field_name="Hz", val="real", ax=ax[ii, 0]
    )
    sim_data.plot_field(
        field_monitor_name=f"wg_{8 * ii}", field_name="Sy", val="abs", ax=ax[ii, 1]
    )
    ax[ii, 0].set_title(f"Monitor {8 * ii}, " + r"$\mathrm{Re}(H_z)$")
    ax[ii, 1].set_title(f"Monitor {8 * ii}, $S_y$")

























# Time-reverse a recorded field (conjugate E, negate-conjugate H) for use as a backward source
def time_reverse_field(field):
    return field.updated_copy(
        Ex=np.conj(field.Ex),
        Ez=np.conj(field.Ez),
        Hx=-np.conj(field.Hx),
        Hz=-np.conj(field.Hz),
    )


# Convert each recorded array-arm field into a custom field source, paired with its
# group-delayed pulse
CFS_list = [
    time_reverse_field(sim_data[f"wg_{n}"]).to_source(
        source_time=delayed_source_time, center=sim_data[f"wg_{n}"].monitor.center
    )
    for n, delayed_source_time in enumerate(delayed_source_time_list)
]











Output Simulation Setup


A broadband flux monitor is placed across each of the nine output waveguides. The simulation combines the output coupler with the 17 custom array-arm sources and the input coupler’s slab/facet field monitors, now recorded at multiple frequencies for the wavelength-dependent plots below; a longer run time is used to accommodate the delayed pulses.






# Broadband flux monitors across each output waveguide
r_mnt_list = Ra * np.cos(theta_out) + 5.0 + lda_max * 1.5

mnts_wg_out = [
    td.FluxMonitor(
        center=(-r * np.sin(np.pi + theta), Ra / 2 + r * np.cos(np.pi + theta), 0),
        size=(
            w_wg / np.abs(np.cos(np.pi + theta)) * box_factor[0],
            0,
            2 * t_LiTaO3 * box_factor[1],
        ),
        freqs=freqs,
        name=f"wg_out_{n}",
    )
    for n, (r, theta) in enumerate(zip(r_mnt_list, theta_out))
]

# Assemble the output-coupler simulation: time-reversed array-arm sources + output monitors
# (reusing the input slab/facet monitors, now recorded at freqs_reduced for the spectrum plots)
sim_out = td.Simulation(
    structures=[*base_layers, Rowland, *wg_out, *wg_gratings],
    size=sim_in.size,
    sources=CFS_list,
    monitors=[
        mnt_field.updated_copy(freqs=freqs_reduced),
        mnt_field_facets.updated_copy(freqs=freqs_reduced),
        *mnts_wg_out,
    ],
    boundary_spec=boundary_spec,
    grid_spec=grid_spec,
    run_time=run_time_out,
)

fig, ax = plt.subplots(tight_layout=True)
sim_out.plot(z=1e-5, ax=ax);

























Run the Output Simulation


The output-coupler simulation is submitted, and the fields and output-port spectra are collected.






# Submit the output-coupler simulation and retrieve the results
sim_data_out = web.run(
    simulation=sim_out,
    task_name="Rowland_output",
    folder_name="StarCoupler",
    path="output/sim_data_out.hdf5",
    verbose=True,
)







19:30:59 UTC Created task 'Rowland_output' with resource_id                     
             'fdve-4c3595b2-03ea-42d2-b1e0-a8051ab7ca3c' and task_type 'FDTD'.  







             View task using web UI at                                          
             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-4c3595b2-03e
             a-42d2-b1e0-a8051ab7ca3c'.                                         







             Task folder: 'default'.                                            


















             Estimated FlexCredit cost: 67.378. This assumes the FDTD solver    
             runs for the full simulation time; if early shutoff is reached, the
             billed cost can be lower. Use 'web.real_cost(task_id)' to get the  
             billed FlexCredit cost after a simulation run.                     







             status = queued                                                    







             To cancel the simulation, use 'web.abort(task_id)' or              
             'web.delete(task_id)' or abort/delete the task in the web UI.      
             Terminating the Python script will not stop the job running on the 
             cloud.                                                             












19:31:13 UTC status = preprocess                                                













19:31:17 UTC starting up solver                                                 







             running solver                                                     
















Output Interference Pattern


The output power flow and \(\mathrm{Re}(H_z)\) near the array facets are plotted at two frequencies, \(f_0-4\Delta f\) and \(f_0+2.5\Delta f\), side by side, illustrating how the interference focus shifts along the output arc between channels.






# Plot the power-flux distribution at two frequencies, showing the focus shift along the output arc
S = inplane_poynting(sim_data_out["field"])
Hz_facets = sim_data_out["field_facets"].Hz.squeeze()

f_tar = [freq0 - 4 * f_spacing, freq0 + 2.5 * f_spacing]
titles = [r"$f = f_0 - 4\Delta f$", r"$f = f_0 +2.5\Delta f$"]

fig, ax = plt.subplots(
    1, 2, figsize=(10, 7), sharex=True, sharey=True, tight_layout=True
)
for aa, ftar, title in zip(ax, f_tar, titles):
    plot_slab_with_facet_inset(
        aa, S, Hz_facets, f=ftar, vmax_frac=1 / 4, output=True, title=title
    )

























False-color Visualization of Rowland On-chip Spectroscopy


Each simulated frequency is assigned a color, and the weighted field intensities are blended into a single false-color image, illustrating how the star coupler spatially separates the spectral channels in the manner of a free-space spectrometer.






S = inplane_poynting(sim_data_out["field"])
S = S / np.max(S)

# Assign each frequency a color and blend the weighted intensities into a false-color image
cmap = plt.cm.jet_r
spectrum = cmap(np.linspace(0, 1, len(freqs_reduced)))[:, :3]

RGB = (S.data.reshape(*S.shape, 1) * spectrum).sum(axis=2)
RGB = np.clip(RGB, a_min=0, a_max=1)

x, y = meshgrid_xy(S)

# Plot the false-color spectrometer image with waveguide and arc outlines overlaid
fig, ax = plt.subplots(figsize=(6, 6))
draw_outline(ax, output=True)
ax.set(
    aspect="equal",
    xlabel=r"$x$ (um)",
    ylabel=r"$y$ (um)",
    xlim=(-sim_in.size[0] / 2, sim_in.size[0] / 2),
    ylim=(-sim_in.size[1] / 2, sim_in.size[1] / 2),
)
ax.pcolormesh(x, y, RGB)

cax = ax.inset_axes([1.02, 0.2, 0.03, 0.6])
norm = mpl.colors.Normalize(vmin=freqs_reduced[0] / 1e12, vmax=freqs_reduced[-1] / 1e12)
sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm)
cbar = plt.colorbar(sm, cax=cax)
cbar.ax.set_ylabel("Frequency (THz)");





















The flux spectrum recorded at each output monitor is plotted, with vertical dashed lines marking the theoretical peak frequencies at 100 GHz spacing centered on \(f_0 = 198.4\) THz.






# Plot the transmission spectrum at each output waveguide, with dashed lines at the
# theoretical peak (channel) frequencies
fig, ax = plt.subplots()
for n in range(n_output_channels):
    flux = -sim_data_out[f"wg_out_{n}"].flux
    ax.plot(
        freqs / 1e12,
        10 * np.log10(flux),
        color=cmap(0.1 + 0.8 * n / (n_output_channels - 1)),
    )
    ax.axvline(
        x=(freq0 + f_spacing * (n - n_output_channels // 2)) / 1e12,
        color="k",
        ls="--",
        lw=0.75,
    )

ax.set(
    xlabel=r"Frequency (THz)",
    ylabel=r"Transmission (dB)",
    xlim=(fmin / 1e12, fmax / 1e12),
    ylim=(-20, 0),
);





















The simulated output spectra peak close to the dashed lines marking the theoretical channel frequencies, confirming that the observed frequency spacing and corresponding output-waveguide positions agree well with the Rowland focusing theory stated above.







          

          

            Example library
          

          
          

          

            

              


                

                  More Examples
                



              Read More
            

            

              

              

                
                
                

                  
                    Waveguide Y junction
                    

                    Waveguide Y junction
                  
                

                
                

                  
                    Waveguide to ring coupling
                    

                    Waveguide to ring coupling
                  
                

                
                

                  
                    Waveguide crossing based on cosine tapers
                    

                    Waveguide crossing based on cosine tapers
                  
                

                
                
              

              

              

                
No matching results found.

              

            

          

          

          
          

          

            
            
              
              Previous Tutorial
            
            
            
            
              Next Tutorial