How many mutually orthogonal circles are possible? How many mutually orthogonal circles is it possible to have? 
It is easy to construct $3$ mutually orthogonal circles, e.g. $3$ circles with radius $1$ and centers at the vertices of an equilateral triangle of side length $\sqrt 2$. But how does one construct 4 mutually orthogonal circles? Apparently it is possible given this link: http://www.jstor.org/stable/1968006?seq=1#page_scan_tab_contents
It isn't clear from that link how to do it, though. In this link http://mathworld.wolfram.com/OrthogonalCircles.html they describe how to construct a circle which is orthogonal to $3$ existing circles, but the method doesn't work for my example above as the radical center of the new circle is inside the $3$ existing circles. And it seems clear that no circle can be orthogonal to another circle if the center of one circle is inside the other.
So I guess I have 2 questions:
1) How can one construct 4 mutually orthogonal circles?
2) What is the maximum (if any) of the number of mutually orthogonal circles?

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.)
