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.
If at every time $t$, $\phi(\mathbf{x},t)$ is a nice enough function that it has a Fourier transform, then $$\phi(\mathbf{x},t)=\int_{-\infty}^{\infty}\frac{d^{3}k}{(2\pi)^{3}}\widetilde{\phi}(\mathbf{k},t)e^{i\mathbf{k}\cdot\mathbf{x}},$$ where $\widetilde{\phi}(\mathbf{k},t)$ is just the Fourier coefficient at that time $t.$
But now you ask that the whole function $\phi(\mathbf{x},t)$ (at every time) be a solution to the Klein-Gordon equation. Which means the function at different times need to have the right derivatives. If at every time there is a Fourier transform, then there are Fourier coefficients at every time. So if the wave evolves a certain way in time, then the Fourier coefficients need have values at different times need to be just right to let he wave have the right temporal derivatives.
OK, so the Klein-Gordon equation is second order so we can find the initial function $\phi(\mathbf{x},t=0)$ and its Fourier transform, call it $\theta(\mathbf{k})$ and we can take the initial temporal derivative $\partial_t \phi(\mathbf{x},t)\big\vert_{t=0}$ and it's Fourier transform, call it $\omega(\mathbf{k}).$ Then we know the initial Fourier coefficients and we know their derivatives and the second derivative is enforced to make it satisfy the Klein-Gordon equation, so
$\widetilde{\phi}(\mathbf{k},t)=\theta(\mathbf{k})\cos(E_kt)
+ \frac{1}{E_k}\omega(\mathbf{k})\sin(E_kt).$
Why? Because it has the right initial values and the right initial temporal derivative and it satisfies the Klein-Gordon equation when $E_k=\sqrt{m^2+k^2}.$
So at each time we have Fourier coefficients of our wave. They are designed so that the functions satisfies the Klein-Gordon equation when the Fourier coefficients evolve in time according to a second order equation.
I may have misunderstood your question. But the idea is that from the initial wave and the initial temporal derivative of the wave you get enough initial conditions to know the initial Fourier coefficients and the initial time rate of change of the Fourier coefficients which is all the freedom you have. The rest is determined that the Fourier coefficients have to evolve temporally a certain way so the wave evolves a certain way temporally.
edit to respond to comments
At each time you get a Fourier transform. And you then ask yourself how those Fourier coefficients depend (in time) on each other. In order for the wave to evolve in time by a second order equation, the Fourier transform needs to evolve in time by a second order equation. But the transform gets to do a pointwise evolution that is simpler and that's why we do it.
When you take the Fourier transform at each time, you get a Fourier transform at each time, so you get a $\widetilde{\phi}(\mathbf{k},t)$ that is itself a function of time. So it has a partial derivative $\partial_t \widetilde{\phi}(\mathbf{k},t)\neq 0.$
edit to further respond to comments
Let's say we want to solve $$\left(\partial_t^2-\nabla^2+m^2\right)\phi=.0$$ We start by noting that $\left(A_ke^{i\mathbf{k}\cdot\mathbf{x}}+B_ke^{-i\mathbf{k}\cdot\mathbf{x}}\right)\cos(t\sqrt{m^2+k^2})$ and $\frac{\left(C_ke^{i\mathbf{k}\cdot\mathbf{x}}+D_ke^{-i\mathbf{k}\cdot\mathbf{x}}\right)}{\sqrt{m^2+k^2}}\sin(t\sqrt{m^2+k^2})$ are both solutions for any $\mathbf{k},$ any $A_k, B_k$ and any $C_k, D_k$.
Great.
And the first solution has zero time rate of change at $t=0$ and the second has zero value at $t=0.$ So if our initial wave was $\phi(\mathbf{x},t=0)=3\cos(3x)$ and had an initial time rate of change equal to $\partial_t\phi(\mathbf{x},t)\big\vert_{t=0}=4\cos(2y)$ then we know exactly what the solution to $\left(\partial_t^2-\nabla^2+m^2\right)\phi=0$ is: $$\phi(\mathbf{x},t)=3\cos(3x)\cos(t\sqrt{m^2+3^2})+\frac{4}{\sqrt{m^2+2^2}}\cos(2y)\sin(t\sqrt{m^2+2^2}).$$
And any finite linear combination of $e^{\pm i\mathbf k \cdot \mathbf x}$ for the initial value $\phi(\mathbf{x},t=0)$ is as easily manageable by having a finite linear combination of solutions like $e^{\pm i\mathbf{k}\cdot\mathbf{x}}\cos(t\sqrt{m^2+k^2}).$ And similarly if the initial time rate of change $\partial_t\phi(\mathbf{x},t)\big\vert_{t=0}$ is a finite linear combination of $e^{\pm i\mathbf k \cdot \mathbf x}$ then we add terms that are a finite linear combination of $\frac{e^{\pm i\mathbf{k}\cdot\mathbf{x}}}{\sqrt{m^2+k^2}}\sin(t\sqrt{m^2+k^2}).$
All we are doing is taking solutions and adding them up in combinations that give us the right initial values and right initial time rate of changes. And it is super easy if the initial conditions happen to be a finite linear combination of sines and cosines.
But wait. What if instead of being a finite linear combination of terms like $e^{\pm i\mathbf k \cdot \mathbf x}$ the initial condition was just a function that has a Fourier Transform? Then you can try the same thing. Write your initial value $\phi(\mathbf{x},t=0)$ and your initial time rate of change $\partial_t\phi(\mathbf{x},t)\big\vert_{t=0}$ as inverse Fourier transforms of sines and cosine spatially. Then replace a spatial $\sin(\mathbf k \cdot \mathbf x)$ with the function $\sin(\mathbf k \cdot \mathbf x)\cos(t\sqrt{m^2+k^2})$ and replace the spatial function $\cos(\mathbf k \cdot \mathbf x)$ with the function $\cos(\mathbf k \cdot \mathbf x)\cos(t\sqrt{m^2+k^2}).$ Why? Because each of those satisfies the Klein-Gordon equation. And so for $t=0$ the inverse Fourier transform of those will be the initial value of the wave $\phi(\mathbf{x},t=0).$ So you are taking a $t=0$ spatial Fourier transform of $\phi(\mathbf{x},t=0)$ then replacing every spatial Fourier component $e^{\pm i \mathbf k \cdot \mathbf x}$ with $e^{\pm i \mathbf k \cdot \mathbf x}\cos(t\sqrt{m^2+k^2}).$ This solves the Klein-Gordon equation, has the right initial values and has $\partial_t\phi(\mathbf{x},t)\big\vert_{t=0}=0.$
Next, take your initial time rate of change $\partial_t\phi(\mathbf{x},t)\big\vert_{t=0}$ as an inverse Fourier transforms of sines and cosine spatially. Then replace a spatial $\sin(\mathbf k \cdot \mathbf x)$ with the function $\frac{\sin(\mathbf k \cdot \mathbf x)}{\sqrt{m^2+k^2}}\sin(t\sqrt{m^2+k^2})$ and replace the spatial function $\cos(\mathbf k \cdot \mathbf x)$ with the function $\frac{\cos(\mathbf k \cdot \mathbf x)}{\sqrt{m^2+k^2}}\sin(t\sqrt{m^2+k^2}).$ Why? Because each of those satisfies the Klein-Gordon equation. And yet for $t=0$ if you take the inverse spatial Fourier transform you get a function that has value zero at $t=0$ and has an initial time rate of change that equals $\partial_t\phi(\mathbf{x},t)\big\vert_{t=0}$ of the wave $\phi(\mathbf{x},t=0).$
Why have two solutions? Because this second one also solves the Klein-Gordon equation and has an initial value of zero but has an initial time rate of change that equals $\partial_t\phi(\mathbf{x},t)\big\vert_{t=0}.$ And the first one solves the Klein-Gordon equation, has an initial value that equals $\phi(\mathbf{x},t=0)$ and has a zero initial time rate of change.
So if you add those two solutions together you get a function that (1) solves the Klein-Gordon equation (2) has the right initial value and (3) has the right initial time rate of change. That's what you wanted all along.
If you understand that when the initial wave was $\phi(\mathbf{x},t=0)=3\cos(3x)$ and had an initial time rate of change equal to $\partial_t\phi(\mathbf{x},t)\big\vert_{t=0}=4\cos(2y)$ then he solution was $\left(\partial_t^2-\nabla^2+m^2\right)\phi=0$ is: $$\phi(\mathbf{x},t)=3\cos(3x)\cos(t\sqrt{m^2+3^2})+\frac{4}{\sqrt{m^2+2^2}}\cos(2y)\sin(t\sqrt{m^2+2^2}).$$ If you understand that, then everything else is the same idea applied to finite and "infinite" linear combinations of $e^{\pm i \mathbf k \cdot \mathbf x}$ for the initial conditions.
Best Answer
A nice way to see how the correlation function behaves is described here where it is shown that the propagator goes as $$C(r)=\frac{1}{2\pi}\log(r)$$ which can also be seen as given by Qmechanics hint. Now the interesting thing is not that it diverges at $r=0$ (this happens even in 4D where $C(r)=1/4 \pi^2 r^2$) but that it also diverges as $r\to \infty$. This is an infrared divergence I had not come across before. The Wikipedia article linked above states that this makes a two dimensional massless scalar field slightly tricky to define mathematically and also that you cannot have spontaneous breaking of a continuous symmetry in two dimensions. Very interesting!