Let $B_n$ be the ball with radius $n$, $K_n=C_c^\infty(B_n)$ with its metrizable topology, $\varphi_n\in K_n$ a function with support contained in $B_{n}$ and $\varphi_n(x)=1$ for $x\in B_{n-1}$. First observe that $$
F_n\colon K_n\times K_n \to K_n $$ is a continuous map, which can be easily seen by the defining seminorms for these metric spaces.
Now let $U$ be a convex neighbourhood of $0$, i.e. $U\cap K_n$ is a convex neighbourhood of $0$ in $K_n$ for each $n$. Inductively for each $n$, you can find a $0$-neighbourhood $V_n$ of $K_n$ such that $$
F[V_n,V_n] \subseteq U\cap K_n $$
(by the continuity of $F_n$) and $$ \varphi_k V_n \subseteq V_k\,\,\,\,\, (1\leq k < n).$$
Set $W_n:=V_n\cap K_{n-1}$ and $W$ as the convex hull of $\bigcup_n W_n$. Observe that for each $n$, $W_n$ is neigbourhood of $0$ in $K_{n-1}$, so $W\cap K_{n-1}\supseteq W_n$ is one too, hence $W$ is a neighbourhood of $0$ in $C_c^\infty(\mathbb{R}^d)$. Now $F[W,W]\subseteq U$ would establish the continuity of $F$.
Let $\psi, \chi\in W$, i.e. $\psi=\alpha_1\psi_1+\cdots + \alpha_m\psi_m$ and $\chi=\beta_1 \chi_1 + \cdots + \beta_m \chi_m$ with $\alpha_i, \beta_i\geq 0$, $\sum \alpha_i = \sum \beta_i =1$ and $\psi_i,\chi_i\in V_i$. As $$
F(\psi,\chi)=\psi\cdot \chi = \sum_{i,j} \alpha_i\beta_j \cdot \psi_i\chi_j $$ and $\sum_{i,j} \alpha_i\beta_j = 1$, it it sufficient to verify $\psi_i\chi_j\in U$. Now if $i=j$,
$$ \psi_i\chi_i = F(\psi_i,\chi_i)\in F[V_i,V_i]\subseteq U\cap K_i \subseteq U.$$
If $i\neq j$, e.g. $i<j$, then $\psi_i\in V_i$ and $\chi_j\in V_j$ and so $$
\psi_i\chi_j = (\psi_i \varphi_i) \chi_j =\psi_i (\varphi_i\chi_j)\in V_i\cdot V_i \subseteq U\cap K_i \subseteq U.$$
There were some errors in the OP. First, the Fourier Transform of the Helmholtz Equation (Equation $(1)$ in the OP) is given by
$$(-k'^2+k^2)\hat u(\vec k')=-1$$
where $\vec k'$ is the $3$-D transform variable with $k'=|\vec k'|$.
Second, it then follows that $$\hat u(\vec k')=\frac{1}{k'^2-k^2}\tag 1$$ In what follows, we propose a way forward to take the inverse Fourier Transform of $\hat u(\vec k')$ as given by $(1)$.
METHODOLGY $(1)$:
One way to handle the singularities at $k'=\pm k$ is to let $k$ have a "small," non-zero imaginary part. We will assume that $\text{Im}(k)<0$.
Then, we can write
$$\begin{align}
u(\vec r)&=\frac{1}{(2\pi)^3}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty} \frac{e^{i\vec k'\cdot \vec r}}{k'^2-k^2}\,dk'_x\,dk'_y\,dk'_z\\\\
&=\frac{1}{(2\pi)^3}\int_0^{2\pi}\int_0^\pi\int_0^\infty \frac{e^{ik'r\cos(\theta)}}{k'^2-k^2}\,k'^2\sin(\theta)\,dk'\,d\theta\,d\phi\\\\
&=\frac{1}{(2\pi)^2}\int_0^\pi\int_0^\infty \frac{e^{ik'r\cos(\theta)}}{k'^2-k^2}\,k'^2\sin(\theta)\,dk'\,d\theta\\\\
&=\frac{1}{(2\pi)^2}\int_0^\infty \frac{k'}{k'^2-k^2}\frac{2\sin(k'r)}{r}\,dk'\\\\
&=\frac{1}{(2\pi)^2}\int_{-\infty}^\infty \frac{k'}{k'^2-k^2}\frac{\sin(k'r)}{r}\,dk'\\\\
&=\frac{e^{-ikr}}{4\pi r}
\end{align}$$
Had we assumed that $\text{Im}(k)>0$, we would have recovered the solution $\frac{e^{ikr}}{4\pi r}$.
Finally, we can let the imaginary part of $k$ approach $0$ and recover the result for real $k$.
Recommended references include (i) "Field Theory of Guided Waves," Robert Collin, IEEE Press, (ii) "Waves and Fields in Inhomogeneous Media," W.C. Chew, Van Nostrand Reinhold, and (iii) Radiation and Scattering of Waves, Felsen and Marcuvitz, IEEE Press.
METHODOLGY $(2)$:
We could solve Equation $(1)$ in the OP without the use of integral transformation. Instead, we write
$$\nabla^2 u(\vec r)+k^2u(\vec r)=0$$
for $\vec r\ne 0$. Exploiting the spherical symmetry of the problem, we have
$$\nabla^2 u+k^2 u=\frac1{r^2}\frac{\partial }{\partial r}\left(r^2\frac{\partial u}{\partial r}\right)+k^2u=0$$
which can be rearranged as
$$\frac{\partial^2}{\partial r^2}(r u)+k^2 (ru)=0 \tag 2$$
Solutions to $(2)$ are trivially seen to be
$$u=C^{\pm}\frac{e^{\pm ikr}}{r}$$
We find the constant $C^\pm$ by enforcing the condition $\int_{|\vec r|\le \epsilon}\nabla^2 u(\vec r)\,dV=-1$.
Applying the Divergence Theorem in the sense of distributions reveals
$$\begin{align}
\oint_{|\vec r|=\epsilon}\left.\left(\frac{\partial u(\vec r)}{\partial r}\right)\right|_{|\vec r|=\epsilon}\,\epsilon^2 \sin(\theta)\,d\theta\,d\phi&=-4\pi C^{\pm}\\\\
&=-1\end{align}$$
whereupon solving for $C^{\pm}$ yields $C^\pm=\frac{1}{4\pi}$.
Best Answer
The usual formula gives $$ \Delta^{-1}(x) = \frac{|x|^{2-d}}{\gamma_{d}(2-d)} $$ as the solution to $\Delta\left[\Delta^{-1}(x)\right]=\delta(x)$ in $d$ dimensions, where $$ \gamma_{d} = \frac{2\pi^{\frac{d}{2}}}{\Gamma\left(\tfrac{d}{2}\right)} $$ is the $d$-dimensional solid angle. Consider now $d=2-\epsilon$ for small $\epsilon$. Expanding the above expression we find $$ \Delta^{-1}(x) = \frac{1}{\epsilon} + \frac{1}{2\pi}\,\log|x| + \cdots $$ where the dots denote either terms that vanish as $\epsilon\to0$ or are independent of $|x|$. Dropping the $x$-independent terms and sending $\epsilon\to0$ we find $$ \Delta^{-1}(x)=\frac{1}{2\pi}\,\log|x| $$ for $d=2$.
Indeed, consider $$ \frac{1}{2\pi}\int_{\mathbb R^2} \log|x|\ \Delta \chi (x)\,dx = \lim_{\epsilon\to0}\frac{1}{2\pi}\int_{\mathbb R^2\setminus B_\epsilon} \log|x|\ \Delta \chi (x)\,dx $$ for any smooth test function $\chi$, where $B_\epsilon$ denotes a disk centered at the origin with radius $\epsilon$. Integrating by parts twice we have: $$\begin{aligned} \lim_{\epsilon\to0} \frac{1}{2\pi}\int_{\mathbb R^2\setminus B_\epsilon} \log|x|\ \Delta \chi (x)\,dx &= 0 - \lim_{\epsilon\to0}\frac{1}{2\pi}\int_{\mathbb R^2\setminus B_\epsilon} \frac{1}{|x|^2} \langle x, \nabla\chi (x)\rangle\,dx\\ &= + \lim_{\epsilon\to0} \frac{1}{2\pi\epsilon}\int_{\partial B_\epsilon}\chi(x)dS(x) + 0 = \chi(0)\,. \end{aligned} $$ Therefore, $$ \Delta^{-1}(x) = \frac{1}{2\pi}\,\log|x|\,. $$