Spherical tools
At first I understood definition of ruler and compass as operations not in space, but instead operations on the surface of the sphere. So these would spherical tools, not spatial, which is a good starting point.
Using spherical circles, you can construct the perpendicular bisector of the segments $AA'$ and $BB'$, just as you would in the plane, by intersecting circles of equal radius around the endpoints. This bisector is a greatcircle which is the locus of all points at equal distance from $A$ and $A'$ resp. $B$ and $B'$. Note that this great circle will intersect the great circle through $A$ and $A'$ in two points, so there are two antipodal “midpoints” to the segment, not only one. But since we are interested in the bisector itself, we don't have to worry about midpoints.
Now that you have these two bisectors, look at their intersections. They as well will intersect in two antipodal points, which are the points where your axis of rotation intersects the sphere. And the angle of rotation is the angler under which these bisectors intersect.
Spatial tools
So what if your ruler and compass really operate in space? Drawing great circles through two given points, the ruler operation on the sphere, would become a compass operation in space, drawing a circle around the center of the sphere through the two given points.
The easiest way to construct the perpendicular bisectors would involve intersecting three spheres: the sphere on which the four points live and two spheres of equal radius around $A$ and $A'$. But your definition does not involve a tool for this.
But at this point I feel like an extension of your toolbox would be required in any case. Your definition of a compass operation is too tight. For example, it doesn't even cater for the case of drawing two circles of equal radius. Furthermore, one often wants to draw circles around a given point but with arbitrary radius; with your definition that would be impossible because I'd already have to know to points on the preimeter of the circle, i.e. two points equal distance away from the center. One possible generalization here would be a circle with a given center, a given point on the perimeter and a line through the center which the circle will intersect. Conceptually you'd construct a sphere around the center through the given perimeter point, and where that sphere intersects the line would be your second perimeter point for the circle operation as you specified it.
Extended spatial compass
So for now I'll assume that there are two distinct compass operations. In both cases you give the center and a line through that center which the circle has to intersect. In the first case you also give a point on the perimeter, whereas in the second case you give the radius (taken as a distance between two other points) and a second line the circle will intersect as well. So in both cases you fix the center, the plane and the radius, without overspecification.
Using these tools, you can work out the plane of the perpendicular bisector of $A$ and $A'$, and then draw that bisector.
- Connect $A$ and $A'$ with a line $a$
- Draw a circle around center $A$ through line $a$ and an arbitrary point $P$ in space, not on $a$ and at a distance form $A$ which is more than half the distance $AA'$
- Draw a circle around $A'$ with radius $AP$ through lines $a$ and $AP$
- Both circles are in the plane $AA'P$, and have sufficiently large radii, so they will intersect in two points $Q_1$ and $Q_2$
- Construct a circle around the center $O$ of the sphere, passing through both lines $OQ_1$ and $OQ_2$ and with the radius of the sphere, i.e. $OA$
- Perform steps $1$ through $5$ for $B,B'$ instead of $A,A'$
This way you can construct the perpendicular bisectors with tools very similar to what you originally asked for.
I'm not entirely sure if you're asking for a closed expression that gives the points, or just a description of these points and am going to assume the latter.
The first step in solving this problem is to change it into an easier one. We may rotate the sphere in such a way the point $P$ becomes the north pole $N=(0,0,1)$. If we call this rotation (which may be expressed as a $3\times3$ matrix) $R$, a point $Q$ is on the circle with radius $r$ around $P$ iff $R(Q)$ is on the circle with radius $r$ around $N$. This is due to the fact that a rotation of the sphere does not change distances on the sphere.
The next step is to consider these particular circles of radius $r$ around $N$. The best way to see this is to draw a sketch of the intersection of the sphere with the $yz$-plane, i.e. a circle of radius 1 around the origin. We can also see to points $P_1, P_2$ on the circle at distance $r$ from $N$. (I'll assume $r$ is small enough for this to make sense) These points should be such that $P_1$ is obtained from $P_2$ by reflecting around the $z$-axis. The spherical coordinate $\theta$ corresponding to these points can now easily be obtained: this angle is given as an angle wrt the $z$-axis, while $r$ is precisely the arc length corresponding to this angle on a circle with radius 1. This leads to the conclusion that $r=\theta$ in this special case.
The next thing to notice, is that all points at distance $r$ from $N$ have the same $z$-coordinate. By projecting either $P_1$ or $P_2$ at this $z$-coordinate (which again can be done insightful in the same picture), one finds this height to be $z_0=\cos(\theta)=\cos(r)$. Moreover, yet again from considering the same picture, the required circle is a circle parallel to the $xy$-plane with center $(0,0,z_0)$ and radius $\sin(\theta)=\sin(r)$.
We may thus conclude that the points at distance $r$ from $N$ are of the form $(\sin(r)\cos(\phi),\sin(r)\sin(\phi),cos(r))$ with $\phi\in[0,2\pi)$. Using the rotation matrix $R$ from before, we thus get the set of points satisfying your original question:
$$\{R^{-1}(\sin(r)\cos(\phi),\sin(r)\sin(\phi),cos(r))\mid \phi\in[0,2\pi)\}.$$
Best Answer
Remember that a rotation matrix will only give a rotation around an axis that passes through the origin. So you cannot just apply a rotation matrix.
First subtract $(x_0,y_0,z_0)$ from the coordinates of any point you want to rotate. Then apply the rotation matrix. Then add $(x_0,y_0,z_0)$ to the rotated coordinates.
With that in mind, the rotation will actually be performed on a sphere whose center is at $(0,0,0)$ (after subtracting $(x_0,y_0,z_0)$ and before adding $(x_0,y_0,z_0)$ back again).
The simplest rotation, conceptually, is to simply rotate the circle about the axis through the two points where the two circles should intersect. This axis is perpendicular to the vector $(x_1-x_0,y_1-y_0,z_1-z_0),$ and also perpendicular to the projection of that vector onto the $x,y$ plane. The two intersection points are therefore $\left(r \cos\left(\sigma + \frac\pi2\right), r \sin\left(\sigma + \frac\pi2\right), 0\right)$ and $\left(r \cos\left(\sigma - \frac\pi2\right), r \sin\left(\sigma - \frac\pi2\right), 0\right),$ where $r$ is the sphere's radius ("ray" in your formulas). Slightly simpler formulas for these points (using the trig identities for angles plus or minus $\frac\pi2$ radians) are $(- r \sin(\sigma), r \cos(\sigma), 0)$ and $(r \sin(\sigma), -r \cos(\sigma), 0)$.
When describing an axis of rotation we typically just take a unit vector along the axis. The factor of $r$ makes the vector that much longer or shorter, so let's just use $(-\sin(\sigma), \cos(\sigma), 0).$ I chose this vector because it's the one on the "far side" of the sphere in your diagram and a right-hand-rule rotation by $\theta$ will give the result you want. If we used the other vector as an axis we'd have to rotate by $-\theta$ or us a left-hand rule.
From Wikipedia, the matrix for a rotation by angle $\theta$ around axis $(u_x,u_y,u_z)$ is $$\begin{bmatrix}\cos \theta +u_{x}^{2}\left(1-\cos \theta \right)&u_{x}u_{y}\left(1-\cos \theta \right)-u_{z}\sin \theta &u_{x}u_{z}\left(1-\cos \theta \right)+u_{y}\sin \theta \\u_{y}u_{x}\left(1-\cos \theta \right)+u_{z}\sin \theta &\cos \theta +u_{y}^{2}\left(1-\cos \theta \right)&u_{y}u_{z}\left(1-\cos \theta \right)-u_{x}\sin \theta \\u_{z}u_{x}\left(1-\cos \theta \right)-u_{y}\sin \theta &u_{z}u_{y}\left(1-\cos \theta \right)+u_{x}\sin \theta &\cos \theta +u_{z}^{2}\left(1-\cos \theta \right)\end{bmatrix}.$$
We want $u_x = -\sin(\sigma),$ $u_y = \cos(\sigma),$ and $u_z = 0.$ Plugging these into the matrix formula, we get $$\begin{bmatrix} \cos\theta + (1-\cos \theta )\sin ^2\sigma & -(1-\cos \theta ) \sin\sigma \cos\sigma & \sin \theta \cos\sigma \\ -\left(1-\cos \theta \right)\sin\sigma \cos\sigma & \cos \theta + (1-\cos \theta)\cos^2\sigma & \sin \theta \sin\sigma \\ -\sin \theta\cos\sigma & -\sin \theta\sin\sigma & \cos \theta \end{bmatrix}.$$
So that's your rotation matrix.
But note that either before or after you apply this rotation, you could rotate the circle in place and it would still be the same circle. When you combine that rotation with the rotation matrix above, you get a rotation about yet another axis which also takes the dotted-line circle to the solid circle.
For example, you could rotate by $\pi$ radians ($180$ degrees) around the vector $\left(\sin\left(\frac\theta2\right)\cos(\sigma), \sin\left(\frac\theta2\right)\sin(\sigma), \cos\left(\frac\theta2\right)\right)$, which is the bisector of the angle between the $z$ axis and $(x_1-x_0,y_1-y_0,z_1-z_0).$ Since we are rotating by $\pi$ radians rather than the angle $\theta,$ we have to substitute $-1$ where the Wikipedia formula has $\cos\theta$ and $0$ where Wikipedia has $\sin\theta.$ That gives a matrix $$\begin{bmatrix} 2u_{x}^{2} - 1 & 2u_{x}u_{y} & 2u_{x}u_{z} \\ 2u_{y}u_{x} & 2u_{y}^{2} - 1 & 2u_{y}u_{z} \\ 2u_{z}u_{x} & 2u_{z}u_{y} & 2u_{z}^{2} - 1 \end{bmatrix}.$$
Set $u_x = \sin\left(\frac\theta2\right)\cos(\sigma),$ $u_y = \sin\left(\frac\theta2\right)\sin(\sigma),$ and $u_z = \cos\left(\frac\theta2\right).$ Also consider the trig identities $2\sin^2\left(\frac\theta2\right) = 1 - \cos\theta$, $2\cos^2\left(\frac\theta2\right) = 1 + \cos\theta,$ and $2\sin\left(\frac\theta2\right) \cos\left(\frac\theta2\right) = \sin\theta,$ we get $$\begin{bmatrix} (1 - \cos\theta)\cos^2\sigma - 1 & (1 - \cos\theta)\sin\sigma\cos\sigma & \sin\theta\cos\sigma \\ (1 - \cos\theta)\sin\sigma\cos\sigma & (1 - \cos\theta)\sin^2\sigma - 1 & \sin\theta\sin\sigma \\ \sin\theta\cos\sigma & \sin\theta\sin\sigma & \cos\theta \end{bmatrix}.$$
A little more manipulation would show that this is very similar to the rotation matrix derived earlier, except that the signs of all entries in the first two columns have been reversed. This is because each individual point of the circle ends up $180$ degrees around the circle from where it does under the other rotation.