Solving Maxwell equations with Fourier transforms

Introduction

The Maxwell equations are the famous set of partial differential equations (PDEs) that govern the behavior of electric and magnetic fields:

$$\begin{array}{rl}\nabla \cdot \mathbf{E}& =\rho /{\epsilon }_0\\ \nabla \cdot \mathbf{B}& =0\\ \nabla \times \mathbf{E}& =-\frac{\partial \mathbf{B}}{\partial t}\\ \nabla \times \mathbf{B}& ={\mu }_0\mathbf{J}+\frac{1}{c^2}\frac{\partial \mathbf{E}}{\partial t}\end{array}$$

When we consider the Maxwell equations in free space, they reduce to the electromagnetic wave equations:

$$\begin{array}{rl}\left(\frac{{\partial }^2}{\partial t^2}-c^2{\nabla }^2\right)\mathbf{E}& =0\\ \left(\frac{{\partial }^2}{\partial t^2}-c^2{\nabla }^2\right)\mathbf{B}& =0\end{array}$$

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,({\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 )$¹:

$$\begin{array}{rl}i\mathbf{k}\cdot \tilde{\mathbf{E}}& =\rho /{\epsilon }_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{array}$$

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{array}{rl}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{array}$$

Taking the cross product of $\mathbf{k}\times \tilde{\mathbf{E}}$ and $\mathbf{k}\times \tilde{\mathbf{B}}$ with $\mathbf{k}$ then gives us²:

$$\begin{array}{rl}\mathbf{k}\times (\mathbf{k}\times \tilde{\mathbf{E}})& =\underset{\mathbf{k}\cdot \tilde{\mathbf{E}}=0}{\underbrace{(\mathbf{k}\cdot \tilde{\mathbf{E}})}}\mathbf{k}-(\mathbf{k}\cdot \mathbf{k})\tilde{\mathbf{E}}\\ & =-k^2\tilde{\mathbf{E}},k^2=\mathbf{k}\cdot \mathbf{k}\\ \mathbf{k}\times (\mathbf{k}\times \tilde{\mathbf{B}})& =\underset{\mathbf{k}\cdot \tilde{\mathbf{B}}=0}{\underbrace{(\mathbf{k}\cdot \tilde{\mathbf{B}})}}\mathbf{k}-(\mathbf{k}\cdot \mathbf{k})\tilde{\mathbf{B}}\\ & =-k^2\tilde{\mathbf{B}}\end{array}$$

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{array}{rl}\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}}\\ \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{array}$$

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³:

$$\begin{array}{r}{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{array}$$

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.

Footnotes

1. From “Maxwell’s equations in Fourier space and coupling between elementary excitations in solids”, Andrade-Neto, Brandão (2023), eqs. 33-36. We simply swapped with $\mathbf{B}={\mu }_0\mathbf{H}$ and $\mathbf{D}={\epsilon }_0\mathbf{E}$ (note that $1/c^2={\mu }_0{\epsilon }_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 $\frac{{\partial }^2}{\partial t^2}\to {\omega }^2$ ↩