# How do we determine refractive index of a photonic crystal?

Normally, refractive index of a medium is defined as

$$n=\frac c{v_\text{p}},$$

where $$c$$ is the speed of light in vacuum, and $$v_\text{p}$$ is the phase speed of light in the medium. Phase speed is defined as

$$v_\text{p}=\frac{\omega}k,$$

where $$\omega$$ is the frequency, and $$k$$ is the wavenumber of light (which is the norm of the wavevector).

But in a photonic crystal we don't have a single well-defined wavevector: we have quasi-wavevector, which is defined up to a vector of the reciprocal lattice. So, how do we then determine the refractive index? Do we simply take the wavevector from the first Brillouin zone and calculate the phase velocity from it? (If yes, how can we justify this?) Or do we actually have multiple refractive indices so that a light beam splits into several beams as if by a diffraction grating?

I would not say that you have ill-defined wave vectors. As you mention, when doing an analysis in the First Brillouin Zone, we commonly obtain the wavenumber modulo $$\pi/(2 a)$$, being $$a$$ the lattice parameter.

For a bilayer material, we have the following dispersion relation

$$\cos(\kappa a) = \cos\left(\frac{\omega a}{2 c_1}\right) \cos\left(\frac{\omega a}{2 c_2}\right) - \frac{c_1^2 + c_2^2}{2 c_1 c_2}\sin\left(\frac{\omega a}{2 c_1}\right) \sin\left(\frac{\omega a}{2 c_2}\right)\, .$$

Plotting the solutions over the first Brillouin zone would obtain the following figure.

If instead, you invert the relationship and consider the shift for each branch you end up with the following figure.

Now, the refractive index is a material parameter and you would need to average to obtain it. Thus, it is a function of frequency and, in general, is anisotropic which leads to a tensor rather than a single scalar:

$$n = [n_{ij}(\omega)]\, .$$

• I agree with your statement about the index of refraction being a material parameter but I believe effective index of refraction is a better description than average index of refraction Jun 10, 2022 at 16:57
• @SteveSaban, when I'm saying "average", I'm referring to obtaining a bulk behavior. Not about computing the "arithmetic mean" of the properties (expected value). Jun 10, 2022 at 17:08

In photonic crystals, the dielectric function $$\varepsilon(r)$$ plays the role of confining potential for the light, so the refractive index $$n = \sqrt{\varepsilon(r)}$$ is not uniform at least along one dimension. And in general, the dispersion relation $$\omega(k)$$ will deviate from the linear dispersion relation, especially near the Brillouin zone boundary. But for a small wavevector, near the zone center, the dispersion relation will be linear, because the light does not "see" the potential variation and thus it only sees an effective average dielectric constant (see the discussion on page 56 on this book, which is freely available to public).

• This is all fine, but it still doesn't address the ambiguity of wavevector vs quasi-wavevector in relation to the phase speed. Jun 10, 2019 at 15:18
• @Ruslan, phase velocity is hard to define exactly because of the ambiguity you mentioned, even when the dispersion relation is approximately linear. You can add reciprocal lattice vector to $k$ and the slope does not change. It's much more helpful to think about group velocity in photonic crystals.
– wcc
Jun 10, 2019 at 15:28
• @Ruslan, your last sentence in OP also hints toward why there is no unique phase velocity in a photonic crystal. There is no unique wavefront to reference and measure the phase, as the multiply scattered waves coherently sum up to make a periodic wavefunction, just like that of electron in a crystalline solid.
– wcc
Jun 10, 2019 at 15:31
• But normal crystals (e.g. diamond) can also be viewed as photonic crystals: they are also periodic, they sum up multiply scattered waves as you say, but still they do have a well-defined phase speed, don't they? Jun 10, 2019 at 15:36
• @Ruslan, or do you mean diamonds have a well-defined dielectric constant in an average sense as I described? Any dielectric material has well-defined group velocity but you can't say the same for phase velocity. And I am not sure diamonds can be viewed as a photonic crystal - the carbon atoms are only angstroms apart, so the light mostly sees an average uniform dielectric material.
– wcc
Jun 10, 2019 at 16:05

In short, the effective refractive index is determined by the radius of curvature of the equi-frequency contour along the dominant direction of propagation. You can do this for the first Brillouin zone, but keep in kind that if you are working in the n-th frequency band, the dominant component of the Bloch mode tend to lie in the n-th Brillouin zone.

Here the effective refractive index is the one that best describes the "free-space" propagation of the optical field. If the equi-frequency contour are perfect circles, then the definitions is omni-directionally valid. However, often these contours are not circular. In that case one would use the radius of curvature along the direction of propagation.

• So, if I have a quasi-monochromatic beam entering from vacuum into a photonic crystal, which point of the Brillouin zone should I take to find the curvature of the equi-frequency contour? It's not quite obvious that the vacuum wavevector will directly correspond to the crystal quasiwavevector (in fact, I think it won't, since we'd generally have a superposition of crystal modes from different points in BZ and possibly different bands). Jul 1, 2020 at 7:45
• It depends on the specific situation, which can vary drastically due to the complex nature of the dispersion diagram (band structure). However, assuming the frequency only cuts through one band, it would determine the equi-frequency contours in that band. Then you need to match the transverse components of the wave vector on either side of the interface. This could be complicated depending on which contour you need to consider. When I have more time, I'll draw a diagram. Hopefully that would make it clear. Jul 1, 2020 at 13:05

## Setup

In this answer we'll use a concrete model as an illustration. It doesn't limit the generality of this discussion, but (hopefully) helps understanding.

In particular, we'll have a $$z$$-homogeneous (thus effectively 2D) non-magnetic $$(\mu=1)$$ crystal lit by a TE-polarized wave. The EM wave will then be described by a single equation for the $$E_z$$ component of the electric field, which we'll denote as simply $$E$$. We'll use the units where vacuum speed of light $$c=1.$$ The equation is then

$$\nabla^2 E=\varepsilon(\vec r)\partial^2_t E.\tag1$$

For a monochromatic wave with frequency $$\omega$$ we can write a time-independent equation:

$$\nabla^2 E=-\varepsilon(\vec r)\omega^2 E.\tag2$$

The concrete relative permittivity function $$\varepsilon(\vec r)$$ we'll take is

$$\varepsilon(\vec r)=1+(\varepsilon_{\text{max}}-1)\left(\sin\left(\frac{\pi x}a\right)\sin\left(\frac{\pi y}a\right)\right)^4.\tag3$$

We set the lattice period $$a=1$$ and peak permittivity $$\varepsilon_{\text{max}}=3.5$$.

The incident wave will come from the left at $$\varphi=10°$$ from the normal.

## Symmetry and conserved quantities

Let's model our vacuum-crystal system the following way. Take the relative permittivity $$\varepsilon(\vec r)$$ from the crystal, aligning a crystal symmetry plane with $$yOz$$ coordinate plane. This crystal has a particular band structure with dispersion relations being $$\omega_n(\vec K),$$ where $$\vec K$$ is the quasi-wavevector, and $$n$$ is the number of band.

Now replace all the values of $$\varepsilon(\vec r)$$ for all $$x<0$$ with $$1$$, the vacuum value. For our example crystal $$(3)$$ the resulting $$\varepsilon(x,y)$$ will look as follows.

The replacement of the LHS of the crystal with vacuum breaks the original translation symmetry of the crystal. But partial symmetry—along $$y$$ and $$z$$ directions—still remains: any translation by a lattice vector lying in $$yOz$$ plane will yield the same $$\varepsilon(\vec r)$$. This means that the corresponding components of $$\vec K$$, i.e. $$K_y$$ and $$K_z$$, are still conserved on transition from the vacuum to the crystal.

For our example crystal, the unit strip (of the translation symmetry in the $$y$$ direction) for the resulting system will look as follows.

## Reflected and transmitted modes

The remaining symmetry that we have lets us represent our full wavefunction as

$$E(\vec r)=u_{K_y,K_z}(\vec r)\exp(i(K_y y + K_z z)),\tag4$$

where $$u_{K_y,K_z}$$ is periodic in $$y$$ and $$z$$ per the Bloch theorem.

This means that to find the transmitted wave we only have to determine the $$K_x$$ component of the quasi-wavevector. To do this, we need to compute the band structure of the crystal—in particular, its cross section at the given values of $$K_y$$ and $$K_z$$.

