simulations, we take $$\bm{u}(0) = (0,0,-1)^{t} \label{eq:sim-ic-1}$$ and specify $\bm{\omega}(0)$ as $$\bm{\omega}(0) = (|\omega_{xy0}|\cos\psi,\  |\omega_{xy0}|\sin\psi,\  -n_{i}) \label{eq:sim-ic-2}$$ in terms of $|\omega_{xy0}|$, $\psi$, and $n_{i}$. According to the simplified dynamics (\[eq:diff-ne-1\])–(\[eq:diff-ne-3\]), the irrelevant mode of oscillation does not affect the dynamics sensitively as long as the relevant mode exists and the initial spin energy is much larger than the initial oscillation energy. Thus we choose $\ket{\omega(0)}= (\omega_{x0}, \omega_{y0})^{t}$ in the direction of the relevant eigenmode, $$\psi = \psi_{p} \text{ for } n_{i}>0,\quad\text{and}\quad\psi = \psi_{r} \text{ for } n_{i}<0, \label{eq:sim-ic-3}$$ where $\psi_{p}$ and $\psi_{r}$ are the angles of the eigenvectors $\ket{\omega_{p}}$ and $\ket{\omega_{r}}$ from the $x$-axis, respectively.

Numerical results are presented in the unit system where $M$, $a$, and $$\tilde{t} \equiv 1/\tilde{\omega} \equiv \sqrt{a/g}$$ as units of mass, length, and time. The size of the time step for the numerical integration is taken to be $0.002\,\tilde{t}$. In numerics, we determine the time for reversal $t_{r}$ by the time at which $n=\bm{\omega}\cdot\bm{u}$ becomes zero for the first time, and they are compared with Garcia-Hubbard formulas (\[eq:trgh-p\]) and (\[eq:trgh-m\]); $n_0$ is determined as $$\frac{\gamma n_{0}^{2}}{2} = \frac{1}{2}(
    \alpha\omega_{x0}^2 + \beta\omega_{y0}^2 + \gamma\omega_{z0}^2),
    \label{eq:def-n0}$$ assuming $I_{\mathrm{eff}} = I_{zz}$ at $t=0$. Here the potential energy $U(\bm{u})$ 