protis

Photonic band gaps in a crystal of dielectric rods

A periodic arrangement of dielectric materials does to light what the periodic potential of atoms in a semiconductor does to electrons: it organizes the allowed frequencies into bands, separated by band gaps, ranges of frequency where no light can propagate through the crystal, whatever its direction. Light at those frequencies is reflected from the crystal, and it can be trapped in a defect or guided along a line defect carved into it.

This notebook computes the band structure of a two-dimensional photonic crystal with protis, a Python implementation of the plane wave expansion method. It follows the square lattice of alumina rods of Joannopoulos et al., Photonic Crystals: Molding the Flow of Light (2008), chapter 5. We check the solver against an exact solution, compute the band diagram and its gap, look at the band surfaces over the whole Brillouin zone in 3D and at the fields of the modes at the gap edges, map how the gap depends on the rod radius, and study the convergence of the method. Finally, we switch to an automatic differentiation backend to compute how the gap responds to a change of permittivity anywhere in the unit cell, the first step of an inverse design.

Two polarizations

We consider a crystal invariant along $z$ and periodic in the $(x, y)$ plane, made of non-magnetic, lossless materials described by a permittivity $\varepsilon(x, y)$. For light propagating in the plane of periodicity, Maxwell’s equations split into two independent polarizations, each described by a single scalar field:

• TM modes have their electric field along the rods, $\mathbf{E} = E_z\,\hat{\mathbf{z}}$, and satisfy

$$-\nabla^2 E_z = k_0^2\, \varepsilon\, E_z,$$

• TE modes have their magnetic field along the rods, $\mathbf{H} = H_z\,\hat{\mathbf{z}}$, and satisfy

$$-\nabla\cdot\left(\frac{1}{\varepsilon}\nabla H_z\right) = k_0^2\, H_z,$$

where $k_0 = \omega / c$. Both are eigenvalue problems: for each wavevector, only a discrete set of frequencies $\omega_n(\mathbf{k})$ is allowed, the bands.

Bloch modes

The permittivity is periodic on a lattice with basis vectors $\mathbf{a}_1$ and $\mathbf{a}_2$. The Bloch theorem says the modes are plane waves modulated by a periodic function, which we expand in a Fourier series:

$$E_z(\mathbf{r}) = e^{i\mathbf{k}\cdot\mathbf{r}} \sum_{\mathbf{G}} e_{\mathbf{G}}\, e^{i\mathbf{G}\cdot\mathbf{r}}, \qquad \mathbf{G} = m\,\mathbf{b}_1 + n\,\mathbf{b}_2,$$

with $\mathbf{a}_i\cdot\mathbf{b}_j = 2\pi\delta_{ij}$. Bloch wavevectors that differ by a reciprocal lattice vector $\mathbf{G}$ describe the same modes, so $\mathbf{k}$ only needs to span one unit cell of the reciprocal lattice, the first Brillouin zone.

Maxwell’s equations have no natural length scale: we measure lengths in units of the lattice period $a$, and frequencies as $\omega a / 2\pi c = a / \lambda$.

import numpy as np
import matplotlib.pyplot as plt
import plotly.graph_objects as go
import plotly.io as pio
import protis as pt

plt.rcParams.update({"figure.dpi": 110, "axes.titlesize": "medium"})
pio.renderers.default = "plotly_mimetype"   # interactive figures, drawn by the page
pi = np.pi

The crystal

A square lattice of alumina rods ($\varepsilon = 8.9$) of radius $r = 0.2\,a$, in air. protis describes a unit cell by its permittivity sampled on a grid, here $2^9 \times 2^9$ points, from which it computes the Fourier coefficients. lattice.circle returns a boolean mask of the grid points inside the rod.

eps_rod = 8.9   # permittivity of the rods (alumina)
r = 0.2         # rod radius (units of a)

lattice = pt.Lattice([[1.0, 0], [0, 1.0]], discretization=2**9)
rod = lattice.circle(center=(0.5, 0.5), radius=r)
epsilon = lattice.ones()
epsilon[rod] = eps_rod

print(f"filling fraction: {rod.mean():.4f} (exact: pi r^2 = {pi * r**2:.4f})")

filling fraction: 0.1249 (exact: pi r^2 = 0.1257)

fig, ax = plt.subplots(figsize=(3.6, 3))
im = pt.plot2d(lattice, lattice.grid, epsilon.real, nper=3, ax=ax, cmap="Greys", vmin=0, vmax=10, show_cell=True, cellstyle="--r")
fig.colorbar(im[0], ax=ax, label=r"$\varepsilon$")
ax.set_xlabel("$x/a$")
ax.set_ylabel("$y/a$")
ax.set_title("permittivity, 3 × 3 unit cells")
plt.tight_layout()
plt.show()

The plane wave expansion

Substituting the Bloch expansion into the wave equation, and expanding the permittivity in a Fourier series $\varepsilon(\mathbf{r}) = \sum_{\mathbf{G}} \hat\varepsilon_{\mathbf{G}}\, e^{i\mathbf{G}\cdot\mathbf{r}}$, turns the product $\varepsilon E_z$ into a convolution of Fourier coefficients. For TM modes:

$$|\mathbf{k} + \mathbf{G}|^2\, e_{\mathbf{G}} = k_0^2 \sum_{\mathbf{G}'} \hat\varepsilon_{\mathbf{G}-\mathbf{G}'}\, e_{\mathbf{G}'}.$$

Keeping $N$ harmonics, this is a generalized eigenvalue problem $A\,\mathbf{e} = k_0^2\, B\,\mathbf{e}$, with $A$ diagonal and $B$ the Toeplitz matrix of the Fourier coefficients of $\varepsilon$. Both are Hermitian and $B$ is positive definite, so the frequencies are real and a single dense eigensolve returns all the bands at once. For TE modes, $A_{\mathbf{G}\mathbf{G}'} = (\mathbf{k}+\mathbf{G})\cdot(\mathbf{k}+\mathbf{G}')\,[\varepsilon^{-1}]_{\mathbf{G}\mathbf{G}'}$, where protis takes the inverse of the Toeplitz matrix of $\varepsilon$ (Ho, Chan and Soukoulis 1990), and $B$ is the identity.

pt.Simulation builds these matrices for nh harmonics (rounded to fill a disc in reciprocal space). We set the Bloch wavevector with sim.k and get the eigenvalues $k_0 = \omega / c$ sorted in increasing order from sim.solve. With $a = 1$, $k_0 / 2\pi$ is the normalized frequency $a/\lambda$.

sim = pt.Simulation(lattice, epsilon=epsilon, nh=100)
print(f"harmonics kept: {sim.nh}")

sim.k = (pi, 0)
freqs = sim.solve("TM", vectors=False) / (2 * pi)
print("first TM frequencies at k = (pi/a, 0):", np.round(freqs[:4].real, 5))

harmonics kept: 97
first TM frequencies at k = (pi/a, 0): [0.27518 0.44327 0.63842 0.77372]

Checking against a uniform medium

Before trusting the bands, we check the solver on a case with a known answer. In a uniform medium of index $n = \sqrt{\varepsilon}$, the modes are plane waves, $\omega = c\,|\mathbf{k}|/n$. Seen through the square lattice, a plane wave with wavevector $\mathbf{k} + \mathbf{G}$ is a Bloch mode with wavevector $\mathbf{k}$, so the bands are the light cone folded back into the Brillouin zone: $\omega_\mathbf{G}(\mathbf{k}) = c\,|\mathbf{k} + \mathbf{G}| / n$. The Fourier series of a constant has a single term, so the plane wave expansion is exact here, for both polarizations.

uniform = pt.Simulation(lattice, epsilon=lattice.ones() * 4.0, nh=100)
uniform.k = (0.3 * pi, 0.1 * pi)
exact = np.sort(np.hypot(uniform.kx, uniform.ky))[:6] / np.sqrt(4.0)

print("  folded light cone:", np.round(exact / (2 * pi), 6))
for polarization in ("TM", "TE"):
    k0 = uniform.solve(polarization, vectors=False)[:6]
    print(f"  {polarization}, max difference: {np.max(abs(k0 - exact)):.1e}")

  folded light cone: [0.079057 0.425735 0.480885 0.53033  0.575543 0.637377]
  TM, max difference: 4.0e-11
  TE, max difference: 4.0e-11

The two agree to about $10^{-10}$. The small difference is deliberate: protis shifts the squared eigenvalues by $10^{-12}\,(2\pi/a)^2$ before taking their square root, which keeps the derivatives finite at $\omega = 0$ when it computes gradients (see the end of this notebook).

The band diagram

A band diagram plots the frequencies along a path through the Brillouin zone. The square lattice is symmetric under the rotations and mirrors of the square, so the bands over the whole zone are determined by those in the irreducible Brillouin zone, the triangle with corners $\Gamma = (0, 0)$, $X = (\pi/a, 0)$ and $M = (\pi/a, \pi/a)$. We follow its edges, $\Gamma \to X \to M \to \Gamma$, with pt.init_bands, solve both polarizations at each point, and draw the result with pt.plot_bands.

Gamma, X, M = (0, 0), (pi, 0), (pi, pi)
sym_points = [Gamma, X, M, Gamma]
nk = 31   # points per segment of the path
kpath = pt.init_bands(sym_points, nk)


def bands(sim, kpath, polarization, nbands=8):
    out = []
    for k in kpath:
        sim.k = tuple(k)
        out.append(sim.solve(polarization, vectors=False)[:nbands].real / (2 * pi))
    return np.array(out)


TM = bands(sim, kpath, "TM")
TE = bands(sim, kpath, "TE")

gap_lo, gap_hi = TM[:, 0].max(), TM[:, 1].min()
ratio = 2 * (gap_hi - gap_lo) / (gap_hi + gap_lo)
print(f"TM gap between bands 1 and 2: {gap_lo:.4f} < a/lambda < {gap_hi:.4f}, gap-midgap ratio {100 * ratio:.1f} %")

fig = plt.figure(figsize=(6.6, 3.8))
lines_TE = pt.plot_bands(sym_points, nk, TE, color="#cf5268")
lines_TM = pt.plot_bands(sym_points, nk, TM, color="#4199b0", xtickslabels=[r"$\Gamma$", "$X$", "$M$", r"$\Gamma$"])
gap = plt.axhspan(gap_lo, gap_hi, color="#4199b0", alpha=0.15, lw=0)
plt.legend([lines_TM[0], lines_TE[0], gap], ["TM modes", "TE modes", "TM gap"], loc="upper left", bbox_to_anchor=(1.01, 1), frameon=False)
plt.ylim(0, 0.8)
plt.ylabel(r"frequency $\omega a / 2\pi c$")
plt.tight_layout()
plt.show()

TM gap between bands 1 and 2: 0.3229 < a/lambda < 0.4433, gap-midgap ratio 31.4 %

At low frequency, both polarizations follow straight lines from $\Gamma$: the wavelength is much larger than the period, and light sees a uniform effective medium. Higher up, the bands bend and split at the edges of the Brillouin zone, where the Bloch wave is Bragg reflected by the lattice. The TM bands open a gap between the first and second bands, from $a/\lambda = 0.323$ to $0.443$, a gap-midgap ratio of $31.4\,\%$, as in Joannopoulos et al. The TE bands have no gap: as we will see, TM gaps favour isolated regions of high permittivity, and TE gaps connected ones.

The bands over the whole Brillouin zone

We followed only the edges of the irreducible zone. A gap is a gap only if no band crosses it for any wavevector, so we check the first three TM bands over the whole Brillouin zone, $-\pi/a \le k_x, k_y \le \pi/a$. Drag to rotate, hover to read the frequencies.

nbz = 41
kxy = np.linspace(-pi, pi, nbz)
surfaces = np.zeros((nbz, nbz, 3))
for i, kx in enumerate(kxy):
    for j, ky in enumerate(kxy):
        sim.k = (kx, ky)
        surfaces[i, j] = sim.solve("TM", vectors=False)[:3].real / (2 * pi)

print(f"band 1: max {surfaces[..., 0].max():.4f} (on the path: {gap_lo:.4f})")
print(f"band 2: min {surfaces[..., 1].min():.4f} (on the path: {gap_hi:.4f})")

f32 = lambda a: np.asarray(a, dtype=np.float32)
kplot = f32(kxy / pi)
colors = ["#2b6f84", "#4199b0", "#8cc4d3"]
fig = go.Figure([
    # surfaces take z[i, j] at (x[j], y[i])
    go.Surface(x=kplot, y=kplot, z=f32(surfaces[..., n].T), colorscale=[[0, c], [1, c]], showscale=False,
               opacity=0.95, name=f"band {n + 1}",
               hovertemplate=f"band {n + 1}<br>a/λ = %{{z:.4f}}<extra></extra>")
    for n, c in enumerate(colors)
] + [
    go.Surface(x=kplot, y=kplot, z=f32(np.full((nbz, nbz), level)), colorscale=[[0, "#cf5268"], [1, "#cf5268"]],
               showscale=False, opacity=0.25, hoverinfo="skip")
    for level in (gap_lo, gap_hi)
])
fig.update_layout(
    title="TM bands 1 to 3 over the Brillouin zone, with the gap edges", height=480,
    margin=dict(l=0, r=0, t=40, b=0), template="simple_white",
    scene=dict(xaxis_title="kx a/π", yaxis_title="ky a/π", zaxis_title="a/λ",
               aspectratio=dict(x=1, y=1, z=0.8), camera=dict(eye=dict(x=1.6, y=-1.6, z=0.9))),
)
fig.show()

band 1: max 0.3229 (on the path: 0.3229)
band 2: min 0.4433 (on the path: 0.4433)

The first band peaks at the corners of the zone ($M$) and the second band dips lowest at the middle of its edges ($X$), so the extrema lie on the path we followed and the gap found from the band diagram is the true one. This is common, but not guaranteed: for less symmetric crystals, a band extremum can fall inside the irreducible zone, and a band diagram along its edges alone would overestimate the gap.

The modes at the gap edges

Why does the gap open? sim.solve with vectors=True (the default) also returns the Fourier coefficients of the modes, and sim.get_mode(n) rebuilds the periodic part of mode $n$ on the grid of one unit cell. Multiplying it by the Bloch phase $e^{i\mathbf{k}\cdot\mathbf{r}}$ gives the field $E_z$, which we draw over $2 \times 2$ cells: unlike the periodic part, the field itself changes from one cell to the next. We look at the first two TM modes at $X$, just below and just above the gap.

The frequencies obey a variational principle: the lowest mode minimizes

$$\frac{\omega^2}{c^2} = \frac{\int |\nabla E_z|^2\,\mathrm{d}\mathbf{r}}{\int \varepsilon\, |E_z|^2\,\mathrm{d}\mathbf{r}},$$

so it puts as much of its electric energy as possible where $\varepsilon$ is large while oscillating as slowly as possible. Each higher mode minimizes the same ratio while staying orthogonal to the modes below it. We measure where the energy goes with the fraction of electric energy density $\varepsilon |E_z|^2$ inside the rods.

sim.k = X
k0, _ = sim.solve("TM")

ncells = 2
x, y = lattice.grid
xs = np.concatenate([x[:, 0] + i for i in range(ncells)])   # grid over ncells x ncells cells
xx, yy = np.meshgrid(xs, xs, indexing="ij")
bloch_phase = np.exp(1j * (X[0] * xx + X[1] * yy))

fig, axes = plt.subplots(1, 2, figsize=(6.6, 3))
for n, (ax, name) in enumerate(zip(axes, ("band 1, below the gap", "band 2, above the gap"))):
    u = sim.get_mode(n)                             # periodic part, one cell
    energy = epsilon.real * abs(u) ** 2             # |E_z| = |u|
    in_rods = energy[rod].sum() / energy.sum()
    Ez = np.tile(u, (ncells, ncells)) * bloch_phase
    Ez = Ez / Ez.flat[np.argmax(abs(Ez))]           # fix the arbitrary phase: real field, maximum 1
    im = ax.pcolormesh(xs, xs, Ez.real.T, cmap="RdBu_r", vmin=-1, vmax=1, shading="nearest", rasterized=True)
    for i in range(ncells):
        for j in range(ncells):
            ax.add_patch(plt.Circle((i + 0.5, j + 0.5), r, fill=False, color="k", lw=0.8, ls=":"))
    ax.set_aspect("equal")
    ax.set_axis_off()
    ax.set_title(f"{name}\n$a/\\lambda = {k0[n].real / (2 * pi):.4f}$, {100 * in_rods:.0f} % in the rods")
    print(f"{name}: {100 * in_rods:.1f} % of the electric energy in the rods")
fig.colorbar(im, ax=axes, shrink=0.8, label="$E_z$")
plt.show()

band 1, below the gap: 83.5 % of the electric energy in the rods
band 2, above the gap: 32.6 % of the electric energy in the rods

At $X$, both modes oscillate in sign from one column of rods to the next, with the same period. The first band concentrates its field, and most of its electric energy, in the rods: it is the dielectric band. The second one must be orthogonal to it, so it moves its nodes onto the rods and puts most of its energy in the air: it is the air band. The large contrast between the two (here $\varepsilon = 8.9$ against $1$) translates into a large difference in frequency, the gap.

This picture also explains why TM gaps favour isolated rods: the field $E_z$ is parallel to the rod walls, so it can sit entirely inside a rod. TE modes have their electric field in the plane, and that field crosses the walls of isolated rods, which makes it hard to concentrate in them. Crystals with connected dielectric veins, such as a slab drilled with air holes, are better at opening TE gaps.

A gap map

How does the gap depend on the rod radius? We repeat the calculation for a range of radii, keeping only the edges of the irreducible zone and the first two TM bands, and trace the gap edges.

radii = np.linspace(0.05, 0.45, 41)
kpath_coarse = pt.init_bands(sym_points, 11)
edges = []
for radius in radii:
    eps_r = lattice.ones()
    eps_r[lattice.circle(center=(0.5, 0.5), radius=radius)] = eps_rod
    TM_r = bands(pt.Simulation(lattice, epsilon=eps_r, nh=100), kpath_coarse, "TM", nbands=2)
    edges.append((TM_r[:, 0].max(), TM_r[:, 1].min()))
lo, hi = np.array(edges).T
ratios = np.where(hi > lo, 2 * (hi - lo) / (hi + lo), 0)
best = np.argmax(ratios)
print(f"gap open for {radii[hi > lo].min():.2f} <= r/a <= {radii[hi > lo].max():.2f}")
print(f"largest gap-midgap ratio: {100 * ratios[best]:.1f} % at r = {radii[best]:.2f} a")

fig, ax = plt.subplots(figsize=(5.2, 3.4))
ax.plot(radii, lo, color="#4199b0", label="top of band 1")
ax.plot(radii, hi, color="#2b6f84", ls="--", label="bottom of band 2")
ax.fill_between(radii, lo, hi, where=hi > lo, color="#4199b0", alpha=0.25, lw=0, label="TM gap")
ax.axvline(r, color="0.6", lw=0.8)
ax.set_xlabel("rod radius $r/a$")
ax.set_ylabel(r"frequency $\omega a / 2\pi c$")
ax.set_xlim(radii[0], radii[-1])
ax.legend(frameon=False)
plt.tight_layout()
plt.show()

gap open for 0.09 <= r/a <= 0.42
largest gap-midgap ratio: 31.6 % at r = 0.21 a

Both edges move down as the rods grow: the modes have more high-permittivity material to sit in, which lowers their frequency. The gap opens once the rods are large enough to scatter light strongly, and closes again when they fill most of the cell, so that the air band, too, has to sit largely in the dielectric. The radius $r = 0.2\,a$ of Joannopoulos et al. sits close to the widest gap relative to its centre frequency.

How many plane waves?

The only approximation is the truncation to $N$ harmonics. We track three gap edges as $N$ grows: the TM bands 1 at $M$ and 2 at $X$, and, for comparison, the TE band 1 at $X$. We take a calculation with about $1000$ harmonics as the reference.

import time

edges_of_interest = {"TM band 1 at M": ("TM", M, 0), "TM band 2 at X": ("TM", X, 1), "TE band 1 at X": ("TE", X, 0)}


def edge_frequencies(nh):
    s = pt.Simulation(lattice, epsilon=epsilon, nh=nh)
    out = []
    for polarization, k, n in edges_of_interest.values():
        s.k = k
        out.append(s.solve(polarization, vectors=False)[n].real / (2 * pi))
    return s.nh, np.array(out)


nh_ref, reference = edge_frequencies(1000)
nhs = [25, 50, 100, 200, 400]
converged = []
for nh in nhs:
    t0 = time.perf_counter()
    n, values = edge_frequencies(nh)
    converged.append((n, values))
    print(f"{n:4d} harmonics: " + "  ".join(f"{v:.5f}" for v in values) + f"   ({time.perf_counter() - t0:.2f} s)")
print(f"{nh_ref:4d} harmonics: " + "  ".join(f"{v:.5f}" for v in reference) + "   (reference)")

fig, ax = plt.subplots(figsize=(5.2, 3.4))
n = np.array([c[0] for c in converged])
errors = abs(np.array([c[1] for c in converged]) - reference) / reference
for (label, (polarization, _, _)), err, marker in zip(edges_of_interest.items(), errors.T, "os^"):
    ax.loglog(n, err, marker + "-", color="#4199b0" if polarization == "TM" else "#cf5268", label=label)
ax.set_xticks(n, [str(v) for v in n], minor=False)
ax.minorticks_off()
ax.set_xlabel("number of harmonics $N$")
ax.set_ylabel("relative error")
ax.legend(frameon=False)
plt.tight_layout()
plt.show()

  21 harmonics: 0.32405  0.44660  0.38822   (0.01 s)
  49 harmonics: 0.32313  0.44352  0.40089   (0.01 s)
  97 harmonics: 0.32295  0.44327  0.40457   (0.01 s)
 197 harmonics: 0.32290  0.44315  0.40869   (0.02 s)
 385 harmonics: 0.32289  0.44311  0.41162   (0.10 s)
 997 harmonics: 0.32288  0.44310  0.41387   (reference)

The TM frequencies converge quickly: with $100$ harmonics, the gap edges are correct to a few $10^{-4}$, plenty for a band diagram. The TE frequency converges much more slowly: it is still $2\,\%$ off with $100$ harmonics, and even the reference moves by $0.5\,\%$ between $400$ and $1000$ harmonics, so its errors here are underestimated. Near the rod walls, the in-plane electric field has a component normal to the wall, which jumps with $\varepsilon$, and the truncated Fourier series of $\varepsilon$ and of its inverse represent such products poorly. This is the same Fourier factorization issue as in the Fourier modal method (Li 1996). The cost of each eigensolve grows as $N^3$, so it pays to check convergence on the quantity of interest rather than adding harmonics blindly.

Gradients by automatic differentiation

protis runs on the numerical backends of nannos: numpy by default, but also JAX or PyTorch, which differentiate through the whole calculation, eigensolver included. We switch to JAX and write the gap-midgap ratio as a function of the permittivity sampled on the grid. jax.grad then returns its derivative with respect to the permittivity of every one of the $2^9 \times 2^9$ grid points, in a single backward pass.

We use the frequency of band 1 at $M$ and of band 2 at $X$ as the gap edges, since we found above that they are the extrema.

import jax
import jax.numpy as jnp

pt.set_backend("jax")   # reloads protis, which now computes with JAX
jax.config.update("jax_enable_x64", True)


def gap_ratio(eps):
    s = pt.Simulation(lattice, epsilon=eps, nh=100)
    s.k = M
    lo = s.solve("TM", vectors=False)[0].real
    s.k = X
    hi = s.solve("TM", vectors=False)[1].real
    return 2 * (hi - lo) / (hi + lo)


eps0 = jnp.asarray(epsilon.real)
t0 = time.perf_counter()
value, sensitivity = jax.value_and_grad(gap_ratio)(eps0)
print(f"gap-midgap ratio: {100 * value:.2f} %, value and gradient in {time.perf_counter() - t0:.1f} s")

# Check: increasing the rod permittivity by delta changes every rod pixel at once.
delta = 1e-3
finite_difference = (gap_ratio(eps0 + delta * rod) - gap_ratio(eps0 - delta * rod)) / (2 * delta)
print(f"d(ratio)/d(eps_rod): automatic differentiation {sensitivity[rod].sum():.6f}, finite differences {finite_difference:.6f}")

gap-midgap ratio: 31.41 %, value and gradient in 4.9 s
d(ratio)/d(eps_rod): automatic differentiation 0.032310, finite differences 0.032310

density = np.asarray(sensitivity) * lattice.discretization[0] * lattice.discretization[1]   # per unit area of the cell
vmax = abs(density).max()

fig, ax = plt.subplots(figsize=(3.8, 3))
im = pt.plot2d(lattice, lattice.grid, density, nper=1, ax=ax, cmap="RdBu_r", vmin=-vmax, vmax=vmax)
ax.add_patch(plt.Circle((0.5, 0.5), r, fill=False, color="k", lw=0.8, ls=":"))
fig.colorbar(im[0], ax=ax, label=r"$\partial (\Delta\omega / \omega_0) / \partial\varepsilon(\mathbf{r})$")
ax.set_axis_off()
ax.set_title("where adding permittivity widens the gap")
plt.tight_layout()
plt.show()

The map reads like the fields above. Adding permittivity where a mode is strong lowers its frequency (to first order, $\delta\omega/\omega = -\tfrac12 \int \delta\varepsilon\, |E|^2 / \int \varepsilon\, |E|^2$). In the rod, where the dielectric band at $M$ is strong, more permittivity lowers the bottom edge and widens the gap. Between the columns of rods, where the air band at $X$ peaks, it lowers the top edge and closes the gap. The map is not symmetric under a quarter turn because the air band at $X$ is not: its partner at $(0, \pi/a)$ has the same frequency, and a complete treatment would include both.

Such a sensitivity map is the starting point of topology optimization: a gradient-based optimizer redistributes the material, step by step, to maximize the gap, with each step costing one forward and one backward pass, whatever the number of grid points.

Things to try

Every cell above is editable. A few ideas:

• Air holes in a dielectric: swap the materials, epsilon = lattice.ones() * 13; epsilon[rod] = 1, with a larger radius. Do the TE modes now open a gap?
• Triangular lattice: use pt.Lattice([[1, 0], [0.5, np.sqrt(3) / 2]], discretization=2**9), put the hole at the centre of the cell, and follow the path $\Gamma \to M \to K \to \Gamma$ with $M = (0, 2\pi/\sqrt3)$ and $K = (2\pi/3, 2\pi/\sqrt3)$. Holes of radius $0.48\,a$ in $\varepsilon = 13$ give a gap for both polarizations at once (Joannopoulos et al., chapter 5).
• Contrast: lower eps_rod and watch the gap shrink. Below which permittivity does it close?
• Optimization: add the gradient to the permittivity in small steps, clipped between $1$ and $8.9$, and watch the gap-midgap ratio grow. Include the mode at $(0, \pi/a)$ in gap_ratio, so the result keeps the symmetry of the square.