This is a complete rewrite of this answer, following an approach which is very different from my first attempt.
Homogeneous parametrization
Suppose you have a point $x=(x_1,x_2,x_3)$ on a circle on the sphere. That means it belongs to the intersection of the sphere and a plane, or in formula, $\lVert x\rVert=1$ and $\langle x,n\rangle=d$. This $n$ is the normal vector of the plane, and $d$ encodes the distance from the origin. If you choose $n\in S^2$ then you can choose $d\in(-1,1)$, so this is your parametrization of circles on the sphere.
Now I'll homogenize things, and move to $\mathbb{RP}^3$. This will allow me to express the stereographic projection in a nice way. For homogeneous coordinates, we represent the unit sphere by the diagonal matrix $S=\operatorname{diag}(1,1,1,-1)$. A point with homogeneous coordinates $x$ lies on that sphere if it satisfies $x^TSx=0$. We also combine the parameters of the plane, writing $p=(n_1,n_2,n_3,-d)$ so a point $x$ lies on the plane $p$ if $\langle x,p\rangle=0$.
Stereographic projection of a point
We'll also need a stereographic projection. For that we'll embed the plane at $z=-1$ and use the point with homogeneous coordinates $N=(0,0,1,1)^T$ as the center of projection. At first we'll map from the plane to the sphere. So suppose now you have a point $(x,y)^T$ in the plane. In $\mathbb{RP}^3$ that point would have coordinates $P=(x,y,-1,1)^T$. To project that onto the sphere, you compute $Q=P^TSP\cdot N - 2N^TSP\cdot P$. This is a a point on the line $PN$ since it is a linear combination of these two. It's distinct from $N$ since it has a non-zero coefficient in front of $P$. And it's on the sphere as you can verify by computing $Q^TSQ$. So the image of $P$ on the sphere is $Q$. In coordinates:
$$P=\begin{pmatrix}x\\y\\-1\\1\end{pmatrix}\mapsto
\begin{pmatrix}4x\\4y\\x^2 + y^2 - 4\\x^2 + y^2 + 4\end{pmatrix}=Q$$
Stereographic projection of a circle
Next I'll project a circle with center $(x,y)$ and radius $r$ in the plane up to the sphere. One way to do this is by projecting three points, say $(x+r,y)$, $(x-r,y)$ and $(x,y+r)$. Projecting these onto the sphere gives three points in space, and three (non-collinear) points in space uniquely define a plane. In this case, the plane in question can be parametrized using
$$p=(4x, 4y, x^2 + y^2 - r^2 - 4, -x^2 - y^2 + r^2 - 4)^T \tag1$$
You can check this by ensuring $\langle Q,p\rangle=0$ for $Q$ being the image of each of the three points I just mentioned. In fact I found that vector by taking these three homogeneous coordinate vectors and computing all $3\times3$ minors, in suitable order and with alternating signs. Any multiple of the vector above will describe the same plane, since those coordinates are homogeneous. The representative above has the special property of $p_3+p_4=-8$. So if you start with any representative $p$, scaling that by $\frac{-8}{p_3+p_4}$ will lead to a representative of the form above.
This allows us to determine the description of the circle in the $z=-1$ plane from the parameters of the plane intersecting the sphere:
$$x=\frac{-2p_1}{p_3+p_4}\qquad y=\frac{-2p_2}{p_3+p_4}
\qquad x^2+y^2-r^2=\frac{-4(p_3-p_4)}{p_3+p_4}$$
Applying the Möbius transformation
Now I want to use this question of mine and this answer I received. One beautiful way of representing circles in the complex plane is using hermitian matrices of the form
$$\begin{pmatrix}x^2+y^2-r^2&x+iy\\x-iy&1\end{pmatrix}$$
or multiples thereof. That approach is completely different from the Apollonian circles you describe in your question, but for this application it's very useful since you can get one such multiple directly from the parameters of your plane:
$$C=\begin{pmatrix}4(p_3-p_4)&2(p_1+ip_2)\\2(p_1-ip_2)&-(p_3+p_4)\end{pmatrix}$$
To this you can apply your Möbius transformation
$$M=\begin{pmatrix}a&b\\c&d\end{pmatrix}$$
using conjugation, as detailed in the answer I mentioned:
$$C'=MCM^*$$
Here $M^*$ denotes the conjugate transpose of $M$.
Back to the sphere
Now we want to project that circle back to the sphere. We could normalize the matrix by dividing it by its lower right entry. Or we could take the matrix as it is, labeling its entries as follows:
$$C'=\begin{pmatrix}s&u+iv\\u-iv&t\end{pmatrix}$$
Then you have
$$x'=\frac{u}{t} \qquad y'=\frac{v}{t} \qquad x'^2+y'^2-r'^2=\frac{s}{t}$$
Plug these back into equation $(1)$:
$$p'\sim\left(4\frac ut, 4\frac vt, \frac st - 4, -\frac st - 4\right)^T$$
Multiply that by $t$ to avoid the divisions:
$$p'=(4u,4v,s-4t,-s-4t)^T$$
If you really want to, you can also normalize this in such a way that $p_1'^2+p_2'^2+p_3'^2=1$, i.e. the normal vector lies on the sphere. But I'd prefer not to do any divisions along the whole way.
Combining everything
If you wanted to, you could combine everything to a single big operation. The step from $p$ to $C$ is linear in $p$. The step from $C$ to $C'$ is quadratic in $M$, but still linear in $C$. The step from $C'$ to $p'$ is linear in $C'$ again. So the whole transformation from $p$ to $p'$ can be written as a $\mathbb R^{4\times 4}$ matrix, with coefficients which are quadratic in the real and imaginary parts of $a$ through $d$, the entries of your Möbius transformation:
$${\tiny\begin{pmatrix}
8b_rc_r+8b_ic_i+8a_rd_r+8a_id_i&
8b_ic_r-8b_rc_i-8a_id_r+8a_rd_i&
16a_rc_r+16a_ic_i-4b_rd_r-4b_id_i&
-16a_rc_r-16a_ic_i-4b_rd_r-4b_id_i\\
8b_ic_r-8b_rc_i+8a_id_r-8a_rd_i&
-8b_rc_r-8b_ic_i+8a_rd_r+8a_id_i&
16a_ic_r-16a_rc_i-4b_id_r+4b_rd_i&
-16a_ic_r+16a_rc_i-4b_id_r+4b_rd_i\\
4a_rb_r+4a_ib_i-16c_rd_r-16c_id_i&
-4a_ib_r+4a_rb_i+16c_id_r-16c_rd_i&
4\lvert a\rvert^2-\lvert b\rvert^2-16\lvert c\rvert^2+4\lvert d\rvert^2&
-4\lvert a\rvert^2-\lvert b\rvert^2+16\lvert c\rvert^2+4\lvert d\rvert^2\\
-4a_rb_r-4a_ib_i-16c_rd_r-16c_id_i&
4a_ib_r-4a_rb_i+16c_id_r-16c_rd_i&
-4\lvert a\rvert^2+\lvert b\rvert^2-16\lvert c\rvert^2+4\lvert d\rvert^2&
4\lvert a\rvert^2+\lvert b\rvert^2+16\lvert c\rvert^2+4\lvert d\rvert^2
\end{pmatrix}}$$
But perhaps this is more readable here if you don't write it as a matrix:
\begin{align*}
p_1'
&= (8b_rc_r+8b_ic_i+8a_rd_r+8a_id_i)p_1 \\
&+ (8b_ic_r-8b_rc_i-8a_id_r+8a_rd_i)p_2 \\
&+ (16a_rc_r+16a_ic_i-4b_rd_r-4b_id_i)p_3 \\
&+ (-16a_rc_r-16a_ic_i-4b_rd_r-4b_id_i)p_4 \\
p_2'
&= (8b_ic_r-8b_rc_i+8a_id_r-8a_rd_i)p_1 \\
&+ (-8b_rc_r-8b_ic_i+8a_rd_r+8a_id_i)p_2 \\
&+ (16a_ic_r-16a_rc_i-4b_id_r+4b_rd_i)p_3 \\
&+ (-16a_ic_r+16a_rc_i-4b_id_r+4b_rd_i)p_4 \\
p_3'
&= (4a_rb_r+4a_ib_i-16c_rd_r-16c_id_i)p_1 \\
&+ (-4a_ib_r+4a_rb_i+16c_id_r-16c_rd_i)p_2 \\
&+ (4\lvert a\rvert^2-\lvert b\rvert^2-16\lvert c\rvert^2+4\lvert d\rvert^2)p_3 \\
&+ (-4\lvert a\rvert^2-\lvert b\rvert^2+16\lvert c\rvert^2+4\lvert d\rvert^2)p_4 \\
p_4'
&= (-4a_rb_r-4a_ib_i-16c_rd_r-16c_id_i)p_1 \\
&+ (4a_ib_r-4a_rb_i+16c_id_r-16c_rd_i)p_2 \\
&+ (-4\lvert a\rvert^2+\lvert b\rvert^2-16\lvert c\rvert^2+4\lvert d\rvert^2)p_3 \\
&+ (4\lvert a\rvert^2+\lvert b\rvert^2+16\lvert c\rvert^2+4\lvert d\rvert^2)p_4
\end{align*}
where $a=a_r+ia_i, b=b_r+ib_i, c=c_r+ic_i, d=d_r+id_i$ are the entries of your Möbius transformation.