Questions on mp.NoisyDrudeSusceptibility

[NatkamonPS]

Hello again,

I tried model thermal radiation by implementing the method described in this paper. In the paper, they added white noise term to the polarization equation to obtain normalized flux spectrum which can be convert into heat flux by multiplying the normalized flux by Planck's distribution in the post processing. As described in meep document, this can be implemented by (i) using a custom source function to input white-noise sources or (ii) using mp.NoisyDrudeSusceptibility/ mp.NoisyLorentzianSusceptibility.

Following the method in the paper, we must calibrate the flux spectrum by comparing the flux(white noise method) and the absorptance (by Kirchoff's law). my results are shown below.

geometry (flat slab with periodic boundary condition in x-direction all the results shown are calculated with sim.k_point = mp.Vector3(0,0,0)

Absorptance from flat slab(using mp.DrudeSusceptibility and the planewave source is mp.GaussianSource normalized by a case with the same source without the slab)

flux spectrum from flat slab(using mp.NoisyDrudeSusceptibility)

Then, the calibration constant can be obtained by comparing the absorptance and flux spectrum.

According to the paper, this constant can be used to convert the flux spectrum(white noise method) to be the emittance from any structure with the same material and discretization parameter. Therefore, I tried to use the same constant obtained from flat slab case to convert to flux spectrum from grating slab to be in the same unit with its absoptance.

The result seems to match with the paper description. However, there are some concerns that I'm still not sure.

(1) Could the calibration constant obtained from flat slab case be used for any cases with the same emitting-material(e.g. parallel flat plates with mp.NoisyDrudeSusceptibility as an emitter and mp.DrudeSusceptibility as a receiver )? and could this constant be used for flux spectrum calculated from other values of sim.k_point?

(2) I had tried to use the mp.NoisyDrudeSusceptibility without any sources specified and got the error RuntimeError: meep: allocate field components (by adding sources) before adding dft objects. Therefore, I specified the source with zero amplitude as shown belows.
How does meep handle this? Is the white-noise term only included in the polarization time-stepping equation in x-direction only(since I specified the source component to be mp.Ex)?

My parameters for mp.NoisyDrudeSusceptibility

plasma_frq = 2.0738e+15 #rad/s
gam0 = 1.3557e+14 #rad/s

frq0 = plasma_frq/(2*math.pi*2.998e+14)#1/um
gam0 = gam0/(2*math.pi*2.998e+14)#1/um
sig0 = (plasma_frq/(2*math.pi*2.998e+14))**2

noisy_susc = [mp.NoisyDrudeSusceptibility(frequency=1, gamma=gam0, sigma=sig0, noise_amp=12*gam0*resolution*resolution/(math.pi*(run_time/0.01)))]

noisy_DSi = mp.Medium(epsilon=11.7, E_susceptibilities=noisy_susc)

source

sources = [mp.Source(
                            mp.CustomSource(src_func=lambda t: 0*np.random.randn()),
                            component=mp.Ex,
                            center=mp.Vector3(0,-(dmid/2)-(2/resolution))
                            )]

(3) As stated in this FAQ and demonstrated in this tutorial, I'm convinced that adding the noise terms into the polarization equation and as current sources in the emitter should be equivalent(since both polarization(P) and current density(J) will be added to the time-stepping equation(Ampère's law) in the same manner). But, my result doesn't seem to be like that(different shape of flux spectrum).

from the same method as above

with mp.DrudeSusceptibility for an emitter(densities of sources described in the figure are for both x and y direction)

for the last two figures the material is modeled by

# for absorptance in 1st figure and fluxes in 2nd figure
susc = [mp.DrudeSusceptibility(frequency=1, gamma=0.3, sigma=10)]
Mat = mp.Medium(epsilon=1, E_susceptibilities=susc)

# for flux in 1st figure
noisy_susc = [mp.NoisyDrudeSusceptibility(frequency=1, gamma=0.3, sigma=10,noise_amp=12*0.3*resolution*resolution/(math.pi*(run_time/0.01)))]
noisy_Mat = mp.Medium(epsilon=1, E_susceptibilities=noisy_susc)

your response would be highly appreciated.
Thank you.