# Material Dispersion

In this example, we will perform a simulation with a **frequency-dependent dielectric** , 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 material-dispersion.ctl. From the dispersion relation , we will compute the numerical via the formula:

We will then compare this with the analytical 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 are plotted below:

We can see that the f=1.1 resonance causes a large change in both the real and imaginary parts of around that frequency. In fact, there is a range of frequencies from 1.1 to 1.2161 where 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 . Nevertheless, it generates a clear peak in the *imaginary* part of , corresponding to a resonant absorption peak.

Now, we'll set up the rest of the simulation. We'll specify a broadband -polarized Gaussian source, create a list of *k* wavevectors that we want to compute 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 via the ratio 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 for the real part of :

The red circles are the computed points from Meep, whereas the blue line is the analytical band diagram from the specified . 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 =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 for a *real* frequency which corresponds to solutions with a *complex* wavevector *k*, whereas Meep is computing at a *complex* for a *real* wavevector *k*. So, the correct comparison is to plug Meep's *complex* into the analytical formula for , 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 look so good, then? The reason is that at real and complex values of are closely related by the analytic properties of . In particular, because is an analytic function on the real- axis, adding a *small* imaginary part to as we are doing here does not change by much. The losses are small for all of the computed *k* points. The change was only significant for the imaginary because the imaginary was small to begin with.