Prove the following beta integral identity I am looking at the following exercise from Grafakos' classical fourier analysis. For $0 < \alpha_1, \alpha_2 < n, \alpha_1 + \alpha_2 > n$, I want to prove that $$\int_{\mathbb{R}^n} \frac{1}{\left| x-t \right|^{\alpha_1} \left| y -t \right|^{\alpha_2}} dt = \pi^{\frac{n}{2}} \frac{\Gamma \left( \frac{n-\alpha_1}{2} \right) \Gamma \left( \frac{n- \alpha_2}{2} \right) \Gamma \left( \frac{\alpha_1 + \alpha_2 -n}{2} \right)}{\Gamma \left( \frac{\alpha_1}{2} \right) \Gamma \left( \frac{\alpha_2}{2} \right) \Gamma \left( n - \frac{\alpha_1 + \alpha_2}{2} \right)}\left| x -y \right|^{n- \alpha_1 - \alpha_2}.$$ 
I've been able to do most of the exercises in Grafakos, but don't even know where to start on this one. 
 A: With the convention
$$\widehat f(p)=\int f(x)\exp(-ix\cdot p)dx$$
I assume you know
$$\widehat{|x|^{-\alpha}}(p)=c_{n-\alpha,n}|p|^{\alpha-n}\text{ where }
c_{\alpha,n}=\pi^{n/2}2^\alpha\frac{\Gamma(\tfrac \alpha 2)}{\Gamma(\tfrac{n-\alpha}2)}$$
for $0<\alpha<n$, and that this needs to be understood in a distributional sense because $|x|^{-\alpha}$ is only locally integrable.
Formal computation - ignoring convergence completely
$$
\begin{align*}
&\int \frac{1}{\left| x-t \right|^{\alpha_1} \left| y -t \right|^{\alpha_2}} dt\\
&=
(2\pi)^{-n}\int c_{n-\alpha_1,n}|p|^{\alpha_1-n}c_{n-\alpha_2,n}|p|^{\alpha_2-n}\exp(-i(x-y)\cdot p) dp\\
&=
(2\pi)^{-n}c_{n-\alpha_1,n}c_{n-\alpha_2,n}\int |p|^{\alpha_1+\alpha_2-2n}\exp(-i(x-y)\cdot p) dp\\
&=
(2\pi)^{-n}c_{n-\alpha_1,n}c_{n-\alpha_2,n}c_{\alpha_1+\alpha_2-n,n}|x-y|^{n-\alpha_1+\alpha_2}\\
&=\pi^{n/2}
\frac{\Gamma \left( \frac{n-\alpha_1}{2} \right) \Gamma \left( \frac{n- \alpha_2}{2} \right) \Gamma \left( \frac{\alpha_1 + \alpha_2 -n}{2} \right)}{\Gamma \left( \frac{\alpha_1}{2} \right) \Gamma \left( \frac{\alpha_2}{2} \right) \Gamma \left( n - \frac{\alpha_1 + \alpha_2}{2} \right)}
|x-y|^{n-\alpha_1-\alpha_2}
\end{align*}
$$
as required.
Making this rigorous seems surprisingly difficult. I think the easiest approach is to apply two Schwarz functions - one just to state a useful identity and another one to prove it.

Lemma
Let $\phi$ be a Schwarz function (thought of as an approximate Dirac delta). Then
  $$\int \frac{1}{\left| x-t \right|^{\alpha_1}} \int \frac{\phi(t-s)}{\left| y -s \right|^{\alpha_2}} ds dt
=
(2\pi)^{-n}c_{n-\alpha_1,n}c_{n-\alpha_2,n}\int \frac{\exp(-i(x-y)\cdot p)\widehat{\phi}(p)}{\left|p\right|^{2n-\alpha_1-\alpha_2}} dp
$$
  Proof.
Multiplying the $\int \frac{\phi(t-s)}{\left| y -s \right|^{\alpha_2}} ds$ by a Schwarz function $\psi(t)$ turns it into a Schwarz function - we've already convolved to remove the singularity at $t=0$. This lets us use the Fourier transform of the distribution $\frac{1}{\left| x-t \right|^{\alpha_1}}$, after which we can use the Fourier transform of $\frac{1}{\left| y -s \right|^{\alpha_2}}$, turning multiplication into convolution and vice versa, to get:
  $$\int \frac{1}{\left| x-t \right|^{\alpha_1}}  \int \psi(t)\frac{\phi(t-s)}{\left| y -s \right|^{\alpha_2}} ds dt\\
=
(2\pi)^{-n}c_{n-\alpha_1,n}c_{n-\alpha_2,n}\int \frac{\exp(-ix\cdot p)}{\left|p\right|^{n-\alpha_1}} \int \hat \psi(p-q) \frac{\exp(-iy\cdot q)\hat \phi(q)}{\left|q\right|^{n-\alpha_2}} dq dp
$$
Now take $\psi$ to be a wide Gaussian so $\psi\to 1$ and $\hat\psi\to\delta_0$. The integrals should converge to the correct values. I think the only subtlety here is the convolution for $\hat\phi$ near the $p=q$ singularity, where we can use Young's inequality on a small ball around $p=0$ to show that the $\int dq$ term is bounded in $L_r$ where $r$ is chosen with $n/\alpha_1<r<n/(n-\alpha_2)$. Then $\frac{1}{\left|p\right|^{n-\alpha_1}}$ is bounded in the dual space $L_{r'}$ with $1/r+1/r'=1$, so the integral over a small ball is bounded by Hölder's inequality by an error term that goes to zero as the radius of the ball goes to zero.

