Suppose we had the Lagrangian:
$$\mathcal{L} = -\frac{1}{4} F^{\mu \nu}F_{\mu \nu} + \overline{\psi} (i \gamma^{\mu}\partial_{\mu} -m)\psi -e\overline{\psi} \gamma_{\mu} \psi A^{\mu} +\frac{1}{2} m_{\gamma}^2 A_{\mu} A^{\mu},$$
which is just the Lagrangian of QED with a massive photon in the $(+,-,-,-)$ sign convention. The lagrangian does not have a local $U(1)$ symmetry but it does have global $U(1)$ symmetry, therefore the Ward identity for global symmetries holds.
Let us consider the correlator:
$\mathcal{\Pi}^{R}_{\alpha \beta} (x) \equiv <T j^{\alpha}(x) j^{\beta}(0)>$,
In Feynman diagrams this correlator is the very much reducible full self-energy of the photon.
By the Ward identity we have:
$\partial^{\alpha}\mathcal{\Pi}^{R}_{\alpha \beta} (x) = -i\delta^4 (x) <\delta j(0)> = 0$,
the last passage is by the invariance under global $U(1)$ of the current. In Fourier space:
$q^{\alpha}\mathcal{\Pi}^{R}_{\alpha \beta} (q) = 0 \; \; \; \; \; (1)$
One can the write the fully dressed photon propagator as:
$G^{\mu \nu} (q) = G_{0}^{\mu \nu} (q) + G_{0}^{\mu \alpha} (q) \Pi^{R}_{\alpha \beta} (q) G_{0}^{\beta \nu} (q) \; \; \; \; \; \; (2) $
The explicit form of the bare photon propagator is:
$G_0^{\mu \nu} (q) = \frac{-i \left( g^{\mu \nu} – \frac{q^{\mu} q^{\nu}}{m_{\gamma}^2} \right)}{q^2-m_{\gamma}^2+i \varepsilon}$.
By contracting with $q_{\mu}$ both sides of equation $(2)$ we get:
$q_{\mu} G^{\mu \nu} (q) = q_{\mu} G_0^{\mu \nu} (q) + \frac{-i \left( 1 – \frac{q^2}{m_{\gamma}^2} \right)}{q^2-m_{\gamma}^2+i \varepsilon} q^{\alpha} \mathcal{\Pi}^{R}_{\alpha \beta} (q) G_0^{\beta \nu} (q) $,
and by $(1)$ the second term on the LHS is $0$, therefore:
$q_{\mu} G^{\mu \nu} (q) = q_{\mu} G_0^{\mu \nu} (q) \; \; \; \; \; (3)$.
But since by definition we have (calling $\mathcal{\Pi}_{\alpha \beta} (q) $ the 1PI photon self-energy):
$G^{\mu \nu} (q) = G^{\mu \nu}_0 (q) + G^{\mu \alpha}_0 (q) \mathcal{\Pi}_{\alpha \beta} (q) G^{\beta \nu} (q)$. By contracting both sides with $q_{\mu}$ we get:
$q_{\mu} G^{\mu \nu} (q) = q_{\mu} G_0^{\mu \nu} (q) + q_{\mu} G^{\mu \alpha}_0 (q) \mathcal{\Pi}_{\alpha \beta} (q) G^{\beta \nu} (q)$,
and by (3) we have:
$q_{\mu} G^{\mu \alpha}_0 (q) \mathcal{\Pi}_{\alpha \beta} (q) G^{\beta \nu} (q) = 0$, so that using the explicit expression of $G_0^{\mu \alpha}$ we get:
$q^{\alpha} \mathcal{\Pi}_{\alpha \beta} (q) G^{\beta \nu} (q) = 0$. Therefore we conclude:
$q^{\alpha} \mathcal{\Pi}_{\alpha \beta} (q) = 0$.
This looks so strange to me. Why should the photon self-energy be transverse even without gauge invariance? Is my math wrong somewhere?
EDIT: Something must be wrong. I tried developing a bit more the math to see where it would lead and I got some upsetting results. If one has $q_{\mu} \mathcal{M}^{\mu…}=0$, where $\mathcal{M}^{\mu}$ is any amplitude with on-shell fermions and one off-shell photon with 4-momentum $q$, then one can disregard the term proportional to $q^{\mu}q^{\nu}$ in the bare propagator. This makes the theory seemingly renormalizable. Also, for the same reason, one can write the full photon propagator as:
$G^{\mu \nu} = \frac{-ig^{\mu \nu}}{(q^2 – m_{\gamma}^2) \left[1-\Pi(q^2)\right]}$, which would mean that the photon mass is protected from renormalization, which is even weirder.
Best Answer
Here is one line of reasoning:
We can incorporate the Proca/massive photon field $A_{\mu}$ into a gauge theory via the Stuckelberg mechanism$^1$ $$ {\cal L}_S~=~-\frac{1}{4}F_{\mu\nu}F^{\mu\nu} -\frac{1}{2}(m_{\gamma} A^{\mu} +\partial^{\mu}\phi) (m_{\gamma} A_{\mu}+\partial_{\mu}\phi).$$ The fermionic matter is governed by the Dirac Lagrangian density $$ \begin{align} {\cal L}_D~=~&\bar{\psi}(i\gamma^{\mu}D_{\mu}-m_e)\psi,\cr D_{\mu} ~=~&\partial_{\mu}-ieA_{\mu} .\end{align} $$
It is possible to choose a $R_{\xi}$-gauge-fixing condition $$ G~=~\partial^{\mu}A_{\mu} + \xi m_{\gamma}\phi$$ so that the Stückelberg scalar $\phi$ becomes free (but massive). In detail, the gauge-fixed Lagrangian density is $$ {\cal L}_{GF}~=~-\frac{G^2}{2\xi},$$ while the Faddeev-Popov (FP) term reads $$ {\cal L}_{FP}~=~\bar{c}\left(\partial^2 -\xi m_{\gamma}^2 \right)c.$$ The total gauge-fixed Lagrangian density becomes$^2$ $$ \begin{align}{\cal L}~=~&{\cal L}_S+{\cal L}_D+{\cal L}_{FP}+{\cal L}_{GF}\cr ~\sim~&\frac{1}{2}A_{\mu}\left(\delta^{\mu}_{\nu} (\partial^2-m_{\gamma}^2)+ (\frac{1}{\xi}-1)\partial^{\mu} \partial_{\nu} \right)A^{\nu}\cr ~+~&\frac{1}{2}\phi(\partial^2-\xi m_{\gamma}^2)\phi +{\cal L}_D+{\cal L}_{FP} . \end{align} $$ Note in particular that the fermionic matter $\psi$ only couples to $A_{\mu}$, not to $\phi.$ In fact the $\bar{\psi}A_{\mu}\psi$ cubic vertex is the only interaction in the entire model! The Faddeev-Popov (FP) ghosts $c$ and $\bar{c}$ also decouple.
Then we can in principle construct the corresponding Ward identities for 1-particle irreducible (1PI) correlator functions, and show renormalizability. Explicit calculations show that the simplest such Ward identity confirms the transversality of the massive photon vacuum polarization tensor/self-energy, cf. OP's title question.
See also e.g. Diego Mazon's Phys.SE answer here.
--
$^1$ This answer uses the $(-,+,+,+)$ sign convention.
$^2$ The $\sim$ symbol means equality up to total derivative terms.