I want a proof or a reference for the identity $$ \int_0^\infty \frac{s^{n-1}}{\Gamma(n)} p_\beta(s,x)\,ds =\frac{x^{n\beta-1}}{\Gamma(\beta n)},\quad x>0, \,n\in\mathbb N, $$ where $x\mapsto p_\beta(s,x)$, $x>0$ is the density of the $\beta$-stable Lévy subordinator at time $s>0$, $\beta\in(0,1)$.
A possible proof is to take the Laplace transforms of both sides above, which are the same, using $\int_0^\infty e^{-\lambda x}p_\beta(s,x)\,dx=e^{-\lambda^\beta s}$, which also proves that the left hand side is finite almost everywhere. If one can prove that the left hand side is also continuous, then uniqueness of Laplace transforms proves the identity.