Solution (source code)

= Solution

At $e_{\max}=\sqrt B$, torque balance gives
$$
\boxed{\sin\phi_0=-\frac{A(1+3B)}{Ca^2(p+q)B^{q/2}}
=-\frac{5A}{Ca^2(5p-q)B^{q/2}}.}
$$
The magnitude of the right-hand side must not exceed one. This fixes the dissipative offset of the <resonant-argument libration> centre.

For $C>0$, put $S=5A/[Ca^2(5p-q)B^{q/2}]$. The two torque-balanced phases are $\phi_0=-\arcsin S$ and $\phi_0=\pi+\arcsin S$, modulo $2\pi$. In the usual weak-perturbation pendulum approximation, $\delta\dot\phi\simeq-3(p+q)n\delta a/(2a)$ and $\delta\dot a\simeq-2(p+q)CaB^{q/2}\cos\phi_0\,\delta\phi$. Restoring libration thus selects $\cos\phi_0<0$, giving
$$
\boxed{\phi_0\simeq\pi+\arcsin S\quad(C>0).}
$$
It approaches $\pi$ when drag is weak compared with resonant torque. For $C<0$, the corresponding restoring branch has $\cos\phi_0>0$ and approaches zero. The sine equation is the direct result independent of this leading-order stability approximation.