Grating coupler

See on github, run on colab, or just follow along with the output below.

[1]:
# get the most recent version of tidy3d
!pip install -q --upgrade tidy3d

# make sure notebook plots inline
%matplotlib inline

# basic imports
import numpy as np
import matplotlib.pylab as plt

# tidy3d imports
import tidy3d as td
import tidy3d.web as web

Problem Setup

In this example, we model a 3D grating coupler in a Silicon on Insulator (SOI) platform.

A basic schematic of the design is shown below. The simulation is about 19um x 4um x 5um with a wavelength of 1.55um and takes about 1 minute to simulate 10,000 time steps.

In the simulation, we inject a modal source into the waveguide and propagate it towards the grating structure. The radiation from the grating coupler is then measured with a near field monitor and we use a far field projection to inspect the angular dependence of the radiation.

fc93ffa278a7437ea82b16c1720dd085

[2]:
# basic parameters (note, all length units are microns)
nm = 1e-3
wavelength = 1550 * nm

# resolution
grids_per_wavelength = 50.0
dl = wavelength / grids_per_wavelength

# waveguide
wg_width = 400 * nm
wg_height = 220 * nm
wg_length = 2 * wavelength

# surrounding
sub_height = 2.0
air_height = 2.0
buffer = 0.5 * wavelength

# coupler
cp_width = 4 * wavelength
cp_length = 8 * wavelength
taper_length = 6 * wavelength
[3]:
# sizes
Lx = buffer + wg_length + taper_length + cp_length
Ly = buffer + cp_width + buffer
Lz = sub_height + wg_height + air_height
sim_size = [Lx, Ly, Lz]

# convenience variables to store center of coupler and waveguide
wg_center_x = +Lx/2 - buffer - (wg_length + taper_length)/2
cp_center_x = -Lx/2 + buffer + cp_length/2
wg_center_z = -Lz/2 + sub_height + wg_height/2
cp_center_z = -Lz/2 + sub_height + wg_height/2

# materials
Air = td.Medium(epsilon=1.0)
Si = td.Medium(epsilon=3.47**2)
SiO2 = td.Medium(epsilon=1.44**2)

# source parameters
freq0 = td.constants.C_0 / wavelength
fwidth = freq0 / 10
run_time = 100 / fwidth

# PML layers
npml = 15

Mode Solve

To determine the pitch of the waveguide for a given design angle, we need to compute the effective index of the waveguide mode being coupled into. For this, we set up a simple simulation of the coupler region and use the mode solver to get the effective index. We will not run this simulation, we just add a ModeMonitor object in order to call the mode solver, sim.compute_modes() below, and get the effective index of the wide-waveguide region.

[4]:
# grating parameters
design_theta_deg = 30
design_theta_rad = np.pi * design_theta_deg / 180
grating_height = 70 * nm

# do a mode solve to get neff of the coupler

sub = td.Box(
    center=[0, 0, -Lz/2],
    size=[td.inf, td.inf, 2 * sub_height],
    material=SiO2,
    name='substrate')

cp = td.Box(
    center=[0, 0, cp_center_z-grating_height/4],
    size=[td.inf, cp_width, wg_height-grating_height/2],
    material=Si,
    name='coupler')

mode_mnt = td.ModeMonitor(
    center=[0, 0, 0],
    size=[0, 8*cp_width, 8*wg_height],
    freqs=freq0)

sim = td.Simulation(
    size=sim_size,
    mesh_step=[dl, dl, dl],
    structures=[sub, cp],
    sources=[],
    monitors=[mode_mnt],
    run_time=run_time,
    pml_layers=[npml, npml, npml])

sim.viz_mat_2D(normal='x');
Initializing simulation...
Mesh step (micron): [3.10e-02, 3.10e-02, 3.10e-02].
Simulation domain in number of grid points: [855, 280, 166].
Total number of computational grid points: 3.97e+07.
Total number of time steps: 96226.
Estimated data size (GB) of monitor monitor: 0.0015.
[4]:
<AxesSubplot:title={'center':'yz-plane at x=0.00e+00um'}, xlabel='y (um)', ylabel='z (um)'>
../_images/examples_GratingCoupler_7_2.png
[5]:
# Compute and visualize the first two modes of the monitor
sim.compute_modes(mode_mnt, Nmodes=1)
sim.viz_modes(mode_mnt)

