Band structure calculation of a photonic crystal slab¶

In this notebook, we simulate a photonic crystal slab consisting of a square lattice of air holes in a dielectric slab. Our goal is to compute the band structure of this photonic crystal slab, as found in:

Shanhui Fan and J. D. Joannopoulos, "Analysis of guided resonances in photonic crystal slabs," Phys. Rev. B 65, 235112 (2002).

To this end, we excite the structure with several PointDipole sources, and we measure the response with several FieldTimeMonitor monitors. We excite modes with a fixed Bloch wavevector by using Bloch boundary conditions. We then use the ResonanceFinder to find the resonant frequencies. By sweeping the Bloch wavevector, we obtain the full band structure of the photonic crystal slab.

See also the api reference for ResonanceFinder here.

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

We will randomly position our sources and monitors. Here, we seed the random number generator to guarantee reproducible results.

Now we define the parameters for the simulation, structure, sources, and monitors.

We take the dipole polarization to be Hz and the symmetry to be (0,0,1) in order to excite only modes which are even with respect to the xy mirror plane.

We define the materials and structures in terms of the above parameters. We take the photonic crystal slab to be lying in the xy plane. Because the photonic crystal slab is periodic in the x and y directions, we only need to simulate a single unit cell, containing the dielectric slab and a single air hole.

We will excite the photonic crystal slab with several PointDipole sources. Each dipole will have a random position and phase.

We create FieldTimeMonitors to record the field as a function of time at several random locations within the photonic crystal slab. Crucially, we start the monitors after the source pulse has decayed.

We will perform 3*Nk different simulations, each with different Bloch boundary conditions, as we sweep the Bloch wavevector over the boundary of the irreducible Brillouin zone. We sweep over three lines, namely $\Gamma X$, $XM$, and $M\Gamma$. We use a PML in the z direction.

Here, we simply define all of the boundary conditions we will use and put them into a single array.

Now we define the simulations we want to run.

Let's check that the structure and source look correct. The source spectrum must fill the entire frequency range of interest.

Now we run the simulations as a Batch.

We set verbose=True to keep track of the status of the jobs in the Batch.

Now that the simulations are complete, we can analyze the data. Let's first look at one of the FieldTimeMonitors to make sure the source has decayed.

We see that the source has mostly decayed by the time we switch on the monitors, and the remaining data shows decay and oscillation due to the resonances inside the system.

Looking at the Fourier transform of this data, we can see resonances at the band frequencies.

We use the ResonanceFinder plugin to find the band frequencies.

We first construct a ResonanceFinder object storing our parameters, and then call run() on our list of FieldTimeData objects. This will add up the signals from all of the monitors before searching for resonances. The ResonanceFinder class has additional methods in case the signal takes another form; see the api reference here.

The run() method returns an xr.Dataset containing the decay rate, Q factor, amplitude, phase, and estimation error for each resonance as a function of frequency.

We see the four resonances from the previous figure. All four have reasonable Q factors, amplitudes, and errors, so they are likely to represent physical resonances. Note that in order to accurately obtain the Q factor for high-Q modes, it may be necessary to run the simulation for a longer time.

Now we are ready to compute the band structure. We run the resonance finder at every Bloch wavevector.

We define a function to filter resonances based on their Q, amplitude, and error.

We plot the band structure with the light line overlaid.

The bandstructure we obtained matches the expected result from the paper. If we were seeing too many resonances, we could change the parameters to our filter_resonances function to eliminate the spurious ones. If we were seeing too few resonances even before filtering, we might have to change the parameters of the ResonanceFinder, for example decreasing rcond or increasing init_num_freqs. If the ResonanceFinder takes too long, we can decrease init_num_freqs. There can also be resonances on the light line associated with Wood's anomaly; we filter those out here based on their small amplitude.