So the inhomogeneous Helmholtz's equation is,
$$\nabla^2 f + \alpha^2 f = -\beta({\mathbf r}).$$
To solve this it helps to first solve the special case $\beta = 0.$ Now you must suppress your first instinct when you do this: your first instinct is to be lazy and assume that $f$ is nice and smooth everywhere. Instead we want to solve this for every $\mathbf r$ except for possibly a singularity at position $\mathbf r_0$. When we do this we are able to, for only a little bit more effort, solve for the case $\beta(\mathbf r) = \delta^3(\mathbf r - \mathbf r_0)$ for the so-called Green's function $f=G(\mathbf r, \mathbf r_0).$ This gives us superpowers because the above differential operator is linear and therefore the solution for general $\beta$ can just be written as $$f(\mathbf r) = \int d^3r_0~\beta(\mathbf r_0)~G(\mathbf r, \mathbf r_0).$$Now how shall we find $G$?
One of the nicest things we can do with this is to operate on the above equation with $\mathcal F_{\mathbf r\to\mathbf k} = \int d^3r~e^{-i{\mathbf k}\cdot{\mathbf r}},$ the 3D Fourier transform. Let me define $G[\mathbf k] = \mathcal F_{\mathbf r\to\mathbf k} G(\mathbf r, \mathbf r_0).$ When we do this we find that we can integrate derivatives by parts so that with suitable decay off at infinity e.g. $$\int dx~e^{-i~k_x~x} ~\partial_x~G = 0-\int dx~e^{-i~k_x~x}~(-i~k_x)~G$$ and thus $\mathcal F\Big(\nabla^2 G(\mathbf r, \mathbf r_0)\Big)=\big((i~k_x)^2 + (i~ k_y)^2 + (i~k_z)^2\big) = -\mathbf k^2G[\mathbf k].$ So in $\mathbf k$-space the above looks like:$$\left(-\mathbf k^2 + \alpha^2\right) G[\mathbf k] = -e^{-i\mathbf k\cdot \mathbf r_0},$$where we've simply done the integral with $\beta = \delta^3(\mathbf r - \mathbf r_0)$ on the right-hand side. Thus if you know the trick of looking at an equation in Fourier space you will simply reflexively write down the inverse-Fourier-transformed solution, that $$G(\mathbf r, \mathbf r_0) = \int \frac{d^3k}{(2\pi)^3} ~ \frac{e^{i~\mathbf k\cdot
(\mathbf r - \mathbf r_0)}}{\mathbf k^2 - \alpha^2}.$$
