First of all a note about cartesian, polar and cylindrical coordinates, in Euclidean space
- in 1D, you have only one coordinate, so no distinction;
- in 2D, you have either Cartesian $(x_1 , x_2)$ or polar $(r ; \alpha)$;
- in 3D, you have Cartesian $(x_1 , x_2 , x_3)$, polar $(r ; \alpha_1,\alpha_2)$ or cylindrical in which two of the cartesian coordinates are replaced
with a couple of 2D polar coordinates;
- in 4D, you have the four cartesian, the "full-polar" (1 radius, 3 angles), and different cylindrical in between, for example
1D c+ 3D spherical, or two 2D polar or 2D c. + 2P polar.
In general it is clear that you can partition a n-D space into subspaces and use for each sub-space any coordinate system.
Let's concentrate on the 4D case. There the equation equation of a unit sphere is
$$
x_{\,1} ^{\,2} + x_{\,2} ^{\,2} + x_{\,3} ^{\,2} + x_{\,4} ^{\,2} = 1\quad \left| {\; - 1 \le x_{\,k} \le 1} \right.
$$
A section at constant $x_4$ is a sphere in the 3D space.
The surface area of a $n-1$-d sphere (embedded into a $n$D space) is
$$
S_{\,n - 1} (r) = {{2\pi ^{\,\,n/2} } \over {\Gamma \left( {n/2} \right)}}r^{\,\,n - 1} \quad \Rightarrow \quad \left\{ \matrix{
S_{\,3} = S_{\,4 - 1} (1) = 2\pi ^{\,\,2} \hfill \cr
S_{\,2} (r) = 4\pi r^{\,\,2} \hfill \cr} \right.
$$
and we get therefrom the values for the spheres mentioned above.
We can think of various repartitions of the points between two different dimensions , as you did and as indicated in
the excellent answer by @ksadri.
But the spiral principle, as far as I understood from the link you cited, is to construct a 1D line of length $L$ which traverses
our 3D surface uniformly, i.e. such that the ratio $dl/dS_3 = L/S_3$ among the length of the line and the area "spanned" by it is constant.
After that we drop the given $N$ points randomly uniformly over $L$.
This can be done in various ways, but the one suggested in the link is to use the equidistribution of the multiples of an irrational
number, that is of the sequence
$$
\left\{ {k\,\alpha } \right\} = \left\{ {\left\lfloor {{l \over L}N} \right\rfloor \,\alpha } \right\}\quad \left| {\;1 \le k \le N} \right.
$$
where $\alpha$ is a chosen irrational number, and the curly brackets indicate the
fractional part.
So let's try and create such a spiral.
We take the spherical coordinates
$$
\left\{ \matrix{
x_{\,1} = \sin \varphi _{\,3} \sin \varphi _{\,2} \sin \varphi _{\,1} \hfill \cr
x_{\,2} = \sin \varphi _{\,3} \sin \varphi _{\,2} \cos \varphi _{\,1} \hfill \cr
x_{\,3} = \sin \varphi _{\,3} \cos \varphi _{\,2} \hfill \cr
x_{\,4} = \cos \varphi _{\,3} \hfill \cr} \right.\quad \left| \matrix{
\;0 \le \varphi _{\,3} ,\varphi _{\,2} \le \pi \hfill \cr
\;0 \le \varphi _{\,1} < 2\pi \hfill \cr} \right.
$$
individuated by the three angles $ \varphi _{\,1}, \varphi _{\,2}, \varphi _{\,3}$.
We have better to parametrize the line in the parameter $t$ and figure it as a time.
With reference to the sketch
we assign to $\varphi _{\,3}$ a constant angular speed $\omega_3 = \pi /T$
$$
\varphi _{\,3} ' = \omega_3 = {\pi \over T}\quad \varphi _{\,3} = \pi {t \over T}\quad \left| {\;0 \le t \le T} \right.
$$
The area of the spherical zone spanned by $\varphi_3$ in $dt$ will be
$$
S_{\,2} (\sin \varphi _{\,3} ) = 4\pi \sin ^2 \varphi _{\,3} \,d\varphi _{\,3} = 4\pi ^2 \sin ^2 \varphi _{\,3} {{dt} \over T}
$$
therefore if we assign to the component of speed in $S_2$ a value proportional to that, we will satisfy the requirement of "uniform coverage"
(the spiral will undergo a number of "turns" proportional to the area).
$$
w_{\,2} = c_2 \sin ^2 \varphi _{\,3}
$$
After that we want again to repart $w_2$ in the same way between
$$
v_{\,2} = \sin \varphi _{\,3} \,\omega _{\,2} \quad v_{\,1} = \sin \varphi _{\,3} \,\sin \varphi _{\,2} \,\omega _{\,1}
$$
which means
$$
\eqalign{
& \left\{ \matrix{
w_{\,2} ^2 = v_{\,2} ^2 + v_{\,1} ^2 \hfill \cr
{{v_{\,1} } \over {v_{\,2} }} \propto {{2\pi \sin \varphi _{\,2} \;d\varphi _{\,2} } \over {d\varphi _{\,2} }} \hfill \cr} \right.\quad \Rightarrow \cr
& \Rightarrow \quad \left\{ \matrix{
{{v_{\,1} } \over {v_{\,2} }} = {{\sin \varphi _{\,2} \,\omega _{\,1} } \over {\,\omega _{\,2} }} = c_{\,1} \sin \varphi _{\,2} \hfill \cr
w_{\,2} ^2 = v_{\,2} ^2 + v_{\,1} ^2 \hfill \cr} \right.\quad \Rightarrow \cr
& \Rightarrow \quad \left\{ \matrix{
{{v_{\,1} } \over {v_{\,2} }} = c_{\,1} \sin \varphi _{\,2} \hfill \cr
c_{\,2} ^{\,2} \sin ^2 \varphi _{\,3} = v_{\,2} ^2 \left( {1 + c_{\,1} ^2 \sin ^2 \varphi _{\,2} } \right) \hfill \cr} \right.\quad \Rightarrow \cr
& \Rightarrow \quad \left\{ \matrix{
\omega _{\,1} = c_{\,1} \,\,\omega _{\,2} \hfill \cr
\omega _{\,2} = {{c_{\,2} \sin \varphi _{\,3} } \over {\sqrt {1 + c_{\,1} ^2 \,\sin ^2 \varphi _{\,2} } }}\quad \, \Rightarrow \quad \varphi _{\,2} ' = {{c_{\,2} \sin \left( {\pi {t \over T}} \right)} \over {\sqrt {1 + c_{\,1} ^2 \,\sin ^2 \varphi _{\,2} } }} \hfill \cr} \right. \cr}
$$
The ode
$$
y(t)'\sqrt {1 + c_{\,1} ^2 \,\sin ^2 y(t)} = g(t)
$$
can be solved by noting that
$$
y'{d \over {dy}}\int_{u = 0}^y {\sqrt {1\, + c_{\,1} ^2 \,\sin ^2 u} \;du} = {d \over {dt}}\int_{u = 0}^y {\sqrt {1\, + c_{\,1} ^2 \,\sin ^2 u} \;du}
$$
thus
$$
{d \over {dt}}\left( {{{\varphi _{\,2} \sqrt {1\, + c_{\,1} ^2 \,\varphi _{\,2} ^2 } } \over 2} + {{\ln \left( {\sqrt {1\, + c_{\,1} ^2 \,\varphi _{\,2} ^2 }
+ c_{\,1} \,\varphi _{\,2} } \right)} \over {2c_{\,1} }}} \right) = c_{\,2} \sin \left( {\pi {t \over T}} \right)
$$
after which it is a matter of calculus, or numerical integration, to find $l(t)$.
