Unexpected reflectance/transmittance results for Ag nanoparticles on Si substrate in Meep (2D)

[rapsanaminanonto-dev]

Hello everyone,

I am simulating silver (Ag) nanoparticles on a silicon (Si) substrate using Meep (Python), but I am not obtaining physically meaningful reflectance, transmittance, or absorption spectra. I would appreciate any feedback on whether my setup or methodology is incorrect.

Simulation setup

Dimension: 2D

Resolution: 60

Cell size: 4 × 8 μm

Boundary condition: PML (1 μm thick)

Polarization: Ez (TM)

Geometry

Silicon block (substrate):

Size: 2 × 0.5 μm

Center: (0, −0.45)

Material: meep.materials.Si

Silver nanoparticles:

Cylinders (2D) with radius R = 0.1 μm

Positioned along x at y = 0

Material: meep.materials.Ag

The geometry consists of 7 Ag cylinders placed above a Si block.

Source and monitors

Source: Gaussian pulse

Center frequency fcen = 3

Frequency width df = 4

Source position: (0, +3)

Flux regions:

Reflection monitor at y = +2.5

Transmission monitor at y = −3

I first perform a normalization run (Si block only) to record incident fields, then a second run with Ag nanoparticles added and subtract the incident fields using load_minus_flux_data.

code :

import matplotlib.pyplot as plt
import numpy as np
import meep as mp
from meep.materials import Ag
from meep.materials import Si

resolution = 60
sx = 4 # size of cell in X direction
sy = 8 # size of cell in Y direction
cell = mp.Vector3(sx, sy, 0)
dpml = 1.0
pml_layers = [mp.PML(dpml)]

pad = 3.5 # padding distance between waveguide and cell edge
w = 1 # width of waveguide
wvg_xcen = 0.5 * (sx - w - 2 * pad) # x center of vert. wvg
wvg_ycen = -0.5 * (sy - w - 2 * pad) # y center of horiz. wvg

geometry = [

    mp.Block(
    size=mp.Vector3(2, .5),
    center=mp.Vector3(0, -.45),
    material=Si 
), 

]

fcen = 3 # pulse center frequency
df = 4 # pulse width (in frequency)

sources = [
mp.Source(
mp.GaussianSource(fcen, fwidth=df),
component=mp.Ez,
center=mp.Vector3(0, +3, 0),
size=mp.Vector3(2*w, 0, 0),
)
]

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

nfreq = 100 # number of frequencies at which to compute flux

reflected flux

refl_fr = mp.FluxRegion(
center=mp.Vector3(0, 2.5, 0),
size=mp.Vector3(2 * w, 0, 0)
)
refl = sim.add_flux(fcen, df, nfreq, refl_fr)

transmitted flux

tran_fr = mp.FluxRegion(
center=mp.Vector3(0, -3, 0),
size=mp.Vector3(2 * w, 0, 0)
)
tran = sim.add_flux(fcen, df, nfreq, tran_fr)

pt = mp.Vector3(0, -2.5)
sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ez, pt, 1e-3))

for normalization run, save flux fields data for reflection plane

straight_refl_data = sim.get_flux_data(refl)

save incident power for transmission plane

straight_tran_flux = mp.get_fluxes(tran)

sim.reset_meep()

R = 0.100 # radius

geometry = [
mp.Cylinder(radius=R, center=mp.Vector3(-.9, 0), material=Ag),
mp.Cylinder(radius=R, center=mp.Vector3(-.6, 0), material=Ag),
mp.Cylinder(radius=R, center=mp.Vector3(-.3, 0), material=Ag),
mp.Cylinder(radius=R, center=mp.Vector3(0, 0), material=Ag),
mp.Cylinder(radius=R, center=mp.Vector3(.3, 0), material=Ag),
mp.Cylinder(radius=R, center=mp.Vector3(.6, 0), material=Ag),
mp.Cylinder(radius=R, center=mp.Vector3(.9, 0), material=Ag),

  mp.Block(
    size=mp.Vector3(2, .5),
    center=mp.Vector3(0, -.45),
    material=Si 
), 

]

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

reflected flux

refl = sim.add_flux(fcen, df, nfreq, refl_fr)

tran_fr = mp.FluxRegion(
center=mp.Vector3(0, -3, 0),
size=mp.Vector3(2 * w, 0, 0)
)
tran = sim.add_flux(fcen, df, nfreq, tran_fr)

for normal run, load negated fields to subtract incident from refl. fields

sim.load_minus_flux_data(refl, straight_refl_data)

pt = mp.Vector3(0, -2.5)
sim.run(until_after_sources=mp.stop_when_fields_decayed(100, mp.Ez, pt, 1e-6))

GET FLUX DATA (THIS WAS MISSING!)

bend_refl_flux = mp.get_fluxes(refl)
bend_tran_flux = mp.get_fluxes(tran)

Visualize the geometry

plt.figure()
sim.plot2D(output_plane=mp.Volume(center=mp.Vector3(), size=mp.Vector3(sx, sy)))
plt.title("2D simulation with Silver Nanoparticles")
plt.savefig("code_6_geometry_2D.png")
plt.show()

Calculate reflectance and transmittance

flux_freqs = mp.get_flux_freqs(refl)
wl = []
Rs = []
Ts = []

Calculate reflectance and transmittance

flux_freqs = mp.get_flux_freqs(refl)
wl = []
Rs = []
Ts = []
for i in range(nfreq):
wl = np.append(wl, 1 / flux_freqs[i])
Rs = np.append(Rs, -bend_refl_flux[i] / straight_tran_flux[i])
A = 1 - Rs[i] - bend_tran_flux[i] / straight_tran_flux[i] # Absorption
Ts = np.append(Ts, 1 - A) # T = 1 - A

if mp.am_master():
plt.figure()
plt.plot(wl, Rs, "bo-", label="reflectance")
plt.plot(wl, Ts, "ro-", label="transmittance (1-A)")
plt.plot(wl, 1 - Rs - Ts, "go-", label="absorption")
plt.axis([0.2, 1.0, 0, 1])
plt.xlabel("wavelength (μm)")
plt.legend(loc="upper right")
plt.savefig("results_RTA_code_6.png")
plt.show()