  Celestial bodies approximated with rigid triaxial ellipsoids in a two-body
system can rotate chaotically due to the time-varying gravitational torque from
the central mass. At small orbital eccentricity values, rotation is short-term
orderly and predictable within the commensurate spin-orbit resonances, while at
eccentricity approaching unity, chaos completely takes over. Here, we present
the full 3D rotational equations of motion around all three principle axes for
triaxial minor planets and two independent methods of numerical solution based
on Euler rotations and quaternion algebra. The domains of chaotic rotation are
numerically investigated over the entire range of eccentricity with a
combination of trial integrations of Euler's equations of motion and the
GALI($k$) method. We quantify the dependence of the order--chaos boundaries on
shape by changing a prolateness parameter, and find that the main 1:1
spin-orbit resonance disappears for specific moderately prolate shapes already
at eccentricities as low as 0.3. The island of short-term stability around the
main 1:1 resonance shrinks with increasing eccentricity at a fixed low degree
of prolateness and completely vanishes at approximately 0.8. This island is
also encroached by chaos on longer time scales indicating longer Lyapunov
exponents. Trajectories in the close vicinity of the 3:2 spin-orbit resonance
become chaotic at smaller eccentricities, but separated enclaves of orderly
rotation emerge at eccentricities as high as 0.8. Initial perturbations of
rotational velocity in latitude away from the exact equilibrium result in a
spectrum of free libration, nutation, and polar wander, which is not well
matched by the linearized analysis omitting the inertial terms.
