Let $X,Y\stackrel{\text{i.i.d.}}{\sim}\mathcal{N}(0,I_d)$, where $I_d$ is the $d$-dimensional identity matrix.
What is the distribution of $\langle X,Y\rangle=X^TY$?
Approach 1:
So far I know that for any $i\in\{1,…,d\}$ the MGF of $X_iY_i$ is:
$$
M_{X_iY_i}(t)=\frac{1}{\sqrt{1-t^2}}.
\qquad\left(=\arcsin'(t)\right)
$$
The derivation of this can be found here. I've verified it – it's correct.
Now, using the product-rule the MGF of $\langle X,Y\rangle$ is:
$$
\begin{align*}
M_{\langle X,Y\rangle}(t)
&=\left(\frac{1}{\sqrt{1-t^2}}\right)^d
\\
&=(1-t^2)^{-\frac{d}{2}}
\\
&=(1-t^2)^{-\frac{d}{2}}
\\
&=(1-t)^{-\frac{d}{2}}(1+t)^{-\frac{d}{2}}
\end{align*}
$$
Now I don't know how to continue from here. Maybe this product might as the product of two typical MGFs. In that case we could express it as the sum of two typical random variables… I don't know how to continue from here.
Approach 2:
Let $\theta$ denote the angle between the two random vectors $X$, $Y$. Then
$$
\begin{align*}
\langle X,Y\rangle
&=
\|X\|\|Y\|\cos(\theta)
\\
&=
\sqrt{\langle X,X \rangle}\sqrt{\langle Y,Y\rangle}\cos(\theta)
\end{align*}
$$
Now, $U:=\langle X,X\rangle\sim \chi^2(d)$ and $V:=\langle Y,Y\rangle\sim\chi^2(d)$ as one can easily verify.
$$
\begin{align*}
&=
\sqrt{U\cdot V}\cos(\theta)
\end{align*}
$$
Now, $\cos(\theta)$ is distributed as the dot-product divided by the norms of the two random vectors $X,Y$. So $\cos(\theta)$ is distributed as the dot-product of the corresponding random unit vectors of $X,Y$. By rotational invariance of the dot-product, we may rotate both random unit vectors (corresponding to $X,Y$), such that one of the unit vectors has only the first coefficient being one. Hence, $\cos(\theta)$ is also distributed as the first-coefficient of a random unit vector. Since the angle between $X$ and $Y$ doesn't affect the outcome of $U$ and $V$ (since they're just the norms of the random vectors; so $\cos(\theta)$ is independent of $U$ and also $V$), we can model the randomness of $\cos(\theta)$ independently (as we've said through the first coefficient of a random unit vector).
We model this random unit vector through sampling a random multivariate gaussian vector $Z\sim\mathcal{N}(0,I_d)$ that we then normalize. Then we just take its first coefficient. Now, w.l.o.g. we may also use Z:=X. So,
$$
W:=
\cos(\theta)\sim\frac{X_1}{\sqrt{\sum_{i=1}^d X_i^2}}.
$$
And finally we have
$$
\begin{align*}
\langle X, Y\rangle
&=\sqrt{U\cdot V}\cdot W
\\
&=\sqrt{\langle X,X\rangle}\sqrt{V}\frac{X_1}{\sqrt{\langle X,X \rangle}}
\\
&=\sqrt{V}X_1
\end{align*}
$$
Now, we have the product of the square root of a chi-squared distributed variable and a standard normal random variable. Do you know if that gives any familiar distribution?
Right now I'm building the product distribution – let's see what comes out.
Approach 3:
It's known that the dot product of two random unit vectors in $\mathbb{R^d}$ follows a $Beta((d-1)/2,(d-1)/2)$ distribution (see this post). So, let
$$
Z\sim Beta((d-1)/2,(d-1)/2)
$$
Then we have that
$$
\langle X,Y \rangle = Z*\|X\|\|Y\|
$$
Where $\|X\|,\|Y\|\stackrel{\text{i.i.d.}}{\sim}\chi^2(d)$. So we have that the dot-product is a product of a Beta-, and two Chi^2-distributed variables.
That might help us to continue from here by building the product distribution.
Sidenote: In the limit $d\to\infty$ we have that it approaches the normal distribution (see here, or Feller II (1971, p. 516, Theorem 2).
Best Answer
The characteristic function of $\langle X,Y\rangle$ when $X$ and $Y$ are independent standard Gaussians in $\mathbb R^d$ is easily seen to be $\phi_d(t)=(1+t^2)^{-d/2}.$ What is sought is a formula for the corresponding density functions $f_d$.
Here is one way of seeing why $\phi_d$ is given by this formula. The inner product $\langle X,Y\rangle$ is the sum of $d$ independent terms, each distributed like $UV$ where $U,V\sim N(0,1)$ are iid. So the desired characteristic function is $E[\exp(itUV)] = E[ E[\exp(itUV)|V)]]$. Here the inner expectation is conditional on $V$, and it evaluates to $\exp(-t^2 V^2/2)$ by the fact that the characteristic function of a univariate standard normal is $t\mapsto\exp(-t^2/2)$. To finish the job, integrate $$\int_{-\infty}^\infty \frac {\exp(-v^2/2)}{\sqrt{2\pi}} \exp( -t^2v^2/2) \,dv,$$ to yield the formula for $\phi_1$. (If you know the norming constant for the $N(0,\sigma^2)$ density, you read off the answer from the case $\sigma^2=1+t^2$.) To finish the job, $\phi_d = \phi_1^d$.
The case $d=1$ has a seemingly simple answer. Starting with the identity $$K_0(x) = \int_0^\infty \frac{\cos(tx)}{\sqrt{t^2+1}}\,dt,$$ found as formula 9.6.21 in Abramowitz and Stegun's Handbook of Mathematical Functions one sees that $$f_1(x) = \frac 1 \pi K_0(|x|),$$ where $K_0$ is a "modified Bessel function of the second kind". (See the DLMF article for useful relevant integrals.)
The $d=2$ case is easier. Here the characteristic function is $\phi_2(t)=1/(1+t^2)$, which corresponds to the density function $f_2(x)=\frac 1 2 \exp(-|x|).$ By repeated convolution one sees that $f_d$ for even $d$ takes the form $f_d(x)=P_d(|x|)\exp(-|x|)$ where $P_d$ is a polynomial.