Transmittance Spectra for Planewave at Normal Incidence of a Binary Grating

The mode-decomposition feature can also be applied to planewaves in homogeneous media with scalar permittivity/permeability (i.e., no anisotropy). This will be demonstrated in this example to compute the diffraction spectrum of a binary phase grating. To compute the diffraction spectrum for a finite-length structure, see Tutorials/Near to Far Field Spectra/Diffraction Spectrum of a Finite Binary Grating.

In this example, we'll define a grating that is periodic in the $y$ direction with periodicity gp and has a rectangular profile of height gh and duty cycle gdc. The grating parameters are gh=0.5 μm, gdc=0.5, and gp=10 μm. There is a semi-infinite substrate of thickness dsub adjacent to the grating. The substrate and grating are glass with a refractive index of 1.5. The surrounding is air/vacuum. Perfectly matched layers (PML) of thickness dpml are used in the $\pm x$ boundaries.

We'll do our analysis first using an ideal quartz that ignores the dispersion of the material. Then we'll include the dispersion to detect any differences.

A pulsed planewave with Ez polarization spanning wavelengths of 0.4 to 0.6 μm is normally incident on the grating from the glass substrate. The eigenmode monitor is placed in the air region. We will use mode decomposition to compute the transmittance — the ratio of the power in the +x direction of the diffracted mode relative to that of the incident planewave — for the first ten diffraction orders.

Two simulations are required: (1) an empty cell of homogeneous glass to obtain the incident power of the source, and (2) the grating structure to obtain the diffraction orders. At the end of the simulation, the wavelength, angle, and transmittance for each diffraction order are computed.

Initially, we'll import our standard libraries:

In [1]:
import meep as mp
import math
import numpy as np
import matplotlib.pyplot as plt
Using MPI version 3.1, 1 processes

We first need to simulate the empty, homogenous glass (fuzed quartz) with a constant refractive index.

In [2]:
fused_quartz = mp.Medium(index=1.5)

resolution = 50        # pixels/μm

dpml = 1.0             # PML thickness
dsub = 3.0             # substrate thickness
dpad = 3.0             # padding between grating and PML
gp = 10.0              # grating period
gh = 0.5               # grating height
gdc = 0.5              # grating duty cycle

sx = dpml+dsub+gh+dpad+dpml
sy = gp

cell_size = mp.Vector3(sx,sy,0)
pml_layers = [mp.PML(thickness=dpml,direction=mp.X)]

wvl_min = 0.4           # min wavelength
wvl_max = 0.6           # max wavelength
fmin = 1/wvl_max        # min frequency
fmax = 1/wvl_min        # max frequency
fcen = 0.5*(fmin+fmax)  # center frequency
df = fmax-fmin          # frequency width

src_pt = mp.Vector3(-0.5*sx+dpml+0.5*dsub)
sources = [mp.Source(mp.GaussianSource(fcen, fwidth=df), component=mp.Ez, center=src_pt, size=mp.Vector3(y=sy))]

k_point = mp.Vector3(0,0,0)

symmetries=[mp.Mirror(mp.Y)]

sim = mp.Simulation(resolution=resolution,
                    cell_size=cell_size,
                    boundary_layers=pml_layers,
                    k_point=k_point,
                    default_material=fused_quartz,
                    sources=sources,
                    symmetries=symmetries)

nfreq = 21
mon_pt = mp.Vector3(0.5*sx-dpml-0.5*dpad)
flux_mon = sim.add_flux(fcen, df, nfreq, mp.FluxRegion(center=mon_pt, size=mp.Vector3(y=sy)))

f = plt.figure(dpi=120)
sim.plot2D(ax=f.gca())
plt.show()
-----------
Initializing structure...

Now, we'll run the simulation and record the fields.

In [3]:
sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ez, mon_pt, 1e-9))
field decay(t = 50.01): 0.10609306658233127 / 0.10609306658233127 = 1.0
field decay(t = 100.01): 8.493199823356459e-20 / 0.10609306658233127 = 8.005424008330923e-19
run 0 finished at t = 100.01 (10001 timesteps)
In [4]:
input_flux = mp.get_fluxes(flux_mon)

Next, we'll simulate the actual grating.

In [5]:
sim.reset_meep()

