MIT Licence

© Alexey A. Shcherbakov, 2024

Lecture 2.2. Photonic crystals. Complete modal basis in 1D case¶

In this lecture we consider an analytical solution for a one-dimensional photonic crystal whose cell consists of two layers of different materials with parameters $\varepsilon_{1,2}$ and $\mu_{1,2}$. This model is a direct analog of the Kronig-Penney model in the electronic theory of solids, but in contrast to the quantum mechanical case, where the periodic system of rectangular barriers is only an approximation of the potential that allows one to obtain an analytical solution, in optics the piecewise constant variation of material parameters of the medium accurately describes existing structures (if we neglect their external boundaries) - Bragg mirrors.

Helmholtz equation¶

Let the boundaries between different materials in a one-dimensional photonic crystal with period $\Lambda$ be parallel to the $XY$ planes of the Cartesian coordinate system, so that the periodicity is given by the direction of the $Z$-axis. Without restriction of generality, the fields can be considered as independent of the coordinate $y$, and Maxwell's equations split for TE and TM polarizations, respectively: \begin{equation}\tag{1} \begin{array}{c} -\dfrac{\partial E_{y}}{\partial z} = i\omega\mu H_{x} \\ \dfrac{\partial E_{y}}{\partial x} = i\omega\mu H_{z} \\ \dfrac{\partial H_{x}}{\partial z}-\dfrac{\partial H_{z}}{\partial x} = -i\omega\varepsilon E_{y} \end{array} \end{equation} and \begin{equation}\tag{2} \begin{array}{c} \dfrac{\partial H_{y}}{\partial z} = i\omega\varepsilon E_{x} \\ \dfrac{\partial H_{y}}{\partial x} = -i\omega\varepsilon E_{z} \\ \dfrac{\partial E_{x}}{\partial z}-\dfrac{\partial E_{z}}{\partial x} = i\omega\mu H_{y} \end{array} \end{equation} The consequent Helmholtz equations can be written in a unified polarization-independent form: \begin{equation}\tag{3} \dfrac{\partial}{\partial z}\left(\dfrac{1}{\eta\left(z\right)}\dfrac{\partial F}{\partial z}\right)+\dfrac{1}{\eta\left(z\right)}\dfrac{\partial^{2}F}{\partial x^{2}}+\omega^{2}\zeta\left(z\right)F=0 \end{equation} where $F^{s}\left(x,z\right) = E_{y}\left(x,z\right)$, $\eta^s=\mu$ for the s-polarization, and $F^p\left(x,y\right) = H_{y}\left(x,y\right)$, $\eta^p=\varepsilon$ for the p-polarization. Additionally we denote the second field components in the layer plane as $G$, so that $G^s = H_x$ и $G^p = E_x$. This component can be derived from the equations (1),(2): $G^{x,p} = (\pm)_{s,p}(i/\omega\eta^{s,p}) dF^{s,p}/dz$. The Helmholtz equation is supplemented with the continuity conditions for the fields $F$ and $G$ at the layer boundaries together with the quasiperiodicity conditions.

Dispersion equation¶

To solve the equation, let us invoke the method of separation of variables. Writing $F=u\left(x\right)v\left(z\right)$, we arrive at two equations for the new functions: \begin{split}\tag{4} \dfrac{1}{u\left(x\right)}\dfrac{d^{2}u\left(x\right)}{dx^{2}} &= -\beta^{2} \\ \dfrac{\eta\left(z\right)}{v\left(z\right)}\dfrac{d}{dz}\left(\dfrac{1}{\eta\left(z\right)}\dfrac{\partial v\left(z\right)}{\partial z}\right) & + \left[ \omega^{2}\varepsilon\left(z\right)\mu\left(z\right)-\beta^{2} \right] = 0 \end{split} Solution to the first equation corresponds to the field translational invariance along the $X$ axis and explicitly is given by the exponent \begin{equation}\tag{5} u\left(x\right)\sim\exp\left(i\beta x\right) \end{equation} with some propagation constant $\beta$. The second operator in the case when the material parameters are real functions of the coordinates is self-adjoint. Hence, its eigenvalues $\beta^2$ are real, and eigenvectors corresponding to these values are mutually orthogonal. The orthogonality relation is written as \begin{equation}\tag{6} \dfrac{1}{\Lambda}\intop_{0}^{\Lambda}\dfrac{v_{m}\left(z\right) v^*_{n}\left(z\right)}{\omega\eta\left(z\right)}dz=\delta_{mn} \end{equation}

