This question is adapted from another question on the 2008 Putnam test which asks specifically for the case when $n = 4$ and $k = 2$. The answer is $\dfrac{1}{2}\sqrt{\dfrac{n}{k}}$ but I am looking for a proof or other justification.
In many of the attempts I have seen for similar problems, the hypercube is scaled by a factor of $2$ and centered at the origin.
Here are some posts regarding the Putnam and a discussion for the case when $k = 2$ for all values of $n$.
I would also be grateful for intuition regarding why the maximum radius increases and decreases with square roots when the dimension of the ball and cube are changed.
Like you said, it's easiest to consider the cube $[-1,1]^n$. I'll try to generalize the argument made in the link you submitted. We can assume by symmetry that the center of the sphere is the origin. Consider a $k-1$ sphere in $k-1$ spherical coordinates $\theta_i$ for $i=1,k-1$ with the range of $\theta_{k-1}$ being $[0,2\pi)$ and each other having range $[0,\pi]$, and a radius of $r$. For any orthonormal set $\{x^i\}_{i=1}^k$ of vectors in $\Bbb{R}^n$ the sphere can be parametrized as
$$S(\theta_1,\ldots,\theta_{k-1}) = r(\cos(\theta_1)x^1 + \sin(\theta_1)\cos(\theta_2)x^2+\cdots+\sin(\theta_1)\sin(\theta_2)\cdots\cos(\theta_{k-1})x^{k-1}\\ +\sin(\theta_1)\cdots\sin(\theta_{k-1})x^{k})$$
Each component of this must lie in $[-1,1]$. Let $x^j_i = \langle e_i,x^j\rangle$ be the $i$th component of $x^j$. The $i$th component of $S$ is then $$S_i = r(\cos(\theta_1)x^1_i + \sin(\theta_1)\cos(\theta_2)x^2_i+\cdots+\sin(\theta_1)\sin(\theta_2)\cdots\cos(\theta_{k-1})x^{k-1}_i\\ +\sin(\theta_1)\cdots\sin(\theta_{k-1})x^{k}_i)$$
This can be written as $r\langle x_i^T, \Theta\rangle$ where $\Theta = (\cos\theta_1,\dots, \sin\theta_1\dots\sin\theta_{k-1})$. We have $||\Theta|| = 1$. Thus by the identity $\langle x,y\rangle = ||x|| ||y|| \cos(\theta) $ and the fact that $\Theta$ ranges over the entire sphere (hence $\theta$ can take the value $0$) we get $r||x|| \le 1 $ in order to fit inside the interval $[-1,1]$ and this gives $$\sum_{j=1}^k (x_i^j)^2 \le \frac{1}{r^2}.$$ Summing over $i$ gives $$\sum_{i=1}^n\sum_{j=1}^k (x_i^j)^2 \le \frac{n}{r^2}$$ and by changing the order or summation and using orthonormality we get $$k\le \frac n {r^2}$$ or, as desired $$r\le \sqrt\frac{n}{k}.$$
Scaling back down to a unit cube gives the bound that you proposed.
EDIT:
I was originally wrong about the way of attaining the bound. Googling the answer I found an article which proved a slightly general result: it's an interesting read.