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 matplotlib.pyplot as plt
import numpy as np
import photonforge as pf
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 a 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{x h_{\text{bend}}}{l_{\text{bend}}} - \frac{h_{\text{bend}}}{2\pi} \sin\left(\frac{2\pi x}{l_{\text{bend}}}\right). $$
The easiest way is using PhotonForge’s layout capabilities. PhotonForge also provides a complete set of tools for photonic design automation, including PDK integration, parametric components, connectivity management, and technology definitions. You can check the many applications of this powerful tool here.
In this example, the waveguide bends are defined using the pf.Path object, which accepts the arguments origin and width. The full trajectory of the bend is defined by assigning a sequence of points to the path, reproducing the desired smooth S-shaped transition.
The straight output waveguide is defined separately as another pf.Path object and created with the segment method, ensuring continuity after the S-bend.
Finally, we use pf.boolean to combine the structures into polygons, and define the Tidy3D geometry as a PolySlab object.
x_start = l_in + l_junction # starting x-position of the bends
dx = l_bend # horizontal length of the S-bend
dy = h_bend # vertical offset of the S-bend
y0 = w13 / 2 - w1 / 2 # initial vertical position (top waveguide center)
# create top S-bend + straight in a single path
bend_up = pf.Path(origin=(x_start, y0), width=w1)
bend_up.s_bend((dx, dy), relative=True)
bend_up.segment((inf_eff, 0), relative=True)
# create bottom by mirroring
bend_down = bend_up.copy().mirror()
# boolean union of all paths into polygons
polygons = pf.boolean([bend_up, bend_down], [], "+")
# convert polygons into Tidy3D PolySlab geometries
geoms = [td.PolySlab(vertices=p.vertices, axis=2, slab_bounds=(0, t)) for p in polygons]
# create Tidy3D structures
wg_bend_1 = td.Structure(geometry=geoms[0], medium=si)
wg_bend_2 = td.Structure(geometry=geoms[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 anyway, we will simply model the whole device without using symmetry. To use symmetry, we just need to set symmetry=(0,-1,1) in the Simulation object initiation. 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()
11:37:00 -03 WARNING: Use the remote mode solver with subpixel averaging for better accuracy through 'tidy3d.web.run(...)' or the deprecated '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)
Loading simulation from local cache. View cached task using web UI at 'https://tidy3d.simulation.cloud/workbench?taskId=fdve-f40dbc85-345 4-4aa6-b46f-96830081d830'.
11:37:03 -03 Estimated FlexCredit cost: 0.113. 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 reasonable so we can run the simulation.
sim_data = job.run(path="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(r"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()
TERMS & CONDITIONS
How the process works:
Try this real URL to see an example: