MIT Licence
© Alexey A. Shcherbakov, 2024
In this lecture we will consider an efficient method for calculating eigen solutions in photonic crystals based on the Fourier decomposition of periodic functions. First, the idea will be illustrated on the simplest example of a one-dimensional crystal, and then general formulas for the two-dimensional and three-dimensional cases will be given.
Let us start with a simple example and consider eigenwaves propagating perpendicular to the layers in a one-dimensional photonic crystal with period $\Lambda$. Let the layers of the crystal be parallel to the $XY$ plane of the Cartesian coordinate system. Then it is sufficient to consider the scalar Helmholtz equation for one of the field components, for example, \begin{equation}\tag{1} \left[ \dfrac{d^{2}}{dz^{2}} + \omega^{2}\varepsilon\left(z\right)\mu_0\right] E_x\left(z\right) = 0 \end{equation} For simplicity all the media are supposed to be non-magnetic.
First we write Bloch's theorem and expand the periodic part of the field amplitude into a Fourier series: \begin{equation}\tag{2} E_x\left(z\right) = e^{ik_{B}z} \varphi\left(z\right) = \sum_{n=-\infty}^{\infty}\varphi_{n}e^{i\left(k_{B}+\frac{2\pi n}{\Lambda}\right)z} \end{equation} Direct and inverse Fourier expansion for dielectric permittivity is \begin{split} & \varepsilon\left(z\right)=\sum_{n=-\infty}^{\infty}\epsilon_{n}e^{2\pi ni\frac{z}{\Lambda}} \\ & \epsilon_{n} = \dfrac{1}{\Lambda}\intop_{-\Lambda/2}^{\Lambda/2}\varepsilon\left(z\right)e^{-2\pi ni\frac{z}{\Lambda}}dz \end{split} As a rule, these Fourier coefficients can be easily found analytically. For example, in the case of two layers on the period \begin{equation*} \varepsilon\left(z\right)=\begin{cases} \varepsilon_{1} & z\in\left[-d_{1}/2;d_{1}/2\right)\\ \varepsilon_{2} & z\in\left[-\Lambda/2;-d_{1}/2\right)\cup\left[d_{1}/2;\Lambda/2\right) \end{cases}\Rightarrow\epsilon_{n}=\begin{cases} \varepsilon_{1}\dfrac{d_{1}}{\Lambda}+\varepsilon_{2}\dfrac{d_{2}}{\Lambda} & n=0\\ \left(\varepsilon_{1}-\varepsilon_{2}\right)\dfrac{\sin\left(\pi nd_{1}/\Lambda\right)}{\pi n} & n\neq0 \end{cases} \end{equation*} Now we substitute the field and dielectric constant decomposition into the Helmholtz equation: \begin{equation}\tag{3} \begin{split} -\sum_{n=-\infty}^{\infty}\left(k_{B}+\dfrac{2\pi n}{\Lambda}\right)^{2} & \varphi_{n}\exp\left[i\left(k_{B}+\dfrac{2\pi n}{\Lambda}\right)z\right] + \\ + \omega^{2}\left\{ \sum_{m=-\infty}^{\infty}\epsilon_{m}\exp\left(2\pi im\dfrac{z}{\Lambda}\right)\right\} & \left\{ \sum_{p=-\infty}^{\infty}\varphi_{p}\exp\left[i\left(k_{B}+\dfrac{2\pi p}{\Lambda}\right)z\right]\right\} =0 \end{split} \end{equation} We take advantage of the orthogonality of the exponential factors by multiplying both parts of the equation by $\exp\left[-i\left(k_{B}+2\pi q/\Lambda\right)z\right]$ and integrating over the period: \begin{equation}\tag{4} \left(k_{B}+\dfrac{2\pi q}{\Lambda}\right)^{2}\varphi_{q}-\omega^{2}\sum_{m=-\infty}^{\infty}\epsilon_{q-m}\varphi_{m} = 0 \end{equation} By restricting the summation to the maximum index $\max\left|n\right| = N$ for the computations, we arrive at the generalized eigenvalue equation. In normalized form, it can be written as \begin{equation}\tag{5} \left(\dfrac{k_{B}\Lambda}{2\pi}+q\right)^{2}\varphi_{q}=\left(\dfrac{k_{0}\Lambda}{2\pi}\right)^{2}\sum_{m=-N}^{N}\epsilon_{q-m}\varphi_{m} \end{equation} The vector $\{ \varphi_n \}$ contains $2N+1$ elements, so that one needs to take into account $4N+1$ terms in the permittivity expansion.
The eigenvalue problem for non-magnetic crystals is commonly written for the magnetic field: \begin{equation}\tag{6} \nabla \times \frac{1}{\varepsilon\left(\boldsymbol{r}\right)} \nabla \times\boldsymbol{H}\left(\boldsymbol{r}\right) = \omega^{2}\mu_{0} \boldsymbol{H}\left(\boldsymbol{r}\right) \end{equation} The Bloch theorem in 3D: \begin{equation}\tag{7} \boldsymbol{H}\left(\boldsymbol{r}\right) = \exp\left(i\boldsymbol{k}_{B}\boldsymbol{r}\right) \boldsymbol{H}_{\boldsymbol{k}}\left(\boldsymbol{r}\right) \end{equation} with periodic function $\boldsymbol{H}_{\boldsymbol{k}}(\boldsymbol{r}) = \boldsymbol{H}_{\boldsymbol{k}}(\boldsymbol{r}+\boldsymbol{R}_{m})$ and Bloch wavevector $\boldsymbol{k}_{B}=\left(k_{Bx}, k_{By},k_{Bz}\right)^{T}$. Fourier expansion of the periodic part of the field and the inverse permeability \begin{equation}\tag{8} \begin{split} 1/\varepsilon\left(\boldsymbol{r}\right) &= \sum_{m}f_{m}\exp\left(i\boldsymbol{G}_{m}\boldsymbol{r}\right) \\ \boldsymbol{H}_{\boldsymbol{k}}\left(\boldsymbol{r}\right) &= \sum_{m}\boldsymbol{h}_{m}\exp\left(i\boldsymbol{G}_{m}\boldsymbol{r}\right) \end{split} \end{equation} where $\boldsymbol{G}_{m}$ is the reciprocal lattice vector, and the summation is done over the 3D index. Substitution of the Fourier expansion into the Helmholtz equations gives \begin{equation}\tag{9} \sum_{p} \sum_{n} f_{p} \left[ \left( \boldsymbol{k}_{n} + \boldsymbol{G}_{p} \right) \times \boldsymbol{k}_{n} \times \boldsymbol{h}_{\boldsymbol{k}n} \right] \exp \left[ i \left( \boldsymbol{k}_{n} + \boldsymbol{G}_{p} \right) \boldsymbol{r} \right] = -\omega^{2} \mu_{0} \sum_{m} \boldsymbol{h}_{m} \exp\left(i\boldsymbol{k}_{m}\boldsymbol{r}\right) \end{equation} As in the one-dimensional case, using the orthogonality of the exponential factors yields an equation on the eigenvalues for the Fourier amplitude vectors of the eigenfields: \begin{equation}\tag{10} \sum_{n} f_{m-n} \left(\boldsymbol{k}_{m} \times \boldsymbol{k}_{n}\times\boldsymbol{h}_{\boldsymbol{k}n}\right)=-\omega^{2}\mu_{0}\boldsymbol{h}_{\boldsymbol{k}m} \end{equation} The convolution and vector multiplications in the left-hand part of the equation can be represented as multiplication of some matrix by a vector. Solving the eigenvalue problem for three-dimensional structures becomes demanding on computational resources, so algorithmic ways to reduce the computational complexity $O(N^3)$ of this problem are applied.
If the dielectric constant of one of the crystal materials depends on the wavelength, the above eigenvalue problem for a fixed Bloch vector turns out to be essentially nonlinear. It is convenient to reformulate it as a problem to find the modulus of the Bloch vector $k_B$ for a fixed frequency $\omega$ and the direction of this vector $\hat{\boldsymbol{k}}_B$. Let us rewrite the equation by substituting the wavevectors in the form $\boldsymbol{k}_{m}=k_{B}\hat{\boldsymbol{k}}_{B}+\boldsymbol{G}_{m}$: \begin{equation}\tag{11} \begin{split} & k_{B}^{2} \left[\sum_{n}f_{m-n}\left(\hat{\boldsymbol{k}}_{B} \times \hat{\boldsymbol{k}}_{B} \times \boldsymbol{h}_{n}\right) \right] + \\ & + k_{B} \left[\sum_{n} f_{m-n} \left( \boldsymbol{G}_{m} \times \hat{\boldsymbol{k}}_{B} \times \boldsymbol{h}_{n} + \hat{\boldsymbol{e}}_{B} \times \boldsymbol{G}_{n} \times \boldsymbol{h}_{n} \right) \right] + \\ & + \sum_{n} f_{m-n} \left(\boldsymbol{G}_{m} \times\boldsymbol{G}_{n} \times \boldsymbol{h}_{n} \right) + \omega^{2} \mu_{0} \boldsymbol{h}_{m} = 0 \end{split} \end{equation} This equation can be rewritten in the form of a quadratic eigenvalue problem \begin{equation}\tag{12} k_{B}^{2}M_{2}\boldsymbol{h} + k_{B}M_{1}\boldsymbol{h} + M_{0}\boldsymbol{h}=0 \end{equation} which can be reduced to a generalized eigenvalue problem of twice the size: \begin{equation}\tag{13} \left(\begin{array}{cc} M_{0} & M_{1}\\ 0 & I \end{array}\right)\left(\begin{array}{c} \boldsymbol{h}_1\\ \boldsymbol{h}_2 \end{array}\right)=k_{B}\left(\begin{array}{cc} 0 & -M_{2}\\ I & 0 \end{array}\right)\left(\begin{array}{c} \boldsymbol{h}_1\\ \boldsymbol{h}_2 \end{array}\right) \end{equation} The latter equation can be solved by standard functions present in all math packages.
In the case of two dimensions, the general relations discussed above are simplified, the Fourier decomposition needs to be performed only in the plane of the variable $\boldsymbol{\rho}=(x,y)^T$, the fields are invariant with respect to translations along the $Z$-axis perpendicular to the crystal plane. The periodic field function can be decomposed by two orthogonal polarization vectors \begin{equation}\tag{14} \boldsymbol{H}_{\boldsymbol{k}} \left(\boldsymbol{\rho}\right) = \hat{\boldsymbol{e}}_{\boldsymbol{k}}^{s} H_{\boldsymbol{k}}^{s} \left(\boldsymbol{\rho}\right) + \hat{\boldsymbol{e}}_{\boldsymbol{k}}^{p} H_{\boldsymbol{k}}^{p} \left(\boldsymbol{\rho}\right) \end{equation} where $\hat{\boldsymbol{e}}_{\boldsymbol{k}}^{p} = \hat{\boldsymbol{e}}_{z}$, $\hat{\boldsymbol{e}}_{\boldsymbol{k}}^{s} = (1/k) \boldsymbol{k} \times \hat{\boldsymbol{e}}_{z}$, so that $\boldsymbol{k} \cdot \hat{\boldsymbol{e}}_{z} = 0$. The Helmholtz equation reduces to \begin{equation}\tag{15} \sum_{n}f_{m-n} \left[ \boldsymbol{k}_{m} \times \boldsymbol{k}_{n} \times \left( \hat{\boldsymbol{e}}_{n}^{s} H_{\boldsymbol{k},n}^{s} + \hat{\boldsymbol{e}}_{n}^{p} H_{\boldsymbol{k},n}^{p} \right) \right] = - \omega^{2} \mu_{0} \left( \hat{\boldsymbol{e}}_{m}^{s} H_{\boldsymbol{k}m}^{s} + \hat{\boldsymbol{e}}_{m}^{p} H_{\boldsymbol{k}m}^{p}\right) \end{equation} which, in turn, splits into two independent equations for the two polarizations: \begin{equation}\tag{16} \begin{split} \sum_{n} f_{m-n} \left|\boldsymbol{k}_{m}\right| \left|\boldsymbol{k}_{n}\right| H_{\boldsymbol{k}n}^{s} = \omega^{2} \mu_{0} H_{\boldsymbol{k}m}^{s} \\ \sum_{n} f_{m-n} \left(\boldsymbol{k}_{m} \cdot \boldsymbol{k}_{n}\right) H_{\boldsymbol{k}n}^{p} = \omega^{2} \mu_{0} H_{\boldsymbol{k}m}^{p} \end{split} \end{equation} They can be solved either with respect to frequency or to the modulus of the wave vector.