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

**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$