# Get the data for the mode corresponding to frequency index 0 and mode index 0
mode = sim.data(mode_mnt)["modes"][0][0]
neff = mode.neff
print(f'\n\neffective index of coupler region, mode 0 = {neff}')


effective index of coupler region, mode 0 = 2.6997185508451422
/home/momchil/miniconda3/envs/tidy_tmp/lib/python3.6/site-packages/IPython/core/pylabtools.py:132: UserWarning: constrained_layout not applied.  At least one axes collapsed to zero width or height.
  fig.canvas.print_figure(bytes_io, **kw)
../_images/examples_GratingCoupler_8_2.png

Create Simulation

Now we set up the grating coupler to simulate in Tidy3D.

[6]:
# gratings
pitch = wavelength / (neff - np.sin(design_theta_rad))
grating_length = pitch / 2.0
num_gratings = int(cp_length / pitch)

sub = td.Box(
    center=[0, 0, -Lz/2],
    size=[td.inf, td.inf, 2 * sub_height],
    material=SiO2,
    name='substrate')

wg = td.Box(
    center=[wg_center_x, 0, wg_center_z],
    size=[buffer + wg_length + taper_length + cp_length/2, wg_width, wg_height],
    material=Si,
    name='waveguide')

cp = td.Box(
    center=[cp_center_x, 0, cp_center_z],
    size=[cp_length, cp_width, wg_height],
    material=Si,
    name='coupler')

tp = td.PolySlab(
    vertices=[
        [cp_center_x + cp_length/2 + taper_length, + wg_width/2],
        [cp_center_x + cp_length/2 + taper_length, - wg_width/2],
        [cp_center_x + cp_length/2, - cp_width/2],
        [cp_center_x + cp_length/2, + cp_width/2]],
    z_cent=wg_center_z,
    z_size=wg_height,
    material=Si,
    name='taper')

grating_left_x = cp_center_x - cp_length/2
gratings = [
    td.Box(
        center=[grating_left_x + (i + 0.5) * pitch, 0, cp_center_z + wg_height/2 - grating_height/2],
        size=[grating_length, cp_width, grating_height],
        material=Air,
        name=f'{i}th_grating')
    for i in range(num_gratings)]
[7]:
mode_source = td.ModeSource(
    td.GaussianPulse(freq0, fwidth, phase=0),
    center=[Lx/2 - buffer, 0, cp_center_z],
    size=[0, 8*wg_width, 8*wg_height],
    direction='backward',
    amplitude=1.0,
    name='modal_source')
[8]:
# distance to near field monitor
nf_offset = 50 * nm

plane_monitor = td.FreqMonitor(
    center=[0, 0, cp_center_z],
    size=[Lx, Ly, 0],
    freqs=freq0,
    name='full_domain_fields')

rad_monitor = td.FreqMonitor(
    center=[0, 0, 0],
    size=[Lx, 0, Lz],
    freqs=freq0,
    name='full_domain_fields')

near_field_monitor = td.FreqMonitor(
    center=[cp_center_x, 0, cp_center_z + wg_height/2 + nf_offset],
    size=[cp_length, cp_width, 0],
    freqs=freq0,
    name='radiated_near_fields')
[9]:
sim = td.Simulation(
    size=sim_size,
    mesh_step=[dl, dl, dl],
    structures=[sub, wg, cp, tp] + gratings,
    sources=[mode_source],
    monitors=[plane_monitor, rad_monitor, near_field_monitor],
    run_time=run_time,
    pml_layers=[npml, npml, npml])
Initializing simulation...
Mesh step (micron): [3.10e-02, 3.10e-02, 3.10e-02].
Simulation domain in number of grid points: [855, 280, 166].
Total number of computational grid points: 3.97e+07.
Total number of time steps: 96226.
Estimated data size (GB) of monitor full_domain_fields: 0.0099.
Estimated data size (GB) of monitor full_domain_fields_1: 0.0054.
Estimated data size (GB) of monitor radiated_near_fields: 0.0038.
[10]:
fig, axes = plt.subplots(1, 3, tight_layout=True, figsize=(14, 3))
for val, pos, ax in zip('xyz', (0, 0, -Lz/2 + sub_height + wg_height/2), axes):
    sim.viz_eps_2D(normal=val, position=pos, ax=ax, cbar=False)
