Implementing the plane wave expansion method in python from scratch. Still a bit rough but we are getting there.
Maxwell's equations in their differential form in space and time are given by: $$ \nabla \times \mathbf{E} = -\mu \frac{\partial \mathbf{H}}{\partial t}, $$ $$ \nabla \times \mathbf{H} = \varepsilon \frac{\partial \mathbf{E}}{\partial t}. $$ These equations describe how electric ($\mathbf{E}$) and magnetic ($\mathbf{H}$) fields interact and evolve over time. While these equations can be solved numerically using methods like the Finite-Difference Time-Domain (FDTD) method, transforming them into the frequency domain can provide significant advantages for certain types of problems.
When we switch to the frequency domain, Maxwell's equations take the form: $$ \nabla \times \mathbf{E} = -j\omega\mu \mathbf{H}, $$ $$ \nabla \times \mathbf{H} = j\omega\varepsilon \mathbf{E}. $$ Here, $\omega$ represents the angular frequency, and $j$ is the imaginary unit. This transformation is particularly useful because it converts the partial differential equations (PDEs) in time into ordinary differential equations (ODEs) in frequency.
In Fourier space, the spatial derivatives are replaced by wave vectors $\vec{k}$, leading to: $$ j\vec{k} \times \mathbf{E} = -\mu \frac{\partial \mathbf{H}}{\partial t}, $$ $$ j\vec{k} \times \mathbf{H} = \varepsilon \frac{\partial \mathbf{E}}{\partial t}. $$ These equations are advantageous because they simplify the analysis of periodic structures, where the fields can be expressed as a superposition of plane waves with different wave vectors.
Let's recall before proceeding that a scalar field $f(x,y,z) = f(\bf{r})$ can be expressed as a Fourier series $$ f(\bf{r}) = a^{pqr}e^{iG_{pqr}^{k} T^l_k r_l} $$ with the Fourier coefficients being $$ 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$.
When expanding the structure in Fourier space, we really mean to take the Discrete Fourier Transform (DFT) of the $\varepsilon(\bf{r})$ scalar field over a grid. In our case, we can consider a 2D grid of $256\times 256$ pixels containing a disc of $\varepsilon_1=12$ and the background being $\varepsilon_2=1$. The DFT is computed with the Fast Fourier Transform algorithm np.fft.fft2. This operation will leave us with the coefficients $a^{pq}$.
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 assume the notation $\varepsilon_r=\varepsilon$. $$ \epsilon^i_{jk}\partial^jE^k = jk_0\varepsilon H_i $$ In this equation, we will expand the structure: $$ \varepsilon_r(r) = a_{pqr}e^{iG_{k}^{pqr} T^l_k r_l}; $$ and the fields in the unit cell: $$ 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}. $$ Note the presence of a Bloch phase factor $e^{-i\beta_kr^k}$ which represents the orientation of our probing plane-wave relative to the crystal. This vector $\beta$ is the $x$-axis variable you see on band diagrams. For conciseness, 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$. The fields are constructed from this set of plane waves. Fourier coefficients $S_{pqr}^i$ and $U_{pqr}^i$ are independent of space and time. In practice, we limit the number of plane waves (which is analytically infinite) to a practical value like $P\times Q = 9\times9$.
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}$. Just here, the summation is written explicit as $k$ goes up to $n$, no more. $$ (-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} $$ The symbol $\mathcal{A}$ is introduced to denote a convolution tensor that will later be implemented by a Quasi-Toeplitz matrix. The code for computing this matrix is given here:
import numpy as np
def convmat(EPS, P, Q=1, R=1):
'''
Convolution matrix for 3D field.
EPS: np.ndarray, shape = (nx, ny, nz)
P, Q, R: Truncated Fourier space size
C: Convolution Matrix size (PQRxPQR)
'''
nx, ny, nz = EPS.shape
nh = P * Q * R
p = np.arange(-P//2+1, P//2+1)
q = np.arange(-Q//2+1, Q//2+1)
r = np.arange(-R//2+1, R//2+1)
# Take the DFT
EPSFT = np.fft.fftshift(np.fft.fftn(EPS)) / nx / ny / nz
p0 = int(np.floor(nx/2))
q0 = int(np.floor(ny/2))
r0 = int(np.floor(nz/2))
C = np.zeros((nh, nh), dtype=np.complex128)
for rrow in range(R):
for qrow in range(Q):
for prow in range(P):
row = rrow * Q * P + qrow * P + prow
for rcol in range(R):
for qcol in range(Q):
for pcol in range(P):
col = rcol * P * Q + qcol * P + pcol
pfft = p[prow] - p[pcol]
qfft = q[qrow] - q[qcol]
rfft = r[rrow] - r[rcol]
C[row, col] = EPSFT[p0+pfft, q0+qfft, r0+rfft]
return C
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 can 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'} $$ These two sets of equations leave us with a rather straightforward numerical solve!
We can write the previous systems in another symbolic form, embracing their matrix nature: $$ \left\{K\times\right\}U = jk_0\mathcal{A}S $$ and $$ \left\{K\times\right\}S = jk_0\mathcal{B}U\tag{1} $$ with $$ \left\{K\times\right\}:= \epsilon_{jk}^ik^j_{pqr} =\begin{bmatrix} 0 & -K_z&K_y\\ K_z & 0&-K_x\\ (-)K_y & K_x&0\\ \end{bmatrix} $$ Note that $K_i$ is just a matrix with the wave vectors components $k^j_{pqr}$ on its diagonal.
If we solve $(1)$ for $U$ $$ U=(-)\frac{1}{k_0}\mathcal{B}^{-1}\left\{K\times\right\}S, $$ and substitute the solution in the first one: $$ \left\{K\times\right\}\mathcal{B}^{-1}\left\{K\times\right\}S = -k^2_0\mathcal{A}S; $$ We obtain a generalized eigenvalue problem which is numerically solvable: $$ Ax=\lambda Bx $$
By solving this problem (using scipy.linalg.eig for instance), one gets access to the eigen frequencies through $k_0=\frac{\omega}{c}$. The problem can be defined for various $\bf{\beta}$s, leading to the construction of a band diagram.
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.
The first set: $$ 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 $$ The second set: $$ 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 $$ Check that this leads to the following decoupled TE? / TM? equations: $$ \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, $$ and, $$ \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}$ yourself.
$$ \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