The geodesics on a sphere are arcs, which are pieces of great circles. Great circles on a sphere may be parametrized by $\cos\theta$ and $\sin\theta$ in two coordinates and $0$ in others for some choice of orthonormal coordinates. Equivalently, by $(\cos\theta)\mathbf{a}+(\sin\theta)\mathbf{b}$ for perpendicular vectors $\mathbf{a}$ and $\mathbf{b}$.
For $S^3$ specifically, if $\mathbf{u}$ is a unit vector and $p$ is any versor (unit quaternion) then $p$ and $p\mathbf{u}$ are orthogonal with respect to the inner product, so $(\cos\theta)p+(\sin\theta)p\mathbf{u}=p\exp(\theta\mathbf{u})$ parametrizes a great circle. If we restrict the domain of $\theta$ to an interval $[0,\phi]$ we get an arc.
The angle between $q_0=p$ and $q_1=p\exp(\theta\mathbf{u})$ is $\theta$. Note we can write
$$ \exp(\theta\mathbf{u}) = q_0^{-1}q_1 $$
and solve $\theta=\|\ln(q_0^{-1}q_1)\|$ (assuming we define $\ln$ by the rule $\ln\exp(\mathbf{v})=\mathbf{v}$ for $\|\mathbf{v}\|<\pi$).
Note that given $q_0=p$ and $q_1$ arbitrary, it is always possible to write $q_1=p\exp(\theta\mathbf{u})$, as that is just a matter of expressing $q_0^{-1}q_1$ in polar form with a convex angle... unless of course $\theta=\pi$, so $q_0$ and $q_1$ are antipodal, in which case $\ln(-1)$ isn't really defined although you could argue $\|\ln(-1)\|$ should be defined as $\pi$ still.
I think this is basic enough to not need a source in a research article, only a swift explanation. Is there another reason you want a source for this fact?