../_images/examples_GratingCoupler_14_0.png
[11]:
ax1, ax2 = sim.viz_source(mode_source)
ax1.set_xlim((0, 0.1e-12))  # note the pulse extends far beyond this time, adjust lims to inspect
plt.show()
../_images/examples_GratingCoupler_15_0.png
[12]:
sim.compute_modes(mode_source, Nmodes=2)
sim.viz_modes(mode_source)

# Mode 0 is the Ey dominant mode, so we choose that
mode_ind = 0
sim.set_mode(mode_source, mode_ind)
Mode set, recommend verifying using viz_modes.
../_images/examples_GratingCoupler_16_1.png

Run Simulation

Run the simulation and plot the field patterns

[13]:
# create a project, upload to our server to run
project = web.run(sim, task_name='grating_coupler')
with open("output/tidy3d.log") as f:
     print(f.read())
Uploading the json file...
Project 'grating_coupler-16b8755264cf2b50' status: success...

Downloading results...
Loading monitor data into simulation.
Applying source normalization to all frequency monitors using source index 0.
Simulation domain Nx, Ny, Nz: [104, 104, 104]
Applied symmetries: [0, 0, 0]
Number of computational grid points: 1.1249e+06.
Using subpixel averaging: True
Number of time steps: 2308
Automatic shutoff factor: 1.00e-05
Time step (s): 8.6662e-17

Compute source modes time (s):     0.0972
Compute monitor modes time (s):    0.1224

Rest of setup time (s):            0.0605

Starting solver...
- Time step     92 / time 7.97e-15s (  4 % done), field decay: 1.00e+00
- Time step    184 / time 1.59e-14s (  8 % done), field decay: 1.00e+00
- Time step    276 / time 2.39e-14s ( 12 % done), field decay: 1.00e+00
- Time step    369 / time 3.20e-14s ( 16 % done), field decay: 1.00e+00
- Time step    461 / time 4.00e-14s ( 20 % done), field decay: 1.00e+00
- Time step    553 / time 4.79e-14s ( 24 % done), field decay: 1.00e+00
- Time step    646 / time 5.60e-14s ( 28 % done), field decay: 1.00e+00
- Time step    738 / time 6.40e-14s ( 32 % done), field decay: 1.00e+00
- Time step    830 / time 7.19e-14s ( 36 % done), field decay: 1.00e+00
- Time step    918 / time 7.96e-14s ( 39 % done), field decay: 1.00e+00
- Time step    923 / time 8.00e-14s ( 40 % done), field decay: 1.00e+00
- Time step   1015 / time 8.80e-14s ( 44 % done), field decay: 1.00e+00
- Time step   1107 / time 9.59e-14s ( 48 % done), field decay: 8.93e-01
- Time step   1200 / time 1.04e-13s ( 52 % done), field decay: 3.56e-01
- Time step   1292 / time 1.12e-13s ( 56 % done), field decay: 1.34e-01
- Time step   1384 / time 1.20e-13s ( 60 % done), field decay: 4.49e-02
- Time step   1477 / time 1.28e-13s ( 64 % done), field decay: 9.47e-03
- Time step   1569 / time 1.36e-13s ( 68 % done), field decay: 1.50e-03
- Time step   1661 / time 1.44e-13s ( 72 % done), field decay: 2.60e-04
- Time step   1754 / time 1.52e-13s ( 76 % done), field decay: 5.43e-05
- Time step   1846 / time 1.60e-13s ( 80 % done), field decay: 1.15e-05
- Time step   1938 / time 1.68e-13s ( 84 % done), field decay: 2.11e-06
Field decay smaller than shutoff factor, exiting solver.

Solver time (s):                   3.9134
Post-processing time (s):          0.0232

[14]:
fig, axes = plt.subplots(3, 1, tight_layout=True, figsize=(10, 8))
for monitor, cpmp, ax in zip([plane_monitor, rad_monitor, near_field_monitor], 'yyy', axes):
    sim.viz_field_2D(monitor, comp=cpmp, ax=ax, cbar=True)
../_images/examples_GratingCoupler_19_0.png

Far Field Projection

Now we use the Near2Far feature of Tidy3D to compute the anglular dependence of the far field scattering based on the near field monitor.

[15]:
# create range of angles to probe (note: polar coordinates, theta = 0 corresponds to vertical (z axis))
num_angles = 1101
thetas = np.linspace(-np.pi/2, np.pi/2, num_angles)

# make a near field to far field projector with the near field monitor data
n2f = td.Near2Far(sim.data(near_field_monitor))

