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.



Composition of Riesz potentials
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$.


