The axis of rotation is, as you’ve already determined, the line through the center of the two fixed spheres, which can be represented parametrically as $(1-\lambda)\mathbf A+\lambda\mathbf B$. This line has direction vector $\mathbf B-\mathbf A$, which is normalized to $\mathbf w=\|\mathbf B-\mathbf A\|$. When a point that’s not on this axis is rotated about the axis, it moves along a circular arc that lies in a plane perpendicular to the rotation axis and is centered on that axis.
To find this rotation center $\mathbf C$ for the center of the moving sphere $\mathbf X$, observe from the triangle of the sphere centers that $\mathbf C$ is the orthogonal projection of $\mathbf X$ onto $\overline{\mathbf A\mathbf B}$. Therefore, $$\mathbf C=\mathbf A+\mathbf w\cdot(\mathbf X-\mathbf A)\mathbf w.$$ Here, $\mathbf w\cdot(\mathbf X-\mathbf A)$ is the distance $d$ that you were solving for. This point could also be found by directly computing the intersection of the rotation axis and its perpendicular plane through $\mathbf X$. One way to do this is to solve the equation $$[(1-\lambda)\mathbf A+\lambda\mathbf B]\cdot\mathbf w=\mathbf X\cdot\mathbf w$$ for $\lambda$, which leads directly to the above result. (This equation comes from the point-normal form of equation for a plane.)
It seems useful to actually find the circle of motion. Let $\mathbf u=\|\mathbf X-\mathbf C\|$ and $\mathbf v=\mathbf w\times\mathbf u$. These two unit vectors are both parallel to the plane of motion and are orthogonal to each other. A parametric description of $\mathbf X$’s motion is then $\mathbf C+R_X\mathbf u\cos\theta+R_X\mathbf v\sin\theta$, a clockwise rotation through $\theta$ when viewed from $\mathbf A$. (I didn’t really need to normalize $\mathbf u$ for this, but I’m anticipating the next paragraph.)
With the center point $\mathbf C$ and the three unit vectors in hand, it’s fairly easy to construct a homogeneous rotation matrix that implements this rotation. In the $uvw$ coordinate system with origin at $\mathbf C$, it’s just a rotation about the $w$-axis, which can be represented in block form as $\small\begin{bmatrix}R_\theta&\mathbf0_2\\\mathbf0_2&I_2\end{bmatrix}$, where $R_\theta$ is the familiar two-dimensional rotation matrix and $\mathbf0_2$ is a $2\times2$ matrix of zeros. The corresponding matrix in the $xyz$ coordinate system is obtained via a change of basis: $$R=\begin{bmatrix}\mathbf u&\mathbf v&\mathbf w&\mathbf C \\ 0&0&0&1\end{bmatrix}\begin{bmatrix}R_\theta&\mathbf0_2\\\mathbf0_2&I_2\end{bmatrix}\begin{bmatrix}\mathbf u&\mathbf v&\mathbf w&\mathbf C \\ 0&0&0&1\end{bmatrix}^{-1}.$$ The change-of-basis matrix is a composition of a rotation and translation, so its inverse can be computed without going through a generic matrix inversion: $$R=\begin{bmatrix}\mathbf u&\mathbf v&\mathbf w&\mathbf C \\ 0&0&0&1\end{bmatrix}\begin{bmatrix}R_\theta&\mathbf0_2\\\mathbf0_2&I_2\end{bmatrix}\begin{bmatrix}\mathbf u^T&0\\\mathbf v^T&0\\\mathbf w^T&0\\\mathbf0^T&1\end{bmatrix}\begin{bmatrix}I_3&-\mathbf C\\\mathbf0&1\end{bmatrix}.$$ N.B.: in the first matrix, the unit vectors are columns, while in the third they are rows; $\mathbf C$ and $-\mathbf C$ are both entered as columns. Note, too, that any point on the rotation axis can be used instead of $\mathbf C$. You will end up with the same matrix. (If you like, expand the first matrix into its rotation and translation components to restore the symmetry of the product.)