Solution (source code)

= Solution

Use signature $(-,+,+,+)$, for which $\Box=-\partial_t^2+\Delta$. With the <Lorenz constraint in Proca theory>, the <Proca equation> is $(\Box-m^2)A^a=-\mu_0J^a$. A static solution therefore satisfies $(\Delta-m^2)A^a=-\mu_0J^a$.

For $R>0$ direct radial differentiation gives $(\Delta-m^2)(e^{-mR}/R)=0$. At the origin its $1/R$ singularity contributes $-4\pi\delta_0$: integrating the Laplacian over a small sphere gives limiting outward flux $-4\pi$, while the mass term's integral tends to zero. Thus $G_m(R)=e^{-mR}/(4\pi R)$ obeys $(\Delta-m^2)G_m=-\delta_0$. Convolution proves the printed <Yukawa potential> solution
$$
\boxed{A^a(x)=\frac{\mu_0}{4\pi}\int\frac{e^{-m|x-x'|}}{|x-x'|}J^a(x')\,d^3x'.}
$$
The integral requires suitable source decay, for example compact support. For conserved static sources, differentiating the convolution and integrating by parts gives $\partial_aA^a=0$ as well; its spatial divergence is the convolution with $\partial_iJ^i=0$.

For the usual nonnegative mass parameter this is the decaying static solution, and $m=0$ gives the Coulomb kernel with Lorenz gauge imposed. If the arbitrary real parameter is negative, the printed kernel still gives a particular solution for compact sources, but grows at infinity; the decaying boundary condition selects $e^{-|m|R}/(4\pi R)$ instead. Static sources also allow homogeneous time-dependent solutions, so this construction establishes a static solution, not that every solution is static.