I've found a somewhat nicer proof than the version on the Visual Insight blog. This new approach doesn't entail a huge conceptual breakthrough, but it does avoid having to deal with 3072 individual cases.

We have some freedom in choosing $R$ and $S$, but this makes no difference to the resulting sets of axes. We will pick:

$$\displaystyle{R = \frac{1}{4} \left(
\begin{array}{ccc}
2 & 1-\sqrt{5} & -1-\sqrt{5} \\
1-\sqrt{5} & 1+\sqrt{5} & -2 \\
1+\sqrt{5} & 2 & -1+\sqrt{5} \end{array} \right)}$$

$$\displaystyle{S = \frac{1}{4} \left(
\begin{array}{ccc}
\sqrt{5}-1 & -2 & -1-\sqrt{5} \\
2 & 1+\sqrt{5} & 1-\sqrt{5} \\
1+\sqrt{5} & 1-\sqrt{5} & 2 \end{array}\right)}$$

All the powers of these matrices can again be written with a denominator of 4 and numerators taken from $\{\pm 2, \pm 1 \pm \sqrt{5}\}$. We can simplify things a bit by working with integer matrices in 6 dimensions. For each of the four powers of $R$ and $S$, we will multiply the matrix by 4 and then write it as a linear map between 6-dimensional spaces with separate components for the rational and irrational parts of each component of the original vector. For example, for the first power of $R$ we get:

$$\displaystyle{R_6 = \left(
\begin{array}{cccccc}
2 & 0 & 1 & -5 & -1 & -5 \\
0 & 2 & -1 & 1 & -1 & -1 \\
1 & -5 & 1 & 5 & -2 & 0 \\
-1 & 1 & 1 & 1 & 0 & -2 \\
1 & 5 & 2 & 0 & -1 & 5 \\
1 & 1 & 0 & 2 & 1 & -1 \end{array}\right)}$$

and for the first power of $S$ we get:

$$\displaystyle{S_6 = \left(
\begin{array}{cccccc}
-1 & 5 & -2 & 0 & -1 & -5 \\
1 & -1 & 0 & -2 & -1 & -1 \\
2 & 0 & 1 & 5 & 1 & -5 \\
0 & 2 & 1 & 1 & -1 & 1 \\
1 & 5 & 1 & -5 & 2 & 0 \\
1 & 1 & -1 & 1 & 0 & 2 \end{array}\right)}$$

Suppose we have some unit vector $v$ of the form:

$$v = (a + b \sqrt{5}, c + d \sqrt{5}, e + f \sqrt{5}) / 2^{N+1}$$

where $a, b, c, d, e, f$ are integers, with at least one of them odd, and $N \ge 1$. We will work with the integer vector:

$$w = (a, b, c, d, e, f)$$

Because $v$ is a unit vector, the components of $w$ will satisfy the conditions:

$$a^2 + c^2 + e^2 + 5(b^2 + d^2 + f^2) = 4^{N+1}$$

and

$$a b + c d + e f = 0$$

If we multiply any vector of this form by each of the eight $6 \times 6$ matrices corresponding to the four powers of $R$ and $S$, then it turns out that precisely one of those eight matrices will yield a result equal to the zero vector modulo 8, i.e. a 6-tuple of integers all divisible by 8.

To prove this, we take the lattice of vectors in $\mathbb{Z}^6$ equal to the zero vector modulo 8, and multiply it by the inverse of each of the eight matrices in turn, to produce eight new lattices: lattices which yield 6-tuples of integers all divisible by 8 when multiplied by the appropriate matrix. In concrete terms, for a given matrix $M_i$, the basis for the associated lattice is given by the row vectors of $L_i = 8 (M_i^{-1})^T$.

The original claim can now be restated as saying that every vector $w$ that meets the conditions described above will belong to the union of the eight lattices $L_i$, but no such vector will belong to the intersection of any two of the $L_i$.

