Fundamental solution for Helmholtz equation in higher dimensions 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?
 A: 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)$.
