Solution (source code)

= Solution

To keep the normalization visible, denote the integral defining the printed constant by $A_d$. Substituting $r=|y|^2/(2s)$ gives, for $y\ne0$,
$$
\int_0^\infty p(s,0,y)\,ds=A_d|y|^{2-d},\qquad
A_d=\frac{\Gamma(d/2-1)}{2\pi^{d/2}}.
$$
This is the <Newtonian potential of the Brownian heat kernel>. The singularity at zero is locally integrable in dimension $d$, since its radial integral behaves like $\int_0^1r^{2-d}r^{d-1}dr=1/2$. All integrands are nonnegative, so <Tonelli's theorem> permits both orders of integration.

Integrating in $s$ first gives
$$
I=A_d\int_{\mathbb R^d}p(t,x,y)|y|^{2-d}\,dy.
$$
Integrating in $y$ first uses the <Gaussian heat kernel> <convolution>. Completing the square in $|y-x|^2/t+|y|^2/s$ gives
$$
\int_{\mathbb R^d}p(t,x,y)p(s,0,y)\,dy
=(2\pi(t+s))^{-d/2}e^{-|x|^2/(2(t+s))}=p(t+s,0,x).
$$
Consequently $I=\int_t^\infty p(r,0,x)dr$, and the actual identity is
$$
\boxed{\int_{\mathbb R^d}p(t,x,y)|y|^{2-d}\,dy
=A_d^{-1}\int_t^\infty p(r,0,x)\,dr.}
$$
The original PDF uses $A_d$ rather than $A_d^{-1}$ as the multiplier while defining its constant to equal $A_d$. That formula is false. For example $A_3=1/(2\pi)$, so the required multiplier in dimension three is $2\pi$. Equivalently, as $t\downarrow0$ at $x\ne0$, the left side tends to $|x|^{2-d}$, whereas the printed right side tends to $A_d^2|x|^{2-d}$. One can retain the printed form only by defining its multiplier as $c_d=A_d^{-1}$.