Plasmonic materials in MEEP

 

The aim of this post is to share my experience in incorporating dielectric function of metals such as gold and silver into MEEP (a free finite difference time domain package) code. The incorporation is not an easy task and can be daunting for the first time user.

Metals such as gold and silver have both Drude and Lorentz components for the dielectric function. There are many forms of Lorentz-Drude expressions in literature with slight notation differences.

I prefer the Lorentz-Drude expression mentioned in Rakic et al., Optical properties of metallic films for vertical-cavity optoelectronic devices, Applied Optics (1998) and is given by:

\epsilon(\omega)=1-\frac{f_1\omega_p'^2}{(\omega^2+i\Gamma_1'\omega)}+\sum_{j=2}^{n}\frac{f_j\omega_p'^2}{(\omega_{o,j}'^2-\omega^2-i\Gamma_j'\omega)}.

I took it as reference as they also give the parameters for metals such as Ag and Au and these models fit decently well to experimental data over a larger frequency range (0.1 eV to 6 eV).

For Silver, the parameters are :

\omega_p = 9.01;% Plasma frequency in eV

f = [0.845 0.065 0.124 0.011 0.840 5.646];% Oscillators’ strengths

\Gamma = [0.048 3.886 0.452 0.065 0.916 2.419];% Damping frequency of each oscillator in eV

\omega_o = [0.000 0.816 4.481 8.185 9.083 20.29];% Resonant frequency of each oscillator in eV

For Gold, the parameters are :

\omega_p = 9.03;% Plasma frequency in eV

f = [0.760 0.024 0.010 0.071 0.601 4.384];% Oscillators’ strengths

\Gamma = [0.053 0.241 0.345 0.870 2.494 2.214];% Damping frequency of each oscillator in eV

\omega_o = [0.000 0.415 0.830 2.969 4.304 13.32];% Resonant frequency of each oscillator in eV

To convert these parameters usable for MEEP codes, one needs to do the following:

1) Normalization:

\omega_p'=\omega_p \times \textrm{norm},

\Gamma'=\Gamma \times \textrm{norm},

\omega_o'=\omega_o \times \textrm{norm},

where \textrm{norm}=\frac{ev2w}{\frac{2\pi c}{a}}, ev2w=\frac{2\pi}{4.135666e-15} .

Essentially, we are converting the units of \omega_p,\Gamma and \omega_o from eV to angular frequency, which is done by multiplying them with ev2w. Further, they need to be normalized to MEEP units, divide by \frac{2\pi c}{a}, where a is the length unit in meters used in your MEEP code and obviously c is 299792458 m/s.

2) Casting Lorentz-Drude model to MEEP’s model

Current version of MEEP (1.1.1), supports dispersive materials through Lorentz model and is given by:

\epsilon(\omega)=\epsilon_{\infty}+\sum_{j=1}^{n}\frac{\sigma_{L,j}\omega_{L,j}^2} {(\omega_{L,j}^2-\omega^2-i\Gamma_{L,j}\omega)}.

The normalized parameters, \omega_p',\Gamma' and \omega_o' from step 1, need to be cast into the MEEP’s dispersive materials model by the following modifications:

1) \omega_{L,1}=1e^{-20} (Any small number should be fine.)

2) \omega_{L,j>1}=\omega_{o,j>1}'

3) \sigma_{L,j}=\frac{f_j\omega_p'^2}{\omega_{L,j}^2}

4) \epsilon_{\infty}=1

5) \Gamma_{L,j}=\Gamma_j'

Here is scheme code snippet for silver with a=100e-9m (100 nm) after normalization and casting them into MEEP’s material model.

(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))
)))

The first polarizability is the Drude component and the rest 5 polarzabilities represent the Lorentz components. I have used this snippet in a MEEP code that calculates the reflection from a thin film of silver (to be consistent, all my units in this MEEP code are also normalized to 100 nm). Later I compared the reflection spectra with an analytical solution. These results are shown below. The MEEP calculation results match very well the analytical solutions. One can also clearly see the Ag interband effects in UV regime. The results also match with an other reference.

If anyone needs the codes that generate the scheme snippet for material polarization and calculation of reflectivity of thin metal films, please shoot me an email.

Reflectivity of Ag thin film (MEEP vs. analytical model)

You can find my project file here. Please read the readme.txt for further instructions.

UPDATE: Aaron Webster used this code for a 1 × 1 × 1 pixel region and epilion was exported with the command meep-fields-analytic-chi1. He showed that the dielectric function seen by meep through using my snippet matches with the experimental dielectric function. You can find his report here. Thank you Aaron.

Majid has written a python code to do the operation described above. Please find it here

