Waveguide Y junction

Power splitters such as Y-junctions are widely used in photonic integrated circuits across different applications. When designing a power splitter, we aim to achieve a flat broadband response, low insertion loss, and compact footprint. At the same time, the design needs to comply with the fabrication resolution and tolerance.

In this example, we demonstrate the modeling of a Y-junction for integrated photonics. The designed device shows an average insertion loss below 0.2 dB in the wavelength range of 1500 nm to 1600 nm. At the same time, it has a small footprint. The junction area is smaller than 2 \(\mu m\) by 2 \(\mu m\), much smaller than the typical power splitters based on multimode interference devices. The design is adapted from Yi Zhang, Shuyu Yang, Andy Eu-Jin Lim, Guo-Qiang Lo, Christophe Galland, Tom Baehr-Jones, and Michael Hochberg, “A compact and low loss Y-junction for submicron silicon waveguide,” Opt. Express 21, 1310-1316 (2013).

For more integrated photonic examples such as the 8-Channel mode and polarization de-multiplexer, the broadband bi-level taper polarization rotator-splitter, and the broadband directional coupler, please visit our examples page.

If you are new to the finite-difference time-domain (FDTD) method, we highly recommend going through our FDTD101 tutorials.

FDTD simulations can diverge due to various reasons. If you run into any simulation divergence issues, please follow the steps outlined in our troubleshooting guide to resolve it.

import numpy as np
import matplotlib.pyplot as plt
import gdstk

import tidy3d as td
import tidy3d.web as web
from tidy3d.plugins.mode import ModeSolver

Simulation Setup

Define simulation wavelength range to be 1.5 \(\mu m\) to 1.6 \(\mu m\).

lda0 = 1.55  # central wavelength
freq0 = td.C_0 / lda0  # central frequency
ldas = np.linspace(1.5, 1.6, 101)  # wavelength range
freqs = td.C_0 / ldas  # frequency range
fwidth = 0.5 * (np.max(freqs) - np.min(freqs))  # width of the source frequency range

In this model, the Y-junction is made of silicon. The top cladding is made of silicon oxide. We will directly use the silicon and oxide media from Tidy3D’s material library. More specifically, we use the data from the widely used Handbook of Optical Constants of Solids by Palik.

si = td.material_library["cSi"]["Palik_Lossless"]

sio2 = td.material_library["SiO2"]["Palik_Lossless"]

The junction is discretized into 13 segments. Each segment is a tapper with the given widths. The optimum design is obtained by optimizing the 13 width parameters using the Particle Swarm Optimization algorithm. For the sake of simplicity, in this notebook, we skip the optimization procedure and only present the optimized result.

t = 0.22  # thickness of the silicon layer

# width of the 13 segments
w1 = 0.5
w2 = 0.5
w3 = 0.6
w4 = 0.7
w5 = 0.9
w6 = 1.26
w7 = 1.4
w8 = 1.4
w9 = 1.4
w10 = 1.4
w11 = 1.31
w12 = 1.2
w13 = 1.2

l_in = 1  # input waveguide length
l_junction = 2  # length of the junction
l_bend = 6  # horizontal length of the waveguide bend
h_bend = 2  # vertical offset of the waveguide bend
l_out = 1  # output waveguide length
inf_eff = 100  # effective infinity

First, define the junction structure by using a PolySlab. The vertices are given by the widths of the segments defined above. If a smooth curve is desirable, one can interpolate the vertices to a finer grid using spline for example.

Before proceeding further to construct other structures, we can use the plot method to inspect the geometry.

x = np.linspace(l_in, l_in + l_junction, 13)  # x coordinates of the top edge vertices
y = np.array(
    [w1, w2, w3, w4, w5, w6, w7, w8, w9, w10, w11, w12, w13]
)  # y coordinates of the top edge vertices

# using concatenate to include bottom edge vertices
x = np.concatenate((x, np.flipud(x)))
y = np.concatenate((y / 2, -np.flipud(y / 2)))

# stacking x and y coordinates to form vertices pairs
vertices = np.transpose(np.vstack((x, y)))

junction = td.Structure(
    geometry=td.PolySlab(vertices=vertices, axis=2, slab_bounds=(0, t)), medium=si
)
junction.plot(z=t / 2)
plt.show()

For the output waveguide bends, we use S bend sine waveguides, which are described by the function

\(y = \frac{xh_{band}}{l_{bend}}-\frac{h_{bend}}{2\pi}sin(\frac{2\pi x}{l_{bend}})\).

Different types of bend can also be used here for similar performance. To define the bends, the most convenient way is to use gdstk. First compute the \(x\) and \(y\) coordinates of the bend, then define a FlexPath and add it to a gds cell. The paths in the cell can be converted to Tidy3D PolySlabs using the from_gds method.