geometry = [mp.Block(material=fused_quartz, size=mp.Vector3(dpml+dsub,mp.inf,mp.inf), center=mp.Vector3(-0.5*sx+0.5*(dpml+dsub))),
            mp.Block(material=fused_quartz, size=mp.Vector3(gh,gdc*gp,mp.inf), center=mp.Vector3(-0.5*sx+dpml+dsub+0.5*gh))]

sim = mp.Simulation(resolution=resolution,
                    cell_size=cell_size,
                    boundary_layers=pml_layers,
                    geometry=geometry,
                    k_point=k_point,
                    sources=sources,
                    symmetries=symmetries)

mode_mon = sim.add_flux(fcen, df, nfreq, mp.FluxRegion(center=mon_pt, size=mp.Vector3(y=sy)))

f2 = plt.figure(dpi=120)
sim.plot2D(ax=f2.gca())
plt.show()
-----------
Initializing structure...
     block, center = (-2.25,0,0)
          size (4,1e+20,1e+20)
          axes (1,0,0), (0,1,0), (0,0,1)
     block, center = (0,0,0)
          size (0.5,5,1e+20)
          axes (1,0,0), (0,1,0), (0,0,1)
In [6]:
sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ez, mon_pt, 1e-9))
field decay(t = 50.01): 0.1031398354415896 / 0.1031398354415896 = 1.0
field decay(t = 100.01): 8.275841626038753e-06 / 0.1031398354415896 = 8.023904237006029e-05
field decay(t = 150.02): 7.578862246277466e-06 / 0.1031398354415896 = 7.348142658778572e-05
field decay(t = 200.03): 2.633198313201328e-06 / 0.1031398354415896 = 2.5530371479917354e-05
field decay(t = 250.04): 1.0595609940381944e-06 / 0.1031398354415896 = 1.0273052981922272e-05
field decay(t = 300.04): 4.182093600425728e-07 / 0.1031398354415896 = 4.054780175399971e-06
field decay(t = 350.05): 1.7897453529966721e-07 / 0.1031398354415896 = 1.7352610127153491e-06
field decay(t = 400.06): 7.323581231104047e-08 / 0.1031398354415896 = 7.100633038387535e-07
field decay(t = 450.07): 2.934107857572782e-08 / 0.1031398354415896 = 2.844786250637789e-07
field decay(t = 500.08): 1.1841535133169714e-08 / 0.1031398354415896 = 1.1481049084934639e-07
field decay(t = 550.08): 4.9984066703627664e-09 / 0.1031398354415896 = 4.8462426267816536e-08
field decay(t = 600.09): 2.35073055720754e-09 / 0.1031398354415896 = 2.2791684194016495e-08
field decay(t = 650.1): 1.1816026525464885e-09 / 0.1031398354415896 = 1.1456317023267471e-08
field decay(t = 700.11): 3.9576427413983696e-10 / 0.1031398354415896 = 3.837162163827255e-09
field decay(t = 750.12): 1.4213834548284144e-10 / 0.1031398354415896 = 1.3781129752076983e-09
field decay(t = 800.13): 8.16132061584665e-11 / 0.1031398354415896 = 7.912869533778332e-10
run 0 finished at t = 800.13 (80013 timesteps)

Now we can compute the diffraction orders as a function of wavelength:

In [7]:
freqs = mp.get_eigenmode_freqs(mode_mon)

nmode = 10
res = sim.get_eigenmode_coefficients(mode_mon, range(1,nmode+1), eig_parity=mp.ODD_Z+mp.EVEN_Y)
coeffs = res.alpha
kdom = res.kdom

mode_wvl = []
mode_angle = []
mode_tran = []

for nm in range(nmode):
  for nf in range(nfreq):
    mode_wvl.append(1/freqs[nf])
    mode_angle.append(math.degrees(math.acos(kdom[nm*nfreq+nf].x/freqs[nf])))
    tran = abs(coeffs[nm,nf,0])**2/input_flux[nf]
    mode_tran.append(0.5*tran if nm != 0 else tran)

tran_max = round(max(mode_tran),1)

