Wrong results for scattering from a gold nanoparticle

[Dev-linkon]

Greetings!!

I want to study the scattering from a gold nanoparticle. I followed the steps given in the section "Mie Scattering of a Lossless Dielectric Sphere" in the MEEP tutorial (https://meep.readthedocs.io/en/latest/Python_Tutorials/Basics/#mie-scattering-of-a-lossless-dielectric-sphere).
The python code is as given below and the result I get is in the image attached. Please help me correct my code to get the correct result.

`import meep as mp
import numpy as np
import matplotlib.pyplot as plt
#---------------------------------------------------------------------------------------------------------#
um_scale = 1.0
eV_um_scale = um_scale/1.23984193
metal_visible_range = mp.FreqRange(min=um_scale/0.8, max=um_scale/0.4)
Au_JC_visible_frq0 = 1/(2.68731967um_scale)
Au_JC_visible_gam0 = 1/(22.43470716um_scale)
Au_JC_visible_sig0 = 1.09954246
Au_JC_visible_frq1 = 1/(0.71551082um_scale)
Au_JC_visible_gam1 = 1/(32.25215811um_scale)
Au_JC_visible_sig1 = 0
Au_JC_visible_susc = [mp.DrudeSusceptibility(frequency=Au_JC_visible_frq0, gamma=Au_JC_visible_gam0, sigma=Au_JC_visible_sig0),
mp.LorentzianSusceptibility(frequency=Au_JC_visible_frq1, gamma=Au_JC_visible_gam1, sigma=Au_JC_visible_sig1)]
Au_JC_visible = mp.Medium(epsilon=5.6902532, E_susceptibilities=Au_JC_visible_susc)
#---------------------------------------------------------------------------------------------------------#

r = 0.025 # radius of PVA
wvl_min = 0.4
wvl_max = 0.8
frq_min = 1/wvl_max
frq_max = 1/wvl_min
frq_cen = 0.5*(frq_min+frq_max)
dfrq = frq_max-frq_min
nfreq = 200

resolution = 250

dpml = 0.25wvl_min
dair = 0.25wvl_min

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

symmetries = [mp.Mirror(mp.Y),
mp.Mirror(mp.Z,phase=-1)]

s = 2*(dpml+dair+r)
cell_size = mp.Vector3(s,s,s)

sources = [mp.Source(mp.GaussianSource(frq_cen,fwidth=dfrq,is_integrated=True),
center=mp.Vector3(-0.5*s+dpml,0,0),
size=mp.Vector3(0,s,s),
component=mp.Ez)]
Courant = 0.2
sim = mp.Simulation(resolution=resolution,
cell_size=cell_size,
boundary_layers=pml_layers,
sources=sources,
k_point=mp.Vector3(),
symmetries=symmetries,
Courant=Courant)

dpad = 0.5dair
box_x1 = sim.add_flux(frq_cen, dfrq, nfreq, mp.FluxRegion(center=mp.Vector3(x=-r-dpad),size=mp.Vector3(0,2(r+dpad),2*(r+dpad))))
box_x2 = sim.add_flux(frq_cen, dfrq, nfreq, mp.FluxRegion(center=mp.Vector3(x=+r+dpad),size=mp.Vector3(0,2*(r+dpad),2*(r+dpad))))
box_y1 = sim.add_flux(frq_cen, dfrq, nfreq, mp.FluxRegion(center=mp.Vector3(y=-r-dpad),size=mp.Vector3(2*(r+dpad),0,2*(r+dpad))))
box_y2 = sim.add_flux(frq_cen, dfrq, nfreq, mp.FluxRegion(center=mp.Vector3(y=+r+dpad),size=mp.Vector3(2*(r+dpad),0,2*(r+dpad))))
box_z1 = sim.add_flux(frq_cen, dfrq, nfreq, mp.FluxRegion(center=mp.Vector3(z=-r-dpad),size=mp.Vector3(2*(r+dpad),2*(r+dpad),0)))
box_z2 = sim.add_flux(frq_cen, dfrq, nfreq, mp.FluxRegion(center=mp.Vector3(z=+r+dpad),size=mp.Vector3(2*(r+dpad),2*(r+dpad),0)))

sim.run(until_after_sources=80)

freqs = mp.get_flux_freqs(box_x1)
box_x1_data = sim.get_flux_data(box_x1)
box_x2_data = sim.get_flux_data(box_x2)
box_y1_data = sim.get_flux_data(box_y1)
box_y2_data = sim.get_flux_data(box_y2)
box_z1_data = sim.get_flux_data(box_z1)
box_z2_data = sim.get_flux_data(box_z2)

box_x1_flux0 = mp.get_fluxes(box_x1)

sim.reset_meep()

Courant = 0.02
``
geometry = [mp.Sphere(material=Au_JC_visible,
center=mp.Vector3(),
radius=r)]

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

box_x1 = sim.add_flux(frq_cen, dfrq, nfreq, mp.FluxRegion(center=mp.Vector3(x=-r-dpad),size=mp.Vector3(0,2*(r+dpad),2*(r+dpad))))
box_x2 = sim.add_flux(frq_cen, dfrq, nfreq, mp.FluxRegion(center=mp.Vector3(x=+r+dpad),size=mp.Vector3(0,2*(r+dpad),2*(r+dpad))))
box_y1 = sim.add_flux(frq_cen, dfrq, nfreq, mp.FluxRegion(center=mp.Vector3(y=-r-dpad),size=mp.Vector3(2*(r+dpad),0,2*(r+dpad))))
box_y2 = sim.add_flux(frq_cen, dfrq, nfreq, mp.FluxRegion(center=mp.Vector3(y=+r+dpad),size=mp.Vector3(2*(r+dpad),0,2*(r+dpad))))
box_z1 = sim.add_flux(frq_cen, dfrq, nfreq, mp.FluxRegion(center=mp.Vector3(z=-r-dpad),size=mp.Vector3(2*(r+dpad),2*(r+dpad),0)))
box_z2 = sim.add_flux(frq_cen, dfrq, nfreq, mp.FluxRegion(center=mp.Vector3(z=+r+dpad),size=mp.Vector3(2*(r+dpad),2*(r+dpad),0)))

sim.load_minus_flux_data(box_x1, box_x1_data)
sim.load_minus_flux_data(box_x2, box_x2_data)
sim.load_minus_flux_data(box_y1, box_y1_data)
sim.load_minus_flux_data(box_y2, box_y2_data)
sim.load_minus_flux_data(box_z1, box_z1_data)
sim.load_minus_flux_data(box_z2, box_z2_data)

sim.run(until_after_sources=80)

box_x1_flux = mp.get_fluxes(box_x1)
box_x2_flux = mp.get_fluxes(box_x2)
box_y1_flux = mp.get_fluxes(box_y1)
box_y2_flux = mp.get_fluxes(box_y2)
box_z1_flux = mp.get_fluxes(box_z1)
box_z2_flux = mp.get_fluxes(box_z2)

scatt_flux = np.array(box_x1_flux)-np.array(box_x2_flux)+np.array(box_y1_flux)-np.array(box_y2_flux)+np.array(box_z1_flux)-np.array(box_z2_flux)
intensity = np.array(box_x1_flux0)/(2r)*2
scatt_cross_section = np.divide(scatt_flux,intensity)
#scatt_eff = scatt_cross_section(-1/(np.pir**2))

wl = []
SC = []
for i in range(nfreq):
wl = np.append(wl, 1000/freqs[i]) # freq to wavelength (in nm)
SC = np.append(SC, scatt_cross_section[i])
loss_data = np.zeros((nfreq, 2))
loss_data[:,0] = wl
loss_data[:,1] = SC
np.savetxt(f"scattering-AuNP-r{r}-res{resolution}-c{Courant}-dpml{dpml}-dair{dair}", loss_data)

if mp.am_master():
plt.figure()
plt.plot(wl,Ls,'g-',linewidth=3.5,label='AuNP')
plt.xticks(fontsize= 14)
plt.yticks(fontsize= 14)
plt.xlabel("wavelength (nm)", weight="bold", size=16, labelpad=8)
plt.ylabel("Scattering efficiency", weight="bold", size=16, labelpad=9)
plt.legend()
plt.grid(color = 'grey', linestyle = '--', linewidth = 0.2)
plt.title('Mie Scattering of AuNP', weight="bold", size=18)
plt.savefig(f'Au_SC-r{r}-res{resolution}-c{Courant}-dpml{dpml}-dair{dair}.png', dpi=220, bbox_inches='tight')`

Aprt from this, I want to know whether the nanoparticles can be simulated in MEEP to see the transmittance, reflectance and loss (=1-refectance-transmittance)? If yes, will I get the peaks for loss at same wavelengths as in the scattering cross-section or extonction cross-section? Please help me with the code as I am also struggling to do this as well.

Thank You!!

[oskooi]

This seems to be a duplicate of #1530.

In general, computing field-derived quantities which are converged in the presence of surface plasmons is challenging due to discretization effects as explained in #1530 (comment).