Finally, we can validate the results for the diffraction spectra of a finite grating via a different approach than computing the far fields: as the (spatial) Fourier transform of the scattered fields. This involves two simulations — one with the grating and the other with just a flat surface — and subtracting the Fourier-transformed fields at a given frequency ω from the two runs to obtain the scattered fields s(y). The Fourier transform of the scattered fields is then computed in post processing: a(ky) = ∫ s(y) exp(ikyy) dy, where |a(ky)|² is the amplitude of the corresponding Fourier component. For a grating with periodicity Λ, we should expect to see peaks in the diffraction spectra at ky=2πm/Λ for m=0, ±1, ±2, ... The total number of diffraction orders is determined by the wavelength as described in Tutorials/Mode Decomposition/Transmittance Spectra for Planewave at Normal Incidence.
The simulation setup is shown in the schematic below. The binary grating has Λ = 1 μm at a wavelength of 0.5 μm via a normally-incident planewave pulse (which must extend into the PML region in order to span the entire width of the cell). The grating structure is terminated with a flat-surface padding in order to give the scattered field space to decay at the edge of the cell.
import meep as mp
import numpy as np
import math
import matplotlib.pyplot as plt
# True: plot the scattered fields in the extended air region adjacent to the grating
# False: plot the diffraction spectra based on a 1d cross section of the scattered fields
field_profile = True
resolution = 50 # pixels/μm
dpml = 1.0 # PML thickness
dsub = 2.0 # substrate thickness
dpad = 1.0 # flat-surface padding
gp = 1.0 # grating periodicity
gh = 0.5 # grating height
gdc = 0.5 # grating duty cycle
num_cells = 5 # number of grating unit cells
# air region thickness adjacent to grating
dair = 10 if field_profile else dpad
wvl = 0.5 # center wavelength
fcen = 1 / wvl # center frequency
k_point = mp.Vector3()
glass = mp.Medium(index=1.5)
pml_layers = [mp.PML(thickness=dpml)]
symmetries = [mp.Mirror(mp.Y)]
sx = dpml + dsub + gh + dair + dpml
sy = dpml + dpad + num_cells * gp + dpad + dpml
cell_size = mp.Vector3(sx, sy)
src_pt = mp.Vector3(-0.5 * sx + dpml + 0.5 * dsub)
sources = [
mp.Source(
mp.GaussianSource(fcen, fwidth=0.2 * fcen),
component=mp.Ez,
center=src_pt,
size=mp.Vector3(y=sy),
)
]
geometry = [
mp.Block(
material=glass,
size=mp.Vector3(dpml + dsub, mp.inf, mp.inf),
center=mp.Vector3(-0.5 * sx + 0.5 * (dpml + dsub)),
)
]
sim = mp.Simulation(
resolution=resolution,
cell_size=cell_size,
boundary_layers=pml_layers,
geometry=geometry,
k_point=k_point,
sources=sources,
symmetries=symmetries,
)
mon_pt = mp.Vector3(0.5 * sx - dpml - 0.5 * dair)
near_fields = sim.add_dft_fields(
[mp.Ez],
fcen,
0,
1,
center=mon_pt,
size=mp.Vector3(dair if field_profile else 0, sy - 2 * dpml),
)
sim.run(until_after_sources=100)
flat_dft = sim.get_dft_array(near_fields, mp.Ez, 0)
sim.reset_meep()
for j in range(num_cells):
geometry.append(
mp.Block(
material=glass,
size=mp.Vector3(gh, gdc * gp, mp.inf),
center=mp.Vector3(
-0.5 * sx + dpml + dsub + 0.5 * gh,
-0.5 * sy + dpml + dpad + (j + 0.5) * gp,
),
)
)
sim = mp.Simulation(
resolution=resolution,
cell_size=cell_size,
boundary_layers=pml_layers,
geometry=geometry,
k_point=k_point,
sources=sources,
symmetries=symmetries,
)
near_fields = sim.add_dft_fields(
[mp.Ez],
fcen,
0,
1,
center=mon_pt,
size=mp.Vector3(dair if field_profile else 0, sy - 2 * dpml),
)
sim.run(until_after_sources=100)
----------- Initializing structure... Halving computational cell along direction y time for choose_chunkdivision = 0.00181007 s Working in 2D dimensions. Computational cell is 14.5 x 9 x 0 with resolution 50 block, center = (-5.75,0,0) size (3,1e+20,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) time for set_epsilon = 0.610417 s ----------- Meep progress: 9.18/125.0 = 7.3% done in 4.0s, 50.5s to go on time step 918 (time=9.18), 0.00435765 s/step Meep progress: 18.490000000000002/125.0 = 14.8% done in 8.0s, 46.1s to go on time step 1849 (time=18.49), 0.00429697 s/step Meep progress: 27.810000000000002/125.0 = 22.2% done in 12.0s, 42.0s to go on time step 2781 (time=27.81), 0.00429364 s/step Meep progress: 37.160000000000004/125.0 = 29.7% done in 16.0s, 37.8s to go on time step 3716 (time=37.16), 0.0042805 s/step Meep progress: 46.44/125.0 = 37.2% done in 20.0s, 33.9s to go on time step 4644 (time=46.44), 0.0043122 s/step Meep progress: 55.75/125.0 = 44.6% done in 24.0s, 29.8s to go on time step 5575 (time=55.75), 0.004299 s/step Meep progress: 65.19/125.0 = 52.2% done in 28.0s, 25.7s to go on time step 6519 (time=65.19), 0.00423808 s/step Meep progress: 74.61/125.0 = 59.7% done in 32.0s, 21.6s to go on time step 7461 (time=74.61), 0.00424806 s/step Meep progress: 83.91/125.0 = 67.1% done in 36.0s, 17.6s to go on time step 8392 (time=83.92), 0.00430039 s/step Meep progress: 93.23/125.0 = 74.6% done in 40.0s, 13.6s to go on time step 9324 (time=93.24), 0.00429399 s/step Meep progress: 102.51/125.0 = 82.0% done in 44.0s, 9.7s to go on time step 10253 (time=102.53), 0.00430959 s/step Meep progress: 111.78/125.0 = 89.4% done in 48.0s, 5.7s to go on time step 11180 (time=111.8), 0.00431786 s/step Meep progress: 121.06/125.0 = 96.8% done in 52.0s, 1.7s to go on time step 12108 (time=121.08), 0.00431251 s/step run 0 finished at t = 125.0 (12500 timesteps) ----------- Initializing structure... Halving computational cell along direction y time for choose_chunkdivision = 0.00184488 s Working in 2D dimensions. Computational cell is 14.5 x 9 x 0 with resolution 50 block, center = (-5.75,0,0) size (3,1e+20,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (-4,-2,0) size (0.5,0.5,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (-4,-1,0) size (0.5,0.5,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (-4,0,0) size (0.5,0.5,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (-4,1,0) size (0.5,0.5,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) block, center = (-4,2,0) size (0.5,0.5,1e+20) axes (1,0,0), (0,1,0), (0,0,1) dielectric constant epsilon diagonal = (2.25,2.25,2.25) time for set_epsilon = 0.655086 s ----------- Meep progress: 9.15/125.0 = 7.3% done in 4.0s, 50.7s to go on time step 915 (time=9.15), 0.00437325 s/step Meep progress: 18.43/125.0 = 14.7% done in 8.0s, 46.3s to go on time step 1844 (time=18.44), 0.00430983 s/step Meep progress: 27.73/125.0 = 22.2% done in 12.0s, 42.1s to go on time step 2774 (time=27.74), 0.00430151 s/step Meep progress: 36.980000000000004/125.0 = 29.6% done in 16.0s, 38.1s to go on time step 3700 (time=37), 0.00432353 s/step Meep progress: 46.13/125.0 = 36.9% done in 20.0s, 34.2s to go on time step 4615 (time=46.15), 0.00437322 s/step Meep progress: 55.32/125.0 = 44.3% done in 24.0s, 30.2s to go on time step 5535 (time=55.35), 0.00435253 s/step Meep progress: 64.53/125.0 = 51.6% done in 28.0s, 26.2s to go on time step 6457 (time=64.57), 0.00434152 s/step Meep progress: 73.79/125.0 = 59.0% done in 32.0s, 22.2s to go on time step 7383 (time=73.83), 0.00432269 s/step Meep progress: 83.04/125.0 = 66.4% done in 36.0s, 18.2s to go on time step 8308 (time=83.08), 0.00432687 s/step Meep progress: 92.32000000000001/125.0 = 73.9% done in 40.0s, 14.2s to go on time step 9237 (time=92.37), 0.00431051 s/step Meep progress: 101.59/125.0 = 81.3% done in 44.0s, 10.1s to go on time step 10164 (time=101.64), 0.00431663 s/step Meep progress: 110.92/125.0 = 88.7% done in 48.0s, 6.1s to go on time step 11097 (time=110.97), 0.00428761 s/step Meep progress: 120.25/125.0 = 96.2% done in 52.0s, 2.1s to go on time step 12030 (time=120.3), 0.00428792 s/step run 0 finished at t = 125.0 (12500 timesteps)
grating_dft = sim.get_dft_array(near_fields, mp.Ez, 0)
scattered_field = grating_dft - flat_dft
scattered_amplitude = np.abs(scattered_field) ** 2
[x, y, z, w] = sim.get_array_metadata(dft_cell=near_fields)
if field_profile:
plt.figure(dpi=150)
plt.pcolormesh(
x,
y,
np.rot90(scattered_amplitude),
cmap="inferno",
shading="gouraud",
vmin=0,
vmax=scattered_amplitude.max(),
)
plt.gca().set_aspect("equal")
plt.xlabel("x (μm)")
plt.ylabel("y (μm)")
# ensure that the height of the colobar matches that of the plot
from mpl_toolkits.axes_grid1 import make_axes_locatable
divider = make_axes_locatable(plt.gca())
cax = divider.append_axes("right", size="5%", pad=0.05)
plt.colorbar(cax=cax)
plt.tight_layout()
plt.show()
else:
ky = np.fft.fftshift(np.fft.fftfreq(len(scattered_field), 1 / resolution))
FT_scattered_field = np.fft.fftshift(np.fft.fft(scattered_field))
plt.figure(dpi=150)
plt.subplots_adjust(hspace=0.3)
plt.subplot(2, 1, 1)
plt.plot(y, scattered_amplitude, "bo-")
plt.xlabel("y (μm)")
plt.ylabel("field amplitude")
plt.subplot(2, 1, 2)
plt.plot(ky, np.abs(FT_scattered_field) ** 2, "ro-")
plt.gca().ticklabel_format(axis="y", style="sci", scilimits=(0, 0))
plt.xlabel(r"wavevector k$_y$, 2π (μm)$^{-1}$")
plt.ylabel("Fourier transform")
plt.gca().set_xlim([-3, 3])
plt.tight_layout(pad=1.0)
plt.show()
Results are shown for two finite gratings with 5 and 20 periods.
The scattered field amplitude profile (the top figure in each of the two sets of results) shows that the fields are nonzero above the grating (which is positioned at the left edge of the figure in the region indicated by the bright spots) and decay to zero away from the grating. The middle figure is the field amplitude along a 1d slice at a distance 1.5 μm above the grating (marked by the dotted green line in the top figure). Note the decaying fields at the edges due to the flat-surface termination. The bottom figure is the Fourier transform of the fields from the 1d slice. As expected, there are only three diffraction orders present at ky=2πm/Λ for m=0, ±1, ±2. These peaks are becoming sharper as the number of grating periods increases.
The sharpness of the peaks directly corresponds to how collimated the diffracted beams are, and in the limit of infinitely many periods the resulting delta-function peaks correspond to diffracted planewaves. (The squared amplitude of each peak is proportional to the power in the corresponding diffraction order.) One can also obtain the collimation of the beams more directly by using Meep's near2far
feature to compute the far-field diffracted waves — this approach is more straightforward, but potentially much more expensive than looking at the Fourier transform of the near field, because one may need a large number of far-field points to resolve the full diffracted beams. In general, there is a tradeoff in computational science between doing direct "numerical experiments" that are conceptually straightforward but often expensive, versus more indirect and tricky calculations that don't directly correspond to laboratory experiments but which can sometimes be vastly more efficient at extracting physical information.
In 3d, the procedure is very similar, but a little more effort is required to disentangle the two polarizations relative to the plane of incidence [the (z,k) plane for each Fourier component k]. For propagation in the $z$ direction, you would Fourier transform both $E_x$ and $E_y$ of the scattered field as a function of k $= (k_x, k_y)$. For each k, you decompose the corresponding E $= (E_x, E_y)$ into the amplitude parallel to k [which gives the p polarization amplitude if you multiply by csc(θ), where cos(θ)=|k|/(nω/c), n is the refractive index of the ambient medium, and ω is the angular frequency] and perpendicular to k [which equals the s polarization amplitude]. Then square these amplitudes to get something proportional to power as above. (Note that this analysis is the same even if the incident wave is at an oblique angle, although the k locations of the diffraction peaks will change.)