Solution (source code)

= Solution

\b[Integral identity and representation.] Define $L=\nabla^2+k^2$ and take $n$ to be the outward unit normal. The <product rule> gives
$$
\nabla\cdot(\phi\nabla\psi-\psi\nabla\phi)=\phi\nabla^2\psi-\psi\nabla^2\phi=\phi L\psi-\psi L\phi.
$$
The $k^2\phi\psi$ terms cancel. Integrating and applying the <divergence theorem> proves the <Helmholtz equation> version of <Green's second identity>:
$$
\boxed{\int_V(\phi L\psi-\psi L\phi)\,dV=\int_S(\phi\partial_n\psi-\psi\partial_n\phi)\,dS.}
$$

Choose the <Dirichlet Green function> with the sign convention
$$
L_xG(x,x')=\delta(x-x'),\qquad G(x,x')=0\quad(x\in S),
$$
where $\delta$ is the three-dimensional <Dirac delta>. Substitute $\psi=G$, use $L\phi=f$, and evaluate the <Dirac delta> integral. The boundary term containing $G$ vanishes and $\phi=g$ on $S$, leaving
$$
\phi(x')-\int_VG(x,x')f(x)\,dV=\int_Sg(x)\partial_nG(x,x')\,dS.
$$
Hence
$$
\boxed{\phi(x')=\int_VG(x,x')f(x)\,dV+\int_Sg(x)\partial_nG(x,x')\,dS.}
$$
The choice of sign in the <Green function> equation matters: it fixes both the source normalization and these signs. For a bounded region this construction presumes that the chosen $k^2$ is not a <Dirichlet eigenvalue> of $-\nabla^2$; at such an eigenvalue a Green inverse for arbitrary data does not exist.

\b[Infinite space.] Put $r=|x-x'|$. Away from the source, spherical symmetry reduces the <Helmholtz equation> to
$$
G''+\frac2rG'+k^2G=0,\qquad (rG)''+k^2(rG)=0.
$$
For $k\ne0$ its general radial solution is
$$
G(r)=\frac{A\cos kr+B\sin kr}{r}.
$$
Integrate $LG=\delta$ over a ball of radius $\varepsilon$ about the source. The <divergence theorem> yields
$$
4\pi\varepsilon^2G'(\varepsilon)+k^2\int_{r<\varepsilon}G\,dV=1.
$$
Since $G=A/r+Bk+O(r)$, the volume term is $O(\varepsilon^2)$ and the flux tends to $-4\pi A$. Therefore $A=-1/(4\pi)$, while $B$ remains free:
$$
\boxed{G(r)=-\frac{\cos kr}{4\pi r}+B\frac{\sin kr}{r},\qquad r>0.}
$$
The $B$ term extends regularly through the source and solves the homogeneous <Helmholtz equation>, so it cannot change the unit <Dirac delta> normalization. A condition at infinity is needed to select it. For example, with time dependence $e^{-i\omega t}$, the <outgoing Green function for the three-dimensional Helmholtz equation> is $-e^{ikr}/(4\pi r)$, corresponding to $B=-i/(4\pi)$.

\b[Zero wavenumber.] For fixed finite $B$, $\sin(kr)/r\to0$ as $k\to0$, and
$$
\boxed{G(r)\longrightarrow-\frac1{4\pi r},\qquad\nabla^2G=\delta.}
$$
This is the decaying <Green function of the Laplacian>. Directly at $k=0$, the general radial solution is $A/r+C$; source normalization fixes $A=-1/(4\pi)$, and decay at infinity removes $C$. If the free parameter depends on $k$ as $B(k)\sim C/k$, the limit instead includes that arbitrary homogeneous constant. Thus the usual Laplace limit holds for a fixed parameter, or after a suitable condition at infinity has been selected.