Distribution and moments of $\frac{X_iX_j}{\sum_{i=1}^n X_i^2}$ when $X_i$'s are i.i.d $N(0,\sigma^2)$ Suppose $X_1,X_2,\ldots,X_n$ are independent $N(0,\sigma^2)$ random variables.
For $i,j\in \{1,2,\ldots,n\}$, consider $$U=\frac{X_iX_j}{\sum_{i=1}^n X_i^2}$$
Provided $n>1$, we know that $U$ has a Beta distribution when $i=j$ :
$$U=\frac{X_i^2/\sigma^2}{\sum_{i=1}^n X_i^2/\sigma^2} \sim \text{Beta}\left(\frac12,\frac{n-1}2\right) \quad,\,i=1,2,\ldots,n$$

What can we say regarding the distribution of $U$ when $i\ne j$? What are the moments of $U$ in this case?

For $n=2$, if we transform $(X_1,X_2)$ to polar coordinates $(R,\Theta)$, then
$$U=\frac{X_1X_2}{X_1^2+X_2^2}=\frac{R^2\cos\Theta\sin\Theta}{R^2}=\frac{\sin(2\Theta)}{2}$$
Since $\Theta$ is uniformly distributed on $(0,2\pi)$, it seems $\sin(2\Theta)$ has an $\text{Arcsine}(-1,1)$ distribution with pdf
$$f(x)=\frac1{\pi \sqrt{1-x^2}}\mathbf 1_{(-1,1)}(x)$$
So that $U$ has pdf
$$f_U(u)=2 f(2u)=\frac2{\pi\sqrt{1-4u^2}}\mathbf1_{\left(-\frac12,\frac12\right)}(u)$$
If $\boldsymbol X=(X_1,X_2,\ldots,X_n)^T$, we can think of $U$ as the product of $i$th and $j$th components of the vector $\frac{\boldsymbol X}{\lVert \boldsymbol X \rVert}$. And we know that $\frac{\boldsymbol X}{\lVert \boldsymbol X \rVert}$ is uniformly distributed on the surface of a unit sphere. I am not sure if this helps in any way.
 A: Note that we may assume $X_i$'s are standard normal. Then
