nannos

Guided resonances in a photonic crystal slab

Drill a square lattice of holes through a thin slab of silicon and it becomes two things at once: a waveguide, which traps light inside the slab by total internal reflection, and a grating, which couples light between the slab and free space. Where both roles meet, light shone on the slab from outside can couple to modes guided inside it, for a narrow band of frequencies. These guided resonances show up as sharp lines in the transmission spectrum, where the slab switches from transparent to perfectly reflecting within a fraction of a percent in frequency.

This notebook computes the transmission spectrum of such a slab with nannos, a Python implementation of the Fourier modal method (also known as rigorous coupled wave analysis, RCWA). It follows the structure studied by Fan and Joannopoulos (2002), checks the solver against an exact solution, looks at the fields of a resonance in 3D, and studies how the result converges with the number of Fourier harmonics.

Periodic structures and Bloch waves

We work with time-harmonic fields, $\mathbf{E}(\mathbf{r}, t) = \mathrm{Re}\left[\mathbf{E}(\mathbf{r})\, e^{-i\omega t}\right]$, which is the convention used by nannos (losses therefore appear as a positive imaginary part of the permittivity). Without sources, Maxwell’s equations read

The permittivity $\varepsilon(x, y)$ is periodic in the plane of the slab, on a lattice with basis vectors $\mathbf{a}_1$ and $\mathbf{a}_2$. The Bloch theorem then says that, for an incident plane wave with in-plane wavevector $\mathbf{k}_\parallel$, the field is a sum of plane waves whose in-plane wavevectors differ by vectors $\mathbf{G}$ of the reciprocal lattice:

with $\mathbf{a}_i\cdot\mathbf{b}_j = 2\pi\delta_{ij}$. Each term $(m, n)$ is a diffraction order. Above and below the slab, an order carries energy away only if $|\mathbf{k}_\parallel + \mathbf{G}| < k_0$, where $k_0 = \omega/c = 2\pi/\lambda$; otherwise it decays exponentially away from the slab.

Maxwell’s equations have no natural length scale, so we measure lengths in units of the lattice period $a$ and frequencies as $a/\lambda = \omega a / 2\pi c$. For a square lattice at normal incidence ($\mathbf{k}_\parallel = 0$), the first orders $(\pm 1, 0)$ and $(0, \pm 1)$ start to propagate in air at $a/\lambda = 1$. Below that, only the specular order $(0, 0)$ leaves the slab: everything that is not transmitted straight through is reflected straight back.

The structure

The slab has a permittivity $\varepsilon = 12$ (close to silicon or GaAs in the near infrared), a thickness $h = 0.5\,a$, and holes of radius $r = 0.2\,a$. It is suspended in air, so the media above (superstrate) and below (substrate) are vacuum.

nannos describes the structure as a stack of layers, each invariant along $z$, sharing one lattice. The superstrate and substrate are semi-infinite and uniform. The patterned layer stores its permittivity on a grid over one unit cell, here $2^9 \times 2^9$ points, from which nannos computes the Fourier coefficients. We draw the hole once with shapely: lattice.geometry_mask finds the grid points inside it, and the 3D views below reuse it.

A 3D view

The unit cell repeats in both directions. nannos draws structures in 3D with pyvista (sim.plot_structure), which needs a running kernel, so for this page we build the solid ourselves with plotly. We cut the holes out of a $3\times3$ block of cells and extrude the resulting polygon over the thickness of the slab: a constrained Delaunay triangulation of the polygon gives the top and bottom faces, and each edge of its outline gives a rectangle of the side walls. Drag to rotate, scroll to zoom.

We place the top surface of the slab at height $0$, with light arriving from above.

Layers, Fourier series and the eigenproblem

Inside each layer, nothing depends on $z$, so we expand the permittivity in a Fourier series over the unit cell,

and substitute it, with the Bloch expansion of the fields, into Maxwell’s equations. A product such as $\varepsilon\mathbf{E}$ becomes a convolution of Fourier coefficients, that is a matrix (a Toeplitz matrix built from the $\hat\varepsilon_{\mathbf{G}}$) acting on the vector of field coefficients. The $z$ dependence then decouples: each layer supports a set of modes, solutions of an algebraic eigenproblem of size $2N$ for $N$ harmonics, which vary as $e^{\pm i\gamma z}$ along the stack. In the uniform superstrate and substrate, these modes are just the plane waves of the diffraction orders.

The layers are then joined by requiring the tangential components of $\mathbf{E}$ and $\mathbf{H}$ to be continuous at each interface. nannos does this with scattering matrices, which stay numerically stable however thick the layers are, and returns the amplitudes of every order in reflection and transmission.

The only approximation is the truncation to a finite set of $N$ harmonics. Here is what the permittivity looks like when rebuilt from the harmonics kept with nh=100, as returned by sim.get_epsilon:

With a few dozen harmonics per direction, the steep edges of the hole are smoothed and ring (the Gibbs phenomenon). This sets how fast the method converges, and we will come back to it at the end.

Transmission at one frequency

The excitation is a plane wave of unit amplitude. Its direction is given by the polar and azimuthal angles $(\theta, \varphi)$ and its polarization by the angle $\psi$: at normal incidence, $\psi = 0$ gives an electric field along $x$. For the square lattice and circular holes the two polarizations are equivalent at normal incidence.

diffraction_efficiencies returns the reflected and transmitted powers $R$ and $T$, summed over the propagating orders and normalized by the incident power. The slab is lossless, so energy conservation requires $R + T = 1$.

Checking against a uniform slab

