How can I construct a basis of $\mathbb{R}^n$ with this property? Given $\theta \in (0,\pi/2]$ I would like to find a basis of unit length vectors $\{v_1, \dots, v_n\}$ for $\mathbb{R}^n$ such that each $v_i$ has angle $\theta$ from the subspace spanned by the remaining vectors. Are there any references for something like this? A (recursive) family of examples would be great but even some algorithmic ideas would be helpful.
For $n = 2$ it is trivial, take any unit length vector and rotate it by $\theta$.
For $n = 3$ I am already struggling. I have an ugly numeric approach that can find me one solution, however I am having a hard time generalizing it to dimensions $n > 3$. (I'm happy to share more details but it's not clever at all).
The properties of such a basis also means the pairwise angles between the $v_i$ will be the same, but for $n > 2$ it is larger than $\theta$ and (I imagine) continues to increase with the dimension. Knowing this angle in advance would also make constructing examples of these bases by hand much easier.
 A: Edit: Here's a construction for $n \ge 3$, generalizing the argument that was previously presented here for $n = 3$. For $n = 3$ the idea is to think of three vectors $v_1, v_2, v_3$ forming an equilateral triangle such that the line $x = y = z$ passes through its center.
We'll take our basis $\{ v_1, \dots v_n \}$ to consist of cyclic permutations of a single vector $v_1$; this guarantees the desired condition by cyclic symmetry, and it remains to find a choice of $v_1$ that lets us clearly control the resulting angle $\theta$. We'll take $v_1 = (a, b, b, b, \dots )$, so that $v_2 = (b, a, b, \dots)$ etc., where $a^2 + (n-1) b^2 = 1$. One way to motivate this choice is that when $a = 1, b = 0$ this recovers an orthonormal basis so $\theta = \frac{\pi}{2}$ whereas when $a = b = \frac{1}{\sqrt{n}}$ all the vectors agree so $\theta = 0$ (of course we don't have a basis in this case but we can take a limit as $a$ and $b$ approach this point).
Let $V = \text{span}(v_2, \dots v_n)$. We want to compute the angle between $v_1$ and $V$, or equivalently the angle between $v_1$ and $\text{proj}_V(v_1)$. By replacing $v_2, \dots v_n$ with $v_2, v_3 - v_2, v_4 - v_2, \dots$, and observing that $v_n - v_2 = (b - a)(0, 1, \dots, -1, \dots, 0)$ we see that $V = \text{span}(v_2) \oplus W$ where
$$\begin{align*} W &= \{ (x_1, \dots x_n) \in \mathbb{R}^n : x_2 + \dots + x_n = 0 \} \\
 &= \text{span}( (1, 0, 0, \dots), (0, 1, 1, 1, \dots) )^{\perp}. \end{align*}.$$
Write $e_1 = (1, 0, 0, \dots)$ and $e_{-1} = (0, 1, 1, \dots)$ to save notation from this point on. Importantly, $v_1 \perp W$, which simplifies our calculations substantially. This is because, recalling that $\text{proj}_W(v_1) = \text{proj}_{W_1}(v_1) + \text{proj}_{W_2}(v_1)$ if $W = W_1 \oplus W_2$ is an orthogonal direct sum, we have
$$\begin{align*} \text{proj}_V(v_1) &= \text{proj}_W(v_1) + \text{proj}_{v_2 - \text{proj}_W(v_2)}(v_1) \\
 &= \text{proj}_{v_2 - \text{proj}_W(v_2)}(v_1). \end{align*}$$
We can even compute $v_2 - \text{proj}_W(v_2)$ straightforwardly because it's equal to $\text{proj}_{W^{\perp}}(v_2)$ where $W^{\perp} = \text{span}(e_1, e_{-1})$ as above. These vectors are orthogonal which lets us compute this projection as
$$\begin{align*} w = \text{proj}_{W^{\perp}}(v_2) &= \frac{\langle v_2, e_1 \rangle}{\langle e_1, e_1 \rangle} e_1 + \frac{\langle v_2, e_{-1} \rangle}{\langle e_{-1}, e_{-1} \rangle} e_{-1} \\
 &= b e_1 + \frac{a + (n-2) b}{n-1} e_{-1} \\
 &= \left( b, \frac{a + (n-2) b}{n-1}, \frac{a + (n-2) b}{n-1}, \dots \right). \end{align*}$$
Finally, we get that $\text{proj}_V(v_1) = \text{proj}_w(v_1)$, so the angle between $v_1$ and $V$ is just the angle between $v_1$ and $w$, which satisfies
$$\boxed{ \begin{align*}  \cos \theta &= \frac{ \langle v_1, w \rangle }{\| w \|} \\
 &= \frac{2ab + (n-2) b^2}{\sqrt{b^2 + \frac{(a + (n-2) b)^2}{n-1} }} \end{align*} }.$$
This doesn't seem to simplify much even after using $a^2 + (n-1) b^2 = 1$, e.g. making a substitution like $a = \cos \phi, b = \frac{1}{\sqrt{n-1}} \sin \phi$. In any case, to show that every $\theta \in (0, \frac{\pi}{2}]$ is obtainable by this construction it suffices to show that this cosine takes every possible value in $[0, 1)$. The minimum value $0$ occurs when $a = 1, b = 0$ which corresponds to the $v_i$ being the standard orthonormal basis of $\mathbb{R}^n$. And the maximum value $1$ occurs when $a = b = \frac{1}{\sqrt{n}}$ which corresponds to each $v_i$ being the same vector. By continuity every intermediate value occurs.
Actually this was already clear without doing this computation but this computation is for actually finding $a, b$ given $\theta$ numerically. For example, for $\theta = \frac{\pi}{3}$ and $n = 3$, plugging the resulting system in $a$ and $b$ into WolframAlpha gives $a \approx 0.959 \dots, b \approx 0.199 \dots$.
Edit #2: There's some question in the comments about how easy it is to actually compute $a, b$ in general given $\theta$. We can argue as follows: given $\theta$, squaring the above identity, using the trig substitution $a = \cos \phi, b = \frac{1}{\sqrt{n-1}} \sin \phi$, and clearing denominators reduces the question to solving a polynomial in $\cos \phi, \sin \phi$ of degree $4$, which actually simplifies to a polynomial in $\cos 2 \phi, \sin 2 \phi$ of degree $2$. Now making the tangent half-angle substitution
$$\cos 2 \phi = \frac{1 - t^2}{1 + t^2}, \sin  2\phi = \frac{2t}{1 + t^2}$$
where $T = \tan \phi$ and clearing denominators reduces to solving a polynomial in $t$ of degree $4$ which should be straightforward to do numerically.
