Task 5: Sim for Parallel Plates

Now I'm just going to see if I can capture near-field enhancement between two plates. I'm going to set up simulations where I vary the distance between the two plates and get the flux between them. I'm just looking for trends here; I don't know how to work out the units yet.
Results:

I implemented SiC plates. One plate has the noisy polarizations. The other doesn't. In 1D there is no enhancement. There is a ~10x enhancement in 2D as can be seen in the figure. 1 unit on the abscissa is 1 micron.
BUT, after careful examination, this result comes from the flux captured from the gaussian pulse that I had put in (The random current function didn't start with exactly zero initial fields). Also, I'm getting some reflections which can be seen from the negative parts of spectrum. Now when thinking about heat transfer, the other side should be 0K and be totally absorbing. The reflections can be explained as follows in heat transfer terms: the other side is gaining energy and re-emitting back to the source.
I'm going to reduce the amplitude of the noise to see if the reflections are reduced relative to the flux of the source.

Task 4: Understanding the Noisy Lorentzian in MEEP

I need to understand the details of the mathematics involved in the noisy lorentzian susceptibility. My top question is how to define a temperature since the inputs to the command is some amplitude, an omega, a gamma, and a sigma.

I haven't understood this yet but I'm going ahead with the next task.

Task 3: Implement random currents in MEEP

To simulate the effect of temperature on charges, I need to define a random current source in MEEP. I've found the development version of MEEP has this feature. It is given as a noisy lorentzian susceptibility. The math is described in the reference below.

D. Chan, M. Soljačić, and J. Joannopoulos, “Direct calculation of thermal emission for three-dimensionally periodic photonic crystal slabs,” Physical Review E, vol. 74, no. 3, pp. 1-9, Sep. 2006.

Task 2: Surface Plasmon Dispersion

Dispersion curves are useful to understand transport. I set up a MEEP simulation of two surfaces of Ag separated by a 100e-9m gap. You can see 3 modes: surface plasmons, gap propagation, and bulk plasmons. ~0.7 is omega_p the plasma frequency.

(set! geometry-lattice (make lattice (size 20 no-size no-size)))
(set! pml-layers (list (make pml (thickness 1))))
(set-param! resolution 30)

(define myAg (make dielectric (epsilon 1)
(polarizations
(make polarizability
(omega 1e-20) (gamma 0.0038715) (sigma 4.4625e+39))
; (make polarizability
; (omega 0.065815) (gamma 0.31343) (sigma 7.9247))
; (make polarizability
; (omega 0.36142) (gamma 0.036456) (sigma 0.50133))
; (make polarizability
; (omega 0.66017) (gamma 0.0052426) (sigma 0.013329))
; (make polarizability
; (omega 0.73259) (gamma 0.07388) (sigma 0.82655))
; (make polarizability
; (omega 1.6365) (gamma 0.19511) (sigma 1.1133))
)))

(set! geometry (list
                (make block (center -3 0 0 0.0) (size 5 infinity infinity)
                      (material myAg ))
    (make block (center 3 0 0 0.0) (size 5 infinity infinity)
                      (material myAg ))
      ))

(set! sources (list
               (make source
                 (src (make gaussian-src (frequency 2.5) (fwidth 5)))
                 (component Ex)
                 (center 0 0 0))
    (make source
                 (src (make gaussian-src (frequency 2.5) (fwidth 5)))
                 (component Ex)
                 (center -.5 0 0))
    (make source
                 (src (make gaussian-src (frequency 2.5) (fwidth 5)))
                 (component Ex)
                 (center .5 0 0)) 
     ))

(define-param kmin 0.0)
(define-param kmax 4.0)
(define-param k-interp 50)
(define kpts (interpolate k-interp (list (vector3 0 kmin 0) (vector3 0 kmax 0))))
(define all-freqs (run-k-points 200 kpts)) ; a list of lists of frequencies  

And here's a python snippet that gives a list of k-vectors their frequencies by giving it the MEEP output as a file object

def getfreqs(meepoutfileobj):
    fd=[]
    for aline in meepoutfileobj:
        if 'freqs:' in aline:
            ld=aline.split(',')
            kpt=ld[2:5];freqs=ld[5:]
            pts=[( [float(an) for an in kpt]  ,float(afreq)) for afreq in freqs]
            fd.extend(pts)
    if fd==[]:return None
    return fd

Task 1: Understanding MEEP's dispersive media

http://juluribk.com/2011/04/27/plasmonic-materials-in-meep/
This page and the linked pdf explain it rather well. I've implemented it in python where the inputs are lorentz-drude model params and the output is MEEPs (lorentz) model.

#Author: Majid al-Dosari
#code that takes in lorentz-drude model parameters and converts them to code 
#usable by meep. just use eVLD2meepcode
#example: eVLD2meepcode('mysilver',100e-9,9,[[.845,.048,0],[.065,3.886,.816]])
#output:
#(define mysilver (make dielectric (epsilon 1) 
#(polarizations 
#(make polarizability 
#( omega 8.06554724676e-22 ) ( gamma 0.00387146267844 ) ( sigma 6.8445e+41 ))
#( omega 0.0658148655336 ) ( gamma 0.313427166009 ) ( sigma 7.9071150519 ))
#)))

#it's an implementation of
#http://juluribk.com/2011/04/27/plasmonic-materials-in-meep/ by Bala Krishna Juluri


#normalization length a
#w omega
#p plasma
#G Gamma
 

def sip(stuff):#stuff in parens
    return '( ' +str(stuff)+ ' )'
def wordandvalue(word,value):
    return sip(word+' '+str(value))

def meepfactor(a):
    return 2*3.14159265359*299792458/a #2pi*c/a

def eV2w(eV):return eV*(2*3.14159265359)/4.135666e-15

def eV2meep(a,eV):return eV2w(eV)/meepfactor(a)

def sigmaa(f,wp,w):
    """f is oscillator strength, wp is plamsa freq"""
    return (f*wp**2)/w**2

def eVLD2meepL(a,wp,fGw0_list):
    oGs_list=[]
    smallno=1e-20; 
    for f,G,w0 in fGw0_list:
        if w0==0:w0=smallno #hack
        oGs_list.append([eV2meep(a,w0)
        ,eV2meep(a,G)
        ,sigmaa(f,eV2meep(a,wp),eV2meep(a,w0))])
    return oGs_list

def eVLD2meepcode(name,a,wp,fGw0_list):
    """
    input:
        name: for code
        a: meep normalization const
        wp: plasma freq in eV
        fGw0_list: nested list of
            f oscillator strengths
            Gamma relaxation in eV
            omega0 in eV,w=0 is handled by converting it to a very small no.
    """
    return writelorentz(name,oGs_list=eVLD2meepL(a,wp,fGw0_list))


def writelorentz(name,oGs_list=None,epsilon=1):#meep is in lorentz
    txt=''
    txt+='(define '+str(name)+' (make dielectric (epsilon '+str(epsilon)+') \n'
    if oGs_list==None or len(oGs_list)==0:
        txt+='))'
        return txt
    txt+='(polarizations \n'
    for omega,Gamma,sigma in oGs_list:        
        txt+='(make polarizability \n'
        txt+=wordandvalue('omega',omega)+' '
        txt+=wordandvalue('gamma',Gamma)+' '
        txt+=wordandvalue('sigma',sigma)+')\n'
    txt+=')))'
    return txt

#so 1. normalize then 2. cast into meep