x_start = l_in + l_junction  # x coordinate of the starting point of the waveguide bends

x = np.linspace(x_start, x_start + l_bend, 100)  # x coordinates of the top edge vertices

y = (
    (x - x_start) * h_bend / l_bend
    - h_bend * np.sin(2 * np.pi * (x - x_start) / l_bend) / (np.pi * 2)
    + w13 / 2
    - w1 / 2
)  # y coordinates of the top edge vertices

# adding the last point to include the straight waveguide at the output
x = np.append(x, inf_eff)
y = np.append(y, y[-1])

# add path to the cell
cell = gdstk.Cell("bends")
cell.add(gdstk.FlexPath(x + 1j * y, w1, layer=1, datatype=0))  # top waveguide bend
cell.add(gdstk.FlexPath(x - 1j * y, w1, layer=1, datatype=0))  # bottom waveguide bend

# define top waveguide bend structure
wg_bend_1 = td.Structure(
    geometry=td.PolySlab.from_gds(
        cell,
        gds_layer=1,
        axis=2,
        slab_bounds=(0, t),
    )[0],
    medium=si,
)

# define bottom waveguide bend structure
wg_bend_2 = td.Structure(
    geometry=td.PolySlab.from_gds(
        cell,
        gds_layer=1,
        axis=2,
        slab_bounds=(0, t),
    )[1],
    medium=si,
)

# plot the top waveguide bend to visualize
ax = wg_bend_1.plot(z=t / 2)
ax.set_xlim(2, 10)
plt.show()

Lastly, define the straight input waveguide using Box.

# straight input waveguide
wg_in = td.Structure(
    geometry=td.Box.from_bounds(rmin=(-inf_eff, -w1 / 2, 0), rmax=(l_in, w1 / 2, t)),
    medium=si,
)

# the entire model is the collection of all structures defined so far
y_junction = [wg_in, junction, wg_bend_1, wg_bend_2]

Define the simulation domain. Here we ensure sufficient buffer spacing in each direction. In general, we want to make sure that the structure is at least half a wavelength away from the domain boundaries unless it goes into the PML.

Lx = l_in + l_junction + l_out + l_bend  # simulation domain size in x direction
Ly = w13 + 2 * h_bend + 1.5 * lda0  # simulation domain size in y direction
Lz = 10 * t  # simulation domain size in z direction
sim_size = (Lx, Ly, Lz)

We will use a ModeSource to excite the input waveguide using the fundamental TE mode.

A ModeMonitor is placed at the top output waveguide to measure the transmission. A FieldMonitor is added to the xy plane to visualize the power flow.

# add a mode source as excitation
mode_spec = td.ModeSpec(num_modes=1, target_neff=3.5)
mode_source = td.ModeSource(
    center=(l_in / 2, 0, t / 2),
    size=(0, 4 * w1, 6 * t),
    source_time=td.GaussianPulse(freq0=freq0, fwidth=fwidth),
    direction="+",
    mode_spec=mode_spec,
    mode_index=0,
)

# add a mode monitor to measure transmission at the output waveguide
mode_monitor = td.ModeMonitor(
    center=(l_in + l_junction + l_bend + l_out / 2, w13 / 2 - w1 / 2 + h_bend, t / 2),
    size=(0, 4 * w1, 6 * t),
    freqs=freqs,
    mode_spec=mode_spec,
    name="mode",
)

# add a field monitor to visualize field distribution at z=t/2
field_monitor = td.FieldMonitor(
    center=(0, 0, t / 2), size=(td.inf, td.inf, 0), freqs=[freq0], name="field"
)

Set up the simulation with the previously defined structures, source, and monitors. All boundaries are set to PML to mimic infinite open space. Since the top and bottom claddings are silicon oxide, we will set the medium of the background to silicon oxide.

In principle, we can impose symmetry to reduce the computational load. Since this model is relatively small and quick to solve, we will simply model the whole device without using symmetry. To learn more about using symmetry, please refer to the dedicated tutorial.

Using an automatic nonuniform grid is the most efficient and convenient. We set min_steps_per_wvl=20 to achieve a very accurate result while still keeping the simulation cost at a minimum.

run_time = 5e-13  # simulation run time

# construct simulation
sim = td.Simulation(
    center=(Lx / 2, 0, 0),
    size=sim_size,
    grid_spec=td.GridSpec.auto(min_steps_per_wvl=20, wavelength=lda0),
    structures=y_junction,
    sources=[mode_source],
    monitors=[mode_monitor, field_monitor],
    run_time=run_time,
    boundary_spec=td.BoundarySpec.all_sides(boundary=td.PML()),
    medium=sio2,
)

