failing to reproduce the reflectance of a thin film on Si

[MahmutRuzi]

Hello, Dear meep developers.
Thanks for freely sharing this nice code, along with the excellent tutorials and documentation.
I am new to Meep and FDTD simulation in general. I tried to simulate the reflectance spectrum of a film (5.5 um thick, n=1.57) on a substrate (1um thick, n=3.41). The 1D simulation as shown in the examples give similar results to analytical results. However, when I tried to do the simulations in 2D, the result does not look like a reflectance spectrum. Can someone help me tell what I was doing wrong or how to improve?

#%% import libraries
import argparse
import meep as mp
import numpy as np
import math
import cmath
import matplotlib.pyplot as plt
from numpy import savez_compressed

#%%
def main(args):
resolution = args.res
theta = args.theta
l = args.struct_length
dsub = args.sub_thickness
dpml = args.pml_thickness
# w = args.width
# h = args.height

#%%
# wavelength range (in micron)
wvl_min = 2.5
wvl_max = 18.2
dair = 1.3*wvl_max
dthin = 5.5 # thin film thickness

sx = 2*dpml + l
sy = 2*dpml + dair + dsub + dthin

cell_size = mp.Vector3(sx,sy,0)

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

# dft-field parameters
frq_min = 1/wvl_max
frq_max = 1/wvl_min
frq_cen = 0.5*(frq_min+frq_max)
dfrq = frq_max-frq_min
nfrq = 150

#%%

# CCW rotation angle (degrees) about Y-axis of PW current source; 0 degrees along -z axis
#theta = 0

# k with correct length (plane of incidence: XZ) 
k = mp.Vector3(math.sin(theta),0,math.cos(theta)).scale(frq_cen)

# wavevector, important for angle dependent simulations
def pw_amp(k, x0):
    def _pw_amp(x):
        return cmath.exp(1j * 2 * math.pi * k.dot(x + x0))
    return _pw_amp
#%%
src_pos = 0.5*sy -dpml-0.1*dair  # position of the source

#%%
sources = [mp.Source(mp.GaussianSource(frq_cen,fwidth=dfrq,is_integrated=False,cutoff=5),component=mp.Ex,
                     center=mp.Vector3(0,src_pos,0),size=mp.Vector3(l,0,0),
                     amp_func=pw_amp(k, mp.Vector3(0,src_pos,0)))]

sim = mp.Simulation(cell_size=cell_size,
                    boundary_layers=pml_layers,
                    sources=sources,
                    k_point= k,
                    dimensions=2,
                    resolution=resolution,
                    default_material=mp.Medium(index=1.0), 
                    split_chunks_evenly=False)

#%%
# set up flux monitor

fluxmonitor_posY1 = 0.5*sy-dpml-0.2*dair

flux_pos1 = sim.add_flux(frq_cen, dfrq, nfrq, mp.FluxRegion(center=mp.Vector3(0,fluxmonitor_posY1,0),
                                                    size=mp.Vector3(l,0,0)))
#%%                  
sim.run(until_after_sources=mp.stop_when_fields_decayed(10, mp.Ex, mp.Vector3(0,fluxmonitor_posY1,0),1e-7))

#%%
empty_flux_pos1 = mp.get_fluxes(flux_pos1)
empty_data_pos1 = sim.get_flux_data(flux_pos1)
freqs = mp.get_flux_freqs(flux_pos1)

print('First Pos1 REFLECTION FLUX')
sim.display_fluxes(flux_pos1)



#%%
sim.reset_meep()

geometry = [mp.Block(material=mp.Medium(index=3.41), size=mp.Vector3(l, dsub), 
                   center=mp.Vector3(0, -0.5*sy + dpml + 0.5*dsub)),
            mp.Block(material=mp.Medium(index=1.57), size=mp.Vector3(l, dthin), 
                               center=mp.Vector3(0, -0.5*sy + dpml + 0.5*dsub+ 0.5*dthin))]

