Short answer
First the actual answer, for the impatient. The second point of intersection $q$ can be computed as
$$\begin{pmatrix}q_x\\q_y\end{pmatrix} =
\begin{pmatrix}p_x\\p_y\end{pmatrix} +
\tau\begin{pmatrix}y_2-y_1\\x_1-x_2\end{pmatrix}$$
with
$$\tau=2\frac{x_p (y_1 - y_2) - x_2 y_1 + x_1 y_2 - (x_1 - x_2) y_p}
{(x_1 - x_2)^2 + (y_1 - y_2)^2}\;.$$
The line joining the points of intersection is perpendicular to the line joining the centers of the circles, so that direction is described by that last vector in the sum. Interestingly, this formula does not make use of the radii $r_1$ and $r_2$ at all. That makes the problem more consistent: no matter what input you plug into this formula, you'll get the second point of intersection for the two circles around the given centers with radii chosen such that they go through $p$. Or in other words, it's the reflection of $p$ in the line that connects the circle centers.
Long derivation
Actually when coming up with the above, I started with the more general setup of Möbius geometry, as already mentioned in the question. So while the short answer above should be enough for anyone dealing with “real” circles, this answer here may be more useful when dealing with situations where circles may degenerate into lines.
If you write circles as $c_i=[x_i:y_i:x_i^2+y_i^2-r_i^2:1]$ or some multiple thereof, then two circles are perpendicular (i.e. have perpendicular tangents at their points of intersection) iff $c_i^T\cdot B\cdot c_j=0$ with
$$B=\begin{pmatrix}-2&&&\\&-2&&\\&&0&1\\&&1&0\end{pmatrix}$$
where $B$ is the matrix of the bilinear form underlying the Möbius geometry. Of course any multiple of this matrix will work just as well. And if you choose a different basis for your circle representations, it is even possible to make this matrix diagonal with only entries $\pm1$ on the diagonal. So we are essentially operating in some Minkowsky space $\mathbb R^{3,1}$ here. But all of this doesn't have to concern you. Suffice it to say that a point would be a circle of radius zero, so your point $p$ has the vector $[x_p:y_p:x_p^2+y_p^2:1]$. A point is perpendicular to a circle if it lies on that circle. Which means a point is also perpendicular to itself; the set of all points can be characterized as $\{p\mid p^TBp=0\}$.
With this knowledge we can formulate that we're looking for a point $q$ which satisfies
$$c_1^TBq = c_2^TBq = q^TBq = 0\quad q\not\in[p]$$
The first two conditions restrict the original four-dimensional space of coordinate vectors to a two-dimensional linear subspace, namely the space of all circles perpendicular to both the given circles. That space can be expressed using a basis, i.e. as the set of all alinear combinations of two circle vectors satisfying these conditions. One of these vectors would be the above representation of $p$ as a point of radius zero, but we still need a second circle with these properties. One simple choice would be the line that intersects both circles perpendicularly, i.e. the line joining the circle centers. You can find that line $h$ by solving
$$\begin{pmatrix}
-2x_1 & -2y_1 & 1 & x_1^2+y_1^2-r_1^2 \\
-2x_2 & -2y_2 & 1 & x_2^2+y_2^2-r_2^2 \\
0 & 0 & 0 & 1
\end{pmatrix}\cdot h=0$$
The rows of this matrix are $c_1^TB$, $c_2^TB$ and $P_\infty^TB$ where $P_\infty=[0:0:1:0]$ represents the point at infinity. A circle passing through that point at infinity is in fact a line. (Actually you might just as well consider $P_\infty$ to be the line or the circle at infinity, since it has properties of all of these.)
Of course there is a one-dimentional space of possible solutions for $h$, namely the different representatives of its homogeneous coordinate vector. To choose one representative implicitely, one can compute the $3\times3$ minors of this matrix, with suitable order and signs.
\begin{align*}
\begin{vmatrix}-2y_1&1&*\\-2y_2&1&*\\0&0&1\end{vmatrix}&=2(y_2-y_1)&
\begin{vmatrix}-2x_1&-2y_1&*\\-2x_2&-2y_2&*\\0&0&1\end{vmatrix}
&=4(x_1y_2-x_2y_1)\\
\begin{vmatrix}-2x_1&1&*\\-2x_2&1&*\\0&0&1\end{vmatrix}&=2(x_2-x_1)&
\begin{vmatrix}-2x_1&-2y_1&1\\-2x_2&-2y_2&1\\0&0&0\end{vmatrix}&=0
\end{align*}
$$h=\begin{bmatrix}
2(y_2-y_1)\\-2(x_2-x_1)\\4(x_1y_2-x_2y_1)\\0
\end{bmatrix}$$
So you alternate the sign for the minors and obtain a vector from the kernel. Now any circle that is perpendicular to the two given circles can be expressed as $\lambda p+\mu h$. But we want to exclude $p$ itself, which is what you'd get for $\mu=0$. So you can restrict yourself to $\mu\neq 0$ or, since the representation is homogeneous, to $\mu=1$. So write $q=h+\lambda p$ and there is no way that $q$ will be a multiple of $p$ unless $h$ is the null vector (which happens if your circles are concentric).
Now we need to solve
$$0=(h+\lambda p)^T\cdot B\cdot(h+\lambda p)=h^TBh + 2\lambda(p^TBh) + \lambda^2(p^TBp)$$
But if $p$ is a point then $p^TBp=0$ and one can solve for
$$\lambda = \frac{h^TBh}{-2p^TBh}$$
Instead of $q=h+\lambda p$ one can also use a multiple of this, so if you want to avoid the division, you may write
$$q = (h^TBh)p - 2(p^TBh)h$$
If you need the $x$ and $y$ coordinates of $q$, you'd divide the first resp. second coordinate of that vector by the last coordinate. Doing so with the formulas above I obtained the formula I wrote at the very beginning.
In some setups, the machinery to compute $h$ given $c_1,c_2,P_\infty$ may already be in place, and likewise the machinery to evaluate $h^TBh$ and $p^TBh$. In those cases, the formula used in this derivation here may be simpler than the one at the beginning.
Shortcomings
Note that the above will fail of both the circles are indeed lines, and the known point of intersection is the point at infinity. If both circles are lines, then $h$ will be $P_\infty$, and if $p$ is that same point, then all coefficients will be zero and $q=[0:0:0:0]$. One solution would be computing a line-line intersection in this case. Another solution would be using any other circle which is perpendicular to both the given ones, but picking an arbitrary element of a two-dimensional space appears to be harder.