Before trusting the spectrum, we check the solver on a case with a known answer: the same slab without holes. A uniform layer only couples each order to itself, so the Fourier modal method is exact here, and the transmission must match the Fabry–Pérot (Airy) formula for a slab of index $n = \sqrt{\varepsilon}$ and thickness $h$ at normal incidence,

Multiple reflections between the two faces interfere, and the slab becomes fully transparent each time its optical thickness $n h$ is a multiple of half a wavelength.

The two agree to machine precision.

The transmission spectrum

Now the patterned slab. We sweep the frequency between $a/\lambda = 0.25$ and $0.6$. The resonances are very narrow, so a uniform grid would miss most of them: nn.adaptive_sampler starts from a coarse grid and adds points where the curve bends sharply.

As a guide, we also plot the Airy formula for a uniform slab with the average permittivity of the patterned one, $\bar\varepsilon = (1 - f_\text{air})\,\varepsilon + f_\text{air}$.

The spectrum has two parts. The smooth background follows the Fabry–Pérot oscillation of a uniform slab with the average permittivity: at these wavelengths, the light crossing the slab mostly sees an effective medium. On top of it sit sharp, asymmetric lines, the guided resonances. Each one corresponds to a mode guided in the slab whose wavevector, folded back by the lattice, lands at $\mathbf{k}_\parallel = 0$, so that the incident wave can excite it. The mode leaks back out through the same grating, and its re-radiated field interferes with the direct transmission through the slab. The result is a Fano line shape, where the transmission drops to exactly zero and climbs back to one within a tiny range of frequency. This spectrum reproduces Fig. 12(c) of Fan and Joannopoulos.

The fields at a guided resonance

We zoom on the first sharp line near $a/\lambda = 0.382$ and locate the frequency of total reflection with a bounded scalar minimization of $T$.

get_field_grid rebuilds the fields in real space from their Fourier coefficients, in a given layer and at given heights $z$ measured from the top of that layer. We compare the electric field intensity $|\mathbf{E}|^2$ in the middle plane of the slab at the resonance and at $a/\lambda = 0.45$, away from any resonance. The incident wave has unit amplitude, so $|\mathbf{E}|^2$ directly measures the enhancement.

At the resonance the peak intensity is more than a hundred times larger than away from it: the guided mode stores energy, and the narrower the resonance, the larger the build-up. The brightest spots sit inside the hole, against its edges along $x$, the direction of the incident electric field. There the field is normal to the edge, and since $\varepsilon E_n$ is continuous, the field on the air side is $12$ times stronger than on the silicon side.

The fields in 3D

To see how the mode sits in the slab, we combine two cuts in the 3D view: the middle plane of the slab, repeated over the $3\times3$ cells, and a vertical cut through the back row of hole centres that extends one period above and below the slab. The slab itself is drawn transparent. The colour shows the field amplitude $|\mathbf{E}|$ rather than the intensity, so that the much weaker field outside the slab remains visible. Hover over a cut to read its value.

The vertical cut shows the signature of a guided mode: the field is confined to the slab and decays evanescently on both sides. Above the slab, the incident and fully reflected waves also form a standing wave, with an amplitude reaching $2$ at its antinodes, weak next to the field inside the slab. Below the slab the field dies out within about half a period: nothing is transmitted.

How many harmonics?

We used nh=100 harmonics so far. How close is the result to the converged one? We track the transmission at $a/\lambda = 0.45$ as the number of harmonics grows, for two formulations of the method.

The truncation matters most where $\varepsilon$ jumps. At the edge of the hole, the tangential component of $\mathbf{E}$ is continuous, while for the normal component it is $\varepsilon E_n$ that is continuous. The two products $\varepsilon\mathbf{E}$ must therefore be truncated differently: with the Toeplitz matrix of $\varepsilon$ for the tangential component (Laurent’s rule) and with the inverse of the Toeplitz matrix of $1/\varepsilon$ for the normal one (the inverse rule, Li 1996). The original formulation applies Laurent’s rule everywhere. The tangent formulation builds a vector field tangent to the hole edges and applies the appropriate rule to each component.

Here the original formulation settles first: from about 100 harmonics on, its transmission changes only in the fourth decimal. The tangent formulation approaches the same value from above, more slowly, and still differs by about 1 % at 300 harmonics. Which formulation converges faster depends on the structure, the polarization and the quantity of interest, so comparing them, as here, is a cheap way to estimate the error. The cost grows quickly, since the eigenproblem scales with the cube of the number of harmonics.

Sharp resonances are more demanding than the smooth background: with the original formulation, the frequency of total reflection found above moves from $0.3818$ to $0.3820$ between 100 and 200 harmonics. That is small compared with the spacing between resonances, but comparable with the width of the narrowest ones.

Things to try

Every cell above is editable. A few ideas:

• Oblique incidence: set angles=(10, 0, 0) in nn.PlaneWave. The resonances shift and split, and modes that a normally incident wave cannot excite by symmetry appear as new lines.
• Hole size: change r and watch the resonances move. Smaller holes give a weaker grating, hence narrower resonances, and the spectrum tends to the Airy curve.
• Losses: give the slab a small absorption, eps_slab = 12 + 0.05j (positive imaginary part with the $e^{-i\omega t}$ convention). The sharpest resonances are the most affected, since light stays there longest.
• Shape: replace the circle by an ellipse, hole = shapely.affinity.scale(hole, 1.5, 0.75), and compare the two polarizations, $\psi = 0$ and $\psi = 90$.
• Higher frequencies: extend the spectrum beyond $a/\lambda = 1$, where diffracted orders start to propagate, and print sim.diffraction_efficiencies(orders=True).