## Introduction
The Maxwell equations are the famous set of partial differential equations (PDEs) that govern the behavior of electric and magnetic fields:
$
\begin{align*}
\nabla \cdot \mathbf{E} &= \rho/\varepsilon_0 \\
\nabla \cdot \mathbf{B} &= 0 \\
\nabla \times \mathbf{E} &= -\dfrac{\partial \mathbf{B}}{\partial t} \\
\nabla \times \mathbf{B} &= \mu_0 \mathbf{J} + \dfrac{1}{c^2} \dfrac{\partial \mathbf{E}}{\partial t}
\end{align*}
$
When we consider the Maxwell equations in free space, they reduce to the electromagnetic wave equations:
$
\begin{align*}
\left(\dfrac{\partial^2}{\partial t^2} - c^2\nabla^2\right)\mathbf{E} &= 0 \\
\left(\dfrac{\partial^2}{\partial t^2} - c^2\nabla^2\right)\mathbf{B} &= 0
\end{align*}
$
If we assume electromagnetic waves that oscillate at a constant frequency $\omega$, the electromagnetic wave equations further simplify to the Helmholtz equation:
$
(\nabla^2 + k^2)\mathbf{E} = 0, \quad (\nabla^2 + k^2)\mathbf{B} = 0
$
Solving the Maxwell equations is necessary for a multitude of uses, most notably telecommunications and optics. However, analytical solutions to the Maxwell equations are only known for relatively simple cases or as approximate solutions for more complicated systems. Otherwise, we must turn to numerical methods. At Project Elara, we primarily use finite-difference and finite-element methods for these purposes, which are very well-known. However, these are not the only methods to solve the Maxwell equations. Another method is the _Fourier transform method_. This method relies on the fact that if we take the Fourier transform of the Maxwell equations, the Maxwell equations turn into a series of _algebraic equations_ relating the Fourier-transformed fields $\tilde{\mathbf{E}}(\mathbf{k}, \omega), \tilde{\mathbf{B}}(\mathbf{k}, \omega)$ and Fourier-transformed current density $\tilde{\mathbf{J}}(\mathbf{k}, \omega)$[^1]:
$
\begin{align*}
i\mathbf{k} \cdot \tilde{\mathbf{E}} &= \rho/\varepsilon_0 \\
\mathbf{k} \cdot \tilde{\mathbf{B}} &= 0 \\
\mathbf{k} \times \tilde{\mathbf{E}} &= \omega \tilde{\mathbf{B}} \\
i\mathbf{k} \times \tilde{\mathbf{B}} &= \mu_0 \tilde{\mathbf{J}}(\mathbf{k}, \omega)- i(\omega/c^2) \tilde{\mathbf{E}}
\end{align*}
$
In the case of the Maxwell equations in free space, it becomes even simpler, since $\rho = \tilde{\mathbf{J}} = 0$, so the equations simplify to:
$
\begin{align*}
i\mathbf{k} \cdot \tilde{\mathbf{E}} &= 0 \\
\mathbf{k} \cdot \tilde{\mathbf{B}} &= 0 \\
\mathbf{k} \times \tilde{\mathbf{E}} &= \omega \tilde{\mathbf{B}} \\
i\mathbf{k} \times \tilde{\mathbf{B}} &= - i(\omega/c^2) \tilde{\mathbf{E}}
\end{align*}
$
Taking the cross product of $\mathbf{k} \times \tilde{\mathbf{E}}$ and $\mathbf{k} \times \tilde{\mathbf{B}}$ with $\mathbf{k}$ then gives us[^2]:
$
\begin{align*}
% for electric field
\mathbf{k} \times(\mathbf{k \times \tilde{\mathbf{E}}}) &= \underbrace{(\mathbf{k} \cdot \tilde{\mathbf{E}})}_{\mathbf{k}\, \cdot\, \tilde{\mathbf{E}}\, =\, 0}\mathbf{k} - (\mathbf{k} \cdot \mathbf{k}) \tilde{\mathbf{E}} \\
&= -k^2 \tilde{\mathbf{E}}, \quad k^2 = \mathbf{k} \cdot \mathbf{k} \\
% for magnetic field
\mathbf{k} \times(\mathbf{k \times \tilde{\mathbf{B}}}) &= \underbrace{(\mathbf{k} \cdot \tilde{\mathbf{B}})}_{\mathbf{k}\, \cdot\, \tilde{\mathbf{B}}\, =\, 0}\mathbf{k} - (\mathbf{k} \cdot \mathbf{k}) \tilde{\mathbf{B}} \\
&= -k^2 \tilde{\mathbf{B}}
\end{align*}
$
But we also know from the Fourier-transformed version of Faraday's law that $\mathbf{k} \times \tilde{\mathbf{E}} = \omega \tilde{\mathbf{B}}$ and a slight rearrangement of Ampere's law $i\mathbf{k} \times \tilde{\mathbf{B}} = - i(\omega/c^2) \tilde{\mathbf{E}}$ gives us $\mathbf{k} \times \tilde{\mathbf{B}} = - (\omega/c^2) \tilde{\mathbf{E}}$. Now repeating the same procedure and substituting these results, we have:
$
\begin{align*}
% electric field
\mathbf{k} \times(\mathbf{k \times \tilde{\mathbf{E}}}) &= \mathbf{k} \times \omega\tilde{\mathbf{B}} \\
&= \omega(\mathbf{k} \times \tilde{\mathbf{B}}) \\
&= \omega [-(\omega/c^2) \tilde{\mathbf{E}}] \\
&= -(\omega^2 /c^2) \tilde{\mathbf{E}} \\
% magnetic field
\mathbf{k} \times(\mathbf{k \times \tilde{\mathbf{B}}}) &= \mathbf{k} \times - (\omega/c^2) \tilde{\mathbf{E}} \\
&= -(\omega/c^2) (\mathbf{k} \times \tilde{\mathbf{E}}) \\
&= -(\omega/c^2) \omega \tilde{\mathbf{B}} \\
&= -(\omega^2/c^2) \tilde{\mathbf{E}}
\end{align*}
$
Combining our results, we have $-k^2 \tilde{\mathbf{E}} = -(\omega^2/c^2) \tilde{\mathbf{E}}$ and $-k^2 \tilde{\mathbf{B}} = -(\omega^2/c^2) \tilde{\mathbf{B}}$. Some rearrangement gives us the Fourier-transformed versions of the **electromagnetic wave equation**[^3]:
$
\begin{align*}
k^2\tilde{\mathbf{E}}(\omega, \mathbf{k}) = (\omega^2/c^2)\tilde{\mathbf{E}}(\omega, \mathbf{k}) \\
k^2\tilde{\mathbf{B}}(\omega, \mathbf{k}) = (\omega^2/c^2)\tilde{\mathbf{B}}(\omega, \mathbf{k})
\end{align*}
$
## Numerical methods using Fourier transforms
The Fourier transform turns the Maxwell equations from partial differential equations in normal $)\mathbf{x}, t$) space to _algebraic equations_ in $(\omega, \mathbf{k})$ space, and algebraic equations, as we know, are much easier to solve than PDEs! The idea is that using methods we already know for solving systems of linear equations (from linear algebra), we can easily solve the Fourier-transformed Maxwell equations, then take an _inverse Fourier transform_ to get the real-valued fields $\mathbf{E}, \mathbf{B}$. This is especially helpful due to the fact that there exist highly-optimized _fast Fourier transform_ libraries that can do forward/inverse Fourier transforms very quickly, and also highly-optimized libraries for solving systems of linear equations, so in the end there is not that much code we need to write ourselves. To summarize, using Fourier transform methods gives us a powerful tool to be able to solve Maxwell's equations, one that could definitely be a useful supplement to our use of finite-difference and finite-element methods.
[^1]: From [_"Maxwell’s equations in Fourier space and coupling between elementary excitations in solids"_](https://doi.org/10.1590/1806-9126-RBEF-2023-0076), Andrade-Neto, Brandão (2023), eqs. 33-36. We simply swapped with $\mathbf{B} = \mu_0 \mathbf{H}$ and $\mathbf{D} = \varepsilon_0 \mathbf{E}$ (note that $1/c^2 = \mu_0 \varepsilon_0$). This result is derived in the paper, and can be found after a Fourier transform of the fields across both the time domain $(t \to \omega)$ and spatial domain $(\mathbf{x} \to \mathbf{k})$
[^2]: This comes from the identity $\mathbf{a} \times (\mathbf{b} \times \mathbf{c}) = (\mathbf{a} \cdot \mathbf{c})\mathbf{b} - (\mathbf{a} \cdot \mathbf{b})\mathbf{c}$.
[^3]: This form might make more sense if we realize that the Fourier transform essentially converts differentiation into multiplication, so $\nabla^2 \to k^2$ and $\dfrac{\partial^2}{\partial t^2} \to \omega^2$