Guided mode expansion of a PhC slab with a horizontal (xy) symmetry plane

In this example we apply legume to calculate the bands of a PhC slab with a horizontal (xy) mirror plane, and we demonstrate how the bands can be separated according to reflection symmetry \(\sigma_{xy}\). This example is related to Sec. 2.3 of the CPC paper.

The structure parameters are those of Ch. 8, Fig. 2(b) of Molding the Flow of Light , and we focus on how to achieve separation with respect to \(\sigma_{xy}\) by choosing the indices of the guided modes in the expansion.

import legume
import numpy as np
import matplotlib.pyplot as plt
import matplotlib

Define an xy-symmetric photonic crystal slab

Here we consider an xy-symmetric PhC slab, i.e., a slab that is symmetric with respect to mirror reflection in the xy plane. We separate the two polarizations as (see Sec. 2.3, Table 1 of the CPC paper):

• \(xy\)-even or \(\sigma_{xy}=+1\) (or TE-like): we must use guided modes \(\alpha=0, 3, 4, 7, ...\) in the basis

• \(xy\)-odd or \(\sigma_{xy}=-1\) (or TM-like): we must use guided modes \(\alpha=1, 2, 5, 6, ...\) in the basis

In PhC literature, \(xy\)-even modes are often called TE-like, while \(xy\)-odd modes are often called TM-like, where the terminology is borrewed from the 2D case. We prefer to avoid this terminology when using legume, because we reserve the words TE or TM for the polarizations of the guided modes of the effective waveguide, which are defined with respect to a vertical mirror plane.

D = 0.6         # slab thickness in units of a
r = 0.3        # hole radius in units of a
eps_c = 1.0     # permittivity of circular hole
eps_b = 12.0    # background permittivity of slab
eps_lower, eps_upper = 1, 1  # permittivities of lower and upper claddings

lattice = legume.Lattice('hexagonal')
phc = legume.PhotCryst(lattice, eps_l=eps_lower, eps_u=eps_upper)
phc.add_layer(d=D, eps_b=eps_b)
phc.layers[-1].add_shape(legume.Circle(eps=eps_c, r=r, x_cent=0., y_cent=0))

gme = legume.GuidedModeExp(phc, gmax=5, truncate_g='abs')
npw = np.shape(gme.gvec)[1]  # number of plane waves in the expansion
print(f'Number of reciprocal lattice vectors in the expansion: npw = {npw}')

Number of reciprocal lattice vectors in the expansion: npw = 61

Separately calculate xy-even and xy-odd modes

Here we run the guided-mode expansion for the two separate symmetries and store the bands.

path = lattice.bz_path(['G', 'M', 'K', 'G'], [30, 20, 30])
numeig, verbose = 20, True

gme.run(kpoints=path['kpoints'], gmode_inds=[0, 3], numeig=numeig, verbose=True, compute_im=False)
freqs_xyeven = gme.freqs
nkappa, nfreq = freqs_xyeven.shape[0], freqs_xyeven.shape[1]
print(f'\n Number of wavevectors = {nkappa}, number of frequencies = {nfreq}\n')

gme.run(kpoints=path['kpoints'], gmode_inds=[1, 2], numeig=numeig, verbose=True, compute_im=False)
freqs_xyodd = gme.freqs

# X0 is the 1D array with wavevectors, X is the 2D array used for the scatter plot
(X0,X) = legume.viz.calculate_x(path["kpoints"],numeig,k_units=True) # Create an array for bands plotting

Running gme k-points: │------------------------------│ 1 of 81Running gme k-points: │██████████████████████████████│ 81 of 81

┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓
┃ Steps in GuidedModeExp: 61 plane waves and 2 guided modes ┃ Time (s) ┃                 % vs total T ┃
┡━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┩
│ Guided modes computation with gmode_compute='exact'       │ 2.398    │ │████████████████----│   82% │
│ Inverse matrix of Fourier-space permittivity              │ 0.002    │ │--------------------│    0% │
│ Matrix diagionalization using the 'eigh' solver           │ 0.086    │ │--------------------│    3% │
│ Creating GME matrix                                       │ 0.393    │ │██------------------│   13% │
├───────────────────────────────────────────────────────────┼──────────┼──────────────────────────────┤
│ Total time for real part of frequencies for 81 k-points   │ 2.914    │ │████████████████████│  100% │
└───────────────────────────────────────────────────────────┴──────────┴──────────────────────────────┘

Skipping imaginary part computation, use run_im() to run it, or compute_rad() to compute the radiative rates of
selected eigenmodes

 Number of wavevectors = 81, number of frequencies = 20

Running gme k-points: │██████████████████████████████│ 81 of 81

┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓
┃ Steps in GuidedModeExp: 61 plane waves and 2 guided modes ┃ Time (s) ┃                 % vs total T ┃
┡━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┩
│ Guided modes computation with gmode_compute='exact'       │ 2.363    │ │████████████████----│   82% │
│ Inverse matrix of Fourier-space permittivity              │ 0.000    │ │--------------------│    0% │
│ Matrix diagionalization using the 'eigh' solver           │ 0.086    │ │--------------------│    3% │
│ Creating GME matrix                                       │ 0.395    │ │██------------------│   14% │
├───────────────────────────────────────────────────────────┼──────────┼──────────────────────────────┤
│ Total time for real part of frequencies for 81 k-points   │ 2.890    │ │████████████████████│  100% │
└───────────────────────────────────────────────────────────┴──────────┴──────────────────────────────┘

