Suppose that $X$ is distributed uniformly in the scaled $n$sphere $\sqrt{n} \mathbf{S}^{n1} \subset \mathbf{R}^n$. Then apparently the distribution of $(X_1, \dots, X_k)$, the first $k < n$ coordinates of $X$ has density $p(x_1, \dots, x_k)$ with respect to Lebesgue measure in $\mathbf{R}^k$, moreover if $r^2 = x_1^2 + \cdots + x_k^2$, then it is proportional to $$ \left(1  \frac{r^2}{n}\right)^{(nk)/2  1}, \quad \text{if}~0 \leq r^2 \leq n, $$ and otherwise is 0. I tried to compute this using the fact that $(X_1, \dots, X_k) \stackrel{\rm d}{=} \sqrt{n} (g_1, \dots, g_k)/\sqrt{g_1^2 + \cdots + g_n^2}$, when $g_i$ are iid standard normal variables, but it was somewhat unclear to me how to use even this representation to compute the density. Can anyone sketch the details for me?
$\newcommand{\R}{\mathbb{R}} \newcommand{\x}{\mathbf{x}} \newcommand{\X}{\mathbf{X}}$ This is to present a formalization of the answer by Carlo Beenakker, without explicit use of the delta function.
We are going to assume that $\X=(X_1,\dots,X_n)$ is uniformly distributed on the unit sphere $\mathbb S^{n1}$, rather than on $\sqrt n\,\mathbb S^{n1}$.
For each real $t\in(0,1)$, define the measures $\mu_t$ and $\nu_t$ over $\R^n$ by the conditions \begin{equation*} \int f\,d\mu_t=\int_{\R^n}d\x f(\x)1_{1t<\x^2\le1} \end{equation*} and \begin{equation*} \int f\,d\nu_t=\frac{\int f\,d\mu_t}{\int d\mu_t} \end{equation*} for all (say) nonnegative continuous functions $f\colon\R^n\to\R$, where $\cdot$ denotes the Euclidean norm. Then $\nu_t$ is a probability measure converging (as $t\downarrow0$) to the Haar measure $h$ on the unit sphere in $\R^n$, in the sense that (say) \begin{equation*} \int f\,d\nu_t\to\int f\,dh \end{equation*} for all nonnegative continuous functions $f\colon\R^n\to\R$.
Take now any function $f\colon\R^n\to\R$ such that
\begin{equation*}
f(\x)=g(\x_{n1})
\end{equation*}
for some nonnegative continuous function $g\colon\R^{n1}\to\R$ and all $\x\in\R^n$, where $\x_j:=(x_1,\dots,x_j)$ for $\x=(x_1,\dots,x_n)\in\R^n$ and $j=1,\dots,n1$.
Then
\begin{align*}
\int f\,d\mu_t&=\int_{\R^{n1}}d\x_{n1}\,g(\x_{n1})\int_\R du\, 1_{1t<\x_{n1}^2+u^2\le1} \\
&=\int_{\R^{n1}}d\x_{n1}\,g(\x_{n1})
(1+o(1))t\,(1\x_{n1}^2)^{1/2}\,1_{\x_{n1}<1}.
\end{align*}
Also,
\begin{equation*}
\int d\mu_t=\int_{\R^n}d\x\, 1_{1t<\x^2\le1}\propto(1+o(1))t,
\end{equation*}
where $\propto$ means an equality up to a constant factor, depending only on $n,k$.
So,
\begin{equation*}
\int f\,dh=\lim_{t\downarrow0}
\int f\,d\nu_t\propto\int_{\R^{n1}}d\x_{n1}\,g(\x_{n1})
(1\x_{n1}^2)^{1/2}\,1_{\x_{n1}<1}.
\end{equation*}
Thus, the joint pdf of $\X_{n1}=(X_1,\dots,X_{n1})$ is given by
\begin{equation*}
p_{n1}(\x_{n1})\propto(1\x_{n1}^2)^{1/2}\,1_{\x_{n1}<1}.
\end{equation*}
Now successively integrating $p_{n1}(\x_{n1})$ ($n1k$ times) in $x_{n1},\dots,x_{k+1}$ and each time using the formula
\begin{equation*}
\int_0^{b^{1/2}}(bu^2)^p du=c_p b^{p+1/2}
\end{equation*}
for real $b>0$ and $p>1$ with $c_p:=\int_0^1(1u^2)^pdu\in(0,\infty)$ (so that $1/2$ is added to the exponent $p$ after such an integration), we see that the joint pdf of $\X_k=(X_1,\dots,X_k)$ is given by
\begin{equation*}
p_k(\x_k)\propto(1\x_k^2)^{(nk)/21}\,1_{\x_k<1},
\end{equation*}
as desired.

