Solution 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)$.