Fundamental solution for Helmholtz equation in higher dimensions

Question

The fundamental solution for Helmholtz equation $(\Delta + k^2) u = -\delta$ is $e^{i k r}/r$ in 3d and $H_0^1(kr)$ in 2d (up to normalization constants). Is there an explicit expression (eventually in terms of special functions) for the fundamental solution in dimension $\geq 4$? How can one derive it?

Answer 1

We first assume that the point source is located at $\vec r'=0$, $\vec r'\in\mathbb{R}^n$. Spherical symmetry implies that the Green (or "Green's") Function, $G(\vec r|\vec r'=0)$, for the Helmholtz Equation can be written

$$\frac{\partial^2G}{\partial r^2}+\frac{n-1}{r}\frac{\partial G}{\partial r}+k^2G=0\tag 1$$

for $\vec r\ne 0$.

---

FINDING A GENERAL SOLUTION TO $(1)$:

Enforcing the substitution $G(\vec r|\vec r'=0)=r^{1-n/2}g(r)$ in $(1)$ reveals

$$r^2g''(r)+rg'(r)+\left((kr)^2-\left(\frac n2-1\right)^2\right)g(r)=0\tag 2$$

which is Bessel's Differential Equation (with argument scaled by $k$). Solutions to $(2)$ are Bessel Functions of order $n/2-1$ and argument $kr$.

Therefore the general solution to $(1)$ can be written in terms of the first and second kind Hankel functions as

$$\bbox[5px,border:2px solid #C0A000]{G(\vec r|\vec r'=0)=Ar^{1-n/2}H^{(2)}_{n/2-1}(kr)+Br^{1-n/2}H^{(1)}_{n/2-1}(kr)}\tag 3$$

In order for the Green Function to represent an outward travelling wave, either $A=0$ or $B=0$. If the time convention is $e^{i\omega t}$ and $\text{Im}(k)<0$ in a medium with losses, then $B=0$. If the time convention is $e^{-i\omega t}$ and $\text{Im}(k)>0$ in a medium with losses, then $A=0$.

---

HOW TO DETERMINE $B$:

To find $A$ or $B$, we apply the condition

$$\lim_{\epsilon\to 0}\oint_{r=\epsilon}\frac{\partial G}{\partial r}\,dS_{n-1}=-1\tag 4$$

where the integration is over the surface of the $n$-sphere of radius $\epsilon$. The surface area of the $n$-sphere of radius $\epsilon$ is $S_{n-1}=\frac{2\pi^{n/2}\epsilon^{n-1}}{\Gamma(n/2)}$.

Using $(3)$ with $A=0$ in $(4)$ yields

$$\lim_{\epsilon\to 0}\left(S_{n-1}B\left((1-n/2)\epsilon^{-n/2}H_{n/2-1}^{(1)}(k\epsilon)+k\epsilon^{1-n/2} H_{n/2-1}^{(1)'}(k\epsilon) \right)\right)=-1 \tag 5$$

---

FINDING $B$ BY EVALUATING THE LIMIT IN $(5)$:

To evaluate $(5)$, we use the small argument asymptotic approximation for the Hankel function and its derivative. These are

$$\begin{align}H_{n/2-1}^{(1)}(k\epsilon)&\sim -i \frac{\Gamma(n/2-1)}{\pi}\left(\frac{2}{k\epsilon}\right)^{n/2-1} \tag 6\\\\ H_{n/2-1}^{(1)'}(k\epsilon)&\sim i \frac{(n/2-1)\Gamma(n/2-1)}{2\pi}\left(\frac{2}{k\epsilon}\right)^{n/2} \tag 7 \end{align}$$

Putting together $(5)$, $(6)$, and $(7)$ we find that

$$\bbox[5px,border:2px solid #C0A000]{B=\frac i4\left(\frac{k}{2\pi}\right)^{n/2-1}}\tag 8$$

---

PUTTING IT ALL TOGETHER:

Using $(8)$ in $(3)$ and setting $A=0$, we find that

$$G(\vec r|\vec r'=0)=\frac i4 \left(\frac{kr}{2\pi}\right)^{n/2-1}H_{n/2-1}^{(1)}(kr)$$

Finally, shifting the point source to $\vec r'$, we obtain the Green Function for the $n$-Dimensional Helmholtz Equation

$$\bbox[5px,border:2px solid #C0A000]{G(\vec r|\vec r')=\frac i4 \left(\frac{k}{2\pi |\vec r-\vec r'|}\right)^{n/2-1}H_{n/2-1}^{(1)}(k|\vec r-\vec r'|)}\tag 9$$

---

NOTES:

For $n=2$, we recover the familiar $2$-Dimensional Green Function, $\frac i4 H_0^{(1)}(k|\vec r-\vec r'|)$, from $(9)$.

And for $n=3$, recalling that the $1/2$ order Hankel Function is related to the Spherical Bessel Function, we recover the familiar $3$-Dimensional Green Function $\frac{e^{-ik|\vec r-\vec r'|}}{4\pi|\vec r-\vec r'|}$ from $(9)$.

Answer 2

@Manuel Cañizares and @Mark Viola: I found some original work in [1, p231] about the solution steps above, with the scaled homogeneous solution of the corresponding Bessel differential equation, including Sommerfeld's radiation condition to select the physical solution, and finding the scaling B. But it appears that also Sommerfeld did not quite explain how to get the scaling criterion $$ \lim_{r\rightarrow0}\int_{\mathbb{S}_{n-1}}\frac{\partial G}{\partial r}\;r^{n-1}\mathrm{d}\Omega=-1. $$ I think Mark Viola might be wrong about the divergence theorem not applying, because $r$ only needs to approach zero. Anyways, I believe it leads us safely there. In general we integrate $$ \int_V\left[ (\Delta+k^2) G = -\delta \right]\;\mathrm{d}V $$ which yields with the divergence theorem $$ \int_S \frac{\partial G}{\partial n} \;\mathrm{d}S +k^2\int_V G\;\mathrm{d}V= -1 $$ for which it is sufficient to show that the limit towards a ball with small radius annihilates the second term $$ k^2\,S_{n-1}\lim_{r\rightarrow0}\int_0^r G\,r^{n-1}\mathrm{d}r=0, $$ since $\lim_{r\rightarrow0}G\propto r^{-\frac{n-1}{2}}$, and so only $S_{n-1}\,\lim_{r\rightarrow0}\frac{\partial G}{\partial r}=-1$ remains.

I am very thankful for Mark Viola's original answer with the solution steps above, and it was a useful hint in my research that made me find the missing pieces of information that I could not find in textbooks I have, or the more original kind of source [1]. This should now also answer the question that Manuel Cañizares faced.

[1] Sommerfeld, Arnold, Partial differential equations in physics. - Translated by Ernst G. Straus, Pure and Applied Mathematics (Academic Press) 1. New York: Academic Press. xi, 335 p. (1949), MR29463, Zbl 0034.35702.

---

Also [2, p. 342] is quite interesting, because it deals with the Fourier method and its resulting complex contour integral. It is interesting because it does not require to transform the equation into a Bessel differential equation, of which a physical solution needs to be selected and scaled right.

While the reference [2] explains the trick for $n=2$, it is difficult to grasp and it is not all clear.

The Fourier approach basically yields, with $\nu=\frac{n-2}{2}$, $$ \begin{align} G&=\frac{1}{(2\pi)^n} \int \frac{ e^{\mathrm{i}\boldsymbol{k}^\intercal\boldsymbol{x}} }{k^2-\frac{\omega^2}{c^2}} \mathrm{d}^n\boldsymbol{k}\\ &=\frac{1}{2\pi\,(2\pi\,r)^\nu} \int_0^\infty \frac{J_\nu(kr)}{k^2-\frac{\omega^2}{c^2}}k^{\nu+1} \mathrm{d}k\\ &= \frac{1}{4\pi\,(2\pi\,r)^\nu} \int_0^\infty\left[ \frac{k^{\nu}J_\nu(kr)}{k-\frac{\omega}{c}} + \frac{k^{\nu}J_\nu(kr)}{k+\frac{\omega}{c}} \right] \mathrm{d}k\\ &=\frac{\lim_{\epsilon\rightarrow0}}{8\pi\,(2\pi\,r)^\nu} \int_0^\infty\left\{ \frac{k^{\nu}[H^{(1)}_\nu(kr)+H^{(2)}_\nu(kr)]}{k-\frac{\omega}{c}-\mathrm{i}\epsilon} + \frac{k^{\nu}[H^{(1)}_\nu(kr)+H^{(2)}_\nu(kr)]}{k+\frac{\omega}{c}+\mathrm{i}\epsilon} \right\}\,\mathrm{d}k \end{align} $$ which can be closed to a contour by an infinte quarter circle from $\infty$ to $\pm\mathrm{i}\infty$ that vanishes exponentially either for $H_\nu^{(1)}(\mathrm{i}kR)\propto e^{-kR}$ and large positive imaginary values or for $H_\nu^{(2)}(-\mathrm{i}kR)\propto e^{-kR}$ and large negative ones. The clue is that both ways back to $0$ on the respective imaginary half axes have opposing signs $(-z)^{\nu}H_\nu^{(2)}(-zr)=-z^{\nu}H_\nu^{(1)}(zr)$ so that the paths $\lim_{R\rightarrow\infty} \left[\int_{\mathrm{i}R}^0A(z)\,\mathrm{d}z+\int_{-\mathrm{i}R}^0B(z)\,\mathrm{d}z\right]=\lim_{R\rightarrow\infty} \int_{\mathrm{i}R}^0 [A(z)-B(-z)]\, \mathrm{d}z$ cancel in pairs, the 1$^\mathrm{st}$ integrand $\frac{z^{\nu}H_\nu^{(1)}(zr)}{z-\frac{\omega}{c} -\mathrm{i}\epsilon}$ with the 4$^\mathrm{th}$ one $-\frac{(-z)^{\nu}H_\nu^{(2)}(-zr)}{-z+\frac{\omega}{c} +\mathrm{i}\epsilon}$ as well as the $2^\mathrm{nd}$ one $-\frac{(-z)^{\nu}H_\nu^{(1)}(-zr)}{-z-\frac{\omega}{c} -\mathrm{i}\epsilon}$ with the $3^\mathrm{rd}$ one $\frac{z^{\nu}H_\nu^{(1)}(zr)}{z+\frac{\omega}{c} +\mathrm{i}\epsilon}$. Therefore the two paths completing the contours $\int_{C_\pm}\mathrm{d}z=\lim_{R\rightarrow\infty}\left[\int_0^R\mathrm{d}z+\int_0^{\pm\frac{\pi}{2}}Re^{\mathrm{i}\phi}\mathrm{d}\phi +\int_{\pm\mathrm{i}R}^0\mathrm{d}z\right]$ in addition to the real semi-axis do not add anything to the value of the integral, but their formal closure permits utilizing Cauchy's integral formula that yields a vanishing result for the singularity $-\frac{\omega}{c}-\mathrm{i}\epsilon$ that is not contained in any of the contours $C_\pm$, while $\frac{\omega}{c}+\mathrm{i}\epsilon$ is either exclusively enclosed by the contour $C_+$ of the 1$^\mathrm{st}$ integral $\frac{\lim_{\epsilon\rightarrow0}}{8\pi\,(2\pi\,r)^\nu}\int_{C_+}\frac{z^\nu\,H_\nu^{(1)}(zr)}{z-\frac{\omega}{c}-\mathrm{i}\epsilon}\,\mathrm{d}z$ for $\epsilon>0$ or by the contour $C_-$ of the 2$^\mathrm{nd}$ integral $\frac{\lim_{\epsilon\rightarrow0}}{8\pi\,(2\pi\,r)^\nu}\int_{C_-}\frac{z^\nu\,H_\nu^{(2)}(zr)}{z-\frac{\omega}{c}-\mathrm{i}\epsilon}\,\mathrm{d}z$ for $\epsilon<0$. The result is causal in the time domain when choosing the right sign for the vanishing, causal loss $\epsilon$, and stable because we took care of converging integrands as required for existence Fourier representation. I actually prefer $\epsilon<0$ for which the 2$^\mathrm{nd}$ integral yields $\frac{-\mathrm{i}}{4}\left(\frac{k}{2\pi\,r}\right)^\nu\,H^{(2)}_\nu(kr)$, out of signal processing conventions that expand into time by $e^{+\mathrm{i}\omega t}$.

Same solution as in the above transformation of the inhomogeneous Helmholtz equation into an inhomogeneous Bessel differential equation, but the required physical conditions are just imposed by inserting vanishing losses for a causal time function, without requiring to know Sommerfeld's radiation condition. (Bernhard Schnizer's Analytical methods in applied physics lecture notes also turned out to be quite helpful)

[2] Sommerfeld, A., Die Greensche Funktion der Schwingungsgleichung,, Deutsche Math.-Ver. 21, 309-353 (1912). JFM 43.0819.01.