Can the Fibonacci lattice be extended to dimensions higher than 3?

So, as you pointed out, the idea with Fibonacci lattice is based on the fact that when distributed uniformly, both the azimuthal angle $\phi$ and latitude sine (affine with $\cos\theta$) are uniform.

To distribute my points on the $d$ dimensional surface of the sphere $\mathbb{S}^d$ embedded in $d+1$ dimensions, first I need to introduce coordinates. Let's use the generalized polar coordinate system. $$(r, \theta^1, \theta^2, \cdots, \theta^d)$$ This is converted to Cartesian coordinates as $$x^{d+1}=r\cos\theta^d$$ $$x^{d}=r\sin\theta^d\cos\theta^{d-1}$$ $$\vdots$$ $$x^{1}=r\sin\theta^d\cdots\sin\theta^1$$ Now to produce points on $\mathbb{S}^d$, we need $r=1$. It remains to come up with some distribution over $\theta$ coordinates that induces a uniform surface measure on the sphere surface. $$\rho_{\vec\theta}(\vec\theta)=\rho_{\theta^d}(\theta^d)\rho_{\theta^{d-1}|\theta^d}(\theta^{d-1})\cdots\rho_{\theta^1|\theta^2\cdots\theta^d}(\theta^1)$$ An interesting feature of our polar coordinates is that fixing the last $k$ angular coordinates, induces the same coordinate system for a $d+1-k$ dimensional space, assuming one uses the redefinition $$r\rightarrow r\sin\theta^d\cdots\sin\theta^{d-k+1}$$ (Check!) But this means that the angles are independently distributed since the above redefinition only rescales the space and therefore leaves the angular distributions unchanged. We find $$\rho_{\vec\theta}(\vec\theta)=\prod_{\alpha=1}^d\rho_\alpha(\theta^\alpha)$$

It is not hard to prove that $$\rho_\alpha=\frac{1}{\sqrt\pi}\frac{\Gamma\big(\frac{\alpha+1}{2}\big)}{\Gamma\big(\frac{\alpha}{2}\big)}\sin^{\alpha-1}$$ Although for $\alpha=1$ this formula generates a uniform angle in $(0,\pi)$. We assume that this issue has been resolved manually.

Denote the C.D.F's with $$Y^\alpha(\theta^\alpha)=\int_0^{\theta^\alpha}\rho_\alpha(u)du$$ They satisfy (as they should) $$Y^\alpha(0)=0,\hspace{5mm}Y^\alpha(\pi)=1$$ So far we have proved that a point $\vec X$ will have uniform distribution on $\mathbb{S}^d$ if $\vec{X}=\vec{X}(1, \vec\theta)$ and $\vec\theta=\vec\theta(\vec Y)$ if and only if $$\vec Y\sim \textrm{Uni}\big([0,1]^d\big)$$

Remark: The map from $Y$ hypercube to the sphere is $\pi-$Lipschitz!

Now this was what you didn't want :)) A random point ensemble on the sphere. To make it Fibonacci-like, we won't be taking random independent points from the hypercube of $Y$'s. Instead for $n\in\big[N\big]$ generate the points $\vec{Y}_n$ with $$Y^d_n=\frac{n}{N+1}$$ $$Y^{d-1}_n=\{na_1\}$$ $$Y^{d-2}_n=\{na_2\}$$ $$\vdots$$ $$Y^1=\{na_{d-1}\}$$ For fixed irrational numbers $a_i$ satisfying $$\frac{a_i}{a_j}\not\in\mathbb{Q}\hspace{4mm}\forall i\neq j$$ A series of numbers that you can always use for $a$'s is $$\sqrt2, \sqrt3, \sqrt5, \sqrt7, \sqrt{11}, \sqrt {13}, \cdots$$

Example: $d=3$ with above $a$'s. $$t=\cos\psi$$ $$z=\sin\psi\cos\theta$$ $$y=\sin\psi\sin\theta\cos\phi$$ $$x=\sin\psi\sin\theta\sin\phi$$

The angles are given by $$\psi_n=f^{-1}\big(\frac{n\pi}{N+1}\big), \hspace{4mm}f(x)\equiv x+\frac{1}{2}\sin2x$$ $$\theta_n=\arccos\big(1-2\{n\sqrt2\}\big)$$ $$\phi_n=2\pi\{n\sqrt3\}$$