Consider the following representation of the 4-dimensional Dirac delta function,
\begin{equation}
\delta^{(4)}(p-p')=\int\frac{d^4y}{(2 \pi)^4}e^{iy\cdot (p-p')},
\end{equation}
where $p,p',y$ are Cartesian coordinate vectors living in 4-dimensional Euclidean space and "$\cdot$" denotes the dot product.
Now, consider switching to 4-dimensional spherical coordinates,
\begin{align}
y_1&=y\cos(\phi_1)\\
y_2&=y\sin(\phi_1)\cos(\phi_2)\\
y_3&=y\sin(\phi_1)\sin(\phi_2)\cos(\phi_3)\\
y_4&=y\sin(\phi_1)\sin(\phi_2)\sin(\phi_3),\\
\end{align}
where $\phi_1,\phi_2$ range over $[0,\pi]$ while $\phi_3$ ranges over $[0,2\pi]$.
For these coordinates, we then have (choosing to orient our axes such that $(p-p')$ lies along the $y_4$-axis)
\begin{align}
\int\frac{d^4y}{(2 \pi)^4}e^{iy\cdot (p-p')} &= \frac{\pi}{(2\pi)^4}\int^\infty_0 dy y^3 \int^{2\pi}_0 d\phi_3 \exp{[iy|p-p'|\cos{\phi_3}]} \\
&=\frac{2\pi^2}{(2\pi)^4}\int^\infty_0 dy y^3J_0[y|p-p'|],
\end{align}
where $J_0$ is a Bessel function of the first kind.
If my manipulations are correct, then this implies the 4-dimensional delta function can be represented as a 1-dimensional integral of the Bessel function,
\begin{equation}
\delta^{(4)}(p-p')=\frac{1}{2(2\pi)^2}\int_0^\infty dy y^3 J_0[y|p-p'|]\tag{*}\label{*}
\end{equation}
Questions: 1) How might one go about verifying the relation \eqref{*}? 2) Suppose one introduced a cutoff on the $y$-integral—i.e., replaced $\int_0^\infty \to \int_0^{y_{\text{cut}}}$, then can one give a precise statement about how "well" this finite integral expression would approximate the dirac delta for a given $y_\text{cut}$?
Edit: (in response to Paul Garrett's answer) Consider the "rectangular function",
$$
\eta_\epsilon(x)\equiv
\begin{cases}
\epsilon^{-1} & -\epsilon/2 < x < \epsilon/2\\
0 & \text{otherwise}
\end{cases},
$$
which represents the delta function via
\begin{align}
\delta[f] &= \lim_{\epsilon\to 0} \int_\mathbb{R} \eta_{\epsilon} (x) f(x) dx \tag{**} \label{**} \\
&= \lim_{\epsilon\to 0} \epsilon^{-1} \int^{\epsilon/2}_{-\epsilon/2} f(x) dx \\
&= \lim_{\epsilon\to 0} \epsilon^{-1} \int^{\epsilon/2}_{-\epsilon/2} \sum^\infty_{n=0} \frac{x^n}{n!} f^{(n)}(0) dx \\
&= \lim_{\epsilon\to 0} \sum^\infty_{n=0} \frac{(\epsilon/2)^{2n}}{(2n+1)!} f^{(2n)}(0),\\
\end{align}
where I have considered only functions $f(x)$ which equal their Maclaurin series on the interval $x:[-\epsilon/2,\epsilon/2]$.
The analogue to a "cutoff" for the rectangular representation \eqref{**} is taking $\epsilon$ to be finite–for which we find
\begin{equation}
\delta [f] – f(0) = \sum^\infty_{n=1} \frac{(\epsilon/2)^{2n}}{(2n+1)!} f^{(2n)}(0).
\end{equation}
Clearly, under the assumptions made here, the difference $\delta[f]-f(0)$ will go to zero in the limit that $\epsilon\to 0$, but it will do so in a way which depends on the behavior of $f(x)$ in a neighborhood of the origin–in particular, the values of its derivative at $x=0$.
Best Answer
In response to the questioner's comment: one standard viewpoint is that distributions (or tempered distributions, similarly) are "weak-dual-topology" limits of test functions (or of Schwartz functions), where the weak-dual topology (also called "weak-star topology) is given by semi-norms $\mu_\varphi(u)=|u(\varphi)|$ for test functions (or Schwartz) $\varphi$, at first for other test functions (or Schwartz $\varphi$) and then on the (quasi-) completion, by continuity. Thus, in this viewpoint, the integral representation of $\delta$ is not a literal integral, but is an extension-by-continuity of literal integrals giving approximations to $\delta$. I think this viewpoint, of having integrals not being literal, but extensions-by-continuity, may be most consonant with a common physics-y viewpoint... and does (to my mind) convey a more natural sense of "generalized functions/distributions", not as "functionals", but as an extended notion of function.
Of course, the second viewpoint, perhaps more standard in mathematics per se, is that distributions are continuous linear functionals, and that the integrals written about don't make sense... again, at least not as literal integrals, but perhaps more... haughtily said? That is, the same expressions are obtained, but if/because one is insisting on literal integrals, implicitly disallowing writing integral signs for extensions-by-continuity, some new symbols must be made up. (As you can tell, I am not deeply moved by this viewpoint.)
Either way, the underlying mathematics and mechanisms are the same. That is, e.g., for each Schwartz function $\varphi$, the integrals against the truncated Fourier inversion expression for $\delta$ go to $\varphi(0)$ as the cut-off goes to infinity.
The non-uniformity issue/feature is that there is (provably) no estimate for this limit-approach uniform in $\varphi$. It depends, as it turns out, on $\varphi$.
Is this to the point?
EDIT: I forgot to make explicit the point that, whether or not the integrals are pointwise/numerically valued, literal-or-not, etc., the operator/functional indicated by the integral (as extended by continuity, say) has certain expected properties, so that the usual manipulations which would be correct for numerical/literal integrals are also correct for the non-literal (extended-by-continuity) maps written as integrals... because they extend those integrals, and have the same properties... apart, perhaps, from non-critical numerical/pointwise evaluability.
This is why many physics-y heuristics succeed, after all: what is written is not what is literally meant (from the most exclusionary viewpoint), and it is "still ok" because things extend by continuity, with the same correct transformation properties of whatever-it-is operator denoted by "integral" (and extended from a literal integral).
Yes, on good days, we have our cake and eat it, too. This is (to my mind) the marvelous-ness of mathematics, both internally and in applications.
True, on a bad day bad things happen which will be inexplicable to the (naive?) optimist. Yet, crazily enough, even those bad things can often be repurposed to interesting ends. That's the second miracle in this business, I think.
Edit-edit: following up on the edit to the original question, yes, in that example, for $f$ represented by its Taylor-Maclaurin series in a fixed interval $[-\varepsilon_o,\varepsilon_o]$, as $\varepsilon\to 0$ we can estimate the rate of approach to $f(0)$ uniformly.
But if we look at the aggregate of the classes of such functions over all $\varepsilon_o>0$, this uniform-ness is lost.
Also, although I'd not advocate making this one's first worry, there is an old theorem of E. Borel that says that an arbitrary sequence of real (or complex) numbers can be the Taylor-Maclaurin coefficients at $0$ of a smooth function (obviously, then, not necessarily represented by the infinite series, but the finite-series-with-error-term still works...). This shows, strikingly, that there's no way to hope for true uniformity.
Also, as a sort of paraphrase of the Borel resuls: in the classification of distributions supported at $0$ (essentially proven by Taylor-Maclaurin expansions), we find that they are only the finite linear combinations of $\delta$ and its derivatives. Thus, there is no sequence $c_n$ such that $\sum_{n\ge 1} c_n\cdot \delta^{(n)}$ converges as a distribution.