In our considered case of two homogeneous layers with thicknesses $d_{1,2}$ on the period ($d_1+d_2=\Lambda$), the equation on the function $v(z)$ in each layer reduces to a second-order equation with constant coefficients \begin{equation}\tag{7} \dfrac{d^{2}v\left(z\right)}{dz^{2}}+\kappa_{1,2}^{2}v\left(z\right)=0 \end{equation} where $\kappa_{1,2}=\sqrt{\omega^{2}\varepsilon_{1,2}\mu_{1,2}-\beta^{2}}$. We express solutions to this equations as \begin{equation}\tag{8} v\left(z\right)=\begin{cases} \mathcal{A}_{1}^{+}e^{i\kappa_{1}z}+\mathcal{A}_{1}^{-}e^{-i\kappa_{1}z} & 0\leq z\leq d_{1}\\ \mathcal{A}_{2}^{+}e^{i\kappa_{2}\left(z-d_{1}\right)}+\mathcal{A}_{2}^{-}e^{-i\kappa_{2}\left(z-d_{1}\right)} & d_{1}\leq z\leq\Lambda \end{cases} \end{equation} Constants coefficients $\mathcal{A}_{1,2}^{\pm}$ are inter-related by the continuity and quasi-periodicity conditions: \begin{split}\tag{9} v\left(d_{1}-0\right) = v\left(d_{1}+0\right) \\ e^{ik_{B}\Lambda}v\left(0\right) = v\left(\Lambda\right) \end{split} The same conditions written for the second field component $G$, yield analogous equations for the function $\left(1/\eta\right)dv/dz$. The resulting system of four equations is \begin{split}\tag{10} \mathcal{A}_{1}^{+}e^{i\kappa_{1}d_{1}}+\mathcal{A}_{1}^{-}e^{-i\kappa_{1}d_{1}} = \mathcal{A}_{2}^{+}+\mathcal{A}_{2}^{-} \\ e^{ik_{B}\Lambda}\left(\mathcal{A}_{1}^{+} + \mathcal{A}_{1}^{-}\right) = \mathcal{A}_{2}^{+}e^{i\kappa_{2}d_{2}}+\mathcal{A}_{2}^{-}e^{-i\kappa_{2}d_{2}} \\ \dfrac{\kappa_{1}}{\omega\eta_{1}}\left(\mathcal{A}_{1}^{+}e^{i\kappa_{1}d_{1}}-\mathcal{A}_{1}^{-}e^{-i\kappa_{1}d_{1}}\right) = \dfrac{\kappa_{2}}{\omega\eta_{2}}\left(\mathcal{A}_{2}^{+}- \mathcal{A}_{2}^{-}\right) \\ e^{ik_{B}\Lambda}\dfrac{\kappa_{1}}{\omega\eta_{1}}\left(\mathcal{A}_{1}^{+}-\mathcal{A}_{1}^{-}\right) = \dfrac{\kappa_{2}}{\omega\eta_{2}}\left(\mathcal{A}_{2}^{+}e^{i\kappa_{2}d_{2}}-\mathcal{A}_{2}^{-}e^{-i\kappa_{2}d_{2}}\right) \end{split} Equating the determinant to zero, after simple algebraic transformations we obtain the dispersion equation: \begin{equation}\tag{11} \cos\left(\kappa_{1}d_{1}\right)\cos\left(\kappa_{2}d_{2}\right)-\dfrac{1}{2}\left(\dfrac{\kappa_{1}}{\eta_{1}}\dfrac{\eta_{2}}{\kappa_{2}}+\dfrac{\eta_{1}}{\kappa_{1}}\dfrac{\kappa_{2}}{\eta_{2}}\right)\sin\left(\kappa_{1}d_{1}\right)\sin\left(\kappa_{2}d_{2}\right) = \cos\left(k_{B}\Lambda\right) \end{equation} which for real $\varepsilon_{1,2}$, $\mu_{1,2}$ has an infinite countable discrete set of solutions $\beta_m^2$. For given crystal parameters, there is a finite number $\beta_m^2 > 0$ corresponding to modes propagating along the $x$ axis and an infinite number $\beta_m^2 < 0$ corresponding to evanescent solutions. This can be illustrated by plotting the dependence of the left and right-hand parts of the equation on $\beta^2$:

