Joint density function of $(X^\top u,X^\top v)$, where $u$ and $v$ are unit-vectors in $R^d$ and $X$ is uniform on the unit sphere in $R^d$. Let $u$ and $v$ be fixed vectors on the unit-sphere $S_{d-1}$ in $\mathbb R^d$ and set $t:=u^\top v$. Let $X$ be uniformly distributed on $S_{d-1}$.

Question. What is the joint probability density function of $(X^\top u,X^\top v)$ ?

Observation

*

*By rotational symmetry, $(X^\top u,X^\top v)$ has the same distribution as $(X_1,Z)$, where $Z:= tX_1 + \sqrt{1-t^2}X_2$. Let $p_1$ (resp. $F_1$) be the marginal pdf (resp. CDF) of $X_1$. Then the joint distribution of $(X_1,Z)$ is
$$
p_{X_1,Z}(x,y) = \frac{\partial F_1(x)}{\partial x}\frac{\partial F_{12}(\frac{y-tx}{\sqrt{1-t^2}}|x)}{\partial y} = \frac{1}{\sqrt{1-t^2}}p_1(x)p_{12}(\frac{y-tx}{\sqrt{1-t^2}}|x),
\tag{*}
$$
where $F_{12}(y|x)$ is the cdf $X_2$ conditioned on $X_1$.

*For example, in the case where $X \sim N(0,I_d)$, instead of uniform on unit-sphere, we have $p_1(x) = \dfrac{1}{\sqrt{2\pi}}\exp(-x^2/2)$ and $p_{12}(y|x) = p_1(y)$ for all $x,y \in \mathbb R$ (by independence of $X_1$ and $X_2$). Thus, we deduce from (*) that
$$
\begin{split}
p_{X_1,Z}(x,y) &= \dfrac{1}{\sqrt{2\pi}}\exp(-x^2/2)\dfrac{1}{\sqrt{2\pi}\sqrt{1-t^2}}\exp(-\frac{(y-tx)^2}{2(1-t^2})\\
&= \dfrac{1}{2\pi\sqrt{1-t^2}}\exp(-\dfrac{1}{1-t^2}(x^2-2txy+y^2)).
\end{split}
$$
 A: It turns out that the following formal representation of the joint density of these two variables is useful for finding an explicit formula for it:
$$f_{X^T u, X^T v}(\alpha, \beta)=\frac{1}{r^{d-1}\mu(S^{d-1})}\int_{\mathbb{R}^d}d^dx~\delta\left(x^Tu-\alpha\right)\delta\left(x^Tv-\beta \right)\delta\left(||x||-r\right)$$
where we assume here that the random vector $X$ is uniformly distributed on a sphere of radius $r$. We also let the vectors $u,v$ be on this sphere as well. Note here that this integral is non-zero only when $|\alpha|,|\beta|\leq r^2$. We consider the integral in this range of variables from now on.
In $d\geq 3$ where the two random variables are not dependent for most choices of $u,v$ and the sample space is non-trivial, we use the rotational invariance to make the delta function integrations easy to do. We change variables $x\to R x$ with the rotation matrix chosen such that $R^Tu=re_1$ and $R^Tv=R^T[(u\cdot v)u+v_{\perp}]=te_1+\sqrt{r^2-t^2}e_2$, where in the second line we decomposed the second vector into a projection on the first + it's complement and we  used one more degree of freedom of the rotation to take the orthogonal complement to   another vector of the standard basis of $\mathbb{R}^d$. Now the integral reads
$$f(\alpha, \beta)=\frac{1}{r^{d}\mu(S^{d-1})}\int_{\mathbb{R}^d}d^dx~\delta\left(x_1-\alpha/r\right)\delta\left(tx_1+\sqrt{1-t^2}x_2-\beta \right)\delta\left(||x||-r\right)$$
Performing the integration in the first variable is trivial and we get
$$f(\alpha, \beta)=\frac{1}{r^{d+1}\mu(S^{d-1})}\int dx_2...dx_d~\delta\left(t\alpha/r+\sqrt{r^2-t^2}x_2-\beta \right)\left[2r\delta\left(x^2_2+...x_d^2-(r^2-\frac{\alpha^2}{r^2})\right)\right]$$
The integration over the second one is not much more difficult, but it adds the constraint $|\beta-t\alpha/r|\leq r\sqrt{r^2-t^2}$ and the integral now is
$$f(\alpha,\beta)=\frac{1}{r^{d+1}\sqrt{1-(t/r)^2}\mu(S^{d-1})}\int dx_3...dx_d~\left[2r\delta\left(x_3^2+...x_d^2-R^{2}_{\text{eff}}(\alpha,\beta)\right)\right]$$
where we defined the effective radius  $R^2_{\text{eff}}(\alpha,\beta):=r^2-\frac{\alpha^2}{r^2}-\frac{(\beta-t\alpha/r)^2}{r^2-t^2}$. Note here that the constraint $R^2_{\text{eff}}\geq 0$ contains all the constraints collected in the evaluation process and therefore it is the only one that needs to be considered. The remaining integral is proportional to the area of a $(d-3)$-dimensional sphere of radius $R_{\text{eff}}$ and thus
$$\int dx_3...dx_d~\left[2r\delta\left(x_3^2+...x_d^2-R^{2}_{\text{eff}}\right)\right]=\frac{r}{R_{\text{eff}}}\mu(S^{d-3})R_{\text{eff}}^{d-3}$$
and finally we can write the joint distribution of the two variables in the form
$$f(\alpha,\beta)=\frac{d-2}{2\pi r^4\sqrt{1-(t/r)^2}}\left(\frac{R_{\text{eff}}(\alpha,\beta)}{r}\right)^{d-4}$$
whenever $R^2_{\text{eff}}\geq 0 $ (which defines the interior of a rotated and translated ellipse) and zero everywhere else. Notably, in $d=4$ the disribution is constant within the ellipse.
