The traditional (linearly) damped spring is modeled via a standard second-order differential equation: $\ddot{\bf x}(t)+a\dot{\bf x}(t)+b{\bf x}(t)=0$; this can be 'recentered' by letting ${\bf x}(t)={\bf y}(t)-{\bf y_0}$. With appropriate coefficients (particularly for critical damping, where the equation takes on a doubled eigenvalue), this can be used as an 'ease-in' to a specific target position $\bf y_0$.
The same notion — easing-in to a target value — obviously makes sense on $SO(3)$, but the equation itself doesn't: if one takes e.g. the usual quaternionic embedding in $\mathbb{R}^4$ then this equation doesn't guarantee uniticity of $\mathbb{x}(t)$, and if instead one tries to work on the sphere proper, it's not clear to me how one can take what amounts to a 'second-order tangent bundle' here. What would the appropriate mathematical representation of a second-order differential equation like this on the sphere (or presumably on any other manifold, though $SO(3)$ and particularly its quaternion representation is the specific case I'm interested in from a computational standpoint) be?
Not sure if this is exactly what you were thinking of, but here's my take:
If we interpret your equation $\ddot x = - a \dot x - b(x - x_0)$ as a force law, then in some sense the closest equation that stays inside $S^3 \subset\mathbb R^4$ should be $$\ddot x = -a \dot x -b (x- x_0) + \lambda(x,\dot x) x\label{1}\tag{1}$$ where $\lambda$ is something like a Lagrange multiplier keeping the solution in $S^3$. Since $$\frac d {dt} |x|^2 = 2x\cdot\dot x,$$ in order to preserve $|x| = 1$ we must preserve $2x \cdot \dot x = 0$, which (differentiating again) is equivalent to $$x\cdot \ddot x = -|\dot x|^2.$$
Dotting $\eqref{1}$ with $x$, this becomes $$-|\dot x|^2 = -a |\dot x|^2 + \lambda(x,\dot x)-b+b x\cdot x_0;$$ so we can solve for $\lambda = (1-a)|\dot x|^2 + b - b x \cdot x_0.$
Substituting back in to $\eqref{1}$ yields the evolution equation $$\ddot x = -a \dot x +(1-a) |\dot x|^2x + b (x_0 - (x \cdot x_0)x).$$
The first two terms can be considered as a modification of the geodesic equation: when $a=0,b=0$ this really is just the geodesic equation for $S^3$, while for $a>0$ it's something like a geodesic equation with damping.
The term $B(x) = b(x_0 - (x \cdot x_0) x)$ is the projection of $x_0$ on to the orthogonal complement of $x$; i.e. the closest vector to $x_0$ that is contained in the tangent space $T_x S^3$. Thus this term really does act as a driving force towards/away from $x_0$, and has the nice property of being tangent to the sphere; i.e. $B|_{S^3}$ is a section of $TS^3$.
Our calculations show that $x$ will stay in $S^3$ so long as the initial data satisfies $x \cdot \dot x = 0$. Since this is true for $(x,\dot x) \in TS^3$, this equation is a second-order ODE on $S^3$ proper, which can be written intrinsically using the covariant acceleration $\nabla_{\dot x} \dot x$ as
$$\nabla_{\dot x} \dot x=-a \dot x + B(x).$$
(Here the connection $\nabla$ is the Levi-Civita connection of the Riemannian structure on $S^3$ inherited from $\mathbb R^4$, so the covariant acceleration is just the projection of $\ddot x$ on to the tangent space of $S^3$.)
Unfortunately, it seems like this doesn't descend to $SO(3) = S^3 /\{ \pm 1 \}$ unless $b = 0$: quaternions $\pm x$ representing the same rotation would fail to do so at later times. You could fix this by replacing $B$ with any vector field on $S^3$ invariant under the antipodal map - I'm not sure what a physically natural choice would be, but a simple way to achieve this is by multiplying our current $B(x)$ by $x \cdot x_0$.