How many mutually orthogonal circles are possible?

Reformulation using Möbius geometry

It is possible to represent a circle with center $(x,y)$ and radius $r$ using a vector

$$\begin{pmatrix}a\\b\\c\\d\end{pmatrix}= \begin{pmatrix} x \\ y \\ (x^2 + y^2 - r^2 - 1)/2 \\ (x^2 + y^2 - r^2 + 1)/2 \end{pmatrix}$$

(Or any multiple thereof, since this is a homogeneous representation.) Two circles are orthogonal if $(x_1-x_2)^2+(y_1-y_2)^2=r_1^2+r_2^2$. This is equivalent to $a_1a_2+b_1b_2+c_1c_2-d_1d_2=0$ since

$$ -2(a_1a_2+b_1b_2+c_1c_2-d_1d_2) = \\ -2x_1x_2 - 2y_1y_2 + 2(x_1^2+y_1^2-r_1^2)(\tfrac14+\tfrac14) + 2(\tfrac14+\tfrac14)(x_2^2+y_2^2-r_2^2) = \\ x_1^2-2x_1x_2+x_2^2+y_1^2-2y_1y_2+y_2^2-r_1^2-r_2^2 = \\ (x_1-x_2)^2+(y_1-y_2)^2-(r_1^2+r_2^2) $$

This representation is sometimes called Möbius geometry. The bilinear form $a_1a_2+b_1b_2+c_1c_2-d_1d_2$ has signature $(+,+,+,-)$. So the whole geometry describes a Minkowski space $\mathbb R^{3,1}$. You might therefore ask the equivalent question of how many mutually orthogonal (with respect to the given bilinear form) vectors there can be in such a space.

Using known results

Wikipedia writes:

An orthonormal basis for Minkowski space necessarily consists of one timelike and three spacelike unit vectors.

Apparently this is a consequence of Sylvester's law of inertia. The distinction between timelike and spacelike unit vectors depends on the sign of the quadratic form corresponding to the bilinear form. If you assume that $x_1=x_2$ and $y_1=y_2$, then the differences $x_1-x_2$ and $y_1-y_2$ in the above formula vanish. The sign therefore only depends on the sign of $r^2$. And since $r^2$ will always be positive for real circles, all your circles are the same, namely spacelike. A timelike vector would correspond to a circle of imaginary radius.

Well, you don't need unit length here. But if you had four spacelike orthogonal vectors, then you could scale them to unit length without changing the sign of the quadratic form, and without changing the circle they describe either. Which means that if you had four real and pairwise orthogonal circles, you could derive a purely space-like basis from those, in contradiction to the statement above.

Notice that light-like vectors, i.e. circles of radius zero, i.e. points, can not be scaled to unit length. But you can't have two distinct points being orthogonal to one another, nor three circles coincide in one point and still be orthogonal to one another. So points as a special case of circles don't solve this either.

Conclusion

So that paper you reference is either investigating properties of something that cannot exist, or they deliberately consider the case of an imaginary radius as well, as your comment suggests that others did.

Example

Applying the above formalism to a regular triangle of edge length $\sqrt2$, or more precisely the regular triangle with corners $(0,\sqrt{2/3}), (\pm\sqrt{1/2},-\sqrt{1/6})$, you get the circles of radius $1$ described by the vectors

$$ \begin{pmatrix}0\\\sqrt{\frac23}\\-\frac23\\\frac13\end{pmatrix} \qquad \begin{pmatrix}-\sqrt{\frac12}\\-\sqrt{\frac16}\\-\frac23\\\frac13\end{pmatrix} \qquad \begin{pmatrix}\sqrt{\frac12}\\-\sqrt{\frac16}\\-\frac23\\\frac13\end{pmatrix} $$

You can check that these are indeed mutually orthogonal, by computing all pair-wise bilinear forms. You can also find a fourth circle orthogonal to these using the following system of equations, where I already applied the matrix of the bilinear form to the vectors (i.e. negated the last coordinate):

$$\begin{pmatrix} 0&\sqrt{\frac23}&-\frac23&-\frac13 \\ -\sqrt{\frac12}&-\sqrt{\frac16}&-\frac23&-\frac13 \\ \sqrt{\frac12}&-\sqrt{\frac16}&-\frac23&-\frac13 \end{pmatrix} \cdot\begin{pmatrix}a\\b\\c\\d\end{pmatrix} =\begin{pmatrix}0\\[2ex]0\\[2ex]0\end{pmatrix}$$

The matrix has rank $3$, so the kernel has dimension $1$, corresponding to a uniquely defined circle. One element of that kernel would be $(0,0,1,-2)$. In order to obtain a representative of the form given above we need to scale this so that $d-c=1$. We get the vector $(0,0,-1/3,2/3)$ from which we can read $-r^2=x^2+y^2-r^2=c+d=1/3$ and therefore $r=\sqrt{-1/3}$. (You could have arrived at that particular solution more easily, e.g. by assuming that $x=y=0$ for reasons of symmetry. But I wanted to show how the machinery works, so I didn't take any shortcuts.)

Suppose that $a,b,c$ are mutually orthogonal circles. Let $X$ be an intersection point of $a,b$. Consider inversion $\phi$ with respect to a circle centered at $X$ and arbitrary radius $r$.

Using basic properties of inversion we get that $\phi(a), \phi(b)$ are lines passing through some point $Y$ and $\phi(c)$ is a circle. Since inversions are conformal, we get that lines $\phi(a), \phi(b)$ are perpendicular and moreover $\phi(c)$ is orthogonal to both $\phi(a), \phi(b)$, which means that $Y$ is the center of $\phi(c)$.

Now if circle $d$ is orthogonal to $a, b$ then using a similar argument as above we conclude that circle $\phi(d)$ is centered at $Y$. But then $\phi(c), \phi(d)$ cannot be orthogonal. Therefore $c,d$ cannot be orthogonal as well.

Therefore there are no four mutually orthogonal circles on a plane.