Composition of Riesz potentials For $0<\alpha<n$ and $n\geq 2$ we define the Riesz potential by
$$
(I_\alpha f)(x) = \frac{1}{\gamma(\alpha)}
\int_{\mathbb{R}^n} \frac{f(y)}{|x-y|^{n-\alpha}}\, dy\, ,
\quad
\text{where}
\quad
\gamma(\alpha)=
\frac{\pi^{\frac{n}{2}}\, 2^\alpha\,\Gamma\left(\frac{\alpha}{2}\right)}
{\Gamma\left(\frac{n-\alpha}{2}\right)}\, .
$$
The following result is well known and regarded as nearly obvious:

Theorem.
If $\alpha,\beta>0$, $\alpha+\beta<n$, then $I_\alpha(I_\beta\varphi)=I_{\alpha+\beta}\varphi$ for
  $\varphi\in\mathscr S_n$.

Formally
\begin{eqnarray*}
(I_\alpha I_\beta\varphi)^\wedge(\xi)
& = &
\left((-\Delta)^{\alpha/2}(-\Delta)^{\beta/2}\varphi)\right)^\wedge(\xi)
= 
(4\pi^2|\xi|^2)^{\alpha/2}(4\pi^2|\xi|^2)^{\beta/2}\hat{\varphi}(\xi)\\
& = &
(4\pi^2|\xi|^2)^{\frac{\alpha+\beta}{2}}\hat{\varphi}(\xi)
=
\Big((-\Delta)^{\frac{\alpha+\beta}{2}}\varphi\Big)^\wedge(\xi)
=(I_{\alpha+\beta}\varphi)^\wedge(\xi)
\end{eqnarray*}
and the result follows by taking the inverse Fourier transform.
Having this argument in mind, Stein on p.118 in his Singular integrals and differentiable properties of functions writes `deduction of this formula offer no difficulties'. But is it really a rigorous proof? The problem is that $I_\alpha\varphi$ is defined as a convolution with the distribution, $I_\alpha\varphi=u_\alpha*\varphi$, where $u_\alpha=\gamma(\alpha)^{-1}|x|^{\alpha-n}$. We know that in that case we can use the formula $(u_\alpha*\varphi)^\wedge=\hat{\varphi}\hat{u}_\alpha$. Unfortunately, for most of the functions $\varphi\in\mathscr S_n$, $u_\alpha*\varphi=I_\alpha\varphi\not\in\mathscr S_n$ so the formula 
$$
(I_\alpha(I_\beta\varphi))^\wedge=(u_\beta*(u_\alpha*\varphi))^\wedge=
(u_\alpha*\varphi)^\wedge \hat{u}_\beta
$$
is not properly justified.
My question is: Do you know a reference to a rigorous and detailed proof of the above fact?
I will provide an answer to my question by showing how I prove it rigorously, but, perhaps my argument is overly complicated and I am simply not able to see obvious things. Also, I could not find a rigorous and elementary proof anywhere.
 A: This is an extended version of my comment.

Let $g_t$ be the Gauss–Weierstrass kernel,
$$
 g_t(x) = \frac{1}{(4 \pi t)^{n/2}} \, e^{-|x|^2/(4t)} .
$$
If $\alpha \in (0, n)$, we have
$$
\begin{aligned}
 \frac{1}{\Gamma(\tfrac{\alpha}{2})} \int_0^\infty g_t(x) t^{\alpha/2 - 1} dt & = \frac{1}{2^n \pi^{n/2} \Gamma(\tfrac{\alpha}{2})} \int_0^\infty t^{-1 - (n - \alpha)/2} e^{-|x|^2 / (4 t)} dt \\
 & = \frac{2^\alpha}{\pi^{n/2} \Gamma(\tfrac{\alpha}{2}) \, |x|^{n - \alpha}} \int_0^\infty s^{-1 + (n - \alpha)/2} e^{-s} ds \\
 & = \frac{\Gamma(\tfrac{n - \alpha}{2})}{2^\alpha \pi^{n/2} \Gamma(\tfrac{\alpha}{2}) \, |x|^{n - \alpha}} = \frac{1}{\gamma(\alpha)} \, \frac{1}{|x|^{n - \alpha}} \, ;
\end{aligned}
\tag{1}
$$
here we substituted $s = |x|^2 / (4 t)$ and used the gamma integral. Therefore, if $\alpha, \beta, \alpha + \beta \in (0, n)$, we have
\begin{align*}
\hspace{5em} & \hspace{-5em} \int_{\mathbb{R}^n} \frac{1}{\gamma(\alpha)} \, \frac{1}{|y|^{n - \alpha}} \times \frac{1}{\gamma(\beta)} \, \frac{1}{|x - y|^{n - \alpha}} \, dy \\
 & = \frac{1}{\Gamma(\tfrac{\alpha}{2}) \Gamma(\tfrac{\beta}{2})} \int_0^\infty \int_0^\infty g_{t + s}(x) t^{\alpha/2 - 1} s^{\beta/2 - 1} dt ds \\
 & = \frac{1}{\Gamma(\tfrac{\alpha}{2}) \Gamma(\tfrac{\beta}{2})} \int_0^\infty \int_0^u g_u(x) v^{\alpha/2 - 1} (u - v)^{\beta/2 - 1} dv du \\
 & = \frac{1}{\Gamma(\tfrac{\alpha}{2}) \Gamma(\tfrac{\beta}{2})} \int_0^\infty \int_0^1 g_u(x) u^{(\alpha + \beta)/2 - 1} r^{\alpha/2 - 1} (1 - r)^{\beta/2 - 1} dr du \\
 & = \frac{1}{\Gamma(\tfrac{\alpha + \beta}{2})} \int_0^\infty g_u(x) u^{(\alpha + \beta)/2 - 1} du = \frac{1}{\gamma(\alpha + \beta)} \, \frac{1}{|x|^{n - \alpha - \beta}} \, .
\end{align*}
Here we used: Fubini, (1) and $g_t * g_s = g_{t + s}$; then substitutions $s = u - t$ and $v = t$; then $v = u r$; then beta integral; finally again (1). Again by Fubini, we conclude that
$$
 I_\alpha(I_\beta f) = I_{\alpha + \beta} f
$$
for all nonnegative functions $f$.
A: 
Theorem.
If $\alpha,\beta>0$, $\alpha+\beta<n$, then $I_\alpha(I_\beta\varphi)=I_{\alpha+\beta}\varphi$ for
  $\varphi\in\mathscr S_n$.

In the proof we will need the following lemma.

Lemma.
If $\alpha,\beta>0$, $\alpha+\beta<n$, then there is a  constant $C_0=C_0(\alpha,\beta,n)$ such that $$
 \int_{\mathbb{R}^n}\frac{dy}{|x-y|^{n-\alpha}|y|^{n-\beta}}=
 \frac{C_0}{|x|^{n-(\alpha+\beta)}}\, . $$

Proof. First you show that the integral is finite for every $x\neq 0$. By rotational symmetry, the integral on the left hand side depends on $|x|$ only. Denoting its value by $f(|x|)$ a simple change of variables (by scaling) show that 
$f(|x|)=|x|^{\alpha+\beta-n}f(1)$ and the result follows.
$\Box$
Proof of the theorem.
The lemma and the Fubini theorem easily implies that
$$
I_\alpha(I_\beta\varphi)(x)=
\frac{C_0\gamma(\alpha+\beta)}{\gamma(\alpha)\gamma(\beta)}
I_{\alpha+\beta}\varphi(x).
$$
The only problem is to show that the constant is actually equal $1$.
To prove this it suffices to verify that 
$I_\alpha(I_\beta\varphi)= I_{\alpha+\beta}\varphi$
for just one non-zero function $\varphi$.
To this end
let $\varphi\in\mathscr S_n$ be such that $\hat{\varphi}=0$ in a neighborhood of $0$. Then
$$
I_\alpha(I_\beta\varphi)=
I_\alpha\Big(\Big(\underbrace{(4\pi^2|\xi|^2)^{-\beta/2}\hat{\varphi}}_{\in\mathscr S_n}\Big)^\vee\Big) =
\left((4\pi^2|\xi|^2)^{-\alpha/2}(4\pi^2|\xi|^2)^{-\beta/2}\hat{\varphi}\right)^\vee=
I_{\alpha+\beta}\varphi.
$$
The proof is complete. $\Box$