Modal fields¶

Let us rewrite the equations on the coefficients $\mathcal{A}_{1,2m}^{\pm}$ via a new constant: \begin{equation}\tag{12} \begin{array}{c} \mathcal{A}_{1m}^{+}=\sqrt{\mathcal{C}_{m}}\left(1-e^{ik_{B}\Lambda}e^{-i\kappa_{1m}d_{1}}e^{i\kappa_{2m}d_{2}}\right)\left(\dfrac{\eta_{2}}{\eta_{1}}\dfrac{\kappa_{1m}}{\kappa_{2m}}-1\right)\\ \mathcal{A}_{1m}^{-}=\sqrt{\mathcal{C}_{m}}\left(1-e^{ik_{B}\Lambda}e^{i\kappa_{1m}d_{1}}e^{i\kappa_{2m}d_{2}}\right)\left(\dfrac{\eta_{2}}{\eta_{1}}\dfrac{\kappa_{1m}}{\kappa_{2m}}+1\right)\\ \mathcal{A}_{2m}^{+}=\dfrac{1}{2}\sqrt{\mathcal{C}_{m}}\left[\mathcal{A}_{1m}^{+}e^{i\kappa_{1m}d_{1}}\left(1+\dfrac{\eta_{2}}{\eta_{1}}\dfrac{\kappa_{1m}}{\kappa_{2m}}\right)+\mathcal{A}_{1m}^{-}e^{-i\kappa_{1m}d_{1}}\left(1-\dfrac{\eta_{2}}{\eta_{1}}\dfrac{\kappa_{1m}}{\kappa_{2m}}\right)\right]\\ \mathcal{A}_{2m}^{-}=\dfrac{1}{2}\sqrt{\mathcal{C}_{m}}\left[\mathcal{A}_{1m}^{+}e^{i\kappa_{1m}d_{1}}\left(1-\dfrac{\eta_{2}}{\eta_{1}}\dfrac{\kappa_{1m}}{\kappa_{2m}}\right)+\mathcal{A}_{1m}^{-}e^{-i\kappa_{1m}d_{1}}\left(1+\dfrac{\eta_{2}}{\eta_{1}}\dfrac{\kappa_{1m}}{\kappa_{2m}}\right)\right] \end{array} \end{equation} This constant can be derived from the orthogonality relation: \begin{split} 1 &= \dfrac{1}{\Lambda}\int_{0}^{\Lambda} \dfrac{v_{m}\left(z\right)v_{m}^{*}\left(z\right)}{\omega\eta\left(z\right)}dz= \\ &=\left\{ \dfrac{1}{\omega\eta_{1}}\dfrac{d_{1}}{\Lambda}\left[\left(\left|\mathcal{A}_{1m}^{+}\right|^{2}e^{-\Im\kappa_{1m}d_{1}}+\left|\mathcal{A}_{1m}^{-}\right|^{2}e^{\Im\kappa_{1m}d_{1}}\right) \dfrac{\sinh\left(\Im\kappa_{1m}d_{1}\right)}{\Im\kappa_{1m}d_{1}} + 2\Re\left(\mathcal{A}_{1m}^{+}\mathcal{A}_{1m}^{-*}e^{i\Re\kappa_{1m}d_{1}}\right)\dfrac{\sin\left(\Re\kappa_{1m}d_{1}\right)}{\Re\kappa_{1m}d_{1}}\right]+\right. \\ &\left. + \dfrac{1}{\omega\eta_{2}}\dfrac{d_{2}}{\Lambda}\left[\left(\left|\mathcal{A}_{2m}^{+}\right|^{2}e^{-\Im\kappa_{2m}d_{2}}+\left|\mathcal{A}_{2m}^{-}\right|^{2}e^{\Im\kappa_{2m}d_{2}}\right) \dfrac{\sinh\left(\Im\kappa_{2m}d_{2}\right)}{\Im\kappa_{2m}d_{2}}+2\Re\left(\mathcal{A}_{2m}^{+}\mathcal{A}_{2m}^{-*}e^{i\Re\kappa_{2m}d_{2}}\right)\dfrac{\sin\left(\Re\kappa_{2m}d_{2}\right)}{\Re\kappa_{2m}d_{2}}\right] \right\} ^{-1} \end{split}

