How to understand why massive fields decay exponentially with distance? I have read and heard numerous times that a vector field in nonabelian gauge theory which has gained mass through the Higgs mechanism decreases in strength exponentially with distance as it propagates. I would be thankful for a brief explanation of the approach and result and, if possible, a reference for a thorough treatment.
My attempt:
I expected to find evidence for this behavior in the massive vector boson propagator which can be written in the unitary gauge 
$$
\frac{g_{\mu \nu} - k_{\mu} k_{\nu}/m^2}{k^2-m^2}.
$$ 
I would like Fourier transform this to position space, so I integrate the following term by term
$$
\int_{-\infty}^{\infty} \frac{g_{\mu \nu} - k_{\mu} k_{\nu}/m^2}{k^2-m^2} e^{2\pi i k_{\mu}x^{\mu}} d^4k.
$$
I have three types of terms to integrate:
$$
\int_{-\infty}^{\infty} \frac{1}{k^2-a^2} e^{2\pi i kx} dk
$$
$$
\int_{-\infty}^{\infty} \frac{k}{k^2-a^2} e^{2\pi i kx} dk
$$
$$
\int_{-\infty}^{\infty} \frac{k^2}{k^2-a^2} e^{2\pi i kx} dk
$$
But I get divergences in these integrals and I am unsure where to go next. 
 A: Mathematica will make fairly short work of these integrals.  Or they can be looked up in a field theory text book.  I believe, for example, that the Fourier transform of $1/(k^2+m^2)$ is discussed in connection with scalar fields and the Klein-Gordon equation in the first few chapters of Peskin and Schroeder, at least in a Lorentzian setting.
Why the Fourier transform of this propagator is related to the Yukawa potential in the first place is a different question.  If I were to try to derive it, I would calculate tree level exchange of a massive vector in field theory.  Then I would try to figure out what potential in the Born approximation in non-relativistic quantum mechanics would lead to the same scattering.  A slightly different answer is here.  
Let's start with the first integral
$$
I = \int d^4k \frac{e^{i k\cdot x}}{k^2 + m^2} \ .
$$
Note I have flipped the sign of the $m^2$ term.  This is because I want to do the integral in Euclidean signature.  I can adapt the $k$-coordinate system so that $x$ points in the polar or "$z$''-direction.  I can then break up the measure factor in spherical coordinates
$$
d^4 k = k^3 dk \, \sin^2 \theta d \theta \, d \Omega_2 \ .
$$
The last $d\Omega_2$ is the measure on an $S^2$ with unit radius.  As nothing depends on these angles, they will integrate to give $4 \pi$.  Let's focus on the $\theta$ integral.  Mathematica tells us
$$
\int_0^\pi e^{i k x \cos \theta} \sin^2 \theta \, d \theta = \frac{\pi}{kx} J_1(kx) \ .
$$
Mathematica will also handle the final $k$-integral
$$
I =  \frac{4 \pi^2 m}{x} K_1(m x) \ .
$$
If we expand out the Bessel-K function for large argument, we find the desired exponential behavior
$$
I = e^{-m x} \left( \pi^{5/2} m^{1/2} \left(\frac{2}{x} \right)^{3/2} + O(x^{-2}) \right) \ .
$$
A trick for the $k_\mu k_\nu$ integral is that it can be obtained by taking $x$ derivatives of the first integral.  But this will give a similar exponentially damped behavior at large $x$.
