WIP Plane wave expansion method
We start as always with Maxwell's equations $$ \nabla \times E = -\mu \frac{\partial H}{\partial t} $$ $$ \nabla \times H = \varepsilon \frac{\partial E}{\partial t} $$
Maxwell's equations are expressed in spatio-temporal Fourier space
$$ \nabla \times E=-j\omega\mu H $$ $$ \nabla \times H=j\omega\varepsilon E $$ $$ j\vec k\times E=-\mu\frac{\partial H}{\partial t} $$ $$ j\vec k\times H=\varepsilon\frac{\partial E}{\partial t} $$
A scalar field $f(x,y,z) = f(r)$ can be expressed as a Fourier series $$ f(r) = a^{pqr}e^{iG_{pqr}^{k} T^l_k r_l} $$ with $$ a^{pqr}=\frac{1}{V}\int\int\int_V f(r) e^{-iG_{pqr}^{k}T_{k}^l r_l} dV $$ where we have $T^l_k$, a matrix formed with the lattice vectors as columns, $r_l$ the position vector and $G_{pqr}$ the expansion coefficients $p, q, r$ in a column vector. For a square lattice we have $T_{lk}=\frac{2\pi}{\Lambda_l}\delta_{lk}$ with $\Lambda_l$ the pitch of the lattice in direction $l$.
We normalize $H$ to $\bar H = -j \sqrt{\frac{\mu_0}{\varepsilon_0}}H$ and let the $\bar H$ notation fall for now (we should undo this operation when computing quantities in MKSA units). We also take the notation $\varepsilon_r=\varepsilon$. $$ \epsilon^i_{jk}\partial^jE^k = jk_0\varepsilon H_i $$ In this equation, we will expand $$ \varepsilon_r(r) = a_{pqr}e^{iG_{k}^{pqr} T^l_k r_l} $$ and $$ E_i(r) = e^{-i\beta_kr^k}S^i_{pqr}e^{iG^{pqr}_{k} T^k_l r^l} = S^i_{pqr}e^{-i(\beta_l-G^{pqr}_{k} T_l^k)r^l} $$ Finally, the argument of the exponential can be written as the wave vector $k^{pqr}=\beta_l-G^{pqr}_{k} T_l^k$. $$ E_i(r) = S^i_{pqr}e^{-ik^{pqr}_lr^l} $$ $$ H_i(r) = U^i_{pqr}e^{-ik^{pqr}_lr^l} $$ New we clearly see the meaning of plane wave expansion with this set of plane waves numbered by $p$, $q$ and $r$.
We finally substitute the expansions into frequency-domain Maxwell's equations. We take for example: $$ \epsilon^i_{jk}\partial^jH^k=k_0\varepsilon_rE_i $$
$$ \epsilon^i_{jk}\partial^jU^k_{pqr}e^{-ik^{pqr}_lr^l}=k_0a_{p'q'r'}e^{iG^{p'q'r'}_{k} T^k_l r^l}S^i_{pqr}e^{-ik^{pqr}_lr^l} $$ This leaves an horrendous equation, hopefully, applying the derivatives is not too hard.
$$ (-i)\epsilon^i_{jk}k_{pqr}^j U^k_{pqr}e^{-ik^{pqr}_lr^l}=k_0a_{p'q'r'}e^{iG^{p'q'r'}_k T^l_k r_l}S^i_{pqr}e^{-ik_{pqr}^lr_l} $$
On the right-hand side, we have sadly a product of two sums. This can be reformulated using Cauchy product $a_n\delta^{nn}b_m\delta^{mm}$ $=\delta_{nn}\sum^n_k a_kb^{n-k}$. Here the summation is explicit as $k$ goes up to $n$. The symbol $\mathcal{A}$ is introduced to denote a convolution tensor that will later be implemented by a Toeplitz matrix. $$ (-i)\epsilon^i_{jk}k_{pqr}^j U^k_{pqr}e^{-ik_{pqr}^lr_l}=k_0a^{(p-p')(q-q')(r-r')}S^i_{p'q'r'}e^{-ik_{pqr}^lr_l} $$ $$ (-i)\epsilon^i_{jk}k_{pqr}^j U^k_{pqr}e^{-ik_{pqr}^lr_l}=k_0\mathcal{A}^{pqr}_{p'q'r'}S_i^{p'q'r'}e^{-ik_{pqr}^lr_l} $$
We now consider that a choice of indices $p$, $q$ and $r$ lead to independent (non degenerate) terms in the sum such that we can consider each term separately in a matrix. This effectively make the problem a system of $PQR$ equations. We also now simplify the exponential term:
$$ \epsilon^i_{jk}k_{pqr}^j U^k_{pqr}=ik_0\mathcal{A}_{pqr}^{p'q'r'}S^i_{p'q'r'} $$ The operations are similar for the other Maxwell equation, leading to a similar equation set. $\mathcal{B}$ being the convolution tensor for the magnetic permeability. $$ \epsilon^i_{jk}k_{pqr}^j S^k_{pqr}=ik_0\mathcal{B}_{pqr}^{p'q'r'}U^i_{p'q'r'} $$
We can write this in another, symbolic form: $$ \left\{K\times\right\}U = jk_0\mathcal{A}S $$
and $$ \left\{K\times\right\}S = jk_0\mathcal{B}U $$ with $$ \left\{K\times\right\}=\begin{bmatrix} 0 & -K_z&K_y\\ K_z & 0&-K_x\\ (-)K_y & K_x&0\\ \end{bmatrix} $$
If we solve for $U$ and substitute in the first one: $$ (-)\frac{1}{k_0}\mathcal{B}^{-1}\left\{K\times\right\}S =U $$ $$ \left\{K\times\right\}\mathcal{B}^{-1}\left\{K\times\right\}S = -k^2_0\mathcal{A}S $$ This is a generalized eigenvalue problem which is numerically solvable. $$ Ax=\lambda Bx $$
By considering only the 2D case ($K_z=0$) in $\left\{K\times\right\}$, we get a simpler system of equations and TE and TM modes are split. Guess what are the TE and TM equations.
$$ K_xU_y-K_yU_x=jk_0\mathcal{A}S_z $$ $$ K_yS_z=jk_0\mathcal{B} U_x $$ $$ (-)K_xS_z=jk_0\mathcal{B} U_y $$
$$ K_xS_y-K_yS_x=jk_0\mathcal{B} U_z $$ $$ K_yU_z=jk_0\mathcal{A}S_x $$ $$ (-)K_xU_z=jk_0\mathcal{A}S_y $$ Now you have left to express the 2D case as a generalized eigenvalue problem. $$ Ax = \lambda Bx $$ This leads to $$ \left(K_x\mathcal{B}^{-1}K_x+K_y\mathcal{B}^{-1}K_y\right)S_z = k_0^2\mathcal{A}S_z $$
$$ \left(K_x\mathcal{A}^{-1}K_x+K_y\mathcal{A}^{-1}K_y\right)U_z = k_0^2\mathcal{B}U_z $$
How to implement $\mathcal{A}$ and $\mathcal{B}$.
$$ \begin{bmatrix} \mathcal{F}\{\varepsilon\}_{(-2(P-1))(-2(Q-1))(-2(R-1))})&&&\dots\\ &\mathcal{F}\{\varepsilon\}_{(p-p')(q-q')(r-r')})&&\\ &&\mathcal{F}\{\varepsilon\}_{000}&\\ \dots &&&\ddots\\ \end{bmatrix} $$
home