tutorial example for add_near2far with periodic boundaries

[oskooi]

• Python tutorial for add_near2far with periodic boundaries (Scheme version will be added in a separate PR)
• revamp Tutorials/Optical Forces to combine simulation and visualization of results into a single script

[stevengj]

This is a problematic example. You are using a distance L to your far-field plane such that L >> nperiods*a. This means that your periodic Green's function isn't actually converged, so you aren't actually seeing the far fields of a periodic structure.

(An infinite periodic structure actually has no spatial separation of the diffracted orders — they are all present at every far-field point.)

What you are simulating is more similar to the radiation of a finite periodic surface, which does indeed send approximate diffracted orders to different spatial points in the far field — essentially it is a directed antenna. However, this simulation is different from that of a truly finite periodic surface because the edges are handled differently, and I don't think what you are doing exactly corresponds to any reasonable physical simulation.

What you are doing could be viewed as an approximation (akin to the locally periodic approximation that we often use for metasurface design) for a finite periodic surface, where there is an error of O(1/nperiods) in the result due to the incorrect edge handling. For nperiods=10 I'm guessing you have a ≈10% difference in the far fields compared to a true finite simulation.

One could adapt this calculation to illustrate the local periodic approximation. For that, however, I think you would also want to do a truly finite simulation (with a 10x bigger computational cell) for comparison, with a diagram showing the finite structure and the far fields (and maybe an actual picture of the simulation results for the finite structure in a giant computational cell).

It would also be good to show that the error goes roughly as O(1/nperiods) by trying nperiods=10 and nperiods=20 and comparing to the corresponding finite structures.

[stevengj]

(There are infinitely many ways to terminate a finite periodic structure, of course, and different choices will have slightly different errors compared to the periodic approximation. In this case, I would suggest something similar to a reasonable experimental geometry: terminate the finite grating with a simple flat/unetched surface that extends into a PML.)

[oskooi]

As suggested, the following script computes the far fields of the binary grating in two ways via: (1) a unit cell with periodic boundaries and add_near2far containing nperiods, and (2) a supercell containing 2*nperiods+1 unit cells with a flat surface termination extending into PML (schematic shown below). The error is computed as the L₂-norm of the difference of the two far-field datasets with each containing 500 (position) x 21 (frequency) points.

# -*- coding: utf-8 -*-                                                                                                                                       

import meep as mp
import math
import numpy as np
from numpy import linalg as LA
import matplotlib.pyplot as plt

resolution = 40        # 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                                                                                                                   

nperiods = 20          # number of unit cells in supercell                                                                                                    

ff_distance = 1e8      # far-field distance from near-field monitor (100 m or 100/(0.5e-6)=2e8 periods)                                                       
ff_angle = 20          # far-field cone angle                                                                                                                 
ff_npts = 500          # number of far-field points                                                                                                           

ff_length = ff_distance*math.tan(math.radians(ff_angle))
ff_res = ff_npts/ff_length

sx = dpml+dsub+gh+dpad+dpml
cell_size = mp.Vector3(sx)

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)]

k_point = mp.Vector3()

glass = mp.Medium(index=1.5)

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

nfreq = 21
n2f_pt = mp.Vector3(0.5*sx-dpml-0.5*dpad)
n2f_obj = sim.add_near2far(fcen, df, nfreq, mp.Near2FarRegion(center=n2f_pt))

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

ff_source = np.absolute(sim.get_farfields(n2f_obj, ff_res, center=mp.Vector3(ff_distance,0.5*ff_length), size=mp.Vector3(y=ff_length))['Ez'][0])**2

sim.reset_meep()

### unit cell with periodic boundaries                                                                                                                        

sy = gp
cell_size = mp.Vector3(sx,sy)
symmetries = [mp.Mirror(mp.Y)]

sources = [mp.Source(mp.GaussianSource(fcen, fwidth=df), 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))),
            mp.Block(material=glass, 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)

n2f_obj = sim.add_near2far(fcen, df, nfreq, mp.Near2FarRegion(center=n2f_pt, size=mp.Vector3(y=sy)), nperiods=nperiods)

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

ff_unitcell = np.absolute(sim.get_farfields(n2f_obj, ff_res, center=mp.Vector3(ff_distance,0.5*ff_length), size=mp.Vector3(y=ff_length))['Ez'][0])**2

sim.reset_meep()