sim.plot(z=0)
plt.show()

To have a better visualization, we can also plot the simulation in 3D.

sim.plot_3d()

Before submitting the simulation to the server, it is a good practice to visualize the mode profile at the ModeSource to ensure we are launching the fundamental TE mode. To do so, we will use the ModeSolver plugin, which solves for the mode profile on your local computer.

mode_solver = ModeSolver(
    simulation=sim,
    plane=td.Box(center=mode_source.center, size=mode_source.size),
    mode_spec=mode_spec,
    freqs=[freq0],
)
mode_data = mode_solver.solve()

12:57:16 Eastern Standard Time WARNING: Use the remote mode solver with subpixel
                               averaging for better accuracy through            
                               'tidy3d.plugins.mode.web.run(...)'.              

Visualize the mode profile. We confirm that we are exciting the waveguide with the fundamental TE mode.

f, (ax1, ax2, ax3) = plt.subplots(1, 3, tight_layout=True, figsize=(10, 3))
abs(mode_data.Ex.isel(mode_index=0)).plot(x="y", y="z", ax=ax1, cmap="magma")
abs(mode_data.Ey.isel(mode_index=0)).plot(x="y", y="z", ax=ax2, cmap="magma")
abs(mode_data.Ez.isel(mode_index=0)).plot(x="y", y="z", ax=ax3, cmap="magma")

ax1.set_title("|Ex(x, y)|")
ax1.set_aspect("equal")
ax2.set_title("|Ey(x, y)|")
ax2.set_aspect("equal")
ax3.set_title("|Ez(x, y)|")
ax3.set_aspect("equal")
plt.show()

Now that we verified all the settings, we are ready to submit the simulation job to the server. Before running the simulation, we can get a cost estimation using estimate_cost. This prevents us from accidentally running large jobs that we set up by mistake. The estimated cost is the maximum cost corresponding to running all the time steps.

job = web.Job(simulation=sim, task_name="y_junction", verbose=True)
estimated_cost = web.estimate_cost(job.task_id)

12:57:17 Eastern Standard Time Created task 'y_junction' with task_id
                               'fdve-5f0892b7-fb7a-4844-b09a-c0e1cd3f1d8e' and
                               task_type 'FDTD'.

                               View task using web UI at
                               'https://tidy3d.simulation.cloud/workbench?taskId
                               =fdve-5f0892b7-fb7a-4844-b09a-c0e1cd3f1d8e'.

12:57:19 Eastern Standard Time Maximum FlexCredit cost: 0.112. Minimum cost
                               depends on task execution details. Use
                               'web.real_cost(task_id)' to get the billed
                               FlexCredit cost after a simulation run.

The cost is reasonaly so we can run the simulation.

sim_data = job.run(path="data/simulation_data.hdf5")

12:57:20 Eastern Standard Time status = queued

12:57:25 Eastern Standard Time status = preprocess

12:57:27 Eastern Standard Time Maximum FlexCredit cost: 0.112. Use
                               'web.real_cost(task_id)' to get the billed
                               FlexCredit cost after a simulation run.

                               starting up solver

                               running solver

                               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.

12:58:12 Eastern Standard Time early shutoff detected at 76%, exiting.

                               status = postprocess

12:58:18 Eastern Standard Time status = success

12:58:19 Eastern Standard Time View simulation result at
                               'https://tidy3d.simulation.cloud/workbench?taskId
                               =fdve-5f0892b7-fb7a-4844-b09a-c0e1cd3f1d8e'.

12:58:20 Eastern Standard Time loading simulation from data/simulation_data.hdf5

Result Visualization

After the simulation is complete, we first inspect the insertion loss. Within this wavelength range, we see that the insertion loss is generally below 0.2 dB.

# extract the transmission data from the mode monitor
amp = sim_data["mode"].amps.sel(mode_index=0, direction="+")
T = np.abs(amp) ** 2 # transmission to the top waveguide
T_total = 2 * T # total transmission at the two output waveguides

plt.plot(ldas, 10 * np.log10(T_total))
plt.xlim(1.5, 1.6)
plt.ylim(-0.5, 0)
plt.xlabel("Wavelength ($\mu m$)")
plt.ylabel("Insertion loss (dB)")
plt.show()

We can also visualize the field distribution. Here we can see the interference in the junction while no visible higher order modes are excited at the output waveguides.

sim_data.plot_field(
    field_monitor_name="field", field_name="E", val="abs^2", f=freq0, vmin=0, vmax=2e3
)
plt.show()