# loop through angles and record the scattered cross section
Ps = np.zeros(num_angles)
for i in range(num_angles):
    Ps[i] = n2f.get_radar_cross_section(thetas[i], 0.0)
[16]:
# plot the angle dependence
fig, ax = plt.subplots(subplot_kw={'projection': 'polar'}, figsize=(5,5))
ax.plot(thetas, Ps, label='measured')
ax.plot([design_theta_rad, design_theta_rad], [0, np.max(Ps)*0.7], 'r--', alpha=0.8, label='design angle')
ax.set_theta_zero_location("N")
ax.set_yticklabels([])
ax.set_title("Scattered Cross-section (arb. units)", va='bottom')
plt.legend()
plt.show()
../_images/examples_GratingCoupler_22_0.png
[17]:
theta_expected = np.arcsin(np.abs(neff - wavelength / pitch))
print(f'expect angle of {(theta_expected * 180 / np.pi):.2f} degrees')
i_max = np.argmax(Ps)
print(f'got maximum angle of {(thetas[i_max] * 180 / np.pi):.2f} degrees')
expect angle of 30.00 degrees
got maximum angle of 28.64 degrees

The agreement between the target angle and the actual emission angle of the coupler is very good. The small difference comes from the fact that the design is very sensitive to the value of the effective index that we use in the coupler region, and that value depends on which waveguide height we pick in that region: the one with the grating comb, or without. In our setup, we used a thickness that is at the mid-point, but this is a heuristic choice which results in the small final mismatch in angles observed here.

Gaussian beam into the coupler

We can also run the coupler in the opposite way, injecting a Gaussian beam from above and monitoring the transmission into the waveguide.

[18]:
gaussian_beam = td.GaussianBeam(
    normal='z',
    center=[-5, 0, 2],
    source_time=td.GaussianPulse(freq0, fwidth),
    angle_theta=np.pi/6,
    angle_phi=np.pi,
    direction='backward',
    waist_radius=2,
    pol_angle=np.pi/2)

mnt_out = td.ModeMonitor(
    center=[Lx/2 - buffer, 0, cp_center_z],
    size=[0, 8*wg_width, 8*wg_height],
    store=['mode_amps', 'flux'],
    freqs=[freq0])

sim = td.Simulation(
    size=sim_size,
    mesh_step=[dl, dl, dl],
    structures=[sub, wg, cp, tp] + gratings,
    sources=[gaussian_beam],
    monitors=[plane_monitor, rad_monitor, mnt_out],
    run_time=run_time,
    pml_layers=[npml, npml, npml])

fig, axes = plt.subplots(1, 3, tight_layout=True, figsize=(14, 3))
for val, pos, ax in zip('xyz', (0, 0, -Lz/2 + sub_height + wg_height/2), axes):
    sim.viz_eps_2D(normal=val, position=pos, ax=ax, cbar=False)
Initializing simulation...
Mesh step (micron): [3.10e-02, 3.10e-02, 3.10e-02].
Simulation domain in number of grid points: [855, 280, 166].
Total number of computational grid points: 3.97e+07.
Total number of time steps: 96226.
Estimated data size (GB) of monitor full_domain_fields: 0.0099.
Estimated data size (GB) of monitor full_domain_fields_1: 0.0054.
Estimated data size (GB) of monitor monitor: 0.0000.
../_images/examples_GratingCoupler_26_1.png
[19]:
project = web.run(sim, task_name='grating_coupler')
Uploading the json file...
Project 'grating_coupler-16b875707ebf5698' status: success...

Downloading results...
Loading monitor data into simulation.
Applying source normalization to all frequency monitors using source index 0.
[20]:
mdata = sim.data(mnt_out)
print("Transmission into waveguide: ", mdata['flux'][0, 0])

fig, ax = plt.subplots(2, 1, figsize=(8, 8))
sim.viz_field_2D(plane_monitor, val='re', comp='y', ax=ax[0], cbar=True)
sim.viz_field_2D(rad_monitor, val='re', comp='y', ax=ax[1], cbar=True);
Transmission into waveguide:  0.1046307257238687
[20]:
<matplotlib.collections.QuadMesh at 0x7ff6b50a18d0>
../_images/examples_GratingCoupler_28_2.png

The coupler has 10% in-coupling efficiency, and we did not put any effort into optimizing it beyond just defining the grating pitch to target the correct angle!