Let's consider a wave incident on the crystal surface. It is described as

$$E_{\text{inc}}(\vec r)=E_0\exp(i\vec k\vec r),\tag5$$

$$\vec k=\omega\begin{pmatrix}\cos(\varphi)\\ \sin(\varphi)\\ 0\end{pmatrix}.\tag6$$

Taking $$\omega=4.45$$ and our values described above, we have

$$\vec k=\begin{pmatrix}4.3823945\\ 0.7727344\\ 0\end{pmatrix}.\tag7$$

The quasi-wavevector of the incident wave will then be

$$\vec K=\begin{pmatrix}-1.9007908\\ 0.7727344\\ 0\end{pmatrix}.\tag7$$

The band structure of the empty lattice (i.e. vacuum) with our example values at the slice of fixed $$K_y=0.7727344$$ will then look as follows.

From the intersection of $$\omega$$ shown by the gray line with the dispersion relation curves we can find the quasi-wavevector of the reflected wave. Only half of all the intersections correspond to the reflection, another half represent waves that propagate towards the crystal (one of which is the incident wave, others have zero amplitude).

Note how there are generally multiple possible reflected waves (e.g. at $$\omega=8$$). This is the result of diffraction on the periodic interface. If we do some calculations, we'll find that angles of this diffraction exactly follow the diffraction grating formula for different diffraction orders $$m$$:

$$\varphi_m^{\text{reflect}}=\pi-\arcsin\left(\sin\varphi-\frac{m\lambda}a\right).\tag8$$

Now, we are actually interested in the transmitted wave, rather than the reflected one. For this we should examine the band structure of the crystal at the same slice of $$K_y.$$

As we can see, generally there may be multiple transmitted waves. In our case there are only two intersections: one corresponds to the wave propagating to the right, and another to the left.

## Direction of propagation of a crystal mode

Although choosing a particular crystal mode will let us find the phase velocity (though it may appear not so easy to identify for very high bands), this will tell nothing about the direction of a beam of light refracted by the crystal. This spatially localized behavior is described by group velocity instead.

Group velocity in a crystal is defined as usual:

$$\vec v_{\text{g}}=\nabla_{\vec K}\omega(\vec K).\tag9$$

Group velocity is also what we should use to choose from a pair of opposite $$K_x$$ to identify the transmitted wave: the transmitted mode is the one whose group velocity is directed into the crystal (i.e. in our example its $$x$$-component should be positive).

## Snell's law and refractive index

Snell's law follows form conservation of $$K_y$$ and $$K_z$$ on transition through the interface. Considering the incident wave's quasi-wavevector in the incidence plane, we can write it as

$$\vec K_{\text{inc}}=K_{\text{inc}}\begin{pmatrix}\cos\varphi_{\text{inc}}\\ \sin\varphi_{\text{inc}}\end{pmatrix}.\tag{10}$$

Similarly for the transmitted wave we have

$$\vec K_{\text{trans}}=K_{\text{trans}}\begin{pmatrix}\cos\varphi_{\text{trans}}\\ \sin\varphi_{\text{trans}}\end{pmatrix}.\tag{11}$$

Since the $$y$$-components of $$(10)$$ and $$(11)$$ are equal, we get:

$$K_{\text{inc}}\sin\varphi_{\text{inc}}=K_{\text{trans}}\sin\varphi_{\text{trans}}.\tag{12}$$

If the material is homogeneous, then its Brillouin zone is infinite, and $$\vec k=\vec K,$$ and we recover the usual Snell's law by dividing both sides of $$(12)$$ by $$k_{\text{inc}}.$$ Otherwise, we have the generalized version of the Snell's law, where the refractive index is defined in terms of quasi-wavenumbers as

$$n=\frac{K_{\text{trans}}}{K_{\text{inc}}}.\tag{13}$$

Don't forget though, that crystals are generally anisotropic, so refractive index doesn't always let one find the angle of refraction of a beam of light. It's only near the center of the Brillouin zone in the lowest band that we can definitely say that $$\omega(\vec K)$$ becomes isotropic, so we can use Snell's law for this.

In some cases refractive index may become negative (our example is one such case), which will lead to an unusual refraction behavior.