sim = mp.Simulation(cell_size=cell_size,
                    boundary_layers=pml_layers,
                    eps_averaging=True,
                    geometry=geometry,
                    sources=sources,
                    k_point= k,
                    dimensions=2,
                    resolution=resolution,
                    default_material=mp.Medium(index=1.0), 
                    split_chunks_evenly=False)

#%%
# add fluxregion

flux2_pos1 = sim.add_flux(frq_cen, dfrq, nfrq, mp.FluxRegion(center=mp.Vector3(0,fluxmonitor_posY1,0),
                                                    size=mp.Vector3(l,0,0)))

#%%
sim.load_minus_flux_data(flux2_pos1, empty_data_pos1)
#%%
sim.run(until_after_sources=mp.stop_when_fields_decayed(10, mp.Ex, mp.Vector3(0,fluxmonitor_posY1,0),1e-7))

print('Second Pos1 REFLECTION FLUX')
sim.display_fluxes(flux2_pos1)
    
second_flux_pos1 =  mp.get_fluxes(flux2_pos1)


#%%

wl = []
Rs = []
for i in range(nfrq):
    wl = np.append(wl, 1/freqs[i])
    Rs = np.append(Rs,- second_flux_pos1[i]/empty_flux_pos1[i])
    savez_compressed('trans-ref-spectrum-Si-film-'+'Res={}'.format(resolution)+
                      '-StructureX={} um'.format(l)+'-pml_thickness={} um'.format(dpml),
                      [wl, Rs])

#%%
if mp.am_master():
    plt.figure()
    plt.plot(wl,Rs,'bo-',label='reflectance-Pos1',markersize=2)
    #plt.plot(wl,1-Rs-Ts,'go-',label='loss')
    #plt.axis([5.0, 10.0, 0, 1])
    plt.xlabel("wavelength um)")
    plt.legend(loc="best")
    #plt.show()
    plt.tight_layout()
    plt.savefig('trans-ref-spectrum-Si-film-'+'Res={}'.format(resolution)+
                      '-StructureX={} um'.format(l)+'-pml_thickness={} um'.format(dpml)+'.tiff',
                format='tiff')
         
    
    plt.figure()
    sim.plot2D(output_plane=mp.Volume(size=mp.Vector3(sx,sy,0)))
    plt.tight_layout()
    plt.savefig('structure.tiff', format='tiff')
    
    exit()
#%%

#%%
if name == 'main':
parser = argparse.ArgumentParser()
parser.add_argument('-res', type=float, default=200, help='resolution (default: 200 pixels/um)')
parser.add_argument('-theta', type=float, default=0, help='angle of incident planewave (default: 0 degrees)')
parser.add_argument('-struct_length', type=float, default=1, help='structure length (default: 1 um)')
parser.add_argument('-sub_thickness', type=float, default=1, help='substrate thickness (default: 1 um)')
parser.add_argument('-pml_thickness', type=float, default=10, help='pml thickness in microns (default: 10 um)')
args = parser.parse_args()
main(args)

[stevengj]

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

This puts an absorber (not sure why you aren't using a PML, which is better?) on all sides, whereas you should only have PML in the vertical (y) direction. In the other (x) direction, you want periodic boundary conditions (to model an infinite flat surface), so you don't want PML and you want to set the k_point=(0,0,0) to make sure it is periodic.

See also the diffraction grating tutorial.

[MahmutRuzi]

@stevengj Thanks for the prompt reply. I originally tried to simulate a metallic 2D grating (1.3 um long and 100 nm thick Au bar) but couldn't get a convergence in a reasonable time. So I changed to dielectric simple structure with known analytical results to figure out how to get a reasonably accurate results within short time. The mp.Absorber is left from there. With respect to the periodicity, I want to simulate an isolated structure. The reason is to explore the effect of the overlap of beam size and structure, and how that affect (if any) the reflectance/transmittance.
I'll checkout the diffraction grating tutorial. Thank you !