Note the use of the keyword parameter argument eig_parity=mp.ODD_Z+mp.EVEN_Y in the call to get_eigenmode_coefficients. This is important for specifying non-degenerate modes in MPB since the k_point is (0,0,0). ODD_Z is for modes with Ez polarization. EVEN_Y is necessary since each diffraction order which is based on a given kx consists of two modes: one going in the +y direction and the other in the -y direction. EVEN_Y forces MPB to compute only the +ky + -ky (cosine) mode. As a result, the total transmittance must be halved in this case to obtain the transmittance for the individual +ky or -ky mode. For ODD_Y, MPB will compute the +ky - -ky (sine) mode but this will have zero power because the source is even. If the $y$ parity is left out, MPB will return a random superposition of the cosine and sine modes. Alternatively, in this example an input planewave with Hz instead of Ez polarization can be used which requires eig_parity=mp.EVEN_Z+mp.ODD_Y as well as an odd mirror symmetry plane in y. Finally, note the use of add_flux instead of add_mode_monitor when using symmetries.

The diffraction spectrum is then plotted and shown in the figure below.

In [8]:
plt.figure(dpi=200)
plt.pcolormesh(np.reshape(mode_wvl,(nmode,nfreq)),
               np.reshape(mode_angle,(nmode,nfreq)),
               np.reshape(mode_tran,(nmode,nfreq)),
               cmap='Blues',
               shading='flat',
               vmin=0,
               vmax=tran_max)
plt.axis([min(mode_wvl), max(mode_wvl), min(mode_angle), max(mode_angle)])
plt.xlabel("wavelength (μm)")
plt.ylabel("diffraction angle (degrees)")
plt.xticks([t for t in np.linspace(wvl_min,wvl_max,3)])
plt.yticks([t for t in range(0,35,5)])
plt.title("transmittance of diffraction orders")
cbar = plt.colorbar()
cbar.set_ticks([t for t in np.arange(0,tran_max+0.1,0.1)])
cbar.set_ticklabels(["{:.1f}".format(t) for t in np.arange(0,tran_max+0.1,0.1)])

Now let's try using a dispersive medium instead of an ideal, constant refractive index:

In [9]:
from meep.materials import fused_quartz

resolution = 50        # pixels/μm

dpml = 1.0             # PML thickness
dsub = 3.0             # substrate thickness
dpad = 3.0             # padding between grating and PML
gp = 10.0              # grating period
gh = 0.5               # grating height
gdc = 0.5              # grating duty cycle

sx = dpml+dsub+gh+dpad+dpml
sy = gp

cell_size = mp.Vector3(sx,sy,0)
pml_layers = [mp.PML(thickness=dpml,direction=mp.X)]

wvl_min = 0.4           # min wavelength
wvl_max = 0.6           # max wavelength
fmin = 1/wvl_max        # min frequency
fmax = 1/wvl_min        # max frequency
fcen = 0.5*(fmin+fmax)  # center frequency
df = fmax-fmin          # frequency width

src_pt = mp.Vector3(-0.5*sx+dpml+0.5*dsub)
sources = [mp.Source(mp.GaussianSource(fcen, fwidth=df), component=mp.Ez, center=src_pt, size=mp.Vector3(y=sy))]

k_point = mp.Vector3(0,0,0)

symmetries=[mp.Mirror(mp.Y)]

sim = mp.Simulation(resolution=resolution,
                    cell_size=cell_size,
                    boundary_layers=pml_layers,
                    k_point=k_point,
                    default_material=fused_quartz,
                    sources=sources,
                    symmetries=symmetries)

nfreq = 21
mon_pt = mp.Vector3(0.5*sx-dpml-0.5*dpad)
flux_mon = sim.add_flux(fcen, df, nfreq, mp.FluxRegion(center=mon_pt, size=mp.Vector3(y=sy)))

sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ez, mon_pt, 1e-9))

input_flux = mp.get_fluxes(flux_mon)

sim.reset_meep()

geometry = [mp.Block(material=fused_quartz, size=mp.Vector3(dpml+dsub,mp.inf,mp.inf), center=mp.Vector3(-0.5*sx+0.5*(dpml+dsub))),
            mp.Block(material=fused_quartz, size=mp.Vector3(gh,gdc*gp,mp.inf), center=mp.Vector3(-0.5*sx+dpml+dsub+0.5*gh))]

sim = mp.Simulation(resolution=resolution,
                    cell_size=cell_size,
                    boundary_layers=pml_layers,
                    geometry=geometry,
                    k_point=k_point,
                    sources=sources,
                    symmetries=symmetries)

mode_mon = sim.add_flux(fcen, df, nfreq, mp.FluxRegion(center=mon_pt, size=mp.Vector3(y=sy)))

sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ez, mon_pt, 1e-9))