Skipping imaginary part computation, use run_im() to run it, or compute_rad() to compute the radiative rates of
selected eigenmodes

Plot the results and comparison

Calculate the cladding light line for the plot (the one corresponding to the higher permittivity, but here the two claddings are identical)

eps_clad = [gme.phc.claddings[0].eps_avg, gme.phc.claddings[-1].eps_avg]
vec_LL = np.sqrt(
         np.square(gme.kpoints[0, :]) +
         np.square(gme.kpoints[1, :])) / 2 / np.pi / np.sqrt(max(eps_clad))

We make the plots explicitly, without using legume.viz.bands, by making a scatter plot of frquency vs wavevector. This requires extracting a 2D array, with the same shape of gme.r in which the wavevector is normalized from 0 to 1.

markersize = 10
conecolor='lightgrey'
odd_color = "#381a61"
even_color = "#E78429"

fig, ax = plt.subplots(1, constrained_layout=True, figsize=(5, 4))

plt.scatter(X, freqs_xyeven,c=even_color, label='xy-even or TE-like', s=markersize)
plt.scatter(X, freqs_xyodd, c=odd_color, label='xy-odd or TM-like', s=markersize)

ax.set_ylim([0, 0.8])
ax.set_xlim([0, 1])
ax.set_xlabel("Wavevector")
ax.set_ylabel("$\omega a/2\pi c$")
plt.legend(loc='lower center')

ax.fill_between(X0, vec_LL,  max(100, vec_LL.max(), gme.freqs[:].max()),
                facecolor=conecolor,        zorder=0)
ax.set_xticks(path['k_indexes'], path['labels'])
ax.xaxis.grid('True')

The second figure is taken from Ch.8 of Molding the Flow of Light. The GME bands are a bit complicated, aren’t they? This is because we are looking at the modes both below and above the light cone. Moreover, the slab is relatively thick and there are second-order guided modes below the light cone.

To better understand the second point, we compare the band dispersion in the PhC with those of the effective waveguide.

# Run the plane-wave expansion for the two separate polarizations and store the bands
only_gmodes = True  # this keyword parameter select to calculate the guided modes of the effective waveguids
verbose = False

gme = legume.GuidedModeExp(phc, gmax=5, truncate_g='abs') # we have to re-iniitialize the gme object to avoid errors
gme.run(kpoints=path['kpoints'], gmode_inds=[0, 3], numeig=numeig, verbose=verbose, compute_im=False, only_gmodes=only_gmodes)
freqs_xyeven_gmodes = gme.freqs
nkappa, nfreq = freqs_xyeven_gmodes.shape[0], freqs_xyeven_gmodes.shape[1]
print(f'\n Number of wavevectors = {nkappa}, number of frequencies = {nfreq}\n')

gme = legume.GuidedModeExp(phc, gmax=5, truncate_g='abs')
gme.run(kpoints=path['kpoints'], gmode_inds=[1, 2], numeig=numeig, verbose=verbose, compute_im=False, only_gmodes=only_gmodes)
freqs_xyodd_gmodes = gme.freqs
(X0,X_gmodes) = legume.viz.calculate_x(path["kpoints"],numeig,k_units=True) # Create an array for bands plotting

 Number of wavevectors = 81, number of frequencies = 20

ymin, ymax = 0, 0.8
markersize = 4


fig, ax = plt.subplots(nrows=2, ncols=2, figsize = (9, 8))

ax[0,0].scatter(X, freqs_xyeven, c=even_color, label='xy-even',  s=markersize)
ax[1,0].scatter(X, freqs_xyodd , c=odd_color, label='xy-odd' , s=markersize)
ax[0,1].scatter(X_gmodes, freqs_xyeven_gmodes, c=even_color, label='xy-even, eff. waveguide',  s=markersize)
ax[1,1].scatter(X_gmodes, freqs_xyodd_gmodes , c=odd_color, label='xy-odd, eff. waveguide' ,  s=markersize)

for a in ax.flatten():
    a.plot(X0, vec_LL, color='grey')
    a.set_xlim([0, 1])
    a.set_ylim([ymin, ymax])
    a.legend()
    a.set_xticks(path['k_indexes'], path['labels'])
    a.xaxis.grid('True')
ax[0,0].set_ylabel("$\omega a/2\pi c$")
ax[1,0].set_ylabel("$\omega a/2\pi c$")
plt.legend(loc='lower center')

<matplotlib.legend.Legend at 0x7c63d669e310>

The final message is: to see whether a photonic bandgap occurs in a PhC slab, we should look at modes below and above the light cone. And we should keep a sufficient number of guided modes in the basis, to be sure to have all photonic bands in the relevant frequency range.

As a further example, we could calculate the photonic dispersion adding the guided modes one by one in the array gmode_inds, and look at the effect of guided modes of first-order, second-order, …

(It is not easy to do fancy graphics with PhC slab bands, they cannot be plotted with lines because of broken curves. If you want better plots, try increasing the number of points - wait patiently - then play with markersize.)