Solution (source code)

= Solution
{c}

After <integration by parts>, the quadratic <Proca action> with an external source is
$$
S[A,J]=\frac12\int d^4x\,A_\mu K^{\mu\nu}A_\nu+\int d^4x\,J^\mu A_\mu,\qquad
K^{\mu\nu}=g^{\mu\nu}(\partial^2+M^2)-\partial^\mu\partial^\nu.
$$
Since $M>0$, this kernel is invertible with the <Feynman i-epsilon prescription>; the massless <gauge redundancy> is absent in this <unitary gauge> quadratic action. Complete the square by $A\mapsto A+K_F^{-1}J$ in the <Gaussian functional integral>. Its source-independent determinant cancels in the normalized <generating functional>:
$$
\boxed{\frac{Z[J]}{Z[0]}=\exp\left[-\frac i2\int d^4x\,d^4y\,J_\mu(x)(K_F^{-1})^{\mu\nu}(x-y)J_\nu(y)\right].}
$$
In <Fourier space>, the kernel is
$$
K^{\mu\nu}(k)=(M^2-k^2)g^{\mu\nu}+k^\mu k^\nu.
$$
For $k^2\ne0$, use the <transverse projector of a vector field> $P_T^{\mu\nu}=g^{\mu\nu}-k^\mu k^\nu/k^2$ and the <longitudinal projector of a vector field> $P_L^{\mu\nu}=k^\mu k^\nu/k^2$. Their mixed-index products obey $P_TP_L=0$, so
$$
K=(M^2-k^2)P_T+M^2P_L,\qquad
K^{-1}=\frac{P_T}{M^2-k^2}+\frac{P_L}{M^2}=-\frac{g^{\mu\nu}-k^\mu k^\nu/M^2}{k^2-M^2}.
$$
The last expression also extends to $k^2=0$; that point only obstructs the intermediate projector decomposition. Two source <functional derivatives> yield the <gauge-boson propagator> $D_F=iK_F^{-1}$:
$$
\boxed{D_F^{\mu\nu}(k)=\frac{-i}{k^2-M^2+i0}\left(g^{\mu\nu}-\frac{k^\mu k^\nu}{M^2}\right),\qquad
\langle T A^\mu(x)A^\nu(y)\rangle=\int\frac{d^4k}{(2\pi)^4}e^{-ik\cdot(x-y)}D_F^{\mu\nu}(k).}
$$
The pole residue is the physical massive-vector <massive vector polarization sum>. The longitudinal term is essential for inverting the quadratic kernel, and makes the naive $M\to0$ limit of this <unitary gauge> <Proca propagator> singular.