\begin{align*}
\mu_k
&:= \mathbf{E}\left[ \biggl( \frac{ X_i X_j }{X_1^2 + \cdots + X_n^2} \biggr)^k \right] \\
&= \frac{1}{(2\pi)^{n/2}} \int_{\mathbb{R}^n} \biggl( \frac{ x_i x_j }{x_1^2 + \cdots + x_n^2} \biggr)^k e^{-(x_1^2 + \cdots + x_n^2)/2} \, \mathrm{d}x_1 \cdots \mathrm{d}x_n \\
&= \frac{1}{(2\pi)^{n/2}} \left( \int_{0}^{\infty} r^{n-1}e^{-r^2/2} \, \mathrm{d}r \right)\left( \int_{\mathbb{S}^{n-1}} \omega_i^k \omega_j^k \, \sigma(\mathrm{d}\omega) \right),
\end{align*}
where $\mathbb{S}^{n-1} = \{ \mathrm{x} \in \mathbb{R}^n : \|\mathrm{x}\| = 1\}$ is the unit sphere in $\mathbb{R}^n$ and $\sigma$ is the surface measure on $\mathbb{S}^{n-1}$. Then by mimicking the proof of the beta function identity, such as in this article, it is not hard to prove that
$$ \int_{\mathbb{S}^{n-1}} \omega_1^{\alpha_1} \cdots \omega_n^{\alpha_n} \, \sigma(\mathrm{d}\omega)
= \begin{cases}
0, & \text{if some $\alpha_i$ is odd,} \\[0.25em]
\dfrac{2\Gamma(\beta_1)\cdots\Gamma(\beta_n)}{\Gamma(\beta_1 + \cdots + \beta_n)}, & \text{if all $\alpha_i$ are even and $\beta_i=\frac{1}{2}(\alpha_i + 1)$}.
\end{cases} $$
Using this and simplifying the formula, we get
$$ \mu_k
= \begin{cases}
0, & \text{if $k$ is odd,} \\[0.25em]
\dfrac{(1 \cdot 3 \cdot 5 \cdots (k-1))^2}{\prod_{l=0}^{k-1} (n+2l)}, & \text{if $k$ is even.}
\end{cases} $$
For instance, we get
$$ \mu_2 = \frac{1}{n(n+2)} \qquad\text{and}\qquad
\mu_4 = \frac{3^2}{n(n+2)(n+4)(n+6)}.$$
A: You can use the independence of $\color{blue}{U=\frac{X_iX_j}{\sum_{i=1}^n X_i^2}}$ and $V=\sum\limits_{i=1}^n X_i^2$ to find at least the mean and variance of $U$ for $i\ne j$. To argue somewhat statistically, $U$ and $V$ are independent by Basu's theorem because $U$ is an ancillary statistic for $\sigma^2$ (i.e., its distribution is free of $\sigma^2$), and $V$ is a complete sufficient statistic for $\sigma^2$.
So,
$$0=E\left[X_iX_j\right]=E\left[UV\right]=E\left[U\right]E\left[V\right]\,,$$
As $E\left[V\right]\ne 0$, you must have $$\color{blue}{E\left[U\right]=0}$$
Again,
$$\operatorname{Var}(X_iX_j)=\operatorname{Var}(UV)=E\left[U^2\right]E\left[V^2\right]$$
But $$\operatorname{Var}(X_iX_j)=E\left[X_i^2X_j^2\right]=\sigma^4$$
So, $$E\left[U^2\right]=\frac{\sigma^4}{E\left[V^2\right]}$$
As $\frac{V}{\sigma^2} \sim \chi^2_n$, you have
$$E\left[\left(\frac{V}{\sigma^2}\right)^2\right]=\operatorname{Var}\left(\frac V{\sigma^2}\right)+\left(E\left(\frac{V}{\sigma^2}\right)\right)^2=2n+n^2$$
so that $E\left[V^2\right]=\sigma^4 n(n+2)$ and
$$\color{blue}{\operatorname{Var}(U)=E\left[U^2\right]=\frac1{n(n+2)}}$$
A: Due to the symmetry  $\sigma_1=...=\sigma_n=\sigma$, it is enough to consider $i=1$ and $j=2$.
While my first thought was n-dimensional spherical coordinates also, you can immediately do the integrals under consideration with a trick. We start with
$$I(t)=\int_{-\infty}^\infty \cdots \int_{-\infty}^\infty {\rm d}x_1 \cdots {\rm d}x_n \, \frac{e^{-{t(x_1^2+...+x_n^2)}}}{(2\pi \sigma^2)^{n/2}} \, \frac{x_1^2x_2^2}{(x_1^2+...+x_n^2)^2} $$
$$
I''(t)=\frac{1}{(2\pi \sigma^2)^{n/2}}\int_{-\infty}^\infty \cdots \int_{-\infty}^\infty {\rm d}x_1 \cdots {\rm d}x_n \, e^{-t(x_1^2+...+x_n^2)} \, {x_1^2x_2^2} \\
= \frac{\left(\pi/t\right)^{n/2-1}}{(2\pi \sigma^2)^{n/2}} \cdot \frac{\pi}{4t^3} = \frac{t^{-n/2-2}}{4(2 \sigma^2)^{n/2}} \, .$$
Then by integrating we get $$I(t)= \frac{t^{-n/2}}{(2 \sigma^2)^{n/2}} \cdot \frac{1}{n(n+2)} + c_1t + c_2$$
and since $I(\infty)=0$ it follows $c_1=c_2=0$ and the result is then obtained by setting $t=\frac{1}{2\sigma^2}$.
Evidently this can be generalized to k-th momenta with $k$ being even, while for odd $k$ the integrand is odd and the integral thus zero.
A: Let us denote this variable by $Y$, then
$$
\begin{align}
Y &:=\frac{X_i X_j}{\sum_{s=1}^n X_s^2} \\
&=\left(\frac{X_i X_j}{\sum_{s=1}^n X_s^2}+\frac{1}{2}\right)-\frac{1}{2} \\
&= \frac{1}{2}\left(\frac{\sum_{i=1}^nX_i}{\sqrt{\sum_{i=1}^nX_i^2}}  \right)^2  -\frac{1}{2} \\
&= \frac{1}{2}\left(\frac{n\bar{X}}{\sqrt{n-1}\cdot S}  \right)^2  -\frac{1}{2}  \\
&= \frac{1}{2}\frac{n}{n-1}\left(\frac{\bar{X}}{S / \sqrt{n}}  \right)^2  -\frac{1}{2} \tag{1} \\
\end{align}
$$
with $\bar{X} :=\frac{1}{n}\left(\sum_{i=1}^nX_i  \right)$ and $S:=\frac{1}{n-1}\left(\sum_{i=1}^nX_i^2  \right)$ as defined in the Student's t-distribution.
From $(1)$ and the definition of the t-distribution, we deduce that
$$Y {\buildrel d \over =}\frac{n}{2(n-1)}Z^2 -\frac{1}{2}$$
with $Z$ follows the standard t-distribution of $n-1$ degrees of freedom.