22 thoughts on “Plasmonic materials in MEEP

  1. Hi,

    Thank you for your exciting post.
    I am very interested in how MEEP deals with metals.
    Could you send me the code?

    Many thanks in advance,

    Khai

  2. Hi Bala,

    I am very interested in this code, too. Could you send me the source code? Do you think, that MEEP is the right software for computing plasmonic nanoparticles? Are there any disadvantages?

    thank you

    Stefan

  3. Thank you for your interesting article.
    I find your parameters very useful, too. Could you send me the code?
    Thank you
    myriad

  4. Hi Bala Khrishna,

    Thanks for this very interesting example. I would like to try out the metal option in MEEP. Could you send me the code too?

    Kind regards,

    Rainer

  5. Hi,

    I have also tried to define Drude-Lorentz model for Au as follows:

    ;*************************************************************************************
    ; Define silver using Drude-Lorentz model
    (define-param c 3e+14)
    (define norm (/ a (* c 4.135666e-15)))
    (define omega_p (* 9.03 norm))

    (define-param omega1 1e-20)
    (define omega2 (* 0.415 norm))
    (define omega3 (* 0.83 norm))
    (define omega4 (* 2.969 norm))
    (define omega5 (* 4.304 norm))
    (define omega6 (* 13.33 norm))

    (define gamma1 (* 0.053 norm))
    (define gamma2 (* 0.241 norm))
    (define gamma3 (* 0.345 norm))
    (define gamma4 (* 0.87 norm))
    (define gamma5 (* 2.494 norm))
    (define gamma6 (* 2.214 norm))

    (define-param f1 0.76)
    (define-param f2 0.024)
    (define-param f3 0.01)
    (define-param f4 0.071)
    (define-param f5 0.601)
    (define-param f6 4.384)

    (define sigma1 (/ (* f1 (* omega_p omega_p)) (* omega1 omega1)))
    (define sigma2 (/ (* f2 (* omega_p omega_p)) (* omega2 omega2)))
    (define sigma3 (/ (* f3 (* omega_p omega_p)) (* omega3 omega3)))
    (define sigma4 (/ (* f4 (* omega_p omega_p)) (* omega4 omega4)))
    (define sigma5 (/ (* f5 (* omega_p omega_p)) (* omega5 omega5)))
    (define sigma6 (/ (* f6 (* omega_p omega_p)) (* omega6 omega6)))

    (define Au (make dielectric (epsilon 1)
    (polarizations
    (make polarizability
    (omega omega1) (gamma gamma1) (sigma sigma1))
    (make polarizability
    (omega omega2) (gamma gamma2) (sigma sigma2))
    (make polarizability
    (omega omega3) (gamma gamma3) (sigma sigma3))
    (make polarizability
    (omega omega4) (gamma gamma4) (sigma sigma4))
    (make polarizability
    (omega omega5) (gamma gamma5) (sigma sigma5))
    (make polarizability
    (omega omega6) (gamma gamma6) (sigma sigma6)))))
    ;*************************************************************************************

    However, in some cases MEEP is hard to find convergent. :-(

    Cheers,

    Khai

  6. Hi,
    I have also tried to define Drude-Lorentz model for Au as follows:

    ;*************************************************************************************
    ; Define gold using Drude-Lorentz model
    (define-param c 3e+14)
    (define norm (/ a (* c 4.135666e-15)))
    (define omega_p (* 9.03 norm))

    (define-param omega1 1e-20)
    (define omega2 (* 0.415 norm))
    (define omega3 (* 0.83 norm))
    (define omega4 (* 2.969 norm))
    (define omega5 (* 4.304 norm))
    (define omega6 (* 13.33 norm))

    (define gamma1 (* 0.053 norm))
    (define gamma2 (* 0.241 norm))
    (define gamma3 (* 0.345 norm))
    (define gamma4 (* 0.87 norm))
    (define gamma5 (* 2.494 norm))
    (define gamma6 (* 2.214 norm))

    (define-param f1 0.76)
    (define-param f2 0.024)
    (define-param f3 0.01)
    (define-param f4 0.071)
    (define-param f5 0.601)
    (define-param f6 4.384)

    (define sigma1 (/ (* f1 (* omega_p omega_p)) (* omega1 omega1)))
    (define sigma2 (/ (* f2 (* omega_p omega_p)) (* omega2 omega2)))
    (define sigma3 (/ (* f3 (* omega_p omega_p)) (* omega3 omega3)))
    (define sigma4 (/ (* f4 (* omega_p omega_p)) (* omega4 omega4)))
    (define sigma5 (/ (* f5 (* omega_p omega_p)) (* omega5 omega5)))
    (define sigma6 (/ (* f6 (* omega_p omega_p)) (* omega6 omega6)))

    (define Au (make dielectric (epsilon 1)
    (polarizations
    (make polarizability
    (omega omega1) (gamma gamma1) (sigma sigma1))
    (make polarizability
    (omega omega2) (gamma gamma2) (sigma sigma2))
    (make polarizability
    (omega omega3) (gamma gamma3) (sigma sigma3))
    (make polarizability
    (omega omega4) (gamma gamma4) (sigma sigma4))
    (make polarizability
    (omega omega5) (gamma gamma5) (sigma sigma5))
    (make polarizability
    (omega omega6) (gamma gamma6) (sigma sigma6)))))
    ;*************************************************************************************

    However, in some cases MEEP is hard to find convergent.
    Cheers,
    Khai

  7. Hi

    I admire you for your creative mind and the way you look at problems.

    I would like to have your code for reaching to this result.

    Thanks a lot.
    Mehdi

  8. Are u referring to expression 3? if that is case, no need to prime the omega in denominator. because that is already in meep units (see expression 2).

  9. Is there any straightforward way to extend this to other materials that have plasmonic properties like Magnesium? Getting the values for oscillator strengths, etc, all the parameter you would need for the Lorentz-Drude model does not seem to be something that is very easy to obtain. Is there a good source or database for these values?

  10. Hi,
    I’ve found from the ctl file for Ag that the Gamma for the lorentz-drude model is 4 times higher than Gamma for the drude model. Here the example from cylinder_silver2D.ctl:
    (define myAg (make dielectric (epsilon 1.0001)
    (polarizations
    (make polarizability
    (omega 1e-20) (gamma 0.0038715) (sigma 4.4625e+39))
    (make polarizability
    )))
    (define myAg_drude1
    (make dielectric (epsilon 8.926)
    (polarizations
    (make polarizability
    (omega 1e-20) (gamma 0.016373) (sigma 8.7309e+39))
    )))
    Any comments

    With regards

Leave a Reply

Your email address will not be published. Required fields are marked *

You may use these HTML tags and attributes: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong>