(Copiedfrom MSE. Offering four bounties over time, I got no response, other than twenty-nine upvotes.)

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 to Could this three-thruster spherical spacecraft de-tumble itself with zero final velocity? on space.SE.

For any two matrices and , define the commutator , the anticommutator , and the Frobenius inner product , where is the trace and is the transpose.

Let be a symmetric positive-definite matrix, and an antisymmetric matrix. ( describes the distribution of mass in the spacecraft, is a vector locating the thruster relative to the centre of mass, is the force produced by the thruster, and is the torque produced by the thruster. More details at Chiral Anomaly's answer to Rigid body dynamics derivation from Newton's laws for higher dimensions. Everything is described in the rotating reference frame.)

The angular velocity , an antisymmetric matrix, changes with time according to Euler's equation

where is a function (continuous, piecewise-constant, or just integrable) describing when and how strongly the thruster is used.

Question: Can and be chosen such that, for any two antisymmetric matrices and , there exist and such that the solution to Euler's equation with initial value has final value ?

To simplify things, we may take to be a combination of Dirac deltas instead of an ordinary function. Then the angular velocity satisfies except at a finite set of times when changes by a positive multiple of .

Here is Euler's equation in terms of the components of the angular velocity. Assume that is diagonal, with components .

(and permute the indices to get 6 equations like this).

If , then has constant sign; is either always non-increasing, or always non-decreasing. So, if the difference between initial and final values of has the wrong sign compared to , then there is no solution. Thus, we must take , and similarly for .

The Frobenius inner product of two antisymmetric matrices has every term duplicated: (since and ). So it's natural to take half of this.

The angular momentum is . Its squared magnitude is ; the derivative of this is . So the angular momentum has constant magnitude whenever . (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 ; the derivative of this is . So the energy is constant whenever .

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:

Doubled energy:

Halved exterior square of angular momentum:

(Its derivative is the exterior product of angular momentum and torque: .)

These three constants define two ellipsoids and another quadratic surface (with signature ) in the 6D space of antisymmetric matrices. When , must remain on the intersection of these three quadratic surfaces. In terms of instead of (with ), those constants are

We can assume . Taking the squared magnitude of , and subtracting or times the doubled energy, we get

both being constant. Notice the signs of the coefficients; the first expression has only one positive term, so (if the expression is ) it defines a cone, or a pair of opposite cones, around the axis in 6D. And the second expression defines a pair of opposite cones around the axis. If the torque is inside one of these cones, and the initial angular momentum is in the same cone, then stays in that cone, and there's no solution with the final angular momentum outside of the cone. Thus, we must take to be outside of all four of these cones. That is, the first expression (with in place of ) must be negative, and the second expression must be positive.

Notice that Euler's equation has a time symmetry: Given a solution , we can construct another solution , so that , and the initial and final values are swapped (and negated): and . Therefore, any angular velocity can be reached from any other, if and only if any angular velocity can be reached from . (We can reverse a path from to to get a path from to , and concatenate that with a path from to , to get a path from to .)

Also, Euler's equation has a scale symmetry: Given a solution , for any we can construct another solution , so that , and the final value is . Therefore, any angular velocity can be reached from , if and only if any angular velocity in a small neighbourhood of can be reached from .

(It looks like these two symmetries are related by , but we should keep a distinction between positive time and negative time.)

This has the form of a quadratic differential equation, where is some bilinear function. (Specifically, for antisymmetric matrices and , define , so Euler's equation is , as long as . Alternatively, define , so Euler's equation is .)

Fix a norm on the space, and find some constant such that for all .

Given initial value , the equation has the power series solution

where the coefficient of is the average of all possible ways of evaluating using the product . For example, the coefficient of is of

(In the -dimensional algebra, the product is associative, and the power series simplifies to .)

This converges absolutely, for small :

And at certain times (when the thruster is used) an impulse may be applied to , with a fixed direction but an arbitrary magnitude , thus: . These impulses, and the time intervals between them, are many variables that we can control. If the number of variables is at least (or the dimension of the space is in), then there is hope of surrounding in an open set, and thus reaching everywhere in the space (according to the previous section).