freqs = mp.get_eigenmode_freqs(mode_mon)

nmode = 10
res = sim.get_eigenmode_coefficients(mode_mon, range(1,nmode+1), eig_parity=mp.ODD_Z+mp.EVEN_Y)
coeffs = res.alpha
kdom = res.kdom

mode_wvl = []
mode_angle = []
mode_tran = []

for nm in range(nmode):
  for nf in range(nfreq):
    mode_wvl.append(1/freqs[nf])
    mode_angle.append(math.degrees(math.acos(kdom[nm*nfreq+nf].x/freqs[nf])))
    tran = abs(coeffs[nm,nf,0])**2/input_flux[nf]
    mode_tran.append(0.5*tran if nm != 0 else tran)

tran_max = round(max(mode_tran),1)

plt.figure(dpi=200)
plt.pcolormesh(np.reshape(mode_wvl,(nmode,nfreq)),
               np.reshape(mode_angle,(nmode,nfreq)),
               np.reshape(mode_tran,(nmode,nfreq)),
               cmap='Blues',
               shading='flat',
               vmin=0,
               vmax=tran_max)
plt.axis([min(mode_wvl), max(mode_wvl), min(mode_angle), max(mode_angle)])
plt.xlabel("wavelength (μm)")
plt.ylabel("diffraction angle (degrees)")
plt.xticks([t for t in np.linspace(wvl_min,wvl_max,3)])
plt.yticks([t for t in range(0,35,5)])
plt.title("transmittance of diffraction orders")
cbar = plt.colorbar()
cbar.set_ticks([t for t in np.arange(0,tran_max+0.1,0.1)])
cbar.set_ticklabels(["{:.1f}".format(t) for t in np.arange(0,tran_max+0.1,0.1)])
-----------
Initializing structure...
field decay(t = 50.01): 0.11139427530016409 / 0.11139427530016409 = 1.0
field decay(t = 100.01): 1.824305440068526e-15 / 0.11139427530016409 = 1.637701250941967e-14
run 0 finished at t = 100.01 (10001 timesteps)
-----------
Initializing structure...
     block, center = (-2.25,0,0)
          size (4,1e+20,1e+20)
          axes (1,0,0), (0,1,0), (0,0,1)
     block, center = (0,0,0)
          size (0.5,5,1e+20)
          axes (1,0,0), (0,1,0), (0,0,1)
field decay(t = 50.01): 0.1082498119209954 / 0.1082498119209954 = 1.0
field decay(t = 100.01): 7.512926070352665e-06 / 0.1082498119209954 = 6.940359467632024e-05
field decay(t = 150.02): 6.732414916367269e-06 / 0.1082498119209954 = 6.219331744687766e-05
field decay(t = 200.03): 2.3712635292765586e-06 / 0.1082498119209954 = 2.1905474819736332e-05
field decay(t = 250.04): 9.494875348809632e-07 / 0.1082498119209954 = 8.771262675023706e-06
field decay(t = 300.04): 3.8220269083178343e-07 / 0.1082498119209954 = 3.530746927400933e-06
field decay(t = 350.05): 1.5689447663324306e-07 / 0.1082498119209954 = 1.4493741268368232e-06
field decay(t = 400.06): 6.492618040375044e-08 / 0.1082498119209954 = 5.997809996301509e-07
field decay(t = 450.07): 2.9338161615346314e-08 / 0.1082498119209954 = 2.710227490903943e-07
field decay(t = 500.08): 1.115652431764312e-08 / 0.1082498119209954 = 1.0306275936798441e-07
field decay(t = 550.08): 4.421515879496221e-09 / 0.1082498119209954 = 4.084548324872104e-08
field decay(t = 600.09): 1.6989786783603636e-09 / 0.1082498119209954 = 1.5694980418075362e-08
field decay(t = 650.1): 8.779199940057175e-10 / 0.1082498119209954 = 8.110129509014344e-09
field decay(t = 700.11): 2.64630974046326e-10 / 0.1082498119209954 = 2.444632183189964e-09
field decay(t = 750.12): 1.5056115102443512e-10 / 0.1082498119209954 = 1.39086755304776e-09
field decay(t = 800.13): 6.39132960275918e-11 / 0.1082498119209954 = 5.904240838241642e-10
run 0 finished at t = 800.13 (80013 timesteps)

Since fuzed quartz isn't highly dispersive in this range, we don't see much of a difference between the two plots.