Transmittance over 1 of gold metasurface

[2213411860]

Recently I'm trying to reproduce the transmittance-wavelength curve from the paper "Frequency Selective Surface Bandpass Filters Applied to Thermophotovoltaic Generators，https://doi.org/10.1063/1.1841894" .

I use gaussian source and k-point (0,0,0). I calculated the flux power with and without metasurface in two runs and divided the first run by the second.
I compared the results(red line) with the paper and lumerical and found it is rather strange that some transmittances are over 1. (or totally wrong)

I would be appreciate if someone can help me out.
Here is my code:
`import matplotlib.pyplot as plt
import numpy as np
import meep as mp
from meep.materials import Au
import sys
l=0.55 # side length of the filter
w=0.2 # width of the cell
a=0.45
b=0.1
dpml = 1
resolution = 60
h1=3 #distance between source and filter
h2=-3 #distance between probe and filter
cell=mp.Vector3(l, l, 8)

#geometry defination
B1= mp.Block(
size=mp.Vector3(a, b, w),
center=mp.Vector3(0, 0, 0),
material=mp.Medium(epsilon=1),
)
B2= mp.Block(
size=mp.Vector3(b, a, w),
center=mp.Vector3(0, 0, 0),
material=mp.Medium(epsilon=1),
)
B3=mp.Block(
size=mp.Vector3(l, l, w),
center=mp.Vector3(0, 0, 0),
material=Au,
)
geometry1=[B3,B2,B1]
geometry2=[B1]

#source defination
w1=0.5
w2=2
fcen = (1/w1-1/w2)/2 # pulse center frequency
df = (1/w1-1/w2) # pulse width (in frequency)
nfreq = 100
k_point = mp.Vector3(0,0,0)
sources = [
mp.Source(
mp.GaussianSource(fcen, fwidth=df),
component=mp.Ex,
center=mp.Vector3(0, 0, h1),
size=mp.Vector3(l, l, 0),
)
]

#pml layer
pml_layers=[mp.PML(thickness=dpml,direction=mp.Z)]

#simulation1
sim = mp.Simulation(
cell_size=cell,
boundary_layers=pml_layers,
k_point=k_point,
sources=sources,
dimensions=3,
resolution=20,
)

tran_fr = mp.FluxRegion(
center=mp.Vector3(0, 0, h2), size=mp.Vector3(l, l, 0)
)
tran = sim.add_flux(fcen,df , nfreq, tran_fr)
pt = mp.Vector3(0, 0, h2) # Reference point

#round 1
sim.run(until_after_sources=mp.stop_when_fields_decayed(20, mp.Ex, pt, 1e-3))

no_cell_tran_flux = mp.get_fluxes(tran)
empty_data = sim.get_flux_data(tran)
sim.reset_meep()

#simulation2
sim = mp.Simulation(
cell_size=cell,
boundary_layers=pml_layers,
geometry=geometry1,
sources=sources,
k_point=k_point,
dimensions=3,
resolution=resolution,
eps_averaging = False
)
tran_fr = mp.FluxRegion(
center=mp.Vector3(0, 0, h2), size=mp.Vector3(1, 1, 0)
)
tran = sim.add_flux(fcen,df , nfreq, tran_fr)
pt = mp.Vector3(0, 0, h2) # Reference point
#round 2
sim.run(until_after_sources=mp.stop_when_fields_decayed(20, mp.Ex, pt, 1e-3))
with_cell_tran_flux = mp.get_fluxes(tran)

#print(with_cell_tran_flux/ no_cell_tran_flux)
op=open("output.txt","w")
sys.stdout=op
for i in range(nfreq):
print(with_cell_tran_flux[i],no_cell_tran_flux[i],with_cell_tran_flux[i] / no_cell_tran_flux[i])
op.close()
`

[stevengj]

Have you checked convergence with resolution? Metals at IR wavelengths often require extraordinarily high resolution in order to resolve the skin depth (which is important to surface-plasmon resonances). I would also be careful of the material model — you might want to fit your own Drude model based on the parameters cited in the paper, for example …

[2213411860]

Thanks for your reply.
I did reset resolution to 120 but only tiny change was seen compared to resolution=60 curve.
I also checked the epsilon in wavelength range 0.5-2um which seems quite normal.

As for Drude model, I think the original one from meep library could be OK for me.
Maybe I should set resolution even higher? I've been struggling for several days. If there is any other possible resolution to this problem, please just let me know.

[2213411860]

By the way, when I ran the simulation, "DFT frequency 0.0 is out of material's range of 0.16131113692089302-4.0327458966810505" was warned, maybe this had something to do with the problem?