Green's function for the inhomogenous Klein-Gordon equation I'm trying to solve the massive Klein-Gordon equation in good old Minkowski space-time: $$(\square + m^2) \phi = \rho(t,\mathbf{x})$$ where $\square = \partial_{\mu} \partial^{\mu} = \partial_{t}^2 - \nabla^2$. So one can use a Green's function approach to find the fundamental solution of the form $$(\square + m^2) \mathscr{G}_{m} = \delta(x^{\mu} - x'^{\mu})$$ where $\mathscr{G}_{m}$ is the familar Klein-Gordon propagator. One then obtains the solution $\phi$ in position space, given as the familiar solution $$\phi(x^{\mu}) = \int d^{4} x' \mathscr{G}_{m}(x^{\mu},x'^{\mu}) \rho(t',\mathbf{x}') \,\,\,\,\,\,\,\,\,\,\,\,\,\, (\star)$$ I was perfectly happy with this until I needed to to implement an actual $\rho$ and perform the integrals. My best bet so far has been to use the Bessel function representation I found (here: http://functions.wolfram.com/Bessel-TypeFunctions/BesselJ/31/02/) (which I've assumed has a generalisation) of the form: $$\mathscr{G}_{m}(t,t',\mathbf{x},\mathbf{x}') = \frac {\theta(t-t')} {2 \pi} \delta\Big( (t-t')^2 - |\mathbf{x} - \mathbf{x'}|^2 \Big) - \frac {m} {2 \pi} \theta(t-t' - |\mathbf{x} - \mathbf{x}'|)  \frac {J_{1}(m \sqrt{(t-t')^2 - |\mathbf{x} - \mathbf{x}'|^2)}} {m\sqrt{(t-t')^2 - |\mathbf{x} - \mathbf{x}'|^2)}}$$ while this is a nice closed form representation, I am still having real difficulty evaluating the integral $(\star)$. I've looked for quite some time in various places for explicit examples of computing the integral, and have so far come up with very little. Mathematica (my computational program of choice) really disdains these Heaviside functions in the integrals, and offers little guidance. The only case I can do so far is $m \mapsto 0$.
Question: Using the representation of $\mathscr{G}_{m}$ given (or another nicer one), how can one actually go about calculating $(\star)$? Has anyone got a reference wherein some explicit example is calculated where $\rho$ goes beyond a simple $\delta$-function? Even something like $\rho = \rho(r,\theta)$ or $\rho = \rho(r)$ would be of great help.
Thanks!
 A: I'm aware this is an old question and may be considered somewhat obsolete, but let me answer it nevertheless - if only for the sake of completeness.
The position space representation of the Klein-Gordon Green function (propagator) clearly looks intimidating. The trick is to do the calculation in momentum space, where the propagator is just a rational function. Before actually doing this let me point out that in the massless case, $m = 0$, and for a static source, $\rho = \rho (\mathbf{x})$, one is actually solving the Poisson equation. If the source is radially symmetric, $\rho = \rho(r)$, as suggested in the question, the solution is the Coulomb potential, $\phi = \phi(r) \sim 1/r$. Allowing for non-vanishing mass, one obtains a Yukawa potential, $\phi \sim \exp(-mr)/r$.
This can be shown explicitly in terms of Fourier integrals. First, transform the field,
$$ \phi(k) = \int d^4 x  \, e^{i k\cdot x} \, \phi (x) \; , $$
and similarly the the source, $\rho \to \rho(k)$. The momentum space solution of the KG equation is then  
$$ \phi(k) = \frac{\rho(k)}{k^2 - m^2 + i \epsilon} \; , $$
with $k^2 = k_0^2 - \mathbf{k}^2$ and the causal $i\epsilon$-prescription. (Change appropriately for retarded/advanced solutions.) Then transform back to position space, 
$$ \phi (x) = \int \frac{d^4 k}{(2\pi)^4}  \, e^{-i k\cdot x} \, \frac{\rho(k)}{k^2 - m^2 + i \epsilon} \; . \quad (**)$$
As an example, take a Gaussian source, $\rho(r) = \rho_0  \exp(-\alpha r^2)$, with 'normalisation' $ \rho_0 = (\alpha/\pi)^{3/2}$. Its Fourier transfrom is $\rho(k) = 2\pi i \, \delta(k^0) \exp (-k^2/4\alpha)$. The $k^0$-integral in $(**)$ is hence trivial, and the $d^3 k$ integration can be done with the usual residue technique writing $\mathbf{k}^2 + m^2 \equiv (\kappa+im)(\kappa-im)$. The result is
$$ \phi(r) = e^{m^2/4\alpha} \, \frac{e^{-mr}}{4\pi r} \; .$$
In the point source limit ($\alpha \to \infty$) we reobtain the standard Yukawa potential. 
For time dependent sources there will be energy transfer (no $\delta(k^0)$), and the integral $(**)$ will normally be harder. Typically, one can do the $k^0$-integration via residues and the remaining one(s) using stationary phase as e.g. in Ch. 2.1 of Peskin and Schroeder. 