The first part is fairly easy to show. We can obtain a basis for the union of the eight lattices by forming a matrix whose rows are the union of all eight bases, and then reducing that $48 \times 6$ matrix to a $6 \times 6$ matrix by putting it in Hermite Normal Form (the equivalent of reduced row-echelon form for integer matrices), and discarding all rows containing only zeroes. We will call that matrix $L_U$. The test for the vector $w$ belonging to the lattice whose basis is given by the rows of $L_U$ is that the vector $(L_U^{-1})^T w$ has all integer coordinates. When we carry through these calculations, we find:

$$(L_U^{-1})^T w = (2a, b-a, 2c, d-c, e-a-c, \frac{a-b+c-d+f-e}{2})$$

Since $a, b, c, d, e, f$ are integers, the only thing remaining to show is that $a-b+c-d+f-e$ will always be an even integer, given the conditions we've placed on $w$.

We have the conditions:

$$a^2 + c^2 + e^2 + 5(b^2 + d^2 + f^2) = 4^{N+1}$$
$$a b + c d + e f = 0$$

It follows that:

$$a^2 + c^2 + e^2 + 5(b^2 + d^2 + f^2) = 0 \mod 4$$
$$a^2 + c^2 + e^2 + 5(b^2 + d^2 + f^2) - 4(b^2 + d^2 + f^2) - 2(a b + c d + e f)= 0 \mod 4$$
$$(a-b)^2 + (c-d)^2 + (e-f)^2 = 0 \mod 4$$

It's not hard to check that the sum of three squares can only be a multiple of 4 if all three of the numbers being squared are even. So we have $a-b, c-d$ and $e-f$ all individually even, so $a-b+c-d+f-e$ will be even, $w$ will belong to the lattice $L_U$, and at least one of the eight matrices multiplied by $w$ will yield a vector whose components are all divisible by 8.

To prove that only one matrix yields such a result for a given $w$, we need to show that the intersection of any pair of distinct lattices $L_i$ and $L_j$ cannot contain any vector $w$ meeting the conditions we've imposed. There are 28 such pairs of lattices. Finding their bases is a bit more involved than finding the basis for a union. First, we need to construct the dual of each lattice. The dual of a lattice $L_i$ is the set of vectors $d$ such that for every $v \in L_i$, the dot product $d \cdot v$ is an integer. Its basis is given by the rows of the matrix:

$$D_i = (L_i L_i^T)^{-1} L_i$$

We obtain a basis for the intersection of two lattices by forming their dual lattices, finding a basis for the union of those duals (by joining their matrices and reducing it to Hermite Normal Form), and then taking the dual of that union.

If we do this for the 28 pairs of lattices, we find that in 16 cases the intersection of the lattices contains only vectors whose coordinates are all even integers. This violates the requirement that at least one coordinate be odd (which we impose in order that the corresponding vector divided by a power of 2 is in lowest terms).

For the remaining 12 pairs of lattices, the requirement that at least one coordinate be odd can be satisfied if and only if one particular element of the lattice basis has an odd coefficient $\ell$ in the sum that decribes the vector $w$ with respect to that basis. But that in turn clashes with the requirement that:

$$a^2 + c^2 + e^2 + 5(b^2 + d^2 + f^2) = 4^{N+1}$$

The contradiction appears if we require the equation to continue to hold modulo 8. In each case all but one of the lattice coefficients vanish, and what we end up with is:

$$4 \ell^2 = 4^{N+1} \mod 8$$

which is impossible for $N\ge 1$ and odd $\ell$.

Because $R$ and $S$ are rotations of order 5, the set of their first 4 powers can also be seen as the set of inverses of their first 4 powers. Because the corresponding integer matrices are multiplied by a factor of 4, a result that is a multiple of 8 corresponds to a factor of 2 in the original matrices. So what we have established is that, given any unit vector over the golden field with a denominator of $2^{N+1}$ for some $N \ge 1$, the inverse of precisely one of the powers of $R$ and $S$ will take us to another unit vector with a denominator of $2^N$. As we repeat this process, we will move back through the tree to ever smaller denominators, eventually terminating with the original cube.

This means that we can reach every unit vector $v$ of this form as one of the cube axes or their opposites, and also that we can only reach it via a single path.