$\begingroup$ This is a wonderfully simple calculation and was illuminating. Thank you for the follow up questions, your persistence, and both answers. $\endgroup$ May 8 '20 at 20:35

3$\begingroup$ @DrewBrady : Thank you for your kind words. However, I feel grateful to Carlo for teaching me something of value, with patience. $\endgroup$ May 8 '20 at 21:03
With $X$ uniformly distributed over the unit $n$sphere, the joint probability distribution of all $n$ elements of $X$ is a Dirac delta function, $$P(X_1,X_2,\ldots X_n)\propto\delta\left(1\sum_{j=1}^n X_j^2\right).\qquad\qquad(1)$$ Now you integrate out elements one by one, to obtain the marginal distribution $P_k$ of $k<n$ elements. The first integration gives $$P_{n1}(X_2,X_3,\ldots X_n)\propto\left(1\sum_{j=2}^n X_j^2\right)^{1/2}\theta\left(1\sum_{j=2}^n X_j^2\right),\qquad(2)$$ with $\theta$ the unit step function. The second integration gives $$P_{n2}(X_3,\ldots X_n)\propto\theta\left(1\sum_{j=3}^n X_j^2\right),$$ the third integration $$P_{n3}(X_4,\ldots X_n)\propto\left(1\sum_{j=2}^n X_j^2\right)^{1/2}\theta\left(1\sum_{j=2}^n X_j^2\right),$$ and so on. Each additional integration increases the power by 1/2, $$P_{k}(X_{nk+1},\ldots X_n)\propto\left(1\sum_{j=nk+1}^n X_j^2\right)^{(nk)/21}\theta\left(1\sum_{j=nk+1}^n X_j^2\right).$$ This is the answer in the OP _{(without rescaling the radius of the $n$sphere, so $r^2/n\mapsto r^2$).}
_{ As requested in the comments, a more detailed exposition of the various steps. }
_{ • First step: the delta function. Denote the surface measure on the unit $n$sphere as $d\Omega$, and $\int d\Omega=A_n$ the surface area. Uniformity of a distribution on the unit $n$sphere means uniformity with measure $d\Omega$. I maintain that the joint probability distribution of the components of the vector ${\mathbf X}=(X_1,X_2,\ldots X_n)$, uniformly distributed on the unit $n$sphere, is given by Eq. (1) with normalization constant $2/A_n$. Let us check this by calculating the expectation value of an arbitrary function $f$ of ${\mathbf X}$. For that purpose I transform to hyperspherical coordinates $r,\phi_1,\phi_2,\ldots\phi_{n1}$, $$\mathbb{E}[f(\mathbf{X})]=\int dX_1\int dX_2\cdots\int dX_n \,f(X_1,X_2,\ldots X_n)P(X_1,X_2,\ldots X_n)$$ $$\qquad=\int_0^\infty r^{n1} dr \int d\Omega\, f(r,\phi_1,\phi_2,\ldots\phi_{n1})\frac{2}{A_n}\delta(1r^2)$$ $$\qquad=\frac{1}{A_n}\int d\Omega\, f(r=1,\phi_1,\phi_2,\ldots\phi_{n1}).$$ In the last step I used that $\int_0^\infty r^{n1}\delta(1r^2)\,dr=1/2$ for $n\geq 2$. }
_{ • Second step: integration of the delta function, to arrive at Eq. (2). From now on I will ignore the normalization constants, these can easily be recovered at the end. Let me abbreviate $\sum_{j=2}^n X_j^2=s_2$. The marginal distribution $P_1(X_2,X_3,\ldots X_n)$ is obtained by definition upon integration of $P(X_1,X_2,X_3,\ldots X_n)$ over $X_1$. I carry out this integration in cartesian coordinates, changing variables to $q=X_1^2$, $$P_1(X_2,X_3,\ldots X_n)\propto \int_{\infty}^\infty dX_1\delta(1s_2X_1^2),$$ $$\qquad=\int_0^\infty\delta(1s_2q)\frac{dq}{\sqrt q}=(1s_2)^{1/2}\theta(1s_2).$$ }
_{ • Third and following steps: The following steps, subsequent integrations of $X_2,X_3,\ldots$ are now immediate consequences of the integral $$\int_0^a(a^2x^2)^p\,dx=c_p a^{1+2p}.$$ }

