# Material Dispersion

In this example, we will perform a simulation with a frequency-dependent dielectric $\varepsilon(\omega)$, corresponding to material dispersion. See Materials for more information on how material dispersion is supported. In particular, we will model a uniform medium of the dispersive material. See also the material-dispersion.ctl file included with Meep. From the dispersion relation $\omega(k)$, we will compute the numerical $\varepsilon(\omega)$ via the formula:

We will then compare this with the analytical $\varepsilon(\omega)$ that we specified.

Since this is a uniform medium, our computational cell can actually be of zero size (i.e. one pixel), where we will use Bloch-periodic boundary conditions to specify the wavevector k.

(set! geometry-lattice (make lattice (size no-size no-size no-size)))
(set-param! resolution 20)


We will then fill all space with a dispersive material:

(set! default-material
(make dielectric (epsilon 2.25)
(E-susceptibilities
(make lorentzian-susceptibility
(frequency 1.1) (gamma 1e-5) (sigma 0.5))
(make lorentzian-susceptibility
(frequency 0.5) (gamma 0.1) (sigma 2e-5))
)))


corresponding to the dielectric function:

The real and imaginary parts of this dielectric function $\varepsilon(\omega)$ are plotted below:

We can see that the f=1.1 resonance causes a large change in both the real and imaginary parts of $\varepsilon$ around that frequency. In fact, there is a range of frequencies from 1.1 to 1.2161 where $\varepsilon$ is negative. In this range, no propagating modes exist — it is actually a kind of electromagnetic band gap associated with polariton resonances in a material. For more information on the physics of such materials, see e.g. chapter 14 of Introduction to Solid State Physics by C. Kittel.

On the other hand, the f=0.5 resonance, because the sigma numerator is so small, causes very little change in the real part of $\varepsilon$. Nevertheless, it generates a clear peak in the imaginary part of $\varepsilon$, corresponding to a resonant absorption peak.

Now, we'll set up the rest of the simulation. We'll specify a broadband $E_z$-polarized Gaussian source, create a list of k wavevectors that we want to compute $\omega(k)$ over, and compute the associated frequencies by using the run-k-points function:

(define-param fcen 1.0)
(define-param df 2.0)
(set! sources (list (make source
(src (make gaussian-src (frequency fcen) (fwidth df)))
(component Ez) (center 0 0 0))))
(define-param kmin 0.3)
(define-param kmax 2.2)
(define-param k-interp 99)
(define kpts (interpolate k-interp (list (vector3 kmin) (vector3 kmax))))
(define all-freqs (run-k-points 200 kpts)) ; a list of lists of frequencies


The run-k-points function returns a list of lists of frequencies — one list of complex frequencies for each k point — which we store in the all-freqs variable. Finally, we want to loop over this list and print out the corresponding $\varepsilon$ via the ratio $(ck/\omega)^2$ as described above. To do this, we will use the Scheme map function, which applies a given function to every element of a list (or lists), and since we have a list of lists we'll actually nest two map functions:

(map (lambda (kx fs)
(map (lambda (f)
(print "eps:, " (real-part f) ", " (imag-part f)
", " (sqr (/ kx f)) "\n"))
fs))
(map vector3-x kpts) all-freqs)


Alternatively we could just read all of the frequencies into Octave/Matlab or a spreadsheet and compute the ratios there. After running the program with

unix% meep material-dispersion.ctl | tee material-dispersion.out


we can then grep for the frequencies and the computed dielectric function, and plot it. First, let's plot the dispersion relation $\omega(k)$ for the real part of $\omega$:

The red circles are the computed points from Meep, whereas the blue line is the analytical band diagram from the specified $\varepsilon(\omega)$. As you can see, we get two bands at each k, separated by a polaritonic gap (shaded yellow). This dispersion relation can be thought of as the interaction (anti-crossing) between the light line of the ambient $\varepsilon$=2.25 material (dashed black line) and the horizontal line corresponding to the phonon resonance.

Similarly, the computed and analytical real parts of the dielectric function are given by:

which shows excellent agreement between the analytical (blue line) and numerical (red circles) calculations. The imaginary part, however, is more subtle:

The blue line is the analytical calculation from above and the red circles are the numerical value from Meep — why is the agreement so poor? There is nothing wrong with Meep, and this is not a numerical error. The problem is simply that we are comparing apples and oranges.

The blue line is the analytical calculation of $\varepsilon(\omega)$ for a real frequency $\omega$ which corresponds to solutions with a complex wavevector k, whereas Meep is computing $\varepsilon$ at a complex $\omega$ for a real wavevector k. So, the correct comparison is to plug Meep's complex $\omega$ into the analytical formula for $\varepsilon(\omega)$, which results in the green lines on the graph that fall almost on top of the red circles.

Why did our comparison of the real part of $\varepsilon$ look so good, then? The reason is that $\varepsilon(\omega)$ at real and complex values of $\omega$ are closely related by the analytic properties of $\varepsilon$. In particular, because $\varepsilon$ is an analytic function on the real-$\omega$ axis, adding a small imaginary part to $\omega$ as we are doing here does not change $\varepsilon$ by much. The losses are small for all of the computed k points. The change was only significant for the imaginary $\varepsilon$ because the imaginary $\varepsilon$ was small to begin with.