Model thermal radiation with white-noise source

[NatkamonPS]

Hello, I've tried to calculate thermal radiation flux between closely spaced infinite parallel Doped-silicon plates in 2-D(with the Drude model) using a white-noise source as described in this tutorial. Sources are placed with a density of 10 per 1um in both x and y directions.

The reference flux spectrum is the emission from flat surface(lower plate) with no upper plate.

The code is as shown below

import meep as mp
import numpy as np
import matplotlib.pyplot as plt
import math
import argparse


parser = argparse.ArgumentParser()
parser.add_argument('-bb', type=float, default=1, help='batch number')
parser.add_argument('-kk', type=float, default=0, help='paralell wavevector')
args = parser.parse_args()
b = args.bb
k = args.kk


um_scale = 1.0
resolution = 50

dpml = 4.0
pad = 4.0
dtop = 4.5
dmid = 0.5
dbot = 4.5

sy = dpml+pad+dtop+dmid+dbot+pad+dpml
sx = 0.5

cell_size = mp.Vector3(sx,sy)
pml_layers = [mp.PML(direction=mp.Y,
                    thickness=dpml)]


fcen = 0.11
df = 0.2
nfreq = 201
ndipoley = 90 
ndipolex = 10 
ntrial = 15 # no. of trial per batch
run_time = 2*nfreq/df


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 

susc = [mp.DrudeSusceptibility(frequency=frq0, gamma=gam0, sigma=1)]
DSi = mp.Medium(epsilon=11.7, E_susceptibilities=susc)

def compute_flux(kx):
    geometry = [
            mp.Block(size=mp.Vector3(mp.inf,dtop,mp.inf),
                        center=mp.Vector3(0,dtop/2+dmid/2,0),
                        material=DSi),
            mp.Block(size=mp.Vector3(mp.inf,dmid,mp.inf),
                        center=mp.Vector3(),
                        material=mp.Medium(epsilon=1)),
            mp.Block(size=mp.Vector3(mp.inf,dbot,mp.inf),
                        center=mp.Vector3(0,-(dmid/2+dbot/2),0),
                        material=DSi)]
    sources = []

    for i in range(ndipolex+1):
        for j in range(ndipoley+1):
            sources.append(mp.Source(mp.CustomSource(src_func=lambda t: np.random.randn()),
                                        component=mp.Ex,
                                        center=mp.Vector3(-(sx)*(-0.5+i/ndipolex),-(dmid/2)-dbot*(j/ndipoley))))

    sim = mp.Simulation(cell_size=cell_size,
                        boundary_layers=pml_layers,
                        eps_averaging = False,
                        k_point=mp.Vector3(kx,0,0),
                        geometry=geometry,
                        sources=sources,
                        resolution=resolution)

    flux_mon1 = sim.add_flux(fcen,
                            df,
                            nfreq,
                            mp.FluxRegion(center=mp.Vector3(),size=mp.Vector3(sx,0,0)))

    sim.run(until=run_time)

    flux1 = mp.get_fluxes(flux_mon1)
    freqs = mp.get_flux_freqs(flux_mon1)

    sim.reset_meep()
    geometry=[mp.Block(size=mp.Vector3(mp.inf,dbot,mp.inf),
        center=mp.Vector3(0,-(dmid/2)-(dbot/2),0),
        material=DSi)]

    sim = mp.Simulation(cell_size=cell_size,
                        boundary_layers=pml_layers,
                        eps_averaging = False,
                        k_point=mp.Vector3(kx,0,0),
                        geometry=geometry,
                        sources=sources,
                        resolution=resolution)

    input = sim.add_flux(fcen,
                        df,
                        nfreq,
                        mp.FluxRegion(center=mp.Vector3(),size=mp.Vector3(sx,0,0)))

    sim.run(until=run_time)

    input_flux = mp.get_fluxes(input)


    return freqs, flux1, input_flux

fluxes = np.zeros((nfreq,ntrial))
incident_fluxes = np.zeros((nfreq,ntrial))
for t in range(ntrial):
    print("############################")
    print("kx = {} : trial no. {} of {}".format(k,t,ntrial))
    print("############################")
    freqs, fluxes[:,t], input_fluxes[:,t] = compute_flux(k)

data = np.zeros((nfreq,ntrial+1))
data[:,0] = freqs
data[:,1:] = fluxes

data_ff = np.zeros((nfreq,ntrial+1))
data_ff[:,0] = freqs
data_ff[:,1:] = input_fluxes

if mp.am_master():
	np.savetxt("/DSi_Ex_kx{}_{}.txt".format(k,b),data,delimiter=",")
	np.savetxt("/DSi_Ex_kx{}_FF_{}.txt".format(k,b),data_ff,delimiter=",")

I've compared results(averaged over 60 trials) with the solution obtained from analytical green's function at each parallel wavevector(kx), which is in the unit of kx. The shape of the flux is similar (same location of peak values)

I'm not sure that :
(1) how to convert the result at each wavevector from meep to the unit that is comparable to the exact solution? I've tried to use the reference spectrum at the same value of kx but the reference spectrum has values very close to 0 in the frequency range of evanescent mode which blow up when used as a normalization spectrum.

(2) In this paper, they normalized the resulting flux by the flux spectrum in the far-field. Is this the same as calculating the emission from the lower plate with no upper plate at each kx and then integrating over all directions?

(3) Since the polarizations are uncorrelated in this case. For TM polarization(Hz, Ex, Ey polarized), I've to calculate the flux using mp.Ex and mp.Ey separately then sum them together in the post-processing. Is this correct?

[alexrod7]

Hello NatkamonPS,

(1) Often what you want is to normalize by the reference spectrum once you have integrated over "kx", in which case there are no issues with the evanescent window. Because the per-k spectrum is not observable outside of the evanescent window (you would only want to consider its contribution in the near-field regime), this should not be an issue. See my point #2 below.

(2) The results of that paper normalize the spectrum as a function of frequency - the flux has already been integrated over "kx". This isnot the same thing as normalizing the per-k spectrum and then integrating over k.

(3) In principle you need three different calculations (Ex, Ey, and Hz sources).

[NatkamonPS]

Dear alexrod7,

Thank you very much for your response.

Does this mean for each calculation (Ex, Ey, Hz sources), I have to calculate (1) flux spectrum in the case of parallel plates at each kx (2) flux spectrum in the case of single flat plate at each kx. Then, normalize the integrated spectrum in (1) by the integrated spectrum in (2)?

Also, in my understanding, the meaning of this normalized flux is the enhancement/suppression of parallel plate geometry compare to the far-field emissivity of the single plate. Is this correct?

[alexrod7]

Hi NatkamonPS,

Does this mean for each calculation (Ex, Ey, Hz sources), I have to calculate (1) flux spectrum in the case of parallel plates at each kx (2) flux spectrum in the case of single flat plate at each kx. Then, normalize the integrated spectrum in (1) by the integrated spectrum in (2)?

Yes

Also, in my understanding, the meaning of this normalized flux is the enhancement/suppression of parallel plate geometry compare to the far-field emissivity of the single plate. Is this correct?

That is correct.