### supercell with PML boundaries                                                                                                                             

num_cells = 2*nperiods+1
sy = dpml+dpad+num_cells*gp+dpad+dpml
cell_size = mp.Vector3(sx,sy)

pml_layers = [mp.PML(thickness=dpml)]

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

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

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)

n2f_obj = sim.add_near2far(fcen, df, nfreq, mp.Near2FarRegion(center=n2f_pt, size=mp.Vector3(y=sy-2*dpml)))

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

ff_supercell = np.absolute(sim.get_farfields(n2f_obj, ff_res, center=mp.Vector3(ff_distance,0.5*ff_length), size=mp.Vector3(y=ff_length))['Ez'][0])**2

norm_err = LA.norm(ff_unitcell-ff_supercell)
print("norm:, {}, {}".format(nperiods,norm_err))

freqs = mp.get_near2far_freqs(n2f_obj)
wvl = np.divide(1,freqs)
ff_lengths = np.linspace(0,ff_length,ff_npts)
angles = [math.degrees(math.atan(f)) for f in ff_lengths/ff_distance]

wvl_slice = 0.5
idx_slice = np.where(np.asarray(freqs) == 1/wvl_slice)[0][0]

enh_unitcell = ff_unitcell/ff_source
enh_supercell = ff_supercell/ff_source

plt.figure(dpi=200)

plt.subplot(2,2,1)
plt.pcolormesh(wvl,angles,enh_unitcell,cmap='Blues',shading='flat')
plt.axis([wvl_min,wvl_max,0,ff_angle])
plt.xlabel("wavelength (μm)")
plt.ylabel("angle (degrees)")
plt.grid(linewidth=0.5,linestyle='--')
plt.xticks([t for t in np.arange(wvl_min,wvl_max+0.1,0.1)])
plt.yticks([t for t in range(0,ff_angle+1,10)])
plt.title("unit cell: diffraction spectra")

plt.subplot(2,2,2)
plt.plot(angles,enh_unitcell[:,idx_slice],'bo-')
plt.xlim(0,ff_angle)
plt.ylim(0)
plt.xticks([t for t in range(0,ff_angle+1,10)])
plt.xlabel("angle (degrees)")
plt.ylabel("relative enhancement")
plt.grid(axis='x',linewidth=0.5,linestyle='--')
plt.title("unit cell: spectra @  λ = {:.1} μm".format(wvl_slice))

plt.subplot(2,2,3)
plt.pcolormesh(wvl,angles,enh_supercell,cmap='Blues',shading='flat')
plt.axis([wvl_min,wvl_max,0,ff_angle])
plt.xlabel("wavelength (μm)")
plt.ylabel("angle (degrees)")
plt.grid(linewidth=0.5,linestyle='--')
plt.xticks([t for t in np.arange(wvl_min,wvl_max+0.1,0.1)])
plt.yticks([t for t in range(0,ff_angle+1,10)])
plt.title("supercell: diffraction spectra")

plt.subplot(2,2,4)
plt.plot(angles,enh_supercell[:,idx_slice],'bo-')
plt.xlim(0,ff_angle)
plt.ylim(0)
plt.xticks([t for t in range(0,ff_angle+1,10)])
plt.xlabel("angle (degrees)")
plt.ylabel("relative enhancement")
plt.grid(axis='x',linewidth=0.5,linestyle='--')
plt.title("supercell: spectra @  λ = {:.1} μm".format(wvl_slice))

plt.tight_layout(pad=0.5)
plt.show()

Results for nperiods=10 and nperiods=20 are shown in the two figures below.

nperiods=10

nperiods=20

While the far-field spectra from the two methods appears quantitatively similar from the plots, the actual error in the fields increases by ~40% when nperiods is doubled from 10 to 20:

norm:, 10, 4.188976191744251e-06
norm:, 20, 5.852924638020601e-06

Doubling nperiods should result in the error being halved but this does not seem to be the case.

For reference, a schematic of the finite structure containing 10 unit cells and a flat-surface termination is shown below.

[oskooi]

There was a bug in the script posted above involving the computation of the error in the far fields in the line norm_err = LA.norm(ff_unitcell-ff_supercell): the L₂-norm of the difference in the far-fields was not being normalized by nperiods. With this correction, the tutorial has been updated and revised to demonstrate the O(1/nperiods) scaling of the error for the finite periodic structure.