Find the rotation matrix that relates two skew lines

I came up with the following problem that could be interesting to some readers.

You are given two skew lines, specified in parametric form as,

$P_1(t) = a_1 + t b_1$

$P_2(t) = a_2 + t b_2$

where $ b_1 $ and $b_2$ are unit vectors. Determine the axis of rotation and rotation matrix that relates the two skew lines, i.e. specify the rotation of the first line that sends it into the second line.

To limit the number of solutions, it is required that the image of $a_1$ be $a_2$. Here's what I have tried.

The image of $a_1$ is $a_2$ and the image of $a_1 + b_1$ is $a_2 + b_2$, so we have two pairs of points. Since the axis of rotation that revolves $a_1$ into $a_2$ must lie in the perpendicular bisecting plane of the segment $a_1 a_2$, we thus have two planes in which the axis lies, so we should be able to find the axis by solving for the intersection of the two planes. This construction gives a unique axis and a unique rotation matrix.

On the other hand, another possibility pointed out in the comments below, is to make the rotation a composition of two rotations, the first rotating $a_1$ into $a_2$ and $b_1$ into $c_1$, then a second rotation that rotates about the point $a_2$, and rotates $c_1$ into $b_2$. This two-step rotation results in an infinite number of rotation matrices.

I don't understand why this is so, although numerical verification confirms that each is a valid rotation. Why did I get a unique rotation matrix using the first method, and an infinite number of possible rotations using the second two-step method ?

EDIT:

I have found out after some thought that the two-step (composition of two rotations) method results in an affine transformation that is not necessarily a pure rotation. And so if we require the affine transformation, that maps $a_1$ into $a_2$ and the first line into the second, to be a pure rotation then this rotation is unique, as determined by the first single-step method. I would appreciate it if someone could shed some light on that.

Given any two distinct skew lines in Euclidean space, they uniquely fix the points on each line which are at the closest distance to the other line. These two points determine a unique line perpendicular to both skew lines and joining them. Consider the point midway between these two points. There is a unique plane containing this point which is perpendicular to the joining line. Project both skew lines onto this plane. There are exactly two lines in this plane which are equidistant from the projected lines and they are perpendicular to each other. Each of these lines determine a rotation axis where the rotation by a half revolution interchanges the two skew lines. Furthermore, if these two lines are translated parallel to themselves elsewhere in the plane they also determine a rotation axis that interchanges the two skew lines but now with rotation angles less than half of a revolution.

This may not be completely correct. I am working on actual implementation with code.

Assuming that $a_1, a_2$ are respectively the closest points in $P_1(t), P_2(t)$ we choose $a_0 = \frac 12(a_1+a_2)$. Now according to the sign of $\sigma=\vec b_1\cdot\vec b_2$ we choose

$$ \cases{ -1\le\sigma<0\to \vec b_0=\frac{\vec b_1-\vec b_2}{\|\vec b_1-\vec b_2\|}\\ 0<\sigma\le1\to \vec b_0=\frac{\vec b_1+\vec b_2}{\|\vec b_1+\vec b_2\|} } $$

after that, defining a rotation $R(\vec v_0,\theta)$ as a rotation of $\theta$ radians around the axis $\vec b_0$ we have

$$ R(\vec v_0,\pi)(a_1-a_0)+t R(\vec v_0,\pi)\vec b_1 + a_0 = a_2 + t \vec b_2 $$

or

$$ \cases{ R(\vec v_0,\pi)(a_1-a_0)+a_0 = a_2\\ R(\vec v_0,\pi)\vec b_1 = \vec b_2 } $$

NOTE

Here $R(\vec v_0,\theta)=I_3+\sin\theta V_0 +(1-\cos\theta)V_0^2$ (Rodrigues rotation formula) and $V_0=\left(\begin{array}{ccc} 0 & -v_{0z} & v_{0y}\\ v_{0z}& 0 & -v_{0x}\\ -v_{0y} & v_{0x} & 0\end{array}\right)$

EDIT

To obtain the closest points $a_1^*, a_2^*$ given the two skew lines $p_1=a_1+t_1\vec b_1, p_2 = a_2+t_2 \vec b_2$ we formulate the generic distance between the two lines as follows

$$ d^2(t_1,t_2) = \|p_1-p_2\|^2 = \|a_1-a_2\|^2+t_1^2+t_2^2+2(a_1-a_2)\cdot\vec b_1 t_1-2(a_1-a_2)\cdot\vec b_2 t_2 -2\vec b_1\cdot\vec b_2 t_1t_2 $$

so the stationary point due to minimum distance is determined by solving the linear system

$$ \left(\begin{array}{cc} 1 & -\vec b_1\cdot\vec b_2\\ -\vec b_1\cdot\vec b_2 & 1 \end{array}\right)\left(\begin{array}{c} t_1\\ t_2 \end{array}\right) = \left(\begin{array}{c} -(a_1-a_2)\cdot \vec b_1\\ (a_1-a_2)\cdot \vec b_2 \end{array}\right) $$

thus obtaining

$$ \cases{ t_1^* = \frac{1}{1-(\vec b_1\cdot\vec b_2)^2}\left(-(a_1-a_2)\cdot\vec b_1 +\vec b_1\cdot\vec b_2(a_1-a_2)\cdot b_2\right)\\ t_2^* = \frac{1}{1-(\vec b_1\cdot\vec b_2)^2}\left(\vec b_1\cdot\vec b_2(a_1-a_2)\cdot b_1+(a_1-a_2)\vec b_2\right) } $$

and then

$$ \cases{ a_1^* = a_1 + t_1^*\vec b_1\\ a_2^* = a_2 + t_2^*\vec b_2 } $$