Marginal density of uniform spherical distribution

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_{n-1}(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_{n-2}(X_3,\ldots X_n)\propto\theta\left(1-\sum_{j=3}^n X_j^2\right),$$ the third integration $$P_{n-3}(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_{n-k+1},\ldots X_n)\propto\left(1-\sum_{j=n-k+1}^n X_j^2\right)^{(n-k)/2-1}\theta\left(1-\sum_{j=n-k+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_{n-1}$, $$\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^{n-1} dr \int d\Omega\, f(r,\phi_1,\phi_2,\ldots\phi_{n-1})\frac{2}{A_n}\delta(1-r^2)$$ $$\qquad=\frac{1}{A_n}\int d\Omega\, f(r=1,\phi_1,\phi_2,\ldots\phi_{n-1}).$$ In the last step I used that $\int_0^\infty r^{n-1}\delta(1-r^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(1-s_2-X_1^2),$$ $$\qquad=\int_0^\infty\delta(1-s_2-q)\frac{dq}{\sqrt q}=(1-s_2)^{-1/2}\theta(1-s_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^2-x^2)^p\,dx=c_p a^{1+2p}.$$}

$\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^{n-1}$, rather than on $\sqrt n\,\mathbb S^{n-1}$.

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_{1-t<|\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_{n-1}) \end{equation*} for some nonnegative continuous function $g\colon\R^{n-1}\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,n-1$. Then
\begin{align*} \int f\,d\mu_t&=\int_{\R^{n-1}}d\x_{n-1}\,g(\x_{n-1})\int_\R du\, 1_{1-t<|\x_{n-1}|^2+u^2\le1} \\ &=\int_{\R^{n-1}}d\x_{n-1}\,g(\x_{n-1}) (1+o(1))t\,(1-|\x_{n-1}|^2)^{-1/2}\,1_{|\x_{n-1}|<1}. \end{align*} Also, \begin{equation*} \int d\mu_t=\int_{\R^n}d\x\, 1_{1-t<|\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^{n-1}}d\x_{n-1}\,g(\x_{n-1}) (1-|\x_{n-1}|^2)^{-1/2}\,1_{|\x_{n-1}|<1}. \end{equation*} Thus, the joint pdf of $\X_{n-1}=(X_1,\dots,X_{n-1})$ is given by \begin{equation*} p_{n-1}(\x_{n-1})\propto(1-|\x_{n-1}|^2)^{-1/2}\,1_{|\x_{n-1}|<1}. \end{equation*} Now successively integrating $p_{n-1}(\x_{n-1})$ ($n-1-k$ times) in $x_{n-1},\dots,x_{k+1}$ and each time using the formula \begin{equation*} \int_0^{b^{1/2}}(b-u^2)^p du=c_p b^{p+1/2} \end{equation*} for real $b>0$ and $p>-1$ with $c_p:=\int_0^1(1-u^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)^{(n-k)/2-1}\,1_{|\x_k|<1}, \end{equation*} as desired.