MIT Licence
© Alexey A. Shcherbakov, 2024
In this lecture, we recall the known definitions, as well as highlight facts that will be used in the course to derive numerical solutions to the wave equation.
For the scalar Helmholtz equation \begin{equation*} (\nabla^{2}+k^2)\psi(\boldsymbol{r}) = f(\boldsymbol{r}) \end{equation*} the equation on the Green's function \begin{equation*} (\nabla^{2} + k^{2})g(\boldsymbol{r} - \boldsymbol{r}') = -\delta(\boldsymbol{r} - \boldsymbol{r}') \end{equation*} Its solution in a homogeneous isotropic medium with Sommerfeld condition at infinity is explicitly written as \begin{equation*} g_0(\boldsymbol{r} - \boldsymbol{r}') = \frac{\exp\left(ik\left|\boldsymbol{r} - \boldsymbol{r}'\right|\right)} {4\pi\left|\boldsymbol{r} - \boldsymbol{r}'\right|} \end{equation*} Then for a given source $f$ the solution is written as an integral \begin{equation*} \psi\left(\boldsymbol{r}\right) = -\int_{V} g_0\left(\boldsymbol{r} - \boldsymbol{r}'\right) f\left(\boldsymbol{r}'\right) d^{3}\boldsymbol{r}' \end{equation*}
For the corresponding vector equation with electric sources \begin{equation*} \nabla\times\nabla\times\boldsymbol{E} - k^{2}\boldsymbol{E} = i\omega\mu\boldsymbol{J}_e \end{equation*} the equation on the dyadic Green's function has the form \begin{equation*} \nabla\times\nabla\times\mathcal{G}(\boldsymbol{r} - \boldsymbol{r}') - k^{2}\mathcal{G}(\boldsymbol{r} - \boldsymbol{r}') = \mathbb{I} \delta(\boldsymbol{r} - \boldsymbol{r}') \end{equation*} Its solution in a homogeneous isotropic medium with Silver-Muller conditions at infinity is expressed through the scalar Green's function \begin{equation*} \mathcal{G}_0(\boldsymbol{r} - \boldsymbol{r}') = \left(\mathbb{I} + \dfrac{1}{k^{2}} \nabla\nabla \right) g(\boldsymbol{r} - \boldsymbol{r}') \end{equation*} Solution of the Helmholtz equation in the whole space \begin{equation*} \boldsymbol{E} = i\omega\mu \int_{V} \mathcal{G}_0(\boldsymbol{r} - \boldsymbol{r}') \boldsymbol{J}_e \left(\boldsymbol{r}'\right) d^{3}\boldsymbol{r}' \end{equation*}
Explicitly, differentiation of the scalar Green's function leads to the expression \begin{equation*} \mathcal{G}_0(\boldsymbol{r} - \boldsymbol{r}') = \dfrac{e^{ikR}}{4\pi R} \left[ \mathbb{I}\left(1+\dfrac{i}{kR}-\dfrac{1}{k^{2}R^{2}}\right)+\hat{\boldsymbol{e}}_r \hat{\boldsymbol{e}}_r^T \left(\dfrac{3}{k^{2}R^{2}} - \dfrac{3i}{kR} - 1\right)\right] \end{equation*} where $R = |\boldsymbol{r}-\boldsymbol{r}'|$. This expression at the point $\boldsymbol{r}$ is proportional to the field of the electric dipole located at the point with coordinates $\boldsymbol{r}'$.
At large distances $r\gg r'$ \begin{equation*} \mathcal{G}_{0}(\boldsymbol{r}, \boldsymbol{r}') \sim \left(\mathbb{I} - \hat{\boldsymbol{e}}_r \hat{\boldsymbol{e}}_r^T\right)\dfrac{e^{ikr}}{4\pi r}e^{-ik\hat{\boldsymbol{r}}\cdot\boldsymbol{r}'} \end{equation*} which is a transverse spherical wave whose amplitude decreases inversely proportional to the distance to the source. By analogy with the electric field, we can write the condition at infinity \begin{equation*} \lim_{R\rightarrow\infty} R \left[ \nabla\times\mathcal{G}_0(\boldsymbol{r},\boldsymbol{r}') - ik_m\hat{\boldsymbol{R}} \times \mathcal{G}_0(\boldsymbol{r},\boldsymbol{r}') \right] = 0 \end{equation*} with $\boldsymbol{R} = \boldsymbol{r}-\boldsymbol{r}'$.
Additionally the Lorentz reciprocity leads to \begin{equation*} \mathcal{G}_{0}\left(\boldsymbol{r},\boldsymbol{r}'\right) = \mathcal{G}_{0}^{T}\left(\boldsymbol{r}',\boldsymbol{r}\right) \end{equation*}
If we divide the potential into homogeneous and inhomogeneous parts \begin{equation*} \varepsilon = \varepsilon_m + \Delta\varepsilon\Rightarrow k^2 = k_m^2 + \Delta (k^2) \end{equation*} substitute in the equation for the Green's function \begin{equation*} \nabla\times\nabla\times\mathcal{G}(\boldsymbol{r} - \boldsymbol{r}') - k_m^{2}\mathcal{G}(\boldsymbol{r} - \boldsymbol{r}') = \mathbb{I} \delta(\boldsymbol{r} - \boldsymbol{r}') + \Delta (k^2)\mathcal{G}(\boldsymbol{r} - \boldsymbol{r}') \end{equation*} and use the volume integral solution with source \begin{equation*} \mathbb{J} = -\frac{i}{\omega\mu_0}\mathbb{I} \delta(\boldsymbol{r} - \boldsymbol{r}') -i\omega\Delta\varepsilon \mathcal{G}(\boldsymbol{r} - \boldsymbol{r}') \end{equation*} we obtain \begin{equation*} \mathcal{G}(\boldsymbol{r},\boldsymbol{r}') = \mathcal{G}_0(\boldsymbol{r},\boldsymbol{r}') + \int_{V} \mathcal{G}_0(\boldsymbol{r},\boldsymbol{r}'') U(\boldsymbol{r}'') \mathcal{G}(\boldsymbol{r}'',\boldsymbol{r}') d^{3}\boldsymbol{r}'' \end{equation*} This is the Dyson equation.
If we consider a finite volume $\Omega\in\mathbb{R}^3$, the integral solutions of the Helmholtz equation can be obtained from the Green's theorem. For some vector function $\boldsymbol{P}$ and tensor function $\mathcal{Q}$, volume $V$ and its bounding surface $S$, this theorem stands that \begin{equation*} \int_{V} \left[ \boldsymbol{P} \cdot \nabla\times\nabla\times\mathcal{Q} - \left( \nabla\times\nabla\times\boldsymbol{P}\right) \cdot \mathcal{Q} \right] dV = -\oint_{S} \left[ \left(\hat{\boldsymbol{n}} \times\nabla\times\boldsymbol{P}\right) \cdot \mathcal{Q} + \left(\hat{\boldsymbol{n}}\times\boldsymbol{P}\right) \cdot \nabla\times\mathcal{Q} \right] dS \end{equation*} In we substitute \begin{equation*} \boldsymbol{P} = \boldsymbol{E},\thinspace \mathcal{Q} = \mathcal{G}_0(\boldsymbol{r} - \boldsymbol{r}') \end{equation*} use the Helmholtz equation for the electric field, the equation on the Green's function, and Faraday's law, we obtain \begin{equation*} \boldsymbol{E}(\boldsymbol{r}') = ikZ \int_{V} \boldsymbol{J} \cdot \mathcal{G}_0 dV - \oint_{S} \left[ ikZ \left(\hat{\boldsymbol{n}} \times \boldsymbol{H}\right) \cdot \mathcal{G}_0 + \left(\hat{\boldsymbol{n}}\times\boldsymbol{E}\right) \cdot \nabla\times\mathcal{G}_0\right] dS \end{equation*} Here, as before, the free space impedance $Z=\sqrt{\mu/\varepsilon}$$ is introduced.
Analogously, for the magnetic field \begin{equation*} \boldsymbol{H}(\boldsymbol{r}') = \int_{V} \boldsymbol{J} \cdot \nabla\times\mathcal{G}_0 dV + \oint_{S} \left[i\dfrac{k_m}{Z}\left(\hat{\boldsymbol{n}}\times\boldsymbol{E}\right) \cdot \mathcal{G}_0 - \left(\hat{\boldsymbol{n}} \times \boldsymbol{H}\right)\cdot\nabla\times\mathcal{G}_0 \right] dS \end{equation*}
It follows from the explicit expression for the free space Green's tensor function that at the point $\boldsymbol{r} = \boldsymbol{r}'$ the function has a pole of the third order, and the volume integral for the field in the source region is calculated in the sense of the principal value. To calculate this integral, we isolate the integral over an infinitesimal volume near the pole and use the charge conservation law: \begin{split} &\nabla\nabla\int_{V} g_0(\boldsymbol{r},\boldsymbol{r}') \boldsymbol{J} (\boldsymbol{r}') d^{3} \boldsymbol{r}' = \\ &= \lim_{V_{\delta}\rightarrow0} \left[ \nabla\nabla\int_{V \backslash V_{\delta}} g_0 (\boldsymbol{r},\boldsymbol{r}') \boldsymbol{J} (\boldsymbol{r}') d^{3} \boldsymbol{r}' + \nabla\nabla\int_{V_{\delta}} g_0(\boldsymbol{r},\boldsymbol{r}') \boldsymbol{J} (\boldsymbol{r}') d^{3} \boldsymbol{r}' \right] = \\ &= \lim_{V_{\delta}\rightarrow0} \left[ \nabla\nabla\int_{V \backslash V_{\delta}} g_0 (\boldsymbol{r},\boldsymbol{r}') \boldsymbol{J} (\boldsymbol{r}') d^{3} \boldsymbol{r}' - \nabla\int_{V_{\delta}} \nabla' g_0(\boldsymbol{r},\boldsymbol{r}') \boldsymbol{J} (\boldsymbol{r}') d^{3} \boldsymbol{r}' \right] = \\ &= \lim_{V_{\delta}\rightarrow0} \left[ \nabla\nabla\int_{V \backslash V_{\delta}} g_0 (\boldsymbol{r},\boldsymbol{r}') \boldsymbol{J} (\boldsymbol{r}') d^{3} \boldsymbol{r}' - \nabla \left( \int_{S_{\delta}} g_0(\boldsymbol{r},\boldsymbol{r}') \boldsymbol{n} \cdot \boldsymbol{J} (\boldsymbol{r}') dS' - i\omega \int_{V_{\delta}} g_{0}(\boldsymbol{r},\boldsymbol{r}') \rho_e(\boldsymbol{r}') d^{3}\boldsymbol{r}' \right) \right] \end{split} Assuming that the charge density $\rho_e(\boldsymbol{r}')$ is a continuous function of coordinates, we obtain that the last summand is zero. The product $\boldsymbol{n} \cdot \boldsymbol{J}$ defines the surface charge on the surface $S_{\delta}$, so the penultimate summand reduces to the product of a constant multiplier depending on the shape of this surface by the current at the point $\boldsymbol{r}$, so that \begin{equation*} \nabla\nabla\int_{V} g_0(\boldsymbol{r},\boldsymbol{r}') \boldsymbol{J} (\boldsymbol{r}') d^{3} \boldsymbol{r}' = \lim_{V_{\delta}\rightarrow0} \nabla\nabla\int_{V \backslash V_{\delta}} g_0 (\boldsymbol{r},\boldsymbol{r}') \boldsymbol{J} (\boldsymbol{r}') d^{3} \boldsymbol{r}' - \mathcal{L} \boldsymbol{J} (\boldsymbol{r}) \end{equation*} The complete field equation is written as \begin{equation*} \boldsymbol{E} = i\omega\mu_{0} P.V.\int_{V} \mathcal{G}_{0} (\boldsymbol{r} - \boldsymbol{r}') \boldsymbol{J} (\boldsymbol{r}') d^{3}\boldsymbol{r}' + \dfrac{1}{i\omega\varepsilon_{m}} \mathcal{L} \cdot \boldsymbol{J} (\boldsymbol{r}) \end{equation*} For example, for a spherical volume \begin{equation*} \mathcal{L}_{sph} = \frac{1}{3} \mathbb{I} \end{equation*}
Consider the equation for the tensor Green's function. The three-dimensional Fourier transform on the variable $\boldsymbol{r}$ gives \begin{equation*} -\boldsymbol{k} \times \boldsymbol{k} \times \mathcal{G}_0 (\boldsymbol{k},\boldsymbol{r}') - k_{0}^{2} \mathcal{G}_0 (\boldsymbol{k},\boldsymbol{r}') = \mathbb{I} e^{-i\boldsymbol{k}\boldsymbol{r}'} \Rightarrow \mathcal{G}_0 (\boldsymbol{k},\boldsymbol{r}') = \dfrac{\mathbb{I} k_{0}^{2} - \boldsymbol{k}\boldsymbol{k}^T}{k_{0}^{2}\left(k^{2}-k_{0}^{2}\right)} e^{-i\boldsymbol{k}\boldsymbol{r}'} \end{equation*} Using the inverse Fourier transform we obtain the spectral decomposition \begin{equation*} \mathcal{G}_0 (\boldsymbol{r},\boldsymbol{r}') = \dfrac{1}{\left(2\pi\right)^{3}} \int_{-\infty}^{\infty} d^{3}\boldsymbol{k} e^{i\boldsymbol{k}(\boldsymbol{r}-\boldsymbol{r}')} \dfrac{\mathbb{I}k_{0}^{2} - \boldsymbol{k}\boldsymbol{k}^T} {k_{0}^{2}(k^{2}-k_{0}^{2})} \end{equation*} This integral diverges in the classical sense - the fraction under the integral tends to a constant at $|\boldsymbol{k}|\rightarrow\infty$. Consider the limit of $|k_z|\rightarrow\infty$ when $|k_{x,y}|<const$ is bounded: \begin{equation*} \dfrac{\mathbb{I}k_{0}^{2} - \boldsymbol{k}\boldsymbol{k}^T} {k_{0}^{2}(k^{2}-k_{0}^{2})} \sim -\hat{\boldsymbol{e}}_z \hat{\boldsymbol{e}}_z^T \end{equation*} We isolate the singular part in the integral, and integrate the regular part by $k_z$ using the Jordan's lemma, calculating the residues at the poles $\pm\sqrt{k_{0}^{2}-k_{x}^{2}-k_{y}^{2}}$: \begin{align*} \mathcal{G}_0(\boldsymbol{r},\boldsymbol{r}') =& \dfrac{1}{\left(2\pi\right)^{3}} \int_{-\infty}^{\infty} d^{3}\boldsymbol{k} e^{i\boldsymbol{k}(\boldsymbol{r}-\boldsymbol{r}')} \left[ \dfrac{\mathbb{I}k_{0}^{2} - \boldsymbol{k}\boldsymbol{k}} {k_{0}^{2}\left(k^{2}-k_{0}^{2}\right)} + \dfrac{\hat{\boldsymbol{e}}_z \hat{\boldsymbol{e}}_z}{k_{0}^{2}} \right] - \dfrac{\hat{\boldsymbol{e}}_z \hat{\boldsymbol{e}}_z}{\left(2\pi\right)^{3}k_{0}^{2}} \int_{-\infty}^{\infty} d^{3}\boldsymbol{k} e^{i\boldsymbol{k}\left(\boldsymbol{r}-\boldsymbol{r}'\right)} = \\ =& \dfrac{i}{8\pi^{2}} \int_{-\infty}^{\infty} dk_xdk_y e^{i\boldsymbol{\varkappa} (\boldsymbol{\rho} - \boldsymbol{\rho}') + ik_{z} \left|z-z'\right|} \dfrac{\mathbb{I}k_{0}^{2} - \boldsymbol{k}_{0}\boldsymbol{k}_{0}}{k_{0}^{2}k_{z}} - \dfrac{\hat{\boldsymbol{e}}_z \hat{\boldsymbol{e}}_z}{k_{0}^{2}}\delta\left(\boldsymbol{r}-\boldsymbol{r}'\right) \end{align*} where $\boldsymbol{\varkappa} = (k_x,k_y)^T$, $\boldsymbol{\rho} = (x,y)^T$. The resulting expression is the expansion of the free space Green's tensor function by plane waves. The integrand contains only transverse plane waves, and the singularity occurs only in the source region. The comoponents of the wave vector satisfy the dispersion equation, and the sign of the third component is chosen based on the radiation conditions at infinity: \begin{equation*} \boldsymbol{k}_{0} = \begin{cases} \hat{\boldsymbol{e}}_x k_{x} + \hat{\boldsymbol{e}}_y k_{y} + \hat{\boldsymbol{e}}_z k_{z} & z-z'>0 \\ \hat{\boldsymbol{e}}_x k_{x} + \hat{\boldsymbol{e}}_y k_{y} - \hat{\boldsymbol{e}}_z k_{z} & z-z'<0 \end{cases},\thinspace\Im m\left[k_{z}\right] \geq 0 \end{equation*}
Previously we obtained the complete basis of vector spherical waves $\boldsymbol{M}$, $\boldsymbol{N}$, $\boldsymbol{L}$, $\boldsymbol{L}$, the first two of which are solenodal and the third is vortex-free. In Cartesian coordinates, an arbitrary field can be decomposed in these functions as \begin{equation*} \boldsymbol{E} = \iiint_{-\infty}^{\infty} d^{3}\boldsymbol{k} \left[ a_M(\boldsymbol{k})\boldsymbol{M}(\boldsymbol{k},\boldsymbol{r}) + a_N(\boldsymbol{k})\boldsymbol{N}(\boldsymbol{k},\boldsymbol{r}) + a_L(\boldsymbol{k})\boldsymbol{L}(\boldsymbol{k},\boldsymbol{r}) \right] \end{equation*} Considering the orthogonality relations, we can write the unit decomposition (it can be checked by direct substitution) \begin{equation*} \mathcal{I}\delta\left(\boldsymbol{r} - \boldsymbol{r}'\right) = \dfrac{1}{\left(2\pi\right)^{3}} \iiint_{-\infty}^{\infty} d^{3}\boldsymbol{k} \left[ \dfrac{\boldsymbol{M}\left(\boldsymbol{k},\boldsymbol{r}\right)\boldsymbol{M}^T\left(-\boldsymbol{k},\boldsymbol{r}'\right) + \boldsymbol{N}\left(\boldsymbol{k},\boldsymbol{r}\right)\boldsymbol{N}^T\left(-\boldsymbol{k},\boldsymbol{r}'\right)}{\varkappa^{2}} + \dfrac{\boldsymbol{L}\left(\boldsymbol{k},\boldsymbol{r}\right)\boldsymbol{L}^T\left(-\boldsymbol{k},\boldsymbol{r}'\right)} {k_0^{2}}\right] \end{equation*} To express the Green's function, let us write for it a decomposition similar to the decomposition of the field by vector waves: \begin{equation*} \mathcal{G}=\iiint_{-\infty}^{\infty} d^{3}\boldsymbol{k} \left[ \boldsymbol{M}\left(\boldsymbol{k},\boldsymbol{r}\right)\boldsymbol{a}_M^T\left(\boldsymbol{k},\boldsymbol{r}'\right) + \boldsymbol{N}\left(\boldsymbol{k},\boldsymbol{r}\right)\boldsymbol{a}_N^T\left(\boldsymbol{k},\boldsymbol{r}'\right) + \boldsymbol{L}\left(\boldsymbol{k},\boldsymbol{r}\right)\boldsymbol{a}_L^T\left(\boldsymbol{k},\boldsymbol{r}'\right) \right] \end{equation*} Substituting this expansion and the unit expansion into the Helmholtz equation for the Green's function, and using the orthogonality relations gives an explicit expression of the coefficients, so that for the Green's dyadic we have \begin{equation*} \mathcal{G} = \dfrac{1}{\left(2\pi\right)^{3}} \iiint_{-\infty}^{\infty} d^{3}\boldsymbol{k} \left[ \dfrac{\boldsymbol{M}\left(\boldsymbol{k},\boldsymbol{r}\right)\boldsymbol{M}^T\left(-\boldsymbol{k},\boldsymbol{r}'\right)+\boldsymbol{N}\left(\boldsymbol{k},\boldsymbol{r}\right)\boldsymbol{N}^T\left(-\boldsymbol{k},\boldsymbol{r}'\right)}{\left(k^{2}-k_{0}^{2}\right)\varkappa^{2}} - \dfrac{\boldsymbol{L}\left(\boldsymbol{k},\boldsymbol{r}\right)\boldsymbol{L}^T\left(-\boldsymbol{k},\boldsymbol{r}'\right)}{k_{0}^{2}k^{2}}\right] \end{equation*} By analogy with the spectral decomposition, singular part extraction and integration with application of Jordan's lemma leads to a plane wave Green's tensor decomposition, equivalent to the one obtained above when considering the spectral decomposition: \begin{equation*} \mathcal{G} = \dfrac{i}{8\pi^{2}}\iint_{-\infty}^{\infty} dk_xdk_y \left[ \dfrac{\boldsymbol{M}\left(\boldsymbol{\varkappa},\pm k_{z},\boldsymbol{r}\right)\boldsymbol{M}^T\left(-\boldsymbol{\varkappa},\mp k_{z},\boldsymbol{r}\right) + \boldsymbol{N}\left(\boldsymbol{\varkappa},\pm k_{z},\boldsymbol{r}\right)\boldsymbol{N}^T\left(-\boldsymbol{\varkappa},\mp k_{z},\boldsymbol{r}\right)}{k_{z}\varkappa^{2}}\right] -\dfrac{\hat{\boldsymbol{e}}_z\hat{\boldsymbol{e}}_z^T}{k_{0}^{2}}\delta\left(\boldsymbol{r}-\boldsymbol{r}'\right) \end{equation*} Here it is also assumed that the wave vector components satisfy the dispersion equation and a condition on the choice of sign before the propagation constant along the $Z$-axis is applied.
Similarly to the case of plane waves, we obtain an expression for the tensor Green's function decomposed by vector spherical waves: \begin{equation*} \mathcal{G}\left(\boldsymbol{r},\boldsymbol{r}'\right) = ik_{0}\sum_{nm}\dfrac{1}{n\left(n+1\right)}\left[\boldsymbol{M}^{\sigma}_{nm}\left(k_{m}\boldsymbol{r}\right)(\hat{\boldsymbol{M}}^{\sigma'}_{n,-m})^T\left(k_{m}\boldsymbol{r}'\right)+\boldsymbol{N}^{\sigma}_{nm}\left(k_{m}\boldsymbol{r}\right)(\hat{\boldsymbol{N}}^{\sigma'}_{n,-m})^T\left(k_{m}\boldsymbol{r}'\right)\right]-\dfrac{\hat{\boldsymbol{e}}_r\hat{\boldsymbol{e}}_r}{k_{0}^{2}}\delta\left(\boldsymbol{r}-\boldsymbol{r}'\right) \end{equation*} Here for two types of solutions $\sigma=1$, $\sigma'=3$ for $r<r'$, and $\sigma=3$, $\sigma'=1$ for $r>r'$.