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^{n-1}$. 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{1-s_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}(1-s_1)^{-3/2},\quad M:=(1-s_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=|U|Qe_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:=(1-s_1)I_k+|U|^2e_1e_1^T=(1-s_1)I_k+s_1e_1e_1^T$, whence $\det M=\det N=(1-s_1)^{k-1}$. So, \begin{equation*} J=s_2^{k/2}(1-s_1)^{-k/2-1}. \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}{1-s_1}\Big\}s_2^{k/2}(1-s_1)^{-k/2-1}. \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^{n-k}}dz_{k+1}\dots dz_n\,f_n(z_1,\dots,z_n) \\ &=(2\pi)^{-n/2}(1-s_1)^{-k/2-1} \\ &\times \int_{\mathbb R^{n-k}}dz_{k+1}\dots dz_n\,s_2^{k/2}\,\exp\Big\{-\frac12\frac{s_2}{1-s_1}\Big\} \\ &\propto(1-s_1)^{(n-k)/2-1}, \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^{(k-n)/2}\Gamma(n/2)}{\pi^{n/2}\Gamma((n-k)/2)}(1-s_1)^{(n-k)/2-1} \end{align*} for $s_1=\sum_1^k z_j^2\in(0,1)$.

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^{n-1}$. 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{1-s_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}(1-s_1)^{-3/2},\quad M:=(1-s_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=|U|Qe_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:=(1-s_1)I_k+|U|^2e_1e_1^T=(1-s_1)I_k+s_1e_1e_1^T$, whence $\det M=\det N=(1-s_1)^{k-1}$. So, \begin{equation*} J=s_2^{k/2}(1-s_1)^{-k/2-1}. \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}{1-s_1}\Big\}s_2^{k/2}(1-s_1)^{-k/2-1}. \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^{n-k}}dz_{k+1}\dots dz_n\,f_n(z_1,\dots,z_n) \\ &=(2\pi)^{-n/2}(1-s_1)^{-k/2-1} \\ &\times \int_{\mathbb R^{n-k}}dz_{k+1}\dots dz_n\,s_2^{k/2}\,\exp\Big\{-\frac12\frac{s_2}{1-s_1}\Big\} \\ &\propto(1-s_1)^{(n-k)/2-1}, \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{\Gamma(n/2)}{\pi^{k/2}\Gamma((n-k)/2)}(1-s_1)^{(n-k)/2-1} \end{align*} for $s_1=\sum_1^k z_j^2\in(0,1)$.

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^{n-1}$. 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{1-s_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}(1-s_1)^{-3/2},\quad M:=(1-s_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=|U|Qe_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:=(1-s_1)I_k+|U|^2e_1e_1^T=(1-s_1)I_k+s_1e_1e_1^T$, whence $\det M=\det N=(1-s_1)^{k-1}$. So, \begin{equation*} J=s_2^{k/2}(1-s_1)^{-k/2-1}. \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}{1-s_1}\Big\}s_2^{k/2}(1-s_1)^{-k/2-1}. \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^{n-k}}dz_{k+1}\dots dz_n\,f_n(z_1,\dots,z_n) \\ &=(2\pi)^{-n/2}(1-s_1)^{-k/2-1} \\ &\times \int_{\mathbb R^{n-k}}dz_{k+1}\dots dz_n\,s_2^{k/2}\,\exp\Big\{-\frac12\frac{s_2}{1-s_1}\Big\} \\ &\propto(1-s_1)^{(n-k)/2-1}, \end{align*} because $s_2=\sum_{k+1}^n z_j^2$. So, we have the desired result.

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^{n-1}$. 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{1-s_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}(1-s_1)^{-3/2},\quad M:=(1-s_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=|U|Qe_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:=(1-s_1)I_k+|U|^2e_1e_1^T=(1-s_1)I_k+s_1e_1e_1^T$, whence $\det M=\det N=(1-s_1)^{k-1}$. So, \begin{equation*} J=s_2^{k/2}(1-s_1)^{-k/2-1}. \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}{1-s_1}\Big\}s_2^{k/2}(1-s_1)^{-k/2-1}. \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^{n-k}}dz_{k+1}\dots dz_n\,f_n(z_1,\dots,z_n) \\ &=(2\pi)^{-n/2}(1-s_1)^{-k/2-1} \\ &\times \int_{\mathbb R^{n-k}}dz_{k+1}\dots dz_n\,s_2^{k/2}\,\exp\Big\{-\frac12\frac{s_2}{1-s_1}\Big\} \\ &\propto(1-s_1)^{(n-k)/2-1}, \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^{(k-n)/2}\Gamma(n/2)}{\pi^{n/2}\Gamma((n-k)/2)}(1-s_1)^{(n-k)/2-1} \end{align*} for $s_1=\sum_1^k z_j^2\in(0,1)$.