Given this lemma we can use the Fourier transform property of $\frac{1}{|x|^{n-\alpha_1-\alpha_2}}$ to get
$$\int \frac{1}{\left| x-t \right|^{\alpha_1}} \int \frac{\phi(t-s)}{\left| y -s \right|^{\alpha_2}} ds dt\\
=
(2\pi)^{-n}c_{n-\alpha_1,n}c_{n-\alpha_2,n}c_{\alpha_1+\alpha_2-n,n}\int \phi(t)|x-y-t|^{n-\alpha_1+\alpha_2}dt$$
Taking $\phi\to\delta_0$ the limits will converge the the values you'd expect, giving the second to last line of the formal computation.
A: @Dap's answer provides a fundamental understanding on the identity. So let me safely ignore the issue of intuition and just supply a brutal force computation.
Before starting our computation, we remark some identities to be utilized:
1. Variant of the gamma function integral. For $\alpha > 0$ and $s > 0$ we have
$$ \frac{1}{s^{\alpha}} = \frac{2}{\Gamma(\frac{\alpha}{2})} \int_{0}^{\infty} \lambda^{\alpha-1}e^{-s^2\lambda^2} \, d\lambda. \tag{1}$$
2. Beta function integral. For $\alpha, \beta > 0$ we have
$$ \int_{0}^{\frac{\pi}{2}} \sin^{\alpha-1}\theta \cos^{\beta-1}\theta \, d\theta = \frac{\Gamma(\frac{\alpha}{2})\Gamma(\frac{\beta}{2})}{\Gamma(\frac{\alpha+\beta}{2})}. \tag{2}$$
Now let us begin our computation. Let $I$ denote the integral in question. Using $\text{(1)}$, we may write
$$ I = \frac{4}{\Gamma(\frac{\alpha_1}{2})\Gamma(\frac{\alpha_2}{2})} \int_{0}^{\infty} d\lambda_1 \int_{0}^{\infty} d\lambda_2 \int_{\mathbb{R}^n} \lambda_1^{\alpha_1-1}\lambda_2^{\alpha_2-1} e^{-\lambda_1^2 |x-t|^2 - \lambda_2^2|y-t|^2}. \tag{*} $$
(Interchanging the order of integration is justified since the integrand is positive; you may invoke Tonelli's theorem.) Now completing the square of the exponent
$$ -\lambda_1^2 |x-t|^2 - \lambda_2^2|y-t|^2
= -(\lambda_1^2 + \lambda_2^2)\left|t - \frac{\lambda_1^2 x+ \lambda_2^2 y}{\lambda_1^2+\lambda_2^2} \right|^2
- \frac{\lambda_1^2 \lambda_2^2}{\lambda_1^2+\lambda_2^2} |x - y|^2 $$
and applying the gaussian integral, $\text{(*)}$ simplifies to
$$ I = \frac{4\pi^{n/2}}{\Gamma(\frac{\alpha_1}{2})\Gamma(\frac{\alpha_2}{2})} \int_{0}^{\infty} d\lambda_1 \int_{0}^{\infty} d\lambda_2 \frac{\lambda_1^{\alpha_1-1}\lambda_2^{\alpha_2-1}}{(\lambda_1^2 + \lambda_2^2)^{n/2}} \exp\left\{ - \frac{\lambda_1^2 \lambda_2^2}{\lambda_1^2+\lambda_2^2} |x - y|^2\right\}. $$
Then applying the polar coordinates change $(\lambda_1, \lambda_2) = r(\cos\theta, \sin\theta)$, this further simplifies to
\begin{align*}
I
&= \frac{4\pi^{n/2}}{\Gamma(\frac{\alpha_1}{2})\Gamma(\frac{\alpha_2}{2})} \int_{0}^{\frac{\pi}{2}} d\theta \left( \int_{0}^{\infty} dr \, r^{\alpha_1+\alpha_2-n-1} e^{- r^2|x-y|^2 \sin^2\theta \cos^2\theta} \right)  \cos^{\alpha_1 -1}\theta \sin^{\alpha_2-1}\theta \\
&\stackrel{(1)}{=} \frac{2\pi^{n/2}}{\Gamma(\frac{\alpha_1}{2})\Gamma(\frac{\alpha_2}{2})} \int_{0}^{\frac{\pi}{2}} d\theta \left( \frac{\Gamma(\frac{\alpha_1+\alpha_2-n}{2})}{|x-y|^{\alpha_1+\alpha_2-n}\sin^{\alpha_1+\alpha_2-n}\theta \cos^{\alpha_1+\alpha_2-n}\theta} \right)  \cos^{\alpha_1 -1}\theta \sin^{\alpha_2-1}\theta\\
&\stackrel{(2)}{=} \pi^{n/2} \frac{\Gamma(\frac{\alpha_1+\alpha_2-n}{2}) \Gamma(\frac{n-\alpha_1}{2})\Gamma(\frac{n-\alpha_2}{2})}{\Gamma(\frac{\alpha_1}{2})\Gamma(\frac{\alpha_2}{2})\Gamma(n-\frac{\alpha_1+\alpha_2}{2})} |x-y|^{n-\alpha_1-\alpha_2}.
\end{align*}