Expansion of an arbitrary field over the eigenfunctions explicitly writes: \begin{split}\tag{13} F &= \sum_{m=1}^{\infty}v_{m}\left(z\right)\left(a_{m}^{+}e^{i\beta_{m}x}+a_{m}^{-}e^{-i\beta_{m}x}\right) \\ G=\left(\mp_{s,p}\right)\dfrac{1}{i\omega\eta}\dfrac{\partial F}{\partial x} &= \left(\mp_{s,p}\right) \sum_{m=1}^{\infty} v_{m}\left(z\right) \dfrac{\beta_{m}}{\omega\eta} \left(a_{m}^{+}e^{i\beta_{m}x}-a_{m}^{-}e^{-i\beta_{m}x}\right) \end{split} with expansion coefficients $a_m^{\pm}$.

Effective medium limit¶

Now we return to the derived dispersion equation and consider the wave propagation in the $\Gamma$-point. Invoking the trigonometric relations for double angle the equation can be factorized: \begin{split} \dfrac{\eta_{1}}{\kappa_{1}}\dfrac{\kappa_{2}}{\eta_{2}}\left[\sin\left(\dfrac{\kappa_{1}d_{1}}{2}\right)\cos\left(\dfrac{\kappa_{2}d_{2}}{2}\right)+\dfrac{\kappa_{1}}{\eta_{1}}\dfrac{\eta_{2}}{\kappa_{2}}\sin\left(\dfrac{\kappa_{2}d_{2}}{2}\right)\cos\left(\dfrac{\kappa_{1}d_{1}}{2}\right)\right]\times \\ \times\left[\cos\left(\dfrac{\kappa_{1}d_{1}}{2}\right)\sin\left(\dfrac{\kappa_{2}d_{2}}{2}\right)+\dfrac{\kappa_{1}}{\eta_{1}}\dfrac{\eta_{2}}{\kappa_{2}}\sin\left(\dfrac{\kappa_{1}d_{1}}{2}\right)\cos\left(\dfrac{\kappa_{2}d_{2}}{2}\right)\right]=0 \end{split} Then, \begin{equation*} \left[\begin{array}{c} \tan\left(\dfrac{\kappa_{1}d_{1}}{2}\right)=-\dfrac{\kappa_{1}}{\eta_{1}}\dfrac{\eta_{2}}{\kappa_{2}}\tan\left(\dfrac{\kappa_{2}d_{2}}{2}\right)\\ \tan\left(\dfrac{\kappa_{2}d_{2}}{2}\right)=-\dfrac{\kappa_{1}}{\eta_{1}}\dfrac{\eta_{2}}{\kappa_{2}}\tan\left(\dfrac{\kappa_{1}d_{1}}{2}\right) \end{array}\right. \end{equation*} One can show that these two types of solutions correspond to symmetric and antisymmetric modes. Taking the limit $k\Lambda\rightarrow 0$ and substituting $\tan\alpha\approx\alpha$, we arrive at: \begin{equation*} \left[\begin{array}{c} \kappa_{1}d_{1}=-\dfrac{\kappa_{1}}{\eta_{1}}\dfrac{\eta_{2}}{\kappa_{2}}\kappa_{2}d_{2}\\ \kappa_{2}d_{2}=-\dfrac{\kappa_{1}}{\eta_{1}}\dfrac{\eta_{2}}{\kappa_{2}}\tan\kappa_{1}d_{1} \end{array}\right.\Rightarrow\left[\begin{array}{c} \dfrac{d_{1}}{d_{2}}=-\dfrac{\eta_{2}}{\eta_{1}}\\ \kappa_{2}^{2}d_{2}\eta_{1}=-\kappa_{1}^{2}d_{1}\eta_{2} \end{array}\right. \end{equation*} The first equation has no solutions. Using the explicit form of $\kappa_{1,2}$ in the second equation one finds that \begin{equation}\tag{14} \left(\omega^{2}\varepsilon_{2}\mu_{2}-\beta^{2}\right)\eta_{1}d_{2}=-\left(\omega^{2}\varepsilon_{1}\mu_{1}-\beta^{2}\right)\eta_{2}d_{1} \end{equation} In the case of the TE polarization and non-magnetic media \begin{equation}\tag{15} \beta^{2}=\omega^{2}\mu_{0}\dfrac{\varepsilon_{1}d_{1}+\varepsilon_{2}d_{2}}{d_{1}+d_{2}} \end{equation} that is the effective dielectric permittivity is a mean permittivity value. In the case of the TM polarization \begin{equation}\tag{16} \beta^{2} = \omega^{2}\mu_{0} \left( \dfrac{\left(1/\varepsilon_{1}\right)d_{1} + \left(1/\varepsilon_{2}\right) d_{2}}{d_{2}+d_{1}} \right)^{-1} \end{equation} which means averaging of the inverse permittivity.

Taking $\beta=0$, which describes the wave propagation in the direction orthogonal to the crystal layers, we expand the dispersion equaion over the small parameter $k_B\Lambda$, which yields a polarization-independent solution \begin{equation}\tag{17} k_{B}^{2}=\omega^{2}\mu_{0}\dfrac{\varepsilon_{1}d_{1}+\varepsilon_{2}d_{2}}{d_{1}+d_{2}} \end{equation} To summarize, in the large wavelength (small frequency) limit there remains only a sinle propagating mode (for each polarization) and the photonic crystal behaves as an effective uniaxial anisotropic medum, which permittivity tensor in proncipal axes reads \begin{equation}\tag{18} \hat{\varepsilon}=\mathrm{diag}\left\{ \dfrac{\varepsilon_{1}d_{1}+\varepsilon_{2}d_{2}}{d_{1}+d_{2}};\thinspace\dfrac{\varepsilon_{1}d_{1}+\varepsilon_{2}d_{2}}{d_{1}+d_{2}};\thinspace\left(\dfrac{\left(1/\varepsilon_{1}\right)d_{1}+\left(1/\varepsilon_{2}\right)d_{2}}{d_{2}+d_{1}}\right)^{-1}\right\} \end{equation}

References¶

1. Botten, I. C., Craig, M. S., McPhedran, R. C., Adams, J. L., & Andrewartha, J. R. The Dielectric Lamellar Diffraction Grating. Optica Acta: International Journal of Optics, 28(3), 413–428 (1981)
2. Tishchenko, A.V. Phenomenological representation of deep and high contrast lamellar gratings by means of the modal method. Opt Quant Electron 37, 309–330 (2005).