$\begingroup$ Hm. What does it mean that this distribution "is a delta function"? $\endgroup$ May 7 '20 at 16:24

$\begingroup$ @ Carlo Beenakker  are you sure you understand what a joint distribution is? $\endgroup$– R WMay 7 '20 at 18:51

$\begingroup$ certainly, let's work this out for $n=2$; we then have in polar coordinates $x_1=r\cos\phi$, $x_2=r\sin\phi$ and we would expect the uniform distribution on the unit circle to be $P(\phi)=\text{constant}$. So let's check that the delta function distribution indeed gives this: take $P(x_1,x_2)\propto\delta(1x_1^2x_2^2)$ and integrate to obtain $P(\phi)=\int_0^\infty P(r\cos\phi,r\sin\phi) rdr\propto\int_0^\infty \delta(1r^2)rdr=\text{constant}$. The multivariate Gaussian is a convenient way to generate the random coordinates, but for this calculation the delta function is easier. $\endgroup$ May 7 '20 at 21:02

$\begingroup$ @ Carlo Beenakker  how does the RHS of your formula know that you are dealing with the uniform distribution on the sphere? $\endgroup$– R WMay 7 '20 at 21:19

1$\begingroup$ Carlo, I'd really like to understand your solution, especially because it promises an improvement over the Gaussian approach. Alas, I don't understand almost any part of your answer. Especially hard for me to understand your stepbystep integration. Can you write down in detail the integrals for at least the first two steps and how you take them? $\endgroup$ May 7 '20 at 23:29
Let $G_1,\dots,G_n$ be iid standard normal random variables. Then the random vector \begin{equation*} (Y_1,\dots,Y_n):=\Big(\frac{G_1}{\sqrt{\sum_1^n G_j^2}},\dots, \frac{G_n}{\sqrt{\sum_1^n G_j^2}}\Big) \end{equation*} is uniformly distributed on the unit sphere $\mathbb S^{n1}$. Let \begin{equation*} \begin{aligned} Z_i&:=Y_i=\frac{G_i}{\sqrt{\sum_1^n G_j^2}}&\text{ if }i\le k,\\ Z_i&:=G_i&\text{ if }i> k. \end{aligned} \end{equation*} We want to find the joint pdf of $(Z_1,\dots,Z_k)$, which is the same as the joint pdf of $(X_1,\dots,X_k)/\sqrt n$.
The vector $(Z_1,\dots,Z_n)$ is obtained from $(G_1,\dots,G_n)$ by the transformation given by \begin{equation*} \begin{aligned} z_i&:=\frac{g_i}{\sqrt{\sum_1^n g_j^2}}&\text{ if }i\le k,\\ z_i&:=g_i&\text{ if }i> k. \end{aligned} \end{equation*} The transformation inverse to this is given by \begin{equation*} \begin{aligned} g_i&:=\sqrt{s_2}\frac{z_i}{\sqrt{1s_1}}&\text{ if }i\le k,\\ g_i&:=z_i&\text{ if }i> k, \end{aligned} \tag{1} \end{equation*} where \begin{equation*} s_1:=\sum_1^k z_j^2,\quad s_2:=\sum_{k+1}^n z_j^2. \end{equation*} The Jacobian determinant of the inverse transformation is \begin{equation*} J=\det(cM)=c^k\det M, \end{equation*} where \begin{equation*} c:=s_2^{1/2}(1s_1)^{3/2},\quad M:=(1s_1)I_k+UU^T, \end{equation*} $I_k$ is the $k\times k$ identity matrix, and $U:=[z_1,\dots,z_k]^T$.
Write $U=UQe_1$, where $U=\sqrt{s_1}$ is the Euclidean norm of $U$, $Q$ is some orthogonal matrix, and $e_1:=(1,0,\dots,0)$. Then it is clear that the matrix $M$ is similar to $N:=(1s_1)I_k+U^2e_1e_1^T=(1s_1)I_k+s_1e_1e_1^T$, whence $\det M=\det N=(1s_1)^{k1}$. So, \begin{equation*} J=s_2^{k/2}(1s_1)^{k/21}. \tag{2} \end{equation*} Also, the joint pdf of $(G_1,\dots,G_n)$ is given by \begin{equation*} (2\pi)^{n/2}\exp\Big\{\frac12\sum_1^n g_j^2\Big\}. \end{equation*} So, in view of (1) and (2), the joint pdf of $(Z_1,\dots,Z_n)$ is given by \begin{equation*} f_n(z_1,\dots,z_n) =(2\pi)^{n/2}\exp\Big\{\frac12\frac{s_2}{1s_1}\Big\}s_2^{k/2}(1s_1)^{k/21}. \end{equation*} So, the joint pdf of $(Z_1,\dots,Z_k)$ is given by \begin{align*} f_k(z_1,\dots,z_k)&=\int_{\mathbb R^{nk}}dz_{k+1}\dots dz_n\,f_n(z_1,\dots,z_n) \\ &=(2\pi)^{n/2}(1s_1)^{k/21} \\ &\times \int_{\mathbb R^{nk}}dz_{k+1}\dots dz_n\,s_2^{k/2}\,\exp\Big\{\frac12\frac{s_2}{1s_1}\Big\} \\ &\propto(1s_1)^{(nk)/21}, \end{align*} because $s_2=\sum_{k+1}^n z_j^2$. So, we have the desired result.
More specifically, \begin{align*} f_k(z_1,\dots,z_k) =\frac{2^{(kn)/2}\Gamma(n/2)}{\pi^{n/2}\Gamma((nk)/2)}(1s_1)^{(nk)/21} \end{align*} for $s_1=\sum_1^k z_j^2\in(0,1)$.
$\newcommand{\R}{\mathbb{R}} \newcommand{\x}{\mathbf{x}} \newcommand{\X}{\mathbf{X}}$ Here is yet another solution, which is partly informal but I think not hard to completely formalize. Its advantage is a strong and hopefully convincing appeal to geometric intuition.
Again, we are going to assume that $(X_1,\dots,X_n)$ is uniformly distributed on the unit sphere $\mathbb S^{n1}$, rather than on $\sqrt n\,\mathbb S^{n1}$. Let $\X_j:=(X_1,\dots,X_j)$ and $\x_j:=(x_1,\dots,x_j)$ for $\x=(x_1,\dots,x_n)\in\R^n$ and $j=1,\dots,n1$. Let $\cdot$ denotes the Euclidean norm.
The main point is that the probability density $p_{n1}(\x_{n1})$ of $\X_{n1}$ at a point $\x_{n1}\in\R^{n1}$ with $\x_{n1}<1$ is proportional to the ratio $r_{n1}(\x_{n1}):=vol_{n1}(dS)/vol_{n1}(dA)$, where $vol_{n1}$ is of course the $(n1)$volume, $dA$ is an infinitesimal neighborhood of the point $\x_{n1}$ in $\R^{n1}$ and $dS$ is the preimage of $dA$ under the projection of the upper hemisphere $\mathbb S^{n1}_+:=\{\x\in\mathbb S^{n1}\colon\x\cdot e_n\ge0\}$ onto the closed unit ball in $\R^{n1}$; this projection is given by $\mathbb S^{n1}_+\ni(\x_{n1},u)\mapsto\x_{n1}$; here $e_n:=(0,\dots,0,1)$ and $\cdot$ denotes the dot product. But \begin{equation} r_{n1}(\x_{n1})=\frac{vol_{n1}(dS)}{vol_{n1}(dA)}=\frac1{\cos\phi}, \end{equation} where $\phi$ is the angle between the hyperplane $\R^{n1}\times\{0\}$ of $\R^n$ and the tangent hyperplane to $\mathbb S^{n1}$ at the point $(\x_{n1},\sqrt{1\x_{n1}^2})\in\mathbb S^{n1}_+$; that is, $\phi$ is the angle between the corresponding normal vectors $e_n$ and $(\x_{n1},\sqrt{1\x_{n1}^2})$ of these two hyperplanes.
Thus, \begin{equation} p_{n1}(\x_{n1})\propto r_{n1}(\x_{n1})\propto\frac1{\cos\phi}=(1\x_{n1}^2)^{1/2}, \end{equation} where $\propto$ means an equality up to a constant factor, depending only on $n,k$.
Now successively integrating $p_{n1}(\x_{n1})$ ($n1k$ times) in $x_{n1},\dots,x_{k+1}$ and each time using the formula \begin{equation*} \int_0^{b^{1/2}}(bu^2)^p du=c_p b^{p+1/2} \end{equation*} for real $b>0$ and $p>1$ with $c_p:=\int_0^1(1u^2)^pdu\in(0,\infty)$ (so that $1/2$ is added to the exponent $p$ after such an integration), we see that the joint pdf of $\X_k=(X_1,\dots,X_k)$ is given by \begin{equation*} p_k(\x_k)\propto(1\x_k^2)^{(nk)/21}\,1_{\x_k<1}, \end{equation*} as desired.
Here is a picture, for $n=3$, showing the upper hemisphere $\mathbb S^{n1}_+$; a small neighborhood of a point $\x_{n1}$ in the projection of $\mathbb S^{n1}_+$ to the horizontal plane $\R^{n1}\times\{0\}$ of $\R^n$; the preimage of that neighborhood under that projection; and the normal vectors of the horizontal plane and the tangent plane to the sphere  with $\phi$ being the angle between these two normal vectors.
I enjoyed thinking about these answers and this is my attempt to put them into (nonrigorous) geometrical terms. Writing the joint density compositionally as
$$p(\mathbf{x}_k \mid \mathbf{x} = 1)p(\mathbf{x}_{nk} \mid \mathbf{x}_k, \mathbf{x} = 1) = p(\mathbf{x} \mid \mathbf{x} = 1) \propto 1,$$
we want to solve for the first term on the left. But since our density is proportional to a constant, this is just
$$p(\mathbf{x}_k \mid \mathbf{x} = 1) \propto \frac{1}{p(\mathbf{x}_{nk} \mid \mathbf{x}_k, \mathbf{x} = 1)}.$$
Accordingly, instead of performing our calculation by integrating out $X_{k+1} \dots X_{n}$, we can think about it in terms of the conditional density for sampling $\mathbf{X}_{nk}$, given $\mathbf{X}_k$ and the norm constraint $\mathbf{X} = 1$, denoted in the denominator above.
I propose a two step procedure. First, draw a point uniformly from within the $nk1$ dimensional ball with radius $$r_{nk} = \sqrt{1  \mathbf{x}_k^2}.$$ Each such vector has density proportional to $$\left(1  \mathbf{x}_k^2\right)^{(nk1)/2}.$$
This corresponds to the second and higher integrations in the previous answers, whereas here we directly use the formula for the volume of a ball.
Next, the $(nk)$th coordinate must satisfy $\mathbf{x} = 1$, which is achieved by any point on the circle with radius $r_{nk}$, a set with measure proportional to $r_{nk}$. Proving this is the first integration in the previous answers.
Putting these two steps together and taking the reciprocal gives
$$p(\mathbf{x}_k \mid \mathbf{x} = 1) \propto \left(1  \mathbf{x}_k^2\right)^{(nk)/2  1}.$$