home

Rigorous coupled wave analysis for displacement-sentitive photonic crystals

How to use my RCWA library to simulate displacement-sensitive photonic crystal slabs

All matters discussed in this tutorial post are extracted from this article written by Wonjoo Suh.

The structure consists in two stacked layers of dielectric constant $\varepsilon=12$ separated by an air gap. The layers are \(0.55 a\) deep each and the air gap depth varies from \(0.55\) to \(1.35 a\) with \(a\) the lattice constant.

The bilayer system with an interlayer air gap h

Our first objective is to reproduce Figure 2. We compute transmission spectra for different interlayer depths.

There are several ways you can simulate the structures discussed in the article but for maximum convenience, I suggest the Crystal interface. We start by importing relevant modules. 

from khepri.crystal import Crystal
from khepri.expansion import Expansion
import numpy as np
from tqdm import tqdm

We define the parameters of the structures as described above, including the range of frequencies over which we want to compute quantities. All values are without units, dimensions are expressed as fractions of the crystal pitch. 

er_layer  = 12    # Dielectric of the holey film
er_interlayer = 1 # Dielectric of the layer spacing 
r   = 0.4         # Radius of the holes
layer_thickness   = 0.55 # Thickness of holey layer
# We try different interlayer spacings
interlayer_thicknesses = [1.35, 1.1, 0.95, 0.85, 0.75, 0.65, 0.55]
freqsn = np.linspace(0.49, 0.59, 100) # The frequency range

We first draw our unit cell structure in a numpy array. The parameter N has no impact on the performance but it has to be high enough to get good Fourier coefficients in the fast Fourier transform algorithm. 

N = 256
eps = np.ones((N, N)) * er_layer # Uniform dielectric
x = np.linspace(-0.5, 0.5, N)
y = np.linspace(-0.5, 0.5, N)
X, Y = np.meshgrid(x, y)
eps[np.sqrt(X**2+Y**2)< r] = 1.0 # Drilling the hole.

Alternatively, one can use the integrated draw module. 

from khepri.draw import Drawing
N = 256
d = Drawing((N,N), er_layer)
d.circle((0,0), r, 1.0) # Drilling the hole
eps = d.canvas()

Now we devise a function to compute the scattering matrix of the system. 

def create_crystal(interlayer_thickness, pw=(9,9)):
  cl = Crystal(pw)
  # The two membranes described by 'pixmap' eps
  cl.add_layer_pixmap("holey_layer", eps, layer_thickness)
  # The air gap between the membranes, uses analytical solver
  cl.add_layer_uniform("interlayer", er_interlayer, interlayer_thickness)
  cl.set_device(["holey_layer", "interlayer", "holey_layer"])
  return cl

Then we do a loop and compute all scattering matrice and their associeted transmissions. 

transmission_map = list()
for interlayer_thickness in interlayer_thicknesses:
  cl = create_crystal()
  for frequency in tqdm(freqsn):
    wavelength = 1 / frequency
    cl = create_crystal(interlayer_thickness)
    cl.set_source(wavelength, te=1, tm=1, theta=0, phi=0)
    cl.solve()
      R, T = cl.poynting_flux_end()
    transmission_map.append(T)
    
transmission_map = np.asarray(transmission_map).reshape(len(interlayer_thicknesses), len(freqsn))

Plotting this transmission_map using matplotlib gives the figure below. 

import matplotlib.pyplot as plt
fig, ax = plt.subplots()
for i, (t, h) in enumerate(zip(transmission_map, interlayer_thicknesses)):
  ax.plot(freqsn, t, label=f"{h=}")
ax.set_xlabel("Frequency [c/a]")
ax.set_ylabel("Transmission")
plt.savefig("tuto.png")

spectrum

home