I am looking at the following scalar product of divergences in $4$ dimensions (in momentum space):
$$\partial_{p_1} \cdot \partial_{p_2} \frac{1}{p_1^2p_2^2}. \tag{1}$$
There is a further integration over $p_1$ and $p_2$ on the full space. I mention it because it means that $p_1$ and $p_2$ can coincide, and an especially interesting case of that is $p_1^2=p_2^2=0$. Now I am wondering if I can extract a Dirac $\delta$-function, since:
$$\partial_p^2 \frac{1}{p^2} = - (2\pi)^2 \delta^4(p). \tag{2}$$
($1/p^2$ is the Green's function of $\partial_p^2$ in $4$ dimensions) and
$$\partial_p^2 \frac{1}{p^4} = \partial_{p_\mu}\partial_{p_\mu} \left( \frac{1}{p^2}\frac{1}{p^2} \right) = \frac{8}{p_2^6} - \frac{8\pi^2}{p_2^2} \delta^{4}(p). \tag{3}$$
(using the product rule to isolate a $\partial_p^2 (1/p^2)$).
It seems to me that at $p_1=p_2=0$, I have a somewhat similar situation, and ideally I would love to write $(1)$ in the following form:
$$\begin{align}\partial_{p_1} \cdot \partial_{p_2} \frac{1}{p_1^2p_2^2} &= \text{(terms valid at $p_1,p_2\neq0$)} \cr &+ \text{(term with $\delta$-function for case $p_1=p_2=0$)}.\end{align} \tag{4}$$
Is that possible? I think not, but it would be so useful in the computation I am doing that it is worth asking.
If we regularize OP's expressions as a smooth function in $C^{\infty}(\mathbb{R}^4)$, in the sense of generalized functions$^1$ $$ \frac{1}{p^2+\varepsilon}~\rightarrow~ {\rm P.V.}\frac{1}{p^2} \quad\text{for}\quad\varepsilon\to 0^+, \tag{2i}$$ then the derivatives $\partial_{\mu} ~:=~\frac{\partial}{\partial p_{\mu}}$ are well-defined: $$ \partial_{\mu}\frac{1}{p^2+\varepsilon}~=~- \frac{2p_{\mu}}{(p^2+\varepsilon)^2}~\rightarrow~ -{\rm P.V.}\frac{2p_{\mu}}{(p^2)^2} \quad\text{for}\quad\varepsilon\to 0^+, \tag{2ii}$$ $$ \partial^2\frac{1}{p^2+\varepsilon}~=~ -\frac{8\varepsilon}{(p^2+\varepsilon)^3} ~\rightarrow~ (2\pi)^2\delta^4({\bf p}) \quad\text{for}\quad\varepsilon\to 0^+. \tag{2iii}$$ $$\begin{align} \partial_{(1)}\cdot \partial_{(2)}\frac{1}{p^2_{(1)}+\varepsilon}\frac{1}{p^2_{(2)}+\varepsilon}~=~& \frac{4p_{(1)}\cdot p_{(2)}}{(p^2_{(1)}+\varepsilon)^2(p^2_{(2)}+\varepsilon)^2}\cr ~\rightarrow~& {\rm P.V.}\frac{4p_{(1)}\cdot p_{(2)}}{(p^2_{(1)})^2(p^2_{(2)})^2} \quad\text{for}\quad\varepsilon\to 0^+ .\end{align} \tag{1/4}$$ In particular, there are no Dirac delta like contributions in eq. (1/4), cf OP's question. On the other hand, $$ \partial_{\mu}\frac{1}{(p^2+\varepsilon)^2}~=~- \frac{4p_{\mu}}{(p^2+\varepsilon)^3} , \tag{3i}$$ $$ \partial^2\frac{1}{(p^2+\varepsilon)^2}~=~ \frac{8p^2-16\varepsilon}{(p^2+\varepsilon)^4}, \tag{3ii}$$ are too singular to make sense as generalized functions.
--
$^1$Here ${\rm P.V.}$ is the Cauchy principal value and $p^2=\sum_{\mu=1}^4 p_{\mu}^2$.