Why do I have to distinguish between $p$ and $p_1$, $p_2$, $p_3$? That is, why can't I assume that all three differentials - $dp_1$, $dp_2$, $dp_3$ - are simply $dp$?

Because $p$ is a point in 3D space, and its components are three totally different real variables $(p_1,p_2,p_3)$. If we assumed $p_1=p_2=p_3$, we'd essentially just be integrating along the diagonal! Does that answer your question?

If you're more familiar with distinguishing your variables by letter, your comment would be like saying "why can't we just use $dx$ in $\iint dx\ dy\ e^{-(x^2+y^2)}$?" They're essentially the same situation, but here we've sort of labeled $x$ as $p_1$, $y$ as $p_2$, etc.

The formula $$Integrate[Integrate[Integrate[1, p], p], p]$$gives $\frac{p^3}{6}$ which is an analytical answer that seems to agree with the rest of the calculations on the page from which I'm quoting this formula. I'm wondering: is this answer is compatible with yours?

If you provide a source, I can tell you why that would be either super nonstandard notation or incorrect! :)

It's not compatible; the integral $\int_\Omega d^3p$, where $\Omega$ is the region we're integrating over, should be the volume of the region $\Omega$. So, if it's a rectangular prism with dimensions $P_1,P_2,P_3$, the integral should be $\int_0^{P_1}\int_0^{P_2}\int_0^{P_3}dp_3\ dp_2\ dp_1$, which is $P_1 P_2 P_3$ (no division by $6$).

The source is Dodelson, Modern Cosmology, pg 61.

Let us continue this discussion in chat.

Actually, there's no latex rendering in chat, so I think it's best to continue it here, unfortunately.

Ok, so, check out equation 3.2. The relation $E^2 = p^2 + m^2$ actually takes $p$ to be a 3-dimensional vector, and $p^2$ is $p\cdot p$. The constraining equation must integrate over the entirety of the space $p$ is in, not over a single real-valued variable thrice!

Now, it's unfortunate I named the components of the vector $p$ $p_1, p_2, p_3$, because Dodelson uses $p_i$ in eq 3.8 to refer to entirely separate 3D vectors. However, each of them are 3D vectors, and thus to integrate over the entire volume of space, we use $d^3$ to refer to that infinitesimal cube for each of them.

Wonderful. We agree up to that point. Where does that leave us regarding a MMa version of the formula?

Integrate[Integrate[Integrate[f[p], p], p], p] is, mathematically, $\int_0^{p'''} dp''\int_0^{p''}dp'\int_0^{p'} dp\ f(p)$, where i've used primes to distinguish different integrand variables. This is different than integrating over three separate variables, i.e. over three components of a vector, since relative to each other, the three components of a vector are constants, and don't get "picked up" by each successive integration.

In other words, Integrate[Integrate[Integrate[f[p], p], p], p] treats p as a single real variable being iteratively integrated over three times, not as three independent variables, which is what we need in order for $p$ to be a vector.

I think we agree. Thank you.

Ah, ok, great! Glad I could help. :)

OK. I think I see the issue. Is there a way to analyze this spherically (3-sphere) such that, instead of p1, p2, p3, we have a single value for p?

The 3-sphere is 3-dimensional, so you will still need three parameters. However, if you're talking about the "mass shell" defined by $E^2=p^2+m^2$, then those three parameters will range over specific values of $E, p_1, p_2, p_3$, just like a parameterization of the 2-sphere ranges over all three of $x,y,z$. The 4D integral with the dirac delta function actually makes your life easier instead of harder here! In general, you could also try to integrate over each 2-sphere slice of the 3-sphere at constant $E$, and then integrate over $E$. Is that maybe the kind of thing you're looking for?

Someone who has more familiarity with the Dodelson derivation answered a related question about the physics behind this integral. physics.stackexchange.com/questions/621902/… He evaluated the integral with spherical coordinates. I'm still analyzing the formula, but it seems to have the form that the author of the book uses. There's a single 'p' (momentum) term that eventually is converted into mass. That's why I'm trying to find an integral that results in a single 'p' term.

I see. That's a pretty big, albeit common, abuse of notation used there: when used in the context of $\int d^3p$, $p$ is a 3D vector, but when used in the context of $\int_0^\infty p^2\ dp$, it refers to the norm of that vector. The other two degrees of freedom are converted into $\theta$ and $\phi$, but we still have 3 degrees of freedom describing a point in momentum-space. In general the form of integration in spherical coordinates given there does indeed work for integrating over all of momentum space in general:

Integrate[Exp[-E1[{p0 Sin[\[Theta]] Cos[\[Phi]], p0 Sin[\[Theta]] Sin[\[Phi], p0 Cos[\[Theta]]}]] p0^2 Sin[\[Theta]], {p0, 0, Infinity}, {\[Theta], 0, Pi}, {\[Phi], 0, 2 Pi}]

Please have a look at the formula you just pasted above (two comments above). I can't get it to compile.

Missing a bracket: Integrate[Exp[-E1[{p0 Sin[\[Theta]] Cos[\[Phi]], p0 Sin[\[Theta]] Sin[\[Phi]], p0 Cos[\[Theta]]}]] p0^2 Sin[\[Theta]], {p0, 0, Infinity}, {\[Theta], 0, Pi}, {\[Phi], 0, 2 Pi}]