It seems preposterous at first glance. I just want to be sure. Even in 3D the behaviour of rotating objects can be surprising (see the Dzhanibekov effect); in 4D it could be more surprising.
A 2D or 3D spacecraft (with no reaction wheels or gimbaling etc.) needs at least two thrusters to control its spin. See my answer on space.SE.
For any two $4\times4$ matrices $X$ and $Y$, define the commutator $[X,Y]=XY-YX$, the anticommutator $\{X,Y\}=XY+YX$, and the Frobenius inner product $\langle X,Y\rangle=\operatorname{tr}(X^\top Y)$, where $\operatorname{tr}$ is the trace and $\top$ is the transpose.
Let $M$ be a symmetric positive-definite $4\times4$ matrix, and $T=\mathbf f\,\mathbf r^\top-\mathbf r\,\mathbf f^\top$ an antisymmetric $4\times4$ matrix. ($M$ describes the distribution of mass in the spacecraft, $\mathbf r$ is a vector locating the thruster relative to the centre of mass, $\mathbf f$ is the force produced by the thruster, and $T$ is the torque produced by the thruster. More details here. Everything is described in the rotating reference frame.)
The angular velocity $\Omega(t)$, an antisymmetric $4\times4$ matrix, changes with time $t$ according to Euler's equation
$$\{M,\Omega'(t)\}+[\Omega(t),\{M,\Omega(t)\}]=f(t)\,T$$
where $f(t)\geq0$ is a function (continuous, piecewise-constant, or just integrable) describing when and how strongly the thruster is used.
Question: Can $M$ and $T$ be chosen such that, for any two antisymmetric $4\times4$ matrices $\Omega_0$ and $\Omega_1$, there exist $t_1>0$ and $f$ such that the solution $\Omega$ to Euler's equation with initial value $\Omega(0)=\Omega_0$ has final value $\Omega(t_1)=\Omega_1$?
To simplify things, we may take $f$ to be a combination of Dirac deltas instead of an ordinary function. Then the angular velocity satisfies $\{M,\Omega'(t)\}+[\Omega(t),\{M,\Omega(t)\}]=0$ except at a finite set of times when the angular momentum $L(t)=\{M,\Omega(t)\}$ changes by a positive multiple of $T$, or equivalently when $\Omega(t)$ changes by a positive multiple of $A=\{M,\}^{-1}(T)$, where $A$ is the angular acceleration produced by the thruster.
Here is Euler's equation in terms of the components $\omega_{ij}$ of the angular velocity. Assume that $M$ is diagonal, with components $m_i>0$.
$$(m_1+m_2)\,\omega_{12}'(t)+(m_2-m_1)\big(\omega_{13}(t)\,\omega_{32}(t)+\omega_{14}(t)\,\omega_{42}(t)\big)=f(t)\,\tau_{12}$$
(And permute the indices to get 6 equations like this.) Equivalently:
$$\omega_{12}'(t)=\frac{m_1-m_2}{m_1+m_2}\big(\omega_{13}(t)\,\omega_{32}(t)+\omega_{14}(t)\,\omega_{42}(t)\big)+f(t)\,\alpha_{12}$$
If $m_1=m_2$, then $\omega_{12}'$ has constant sign; $\omega_{12}$ is either always non-increasing, or always non-decreasing. So, if the difference between initial and final values of $\omega_{12}$ has the wrong sign compared to $\tau_{12}$, then there is no solution. Thus, we must take $m_1\neq m_2$, and similarly $m_i\neq m_j$ for $i\neq j$.
The Frobenius inner product of two antisymmetric matrices has every term duplicated: $\sum_{i,j}x_{ij}y_{ij}=2\sum_{i<j}x_{ij}y_{ij}$ (since $x_{ji}=-x_{ij}$ and $y_{ji}=-y_{ij}$). So it's natural to take half of this.
The angular momentum's squared magnitude is $\tfrac12\langle L(t),L(t)\rangle$; the derivative of this is $\langle L(t),f(t)T\rangle$. So the angular momentum has constant magnitude whenever $f(t)=0$. (The angular momentum itself would be constant in an inertial reference frame, but here we're using a rotating reference frame.)
The rotational energy is $\tfrac14\langle\Omega(t),L(t)\rangle$; the derivative of this is $\tfrac12\langle\Omega(t),f(t)T\rangle$. So the energy is constant whenever $f(t)=0$.
All of that applies in any dimension. But in 4D there's another constant, the angular momentum bivector's exterior square, or equivalently its inner product with its Hodge dual.
Here are those constants in component form. Squared magnitude of angular momentum:
$$(m_1+m_2)^2\omega_{12}^2+(m_1+m_3)^2\omega_{13}^2+(m_2+m_3)^2\omega_{23}^2\\+(m_1+m_4)^2\omega_{14}^2+(m_2+m_4)^2\omega_{24}^2+(m_3+m_4)^2\omega_{34}^2$$
Doubled energy:
$$(m_1+m_2)\omega_{12}^2+(m_1+m_3)\omega_{13}^2+(m_2+m_3)\omega_{23}^2\\+(m_1+m_4)\omega_{14}^2+(m_2+m_4)\omega_{24}^2+(m_3+m_4)\omega_{34}^2$$
Halved exterior square of angular momentum:
$$(m_1+m_2)(m_3+m_4)\omega_{12}\omega_{34}-(m_1+m_3)(m_2+m_4)\omega_{13}\omega_{24}+(m_2+m_3)(m_1+m_4)\omega_{23}\omega_{14}$$
(Its derivative is the exterior product of angular momentum and torque: $(m_1+m_2)\omega_{12}\,f\,\tau_{34}+\cdots$ .)
These three constants define two ellipsoids and another quadratic surface (with signature $+^3-^3$) in the 6D space of antisymmetric matrices. When $f=0$, $\Omega$ must remain on the intersection of these three quadratic surfaces.
In fact there are more than three, as I recently found! And these aren't specific to 4D. The general form is $Q_i=\sum_{j\neq i}\frac{m_i+m_j}{m_i-m_j}\omega_{ij}^2$; check that $Q_i'=0$, and also note $\sum_iQ_i=0$. Then define $P_i=\sum_{j\leq i}Q_j$; these are positive-semidefinite, on the assumption (WLOG) that $m_i$ are sorted, $m_1>m_2>m_3>m_4>0$. The doubled energy is $\sum_i(m_i-m_{i+1})P_i$, and the momentum's squared magnitude is $\sum_i(m_i^2-m_{i+1}^2)P_i$. But what we should focus on are the quadratic forms $N_i=P_{i-1}-P_i+P_{i+1}$ which are not positive-semidefinite but have a single negative term:
$$N_1=\frac{m_2+m_3}{m_2-m_3}\omega_{23}^2+\frac{m_2+m_4}{m_2-m_4}\omega_{24}^2-\frac{m_1+m_2}{m_1-m_2}\omega_{12}^2$$ $$N_2=\frac{m_1+m_2}{m_1-m_2}\omega_{12}^2+\frac{m_1+m_4}{m_1-m_4}\omega_{14}^2+\frac{m_3+m_4}{m_3-m_4}\omega_{34}^2-\frac{m_2+m_3}{m_2-m_3}\omega_{23}^2$$ $$N_3=\frac{m_1+m_3}{m_1-m_3}\omega_{13}^2+\frac{m_2+m_3}{m_2-m_3}\omega_{23}^2-\frac{m_3+m_4}{m_3-m_4}\omega_{34}^2$$
These are all constant. The first expression (if it's equated to $0$) defines a cone, or a pair of opposite cones, around the $\omega_{12}$ axis in 6D. The second expression defines a pair of opposite cones around the $\omega_{23}$ axis; and the third, around the $\omega_{34}$ axis. If the applied acceleration $A$ is inside one of these cones, and the initial angular velocity $\Omega_0$ is in the same cone, then $\Omega$ stays in that cone, and there's no solution with the final angular velocity $\Omega_1$ outside of the cone. Thus, we must take $A$ to be outside of all six of these cones. That is, all three expressions (with $\alpha_{ij}$ in place of $\omega_{ij}$) must be positive.
Notice that Euler's equation has a time symmetry: Given a solution $\Omega(t)$, we can construct another solution $\tilde\Omega(t)=-\Omega(t_1-t)$, so that $\{M,\tilde\Omega'(t)\}+[\tilde\Omega(t),\{M,\tilde\Omega(t)\}]=f(t_1-t)\,T$, and the initial and final values are swapped (and negated): $\tilde\Omega(0)=-\Omega_1$ and $\tilde\Omega(t_1)=-\Omega_0$. Therefore, any angular velocity can be reached from any other, if and only if any angular velocity can be reached from $0$. (We can reverse a path from $0$ to $-\Omega_0$ to get a path from $\Omega_0$ to $0$, and concatenate that with a path from $0$ to $\Omega_1$, to get a path from $\Omega_0$ to $\Omega_1$.)
Also, Euler's equation has a scale symmetry: Given a solution $\Omega(t)$, for any $k>0$ we can construct another solution $\tilde\Omega(t)=k\,\Omega(kt)$, so that $\{M,\tilde\Omega'(t)\}+[\tilde\Omega(t),\{M,\tilde\Omega(t)\}]=k^2f(kt)\,T$, and the final value is $\tilde\Omega(t_1/k)=k\,\Omega(t_1)=k\,\Omega_1$. Therefore, any angular velocity can be reached from $0$, if and only if any angular velocity in a small neighbourhood of $0$ can be reached from $0$.
(It looks like these two symmetries are related by $k=-1$, but we should keep a distinction between positive time and negative time.)
This has the form of a quadratic differential equation, $\mathbf x'(t)=\mathbf x(t)\odot\mathbf x(t)$ where $\odot$ is some bilinear function. (Specifically, for antisymmetric matrices $X$ and $Y$, define $X\odot Y=-\{M,\}^{-1}([X,\{M,Y\}])$, so Euler's equation is $\Omega'=\Omega\odot\Omega$, as long as $f=0$. Alternatively, define $X\odot Y=-[\{M,\}^{-1}(X),Y]$, so Euler's equation is $L'=L\odot L$.)
(I ought to move this power series stuff to a different page, since it's applicable far beyond this particular Question, and anyway it's bloating this page.)
Fix a norm on the space, and find some constant $\lVert\odot\rVert>0$ such that $\lVert\mathbf x\odot\mathbf y\rVert\leq\lVert\odot\rVert\lVert\mathbf x\rVert\lVert\mathbf y\rVert$ for all $\mathbf x,\mathbf y$.
Given initial value $\mathbf x(0)=\mathbf a$, the equation $\mathbf x'=\mathbf x\odot\mathbf x$ has the power series solution
$$\mathbf x(t)=\sum_{n=0}^\infty t^n\frac{\sum^{n!}\mathbf a^{n+1}}{\sum^{n!}1}$$
where the coefficient of $t^n$ is the average of all possible ways of evaluating $\mathbf a^{n+1}$ using the product $\odot$. For example, the coefficient of $t^3$ is $\tfrac16$ of
$$\sum^6\mathbf a^4=((\mathbf a\mathbf a)\mathbf a)\mathbf a+(\mathbf a(\mathbf a\mathbf a))\mathbf a+2(\mathbf a\mathbf a)(\mathbf a\mathbf a)+\mathbf a((\mathbf a\mathbf a)\mathbf a)+\mathbf a(\mathbf a(\mathbf a\mathbf a))$$ $$\newcommand{\aaaa}[3]{\mathbf a\underset{#1}\odot\mathbf a\underset{#2}\odot\mathbf a\underset{#3}\odot\mathbf a} =\aaaa{1}{2}{3}+\aaaa{2}{1}{3}\begin{matrix}{}+\aaaa{1}{3}{2} \\ {}+\aaaa{2}{3}{1}\end{matrix}+\aaaa{3}{1}{2}+\aaaa{3}{2}{1}.$$
The terms in the series can be obtained recursively:
$$\mathbf u_n=\frac{t^n}{n!}\sum^{n!}\mathbf a^{n+1}$$ $$=\frac tn\sum_{k=0}^{n-1}\mathbf u_{n-1-k}\odot\mathbf u_k.$$
We also have $\lVert\mathbf u_n\rVert\leq|t|^n\lVert\odot\rVert^n\lVert\mathbf a\rVert^{n+1}$, which ensures absolute convergence, for small $t$:
$$\sum_{n=0}^\infty\left\lVert t^n\frac{\sum^{n!}\mathbf a^{n+1}}{\sum^{n!}1}\right\rVert\leq\sum_{n=0}^\infty|t|^n\lVert\odot\rVert^n\lVert\mathbf a\rVert^{n+1}=\frac{\lVert\mathbf a\rVert}{1-|t|\lVert\odot\rVert\lVert\mathbf a\rVert},$$ $$|t|<\frac{1}{\lVert\odot\rVert\lVert\mathbf a\rVert}.$$
And at certain times (when the thruster is used) an impulse may be applied to $\mathbf x$, with a fixed direction $\mathbf b$ but an arbitrary magnitude $c>0$, thus: $\mathbf x(t+)=\mathbf x(t-)+c\mathbf b$. These impulses, and the time intervals between them, are many variables that we can control. If the number of variables is at least $6$ (or the dimension of the space $\mathbf x$ is in), then there is hope of surrounding $0$ in an open set, and thus reaching everywhere in the space (according to the previous section).
$$\mathbf x(0+)=0+c_0\mathbf b$$ $$\mathbf x(t_1-)=\sum_{n=0}^\infty\frac{t_1^n}{n!}\sum^{n!}\mathbf x(0+)^{n+1}$$ $$\mathbf x(t_1+)=\mathbf x(t_1-)+c_1\mathbf b$$ $$\mathbf x(t_1+t_2-)=\sum_{n=0}^\infty\frac{t_2^n}{n!}\sum^{n!}\mathbf x(t_1+)^{n+1}$$ $$\mathbf x(t_1+t_2+)=\mathbf x(t_1+t_2-)+c_2\mathbf b$$ $$\mathbf x(t_1+t_2+t_3-)=\sum_{n=0}^\infty\frac{t_3^n}{n!}\sum^{n!}\mathbf x(t_1+t_2+)^{n+1}$$ $$c_0,t_1,c_1,t_2,c_2,t_3\geq0$$ $$\mathbf x(t_1+t_2+t_3-)\approx0\quad?$$
