Photonic crystal cavity

Run this notebook in your browser using Binder.

In this notebook, we will simulate the commonly used L3 photonic crystal cavity composed of three missing holes in a hexagonal lattice of holes in a silicon slab.

# standard python imports
import numpy as np
import matplotlib.pyplot as plt

# tidy3D import
import tidy3d as td
from tidy3d import web

Coarse simulation

We will first run a broadband simulation to examine the spectrum, and zone in on the fundamental mode of the cavity. We start with defining some general parameters. We will use a fairly low spatial resolution for this initaly simulation. It’s worth remembering that the PML extend beyond the simulation domain, so we don’t need to worry about them covering some of the PhC holes. The one thing we have to remember is to extend the slab through the PML.

In structures with (quasi)-periodicity, that is to say with a well-defined notion of a unit cell, it is usually best to use a grid that is commensurate with the periodicity. This is why here we use a uniform grid in x and y, with a different step size to account for the different periodicity of the PhC lattice in these directions. In z, we use an automatic nonuniform mesh which conforms to the slab thickness and is finer in the silicon region.

# Number of PhC periods in x and y directions
Nx, Ny = 16, 12

# Lattice constant of the PhC in micron
alattice = 0.4

# Regular PhC lattice parameters
ra = 0.25 * alattice # hole radius
d_slab = 0.22        # slab thickness
n_slab = 3.48        # refractive index of the slab

# Materials - air and silicon
air = td.Medium()
si = td.Medium(permittivity=n_slab**2)

# Mesh step in x, y, z, in micron
steps_per_unit_length = 12
grid_spec = td.GridSpec(
    grid_x=td.UniformGrid(dl=alattice / steps_per_unit_length),
    grid_y=td.UniformGrid(dl=alattice / steps_per_unit_length * np.sqrt(3) / 2),
    grid_z=td.AutoGrid(min_steps_per_wvl=steps_per_unit_length)
)

# Central frequency around which we'll look for the cavity mode (Hz)
freq0 = 2e14

# Source bandwidth (Hz)
fwidth = 4e13

# Simulation run time (s)
run_time = 20/fwidth

# Simulation domain size (micron)
sim_size = [(Nx+2)*alattice, ((Ny+1)*alattice)*np.sqrt(3)/2, 4]

Next, we define the positions of the holes that make the photonic crystal structure.

# Define x and y positions in one quadrant of the simulation domain
xp, yp = [], []
nx, ny = Nx//2 + 1, Ny//2 + 1
for iy in range(ny):
    for ix in range(nx):
        xp.append(ix + (iy%2)*0.5)
        yp.append(iy*np.sqrt(3)/2)

# Remove the first two holes to make the L3 defect
xp = xp[2:]
yp = yp[2:]

# Append holes for the other three quadrants
xf, yf = [], []
for x, y in zip(xp, yp):
    xf += [x, x, -x]
    yf += [y, -y, y]
    if x > 0 and y > 0:
        xf += [-x]
        yf += [-y]

Initialize all structures.

slab = td.Structure(
    geometry=td.Box(center=[0, 0, 0], size=[td.inf, td.inf, d_slab]),
    medium=si
)

holes_geo = []
for x, y in zip(xf, yf):
    holes_geo.append(
        td.Cylinder(
            center = (np.array([x, y, 0])*alattice).tolist(),
            axis = 2,
            radius = ra,
            length = d_slab
        )
    )

holes = td.Structure(
    geometry=td.GeometryGroup(geometries=holes_geo),
    medium=air
)

Initialize the source. We are looking for the fundamental mode of the L3 cavity, so we use a y-polarized source at the center of the cavity.

source = td.PointDipole(
    center=(0, 0, 0),
    source_time=td.GaussianPulse(
        freq0=freq0,
        fwidth=fwidth),
    polarization='Ey')

source.source_time.plot(np.linspace(0, run_time, 1001))
plt.show()

<Figure size 432x288 with 1 Axes>

Finally, we also place a time monitor in the same location as the source. We set the time monitor starting time to be after the source decay, such that we can exclude the source signature from the recorded spectrum.

t_start = 1e-13
tmonitor = td.FieldTimeMonitor(center=[0, 0, 0], size=[0, 0, 0], start=t_start, name='field')

Initialize the simulation and visualize the structure. By default, Tidy3D will warn you if you have structures too close to the PML, as this can cause instability in the simulation. In photonic crystals this is sometimes inevitable, however, and it is OK in this case because the fields of the cavity mode are strongly localized around the center of the simulation domain.

# Suppress warnings for some of the holes being too close to the PML
td.config.logging_level = 'error'

sim = td.Simulation(
    size=sim_size,
    grid_spec=grid_spec,
    structures=[slab, holes],
    sources=[source],
    monitors=[tmonitor],
    run_time=run_time,
    boundary_spec=td.BoundarySpec.all_sides(boundary=td.PML())
)

fig, ax = plt.subplots(1, 2, figsize=(12, 4))
sim.plot_eps(z=0, ax=ax[0]);
sim.plot_eps(x=0, ax=ax[1]);

<AxesSubplot:title={'center':'cross section at x=0.00'}, xlabel='y', ylabel='z'>

<Figure size 864x288 with 4 Axes>

Run simulation and examine the spectrum

Now that the simulation is constructed, we can run it using the web API of Tidy3D. First, we submit the project.

# Submit a project to the cluster
job = web.Job(simulation=sim, task_name='L3 low res')
job.start()

And we can continuously monitor the status until the run is succsessful.

job.monitor()

Once the run is successful, we can download the results and load them in the td.SimulationData object.

sim_data = job.load(path='data/sim_data.hdf5')

We finally plot the time dependence of the field in the center of the cavity, and the spectrum computed using a Fourier transform of that field. For the latter, we use the in-built dft_spectrum function.

# Get data from the TimeMonitor
tdata = sim_data['field']

time_series = tdata.Ey.interp(x=0, y=0, z=0)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 3))

# Plot time dependence
time_series.plot(ax=ax1)
# ax[0].set_xlabel("Time [s]")
# ax[0].set_ylabel("Electric field [a.u.]");
# ax[0].set_title("Ey vs. time")

# Make frequency mesh and plot spectrum
dt = sim_data.simulation.dt
fmesh = np.linspace(1.6e14, 2.5e14, 101)

dft_matrix = np.exp(2j * np.pi * fmesh[:, None] * time_series.t.values) / np.sqrt(2 * np.pi)
spectrum = dt * dft_matrix @ np.real(time_series.values)

ax2.plot(fmesh, np.abs(spectrum))
ax2.set_xlim(1.7e14, 2.5e14)
ax2.set_xlabel("Frequency [Hz]")
ax2.set_ylabel("Electric field [a.u.]");
ax2.set_title("Spectrum");

Text(0.5, 1.0, 'Spectrum')

<Figure size 720x216 with 2 Axes>

We see a big peak close to f = 195THz, which is most likely what we are looking for, because a) the fundamental mode is the longest-lived and b) we use a y-polarized source at the center of the simulation domain, which does not excite some of the other modes. Next, we refine the simulation and compute the field profile of this fundamental mode.

Refine simulation, apply symmetries, get mode profile

Now that we’ve seen a clear resonant peak, we can increase the resolution of the simulation to obtain more accurate results, and to get a high-resolution image of the cavity mode. We center the source frequency close to the peak of the spectrum above, and decrease the bandwidth to exclude any other modes. We will also double the spatial resolution, making it 20 pixels per lattice period. Finally, we will also incorporate symmetries to speed up the computation.

# New target frequency based on spectrum above
freq0 = 1.95e14

# Narrow-bandwidth source
source = td.PointDipole(
    center=(0, 0, 0),
    source_time=td.GaussianPulse(
        freq0=freq0,
        fwidth=fwidth/5),
    polarization='Ey')

# Also increase the run time a bit
run_time = 50/fwidth

# 20 pixels per lattice period
steps_per_unit_length = 20
grid_spec = td.GridSpec(
    grid_x=td.UniformGrid(dl=alattice / steps_per_unit_length),
    grid_y=td.UniformGrid(dl=alattice / steps_per_unit_length * np.sqrt(3) / 2),
    grid_z=td.AutoGrid(min_steps_per_wvl=steps_per_unit_length)
)

We can use both a time and a frequency monitor to obtain the field profile, each coming with advantages and disadvantages. The frequency monitor captures accurately the frequency-domain field, but that includes the source signature. On the other hand, examining the time-domain field can capture the “eigenmode” of the system, but only if all the other modes have decayed. This is, to a very large extent, the case in our simulation, so as we’ll see the second approach works very well.

NB: An important thing to note is that a 2D time monitor can result in a very large amount of data. Because of this, we will only record the fields near the last time step, setting start = run_time in the FieldTimeMonitor.

# Time and frequency monitors
tmonitor = td.FieldTimeMonitor(
    center=[0, 0, 0],
    size=[4, 2*np.sqrt(3), 0],
    start=run_time,
    name='final_time')

fmonitor = td.FieldMonitor(
    center=[0, 0, 0],
    size=[4, 2*np.sqrt(3), 0],
    freqs=[freq0],
    name='field')

We initialize the simulation with reflection symmetries defined with respect to the x-, y-, and z-planes. Note that the eigenvalue of the symmetry (plus or minus one) has to be carefully determined, taking into account the vectorial nature of the electric field (and the pseudo-vector nature of the magnetic field). As an extra hint, positive symmetry is equivalent to a PMC plane, where the normal E-field component vanishes, while negative symmetry is equivalent to a PEC plane, where the parallel components of the E-field vanish. The symmetry values can be determined by thinking about a y-polarized electric dipole at the origin: (1, -1, 1).

# Initialize simulation
sim = td.Simulation(
    size=sim_size,
    grid_spec=grid_spec,
    structures=[slab, holes],
    sources=[source],
    monitors=[tmonitor, fmonitor],
    run_time=run_time,
    boundary_spec=td.BoundarySpec.all_sides(boundary=td.PML()),
    symmetry=(1, -1, 1),
)

fig, ax = plt.subplots(1, figsize=(8, 6))
sim.plot_eps(z=0, ax=ax);

<AxesSubplot:title={'center':'cross section at z=0.00'}, xlabel='x', ylabel='y'>

<Figure size 576x432 with 2 Axes>

We run the simulation as above.

sim_data = web.run(sim, task_name='L3 high res', path='data/sim_data.hdf5')

Finally, we plot the field recorded by the frequency monitor, with a rescaled colorbar in order to suppress the strongly dominant feature of the source in the center. On the other hand, the field stored in the time monitor reveals the eigenmode of the cavity.

fig, ax = plt.subplots(1, 2, figsize=(12, 4))
sim_data.plot_field('final_time', 'Ey', val='abs', z=0, time=run_time, ax=ax[0])
sim_data.plot_field('field', 'Ey', val='abs', z=0, freq=freq0, ax=ax[1]);

<AxesSubplot:title={'center':'cross section at z=0.00'}, xlabel='x', ylabel='y'>

<Figure size 864x288 with 4 Axes>