$1/r^2$ force and "photon decay" Background: My question is about the interpretation of the $1/r^2$ force in terms of the fundamental processes of the underlying QFT. We know that if the photon had mass $m$, then we may have a "massive" version of QED that is still renormalizable but has no gauge freedom (the photon mass term kills the gauge symmetry explicitly).
Assuming such a theory, we would observe a modification of the Coulomb law between two charges, namely the Yukawa potential $~e^{-mr}/r$ instead of the usual $~1/r$.
Massless photon case: Now, the $1/r$ potential gives rise to a force that goes as $~1/r^2$. This force decreases as the area of the sphere of radius $r$, indicating that the force $~1/r^2$ is a purely geometrical fact of 3 spatial dimensions: intuitively, we emit $N$ photons from a point and the force decreases with the surface density of the photons on the sphere of radius $r$, see this answer. This works if $N$ is constant, i.e. there should be no diagrams that look like a "photon decay".
Massive photon case: Let's come back to the massive case. It looks like $N$ should decrease with $r$, in such a way as to reproduce the Yukawa potential. However, massive and massless QED look like having the same vertex diagrams (but different propagators).
Question: If the difference between the two theories is only in the photon's propagator, does it mean that the intuitive interpretation linking $1/r^2$ and the conservation of $N$ is misleading, or even wrong?
Is it possible to make the above intuitive argument more precise by considering the different behaviour of the photon number operators in the massless and massive cases?
Note: see this answer for the standard derivation of the $1/r$ electric potential. This question is about the interpretation in terms of photon exchange, not about this well-known calculation.
 A: *

*In a nutshell, to recover the static Yukawa potential between 2 points in 3 spatial dimensions, one should consider the 2-point correlation function/Greens function$^1$
$$G(\vec{r})
~=~\frac{e^{-mr}}{4\pi r}, \qquad (-\vec{\nabla}^2+m^2)G(\vec{r})~=~\delta^3(\vec{r}),\tag{1}
$$
in the static picture, cf. e.g. Wikipedia & this Phys.SE post.


*Example: For a scalar field in 3+0D, the partition function is
$$\begin{align} e^{-\beta F[J]}~=~~& Z[J]\cr
~=~~&\int \! {\cal D}\phi~\exp\left\{\beta \int\! d^3r  \left(-\frac{1}{2}\phi(-\vec{\nabla}^2+m^2)\phi +J\phi\right)\right\} \cr
~\stackrel{\begin{array}{c}\text{Gauss.}\cr\text{int.}\end{array}}{\sim}&
\exp\left\{\frac{\beta}{2}\int\! d^3r\int\! d^3r^{\prime}J(\vec{r})G(\vec{r}-\vec{r}^{\prime})J(\vec{r}^{\prime}) \right\}.
\end{align} \tag{2}$$
Given 2 Dirac delta function point-sources
$$ J(\vec{r})~=~q_1\delta^3(\vec{r}\!-\!\vec{r}_1) +q_2\delta^3(\vec{r}\!-\!\vec{r}_2),\tag{3}$$
the free energy (excluding self-interactions) is the Yukawa potential energy
$$ F[J]~\stackrel{(2)+(3)}{=}~-q_1q_2G(\vec{r}_1\!-\!\vec{r}_2).\tag{4}$$


*For what it's worth, in the "exchange particle decay picture" (that OP asks about) the mass $m$ can be interpreted as an attenuation coefficient of virtual exchange particles. In terms of Feynman diagrams, the mass-term in the action (2) can be seen as a 2-vertex. The only connected Feynman diagrams are necklaces where the mass 2-vertices play the role of pearls.
--
$^1$ Here we work in units where $\hbar=1=c$.
A: I find the "intuitive number" argument to be very misleading. One argument against it is... what is $N$? In fact the state of the electrostatic field around a point particle is a coherent state, which is a superposition of states with different numbers of particles. So it's not really possible to say that there are a definite number of photons "spreading out" from a central source. Additionally, as you point out, I don't see a good reason why the intuitive argument would fail for massive photons. Third, in the particle model, for finite $N$, you would expect the law to break down for large enough distance from the source -- you would expect to have "shot noise" in the force as the "rate" of photons arriving in a given spherical area decreased because the photons will have become spread very thin over a large spherical area; this has never been observed. Finally, the photons that appear in this kind of "static" force are actually virtual particles, not real particles, which is to say that they aren't really "propagating particles" but more like a colorful set of words to describe mathematical terms in a perturbative expansion, which don't need to obey things like continuity equations that hold on shell (e.g., when the classical equations of motion are satisfied).
The straightforward explanation for the $1/r^2$ vs Yukawa suppression, is simply that these force laws follow from the Green's functions (aka propagators) for the massless and massive Poisson equations, respectively.
If you want to look for a deeper explanation as to why massive and massless particles behave differently, in my view, the $1/r^2$ law is really a consequence of Gauss's Law
\begin{equation}
\nabla \cdot \vec{E} = \frac{\rho}{\epsilon_0} \implies E \sim \frac{1}{r^2}\ \ {\rm (point\ charge)}
\end{equation}
and in fact Gauss's law is deeply tied to the masslessness of the photon. In the Hamiltonian formulation of electrodynamics, Gauss's law is a first class constraint associated with the $U(1)$ gauge invariance of the photon. For the massive case, because the gauge invariance is broken, this first class constraint becomes two second class constraints, and Gauss's law no longer holds in its original form. In turn, this means that there is nothing that forces the field to fall off as $1/r^2$, and allows for the Yukawa solution to exist.
You may wonder about scalar fields, which do not have gauge invariance, but still have a difference between a $1/r^2$ behavior exhibited by massless scalar fields, and a Yukawa suppression for massive scalar fields. In fact there is also a symmetry explanation here. A massless scalar field has a shift symmetry under $\phi \rightarrow \phi + c$, which guarantees that the equations of motion can be written in the form of a total divergence. This total divergence plays the role of Gauss's law. A massive scalar breaks the shift symmetry and correspondingly has a different behavior.

A thought struck me while walking. I think this question comes from a similar place as the aether theory of the late 1800s -- an attempt to explain features of (what we now call) relativistic field theory with mechanical models. In your case you want to explain the inverse square law with particles spreading out in space, in the aether case people wanted to explain the medium in which electromagnetic waves are supposedly traveling. I think the modern response to both the "intuitive particle" argument and the aether idea are the same: we have found that mechanical models are both unnecessary and often misleading in describing relativistic field theories. What really matters are the equations satisfied by the fields. While these equations can be used to describe mechanical systems (for example Gauss's law can be used to describe some systems in fluid mechanics), the laws themselves are more general than the specific systems they can be applied to, and there is no need for there to be a mechanical model underlying these equations when they are applied in field theory. The structure matters more than the bricks out of which it is built.
A: The electrostatic potential in QED is represented by virtual photons, which in this case have [no] zero frequency and do not propagate. These photons do not transport energy. They're just Fourier components of the potential. A model of photon decay or density does it herefore not apply.
A: I think your question can really only be answered by observations. If the observed potential appears to be different from $1/r$ though, this would mean that you would not be looking at a pure interaction between two particles but that some 'screening' of the potential due to other effects occurs.
