There has been a growing interest in using silicon photonic phased arrays for applications in light detection and ranging (LIDAR) and free-space communication. For these technologies to be effectively implemented, it is crucial to have large optical apertures. These apertures are necessary to maintain a small diffraction angle and ensure a large area for receiving signals. In the context of one-dimensional arrays, there has been extensive exploration of long waveguide grating antennas (WGAs). WGAs are created by introducing periodic perturbations in an integrated waveguide and have been a popular choice for achieving large apertures, as their length can be extended to several millimeters.
In this notebook, we demonstrate a unidirectional silicon nitride WGA with a dual-layer design based on the work Manan Raval, Christopher V. Poulton, and Michael R. Watts, "Unidirectional waveguide grating antennas with uniform emission for optical phased arrays," Opt. Lett. 42, 2563-2566 (2017) DOI: 10.1364/OL.42.002563.
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.
import matplotlib.pyplot as plt
import numpy as np
import tidy3d as td
import tidy3d.web as web
from tidy3d.plugins.mode import ModeSolver
from tidy3d.plugins.mode.web import run as run_mode_solver
Design Optimization¶
In the first section of the notebook, we will design and optimize the WGA to ensure maximum unidirectional emission. This includes performing mode solving to calculate the effective indices of different waveguide sections and parameter sweeping the offset length between the two waveguides.
The WGA is designed to work at 1550 nm. Since FDTD simulation typically requires launching a broadband light source, we defined the central frequency of the pulse to be td.C_0/(1.55 μm), where td.C_0 is the speed of light in μm/s, and the frequency width to be 10% of the central frequency.
lda0 = 1.55 # central wavelength
freq0 = td.C_0 / lda0 # central frequency
fwidth = freq0 / 10 # width of the source frequency range
The waveguide is made of silicon nitride and the cladding is made of silicon oxide. Both materials can be imported from the material library directly. Note that the refractive index of the materials might not be identical to what is represented in the reference. Therefore, the simulation results are not quantitatively comparable.
SiN = td.material_library["Si3N4"]["Luke2015PMLStable"]
SiO2 = td.material_library["SiO2"]["Palik_Lossless"]
Define geometric parameters including nitride layer thickness, gap size between the layers, and so on. In this particular case, we only model 10 periods of the gratings, similar to what is demonstrated in the reference. Since Tidy3D can handle much larger simulations, a much longer WGA can be simulated in principle.
t = 0.2 # silicon nitride layer thickness
w = 1 # waveguide width of the full width regions
gap = 0.1 # gap size between the layers
L_f = 0.5 # length of the full width regions
N = 10 # number of periods
Mode Solving to Determine the Effective Indices¶
The WGA is engineered to emit vertically at a specific design wavelength of 1550 nm. Consequently, the dimensions of the four waveguide segments must be selected to ensure that the total optical path length for one cycle equals exactly one wavelength. That is,
$$\lambda_0 = (L_f-L_0)n_{full} + (L_p(p)-L_0)n_s(p) + 2L_0n_a(p),$$
where $L_f$ is the length of the full width regions, $L_0$ is the offset in the direction of propagation between the two waveguides, $L_p$ is the length of the perturbed regions, $p$ is the perturbation value, $n_{full}$, $n_s$ and $n_a$ are the effective indices of different cross sections as indicated below as well as in Fig.2(a) of the reference.
These effective indices can be calculated from mode solving. To do so, we define a function cal_n_eff that takes two arguments p1 and p2, which represent the perturbations of the top waveguide and button waveguide, perform the mode solving, and return the effective index.
def cal_n_eff(p1, p2):
# define the top waveguide
top_layer = td.Structure(
geometry=td.Box(center=(0, 0, gap / 2 + t / 2), size=(td.inf, w - 2 * p1, t)), medium=SiN
)
# define the bottom waveguide
bottom_layer = td.Structure(
geometry=td.Box(center=(0, 0, -gap / 2 - t / 2), size=(td.inf, w - 2 * p2, t)), medium=SiN
)
# define a simulation for mode solving
sim = td.Simulation(
size=(1, 3 * w, 10 * t),
grid_spec=td.GridSpec.auto(min_steps_per_wvl=50, wavelength=lda0),
structures=[top_layer, bottom_layer],
run_time=1e-12,
boundary_spec=td.BoundarySpec.all_sides(boundary=td.PML()),
medium=SiO2,
)
# define a mode solver
mode_spec = td.ModeSpec(num_modes=1, target_neff=2)
mode_solver = ModeSolver(
simulation=sim,
plane=td.Box(size=(0, td.inf, td.inf)),
mode_spec=mode_spec,
freqs=[freq0],
)
# solver for the effective index
mode_data = run_mode_solver(mode_solver, verbose=False)
n_eff = mode_data.n_eff.values[0][0]
return n_eff
$n_{full}$ is calculated by p1=p2=0. $n_a$ is calculated by p1=p, p2=0. And $n_s$ is calculated by p1=p2=p. Here we set $p$ to be 100 nm.
p = 0.1 # width perturbation
# calculate the effective indices
n_full = cal_n_eff(0, 0)
n_a = cal_n_eff(p, 0)
n_s = cal_n_eff(p, p)
Parameter Sweep to Determine the Offset¶
A unidirectional emitter can be achieved by utilizing two scattering elements spaced a quarter wavelength apart. This arrangement generates constructive interference in one vertical direction while causing destructive interference in the opposite direction. In this example, the constant thicknesses of the silicon nitride layers and the separation gap hinder the achievement of full constructive interference for upward emission. Therefore, to find the optimal $L_0$, one needs to perform a parameter sweep to find the $L_0$ value such that the downward emission is minimized.
To do so, we define a helper function make_sim that takes a given $L_0$ and returns a simulation.
def make_sim(L_0):
# calculate L_p
L_p = (lda0 + L_0 * (n_full + n_s - 2 * n_a) - L_f * n_full) / n_s
period = L_p + L_f # periodic the grating
# define the grating structures
structures = []
for i in range(N):
structures.append(
td.Structure(
geometry=td.Box(
center=(L_p / 2 + i * period, 0, -gap / 2 - t / 2), size=(L_p, w - 2 * p, t)
),
medium=SiN,
)
)
structures.append(
td.Structure(
geometry=td.Box(
center=(L_p / 2 + L_0 + i * period, 0, gap / 2 + t / 2),
size=(L_p, w - 2 * p, t),
),
medium=SiN,
)
)
structures.append(
td.Structure(
geometry=td.Box(
center=(L_p + L_f / 2 + i * period, 0, -gap / 2 - t / 2), size=(L_f, w, t)
),
medium=SiN,
)
)
structures.append(
td.Structure(
geometry=td.Box(
center=(L_p + L_f / 2 + L_0 + i * period, 0, gap / 2 + t / 2), size=(L_f, w, t)
),
medium=SiN,
)
)
# create the input/output straight waveguide structures
inf_eff = 1e3 # effective infinity
structures.append(
td.Structure(
geometry=td.Box.from_bounds(
rmin=(-inf_eff, -w / 2, -gap / 2 - t), rmax=(0, w / 2, -gap / 2)
),
medium=SiN,
)
)
structures.append(
td.Structure(
geometry=td.Box.from_bounds(
rmin=(-inf_eff, -w / 2, gap / 2), rmax=(L_0, w / 2, gap / 2 + t)
),
medium=SiN,
)
)
structures.append(
td.Structure(
geometry=td.Box.from_bounds(
rmin=(L_p + L_f / 2 + (N - 1) * period, -w / 2, -gap / 2 - t),
rmax=(inf_eff, w / 2, -gap / 2),
),
medium=SiN,
)
)
structures.append(
td.Structure(
geometry=td.Box.from_bounds(
rmin=(L_p + L_f / 2 + L_0 + (N - 1) * period, -w / 2, gap / 2),
rmax=(inf_eff, w / 2, gap / 2 + t),
),
medium=SiN,
)
)
# define a mode source to launch the fundamental mode
mode_spec = td.ModeSpec(num_modes=1, target_neff=2)
mode_source = td.ModeSource(
center=(-lda0, 0, 0),
size=(0, 3 * w, 10 * t),
source_time=td.GaussianPulse(freq0=freq0, fwidth=fwidth),
direction="+",
mode_spec=mode_spec,
mode_index=0,
)
# define a flux monitor to measure downward radiation
flux_down = td.FluxMonitor(
center=(0, 0, -8 * t), size=(td.inf, td.inf, 0), freqs=[freq0], name="flux"
)
run_time = 5e-13 # simulation run time
# define simulation
sim = td.Simulation(
center=(period * (N + 1) // 2, 0, 0),
size=(2 * N * period, 4 * w, 20 * t),
grid_spec=td.GridSpec.auto(min_steps_per_wvl=30, wavelength=lda0),
structures=structures,
sources=[mode_source],
monitors=[flux_down],
run_time=run_time,
boundary_spec=td.BoundarySpec.all_sides(boundary=td.PML()),
medium=SiO2,
symmetry=(0, -1, 0),
)
return sim
To ensure the simulation setup is correct, we create a test simulation and visualize it.
sim = make_sim(0.25)
sim.plot_3d()
We perform a parameter sweep for $L_0$ from 100 nm to 250 nm to determine the optimal offset length.
L_0_list = np.linspace(0.1, 0.25, 16) # range of the parameter sweep
# create simulations
sims = {f"L_0={L_0:.2f}": make_sim(L_0) for L_0 in L_0_list}
# create a batch and run it
batch = web.Batch(simulations=sims, verbose=True)
batch_results = batch.run(path_dir="data")
12:35:45 CEST Started working on Batch containing 16 tasks.
12:36:02 CEST Maximum FlexCredit cost: 1.227 for the whole batch.
Use 'Batch.real_cost()' to get the billed FlexCredit cost after the Batch has completed.
12:36:13 CEST Batch complete.
Extract the downward radiation flux from the sweep.
downward_radiation = []
for L_0 in L_0_list:
sim_data = batch_results[f"L_0={L_0:.2f}"]
downward_radiation.append(np.abs(sim_data["flux"].flux))
Plot the downward radiation flux as a function of $L_0$. Then determine the optimal $L_0$.
plt.plot(L_0_list, downward_radiation, linewidth=3, c="red")
plt.xlabel("Offset (μm)")
plt.ylabel("Downward radiation (%)")
plt.show()
# find the offset that results in minimal downward radiation
i = np.argmin(downward_radiation)
L_0_opt = L_0_list[i]
print(f"The optimal offset is {L_0_opt * 1e3} nm")
The optimal offset is 170.0 nm
Radiation Pattern of the Optimal Design¶
To simulate the optimal design, we also define two FieldProjectionAngleMonitor objects to compute the far-field radiation pattern upward and downward.
sim_opt = make_sim(L_0_opt) # create simulation for the optimal design
# define field projection angle ranges
theta_array = np.linspace(0, np.deg2rad(60), 100) # polar angle
phi_array = np.linspace(0, 2 * np.pi, 200) # azimuthal angle
# define the top field projection monitor
monitor_top = td.FieldProjectionAngleMonitor(
name="top",
center=(0, 0, 8 * t),
size=(td.inf, td.inf, 0),
freqs=[freq0],
theta=theta_array,
phi=phi_array,
)
# define the bottom field projection monitor
monitor_bottom = td.FieldProjectionAngleMonitor(
name="bottom",
normal_dir="-",
center=(0, 0, -8 * t),
size=(td.inf, td.inf, 0),
freqs=[freq0],
theta=np.pi - theta_array,
phi=phi_array,
)
# add the field projection monitors to the simulation
sim_opt = sim_opt.copy(update={"monitors": [monitor_top, monitor_bottom]})
# run the simulation
sim_data_opt = web.run(simulation=sim_opt, task_name="optimized wga")
12:36:33 CEST Created task 'optimized wga' with task_id 'fdve-a3953bd2-7d28-457f-8936-d45738c2f392' and task_type 'FDTD'.
View task using web UI at 'https://tidy3d.simulation.cloud/workbench?taskId=fdve-a3953bd2-7d 28-457f-8936-d45738c2f392'.
Task folder: 'default'.
12:36:35 CEST Maximum FlexCredit cost: 0.081. Minimum cost depends on task execution details. Use 'web.real_cost(task_id)' to get the billed FlexCredit cost after a simulation run.
12:36:36 CEST status = success
12:36:40 CEST loading simulation from simulation_data.hdf5
Plot the normalized upward radiation pattern.
P = sim_data_opt["top"].power.squeeze(drop=True) # extract upward radiated power
P_max = np.max(P)
P_norm = P / P_max
# plot the radiation pattern in polar coordinates
fig = plt.figure()
ax = fig.add_subplot(111, polar=True)
c = ax.pcolor(
phi_array, np.rad2deg(theta_array), P_norm, shading="auto", cmap="hot", vmin=0, vmax=1
)
cb = fig.colorbar(c, ax=ax)
cb.set_label("Intensity")
_ = plt.setp(ax.get_yticklabels(), color="white")
Plot the normalized downward radiation pattern on the same color scale. In comparison, we can see that the downward radiation is nearly invisible since it is orders of magnitude weaker than the upward radiation as desired.
P = sim_data_opt["bottom"].power.squeeze(drop=True) # extract downward radiated power
P_norm = P / P_max
# plot the radiation pattern in polar coordinates
fig = plt.figure()
ax = fig.add_subplot(111, polar=True)
c = ax.pcolor(
phi_array, np.rad2deg(theta_array), P_norm, shading="auto", cmap="hot", vmin=0, vmax=1
)
cb = fig.colorbar(c, ax=ax)
cb.set_label("Intensity")
_ = plt.setp(ax.get_yticklabels(), color="white")
TERMS & CONDITIONS
How the process works:
Try this real URL to see an example: