Finding all combinations of 2 vectors on the 2-sphere generating a point on the 2-radii 3-ball

Supose we have the 3-ball of radii 2, $\mathbb{B}^3_{r=2}=\{ \mathbf{x}\in\mathbb{R}^3: |\mathbf{x}| \leq 2 \}$, and the 2-sphere $\mathbb{S}^2=\{ \mathbf{x}\in\mathbb{R}^3: |\mathbf{x}| = 1 \}$. We know that $\{ \mathbf{w}=\mathbf{v} - \mathbf{u}: \mathbf{v}\in \mathbb{S}^2,\mathbf{u}\in \mathbb{S}^2 \} = \mathbb{B}^3_{r=2}$ thus every point on $\mathbb{B}^3_{r=2}$ can be decomposed into two point on $\mathbb{S}^2$.

Let us say we introduce a parameterization of $\mathbb{S}^2$, $\mathbf{k}(\theta,\phi) = \begin{pmatrix} \cos(\theta)\sin(\phi) \\ \sin(\theta)\sin(\phi) \\ \cos(\phi) \end{pmatrix}$ so that $\{\mathbf{k}(\theta,\phi) \in \mathbb{R}^3 \;\forall\; (\theta,\phi) \in [0,2\pi] \times [0,\pi]\} = \mathbb{S}^2$.

Given a point $\mathbf{b}_0 \in \mathbb{B}^3_{r=2}$, how do we find all combinations of $(\theta_1,\phi_1)$ and $(\theta_2,\phi_2)$ so that $\mathbf{k}(\theta_1,\phi_1) - \mathbf{k}(\theta_2,\phi_2) = \mathbf{b}_0$?

My first try was simply solving the equations as (3 equations, 4 variables, set 1 variable as free): $$\phi_2 = t$$ $$\cos(\phi_1) - \cos(\phi_2) = b_{0z} \Leftrightarrow \phi_1 = \cos^{-1}(b_{0z} + \cos(t))$$ $$\cos(\theta_1)\sin(\phi_1(t)) - \cos(\theta_2)\sin(t) = b_{0x} \Leftrightarrow \theta_1 = \cos^{-1} \left ( \frac{b_{0x} + \cos(\theta_2)\sin(t)}{\sin(\phi_1(t))} \right )$$ $$\sin(\theta_1)\sin(\phi_1(t)) - \sin(\theta_2)\sin(t) = b_{0y} \Leftrightarrow \\ \Leftrightarrow \sin(\cos^{-1} \left ( \frac{b_{0x} + \cos(\theta_2)\sin(t)}{\sin(\phi_1(t))} \right ))\sin(\phi_1(t)) - \sin(\theta_2)\sin(t) = b_{0y} \Leftrightarrow \\ \Leftrightarrow \sqrt{ 1 - \left ( \frac{b_{0x} + \cos(\theta_2)\sin(t)}{\sin(\phi_1(t))} \right )^2}\sin(\phi_1(t)) - \sin(\theta_2)\sin(t) = b_{0y} \Leftrightarrow \\ \Leftrightarrow \sin^2(\phi_1(t)) - \left ( \frac{b_{0x} + \cos(\theta_2)\sin(t)}{\sin(\phi_1(t))} \right )^2\sin^2(\phi_1(t)) = (b_{0y}+ \sin(\theta_2)\sin(t))^2 \Leftrightarrow \\ \Leftrightarrow ... \Leftrightarrow \\ \Leftrightarrow \theta_2 = \sin^{-1} \left (\frac{ \sin^2(\phi_1(t)) - b_{0y}^2 -b_{0x}^2 -\sin^2(t)}{2\sqrt{b_{0y}^2 + b_{0x}^2}\sin(t) } \right ) - atan2(b_{0y},b_{0x})$$

And then i have a closed form expression for $\theta_1(t),\phi_1(t),\theta_2(t),\phi_2(t)$ but when i try this formula numerically it does not work as intended. Can anyone identify either my mistake or give a path (reference to a solution or the solution itself) to the answer?

Surly this question has been solved somewhere before but i cannot find it in any of the books on my shelf nor trough google :P

The solution was:

Rotate $\mathbf{b}_0$ to lie along the z-axis, using its spherical angles $(\theta_b,\phi_b)$, then use the triangle formed in the common plane of the vectors $\mathbf{k}_1,\mathbf{k}_2$.

The triangle will have sides $|\mathbf{b}_0|,1,1$, thus from the law of cosines we know that the zenith angle of the first spherical vector $\mathbf{k}_1$ is

$$\phi_1 = \cos^{-1} \left ( \frac{|\mathbf{b}_0|}{2} \right )$$

We then find the corresponding angle for $-\mathbf{k}_2$ again trough the law of cosines and the triangle,

$$\phi_2 = \pi - \cos^{-1} \left ( \frac{|\mathbf{b}_0|}{2} \right ) \\ \theta_2 = \pi + (\theta_1+\pi) = \theta_1$$

The extra $\pi$ term comes into play as we are looking for the spherical coordinate for $\mathbf{k}_2$ but have used $-\mathbf{k}_2$ in the triangle.

We also see that the only constraint on the common plane is that it contains the $z$-axis. Thus any arbitrary rotation around this axis will produce a valid "triangle" and so our rotated parametric solution is:

$$\theta_1(t) = t, \phi_1(t) = \cos^{-1} \left ( \frac{|\mathbf{b}_0|}{2} \right ) , \\ \theta_2(t) = t, \phi_2(t) =\pi - \cos^{-1} \left ( \frac{|\mathbf{b}_0|}{2} \right )$$

These angles are however only for a vector around the z-axis. To rotate back to the original $\mathbf{b}_0$ we use the Bloch Sphere rotation.

We do this by mapping our spherical coordinates to $\mathbb{C}^2$ by

$$Q(\theta, \phi) = \begin{pmatrix} \cos\left( \frac{\phi}{2} \right ) \\ e^{i\theta}\sin \left( \frac{\phi}{2} \right )\end{pmatrix}.$$

Where the inverse is

$$Q^{-1}\left( \begin{pmatrix} z_1 \\ z_2 \end{pmatrix}\right) = \left( 2\textrm{atan} \frac{|z_2|}{|z_1|}, \textrm{arg}(z_2)-\textrm{arg}(z_1)\right).$$

Then a rotation of $\varphi$ around a axis in the Cartesian coordinates rotates the complex vector by a matrix $R_n(\varphi)$. So to rotate back to $(\theta_b,\phi_b)$ we form

$$(\theta(t),\phi(t)) = Q^{-1}\left( R_z(\theta_b)R_y(\phi_b) Q(\theta'(t),\phi'(t)) \right)$$

where

$$R_y(\varphi) = \begin{pmatrix} \cos\frac{\varphi}{2} & -\sin\frac{\varphi}{2} \\ \sin\frac{\varphi}{2} & \cos\frac{\varphi}{2} \end{pmatrix}, \\ R_z(\varphi) = \begin{pmatrix} e^{-i\frac{\varphi}{2}} & 0 \\ 0 & e^{i\frac{\varphi}{2}} \end{pmatrix}.$$

Given a vector $\mathbf{b}_0 \in \mathbb{B}^3_{r=2}$, all decomposition's into the difference between two vectors from the 2-sphere is given by
$$(\theta_1(t),\phi_1(t)) = Q^{-1}\left( R_z(\theta_b)R_y(\phi_b) Q(\theta_1'(t),\phi_1'(t)) \right), \\ (\theta_2(t),\phi_2(t)) = Q^{-1}\left( R_z(\theta_b)R_y(\phi_b) Q(\theta_2'(t),\phi_2'(t)) \right), \\ \forall\; t\in[0,2\pi]$$
$$\theta_1'(t) = t, \\ \phi_1'(t) = \cos^{-1} \left ( \frac{|\mathbf{b}_0|}{2} \right ) , \\ \theta_2'(t) = t, \\ \phi_2'(t) =\pi - \cos^{-1} \left ( \frac{|\mathbf{b}_0|}{2} \right ), \\ \theta_b = \mathrm{atan2}(b_{0y},b_{0x}), \\ \theta_b = \cos^{-1}(\frac{b_{0z}}{|\mathbf{b}_0|}),$$