Limit Schottky Groups — Steve Trettel

Cusps, limit sets, and the boundary between discrete and indiscrete.

The Phenomenon

When we iterate a classical Schottky group—one built from disjoint circle pairs—the limit set forms a fractal dust scattered across the plane. But something remarkable happens as we bring circles from different pairs closer together: the moment they become tangent, the scattered dust coalesces into a single closed fractal curve.

This is the limit Schottky configuration, and it marks a phase transition in the geometry. The group remains free (the defining property of Schottky groups), but certain words in the generators become parabolic. These parabolic elements create cusps in the quotient surface, and the limit set threads through them in an intricate closed loop.

The limit set transforms from fractal dust to a closed curve at tangency

Why Tangency Changes Everything

To understand what’s happening, recall that a Schottky group of rank 2 is generated by two Möbius transformations $A$ and $B$ , each pairing circles:

$A$ maps $\mathrm{ext}(C_A) \to \mathrm{int}(C_{A'})$
$B$ maps $\mathrm{ext}(C_B) \to \mathrm{int}(C_{B'})$

When all four circles $C_A, C_{A'}, C_B, C_{B'}$ are disjoint, every non-trivial element of the group is loxodromic: it has two fixed points and spirals between them. The limit set is the closure of all these fixed points—a Cantor set of dimension between 0 and 2.

Now bring circles from opposite pairs into tangency. Say $C_A$ touches $C_B$ at a point $p$ . The composite transformation $B^{-1}A$ maps:

\mathrm{ext}(C_A) \xrightarrow{A} \mathrm{int}(C_{A'}) \to \text{[iterations]} \to \mathrm{ext}(C_B) \xrightarrow{B^{-1}} \mathrm{int}(C_A)

At tangency, the point $p$ lies on both $C_A$ and $C_B$ , placing it simultaneously in $\mathrm{ext}(C_A)$ and $\mathrm{ext}(C_B)$ . This means $B^{-1}A$ acts on $p$ and—here’s the key—it fixes $p$ as its unique fixed point. A Möbius transformation with a single fixed point is parabolic.

The geometric picture: $B^{-1}A$ is like a translation along the tangent direction at $p$ , but in the quotient space it creates a cusp—a funnel narrowing to a point that never quite closes.

The Cyclic Tangency Configuration

The most symmetric limit Schottky group has four circles arranged in a square, each tangent to its neighbors in cyclic order:

C_A \frown C_B \frown C_{A'} \frown C_{B'} \frown C_A

Label the tangency points:

$p_1 = C_A \cap C_B$
$p_2 = C_B \cap C_{A'}$
$p_3 = C_{A'} \cap C_{B'}$
$p_4 = C_{B'} \cap C_A$

Then we have four parabolic words:

B^{-1}A \text{ fixes } p_1, \quad A'^{-1}B \text{ fixes } p_2, \quad B'^{-1}A' \text{ fixes } p_3, \quad A^{-1}B' \text{ fixes } p_4

These four cusps determine the topology of the quotient surface: it’s a sphere with four punctures, or equivalently a genus-2 surface with four boundary components removed. The limit set is the preimage of a closed curve on this surface—it visits each cusp infinitely often in an intricate fractal pattern.

From Geometry to Möbius Maps

To implement this in GLSL, we need explicit formulas. A Möbius transformation pairing circles of radii $r_1, r_2$ centered at $c_1, c_2$ takes the form:

f(z) = c_2 + \frac{r_1 r_2}{z - c_1}

This maps the exterior of the first circle to the interior of the second. But there’s a subtlety: where in the target circle does a given point land?

For a loxodromic pairing, we have a free parameter—a rotation inside the target disk—that controls the translation length of the transformation. For the limit case, this rotation must be chosen so that the tangency points match up correctly.

The Rotation Constraint

Suppose $C_A$ and $C_B$ are tangent at $p_1$ , and we want $A$ to map $C_A$ to $C_{A'}$ such that $A(p_1)$ lands at the tangency $p_2 = C_{A'} \cap C_B$ .

The base Möbius map $f(z) = c_{A'} + \frac{r_A^2}{z - c_A}$ sends $p_1$ somewhere on $C_{A'}$ —but generically not to $p_2$ . We need to post-compose with a rotation $R_\theta$ (applied in normalized coordinates of $C_{A'}$ ) to move $f(p_1)$ onto $p_2$ .

The required angle is:

\theta = \arg\left(\frac{p_2 - c_{A'}}{r_{A'}}\right) - \arg\left(\frac{f(p_1) - c_{A'}}{r_{A'}}\right)

This is what computeRotationForTangency computes—it’s the angle that makes the parabolic word $B^{-1}A$ actually fix the tangency point.

float computeRotationForTangency(
    vec2 c1, float r1,    // source circle
    vec2 c2, float r2,    // target circle  
    vec2 sourceTan,       // tangency point to map from
    vec2 targetTan        // tangency point to map to
){
    float k = r1 * r2;
    // Where does base Möbius send the source tangency?
    vec2 base = c2 + complexDiv(vec2(k, 0.0), sourceTan - c1);
    
    // Normalize both points to unit circle
    vec2 uBase   = (base      - c2) / r2;
    vec2 uTarget = (targetTan - c2) / r2;
    
    // Angle between them is the required rotation
    return atan(uTarget.y, uTarget.x) - atan(uBase.y, uBase.x);
}

Implementing the Square Configuration

For four circles tangent in a square, we can set up a beautiful symmetric configuration. Place circles at the four corners, each of radius $r$ tangent to its two neighbors:

void setupGenerators(out Generator genA, out Generator genB, float iTime) {
    // Time-varying radii for animation
    float rA = 0.5 + 0.12 * sin(0.8 * iTime);
    float rB = 0.5 + 0.10 * cos(0.6 * iTime);
    
    // Four circles at square corners
    vec2 cA  = vec2(-rA, -rA);    // bottom-left
    vec2 cB  = vec2( rB, -rA);    // bottom-right  
    vec2 cAp = vec2( rB,  rB);    // top-right
    vec2 cBp = vec2(-rA,  rB);    // top-left
    
    // Four tangency points
    vec2 pAB   = vec2(0.0, -rA);  // bottom edge
    vec2 pBAp  = vec2( rB, 0.0);  // right edge
    vec2 pApBp = vec2(0.0,  rB);  // top edge
    // pBpA = vec2(-rA, 0.0)      // left edge (unused)
    
    // Generator A: bottom-left → top-right
    genA.interior = Circle(cA, rA);
    genA.exterior = Circle(cAp, rA);
    genA.rotation = computeRotationForTangency(
        cA, rA, cAp, rA,
        pAB,    // start at bottom tangency
        pBAp    // end at right tangency
    );
    genA.color = vec3(0.95, 0.2, 0.2);
    
    // Generator B: bottom-right → top-left
    genB.interior = Circle(cB, rB);
    genB.exterior = Circle(cBp, rB);
    genB.rotation = computeRotationForTangency(
        cB, rB, cBp, rB,
        pBAp,   // start at right tangency
        pApBp   // end at top tangency  
    );
    genB.color = vec3(0.2, 0.3, 1.0);
}

Why This Works

The key insight is that we’ve constrained the rotation angles to make the composition $B^{-1}A$ send:

p_1 \xrightarrow{A} p_2 \xrightarrow{B^{-1}} p_1

Since $A$ maps the exterior of $C_A$ to the interior of $C_{A'}$ , and the tangency $p_1$ sits right on the boundary, the point $A(p_1) = p_2$ is technically on the boundary of $C_{A'}$ . Similarly for $B^{-1}$ . This creates the parabolic fixed point.

In practice, numerical iteration will spiral closer and closer to these cusps, creating the characteristic fractal tendrils that define the limit set.

The Schottky Transformations with Rotation

The complete generator includes three steps:

Base Möbius map: $z \mapsto c_2 + \frac{r_1 r_2}{z - c_1}$
Normalize to unit disk: $(w - c_2)/r_2$
Rotate: multiply by $e^{i\theta}$

vec2 schottkyMap(vec2 z, Generator gen, inout float scale) {
    vec2 c1 = gen.interior.center;
    float r1 = gen.interior.radius;
    vec2 c2 = gen.exterior.center;
    float r2 = gen.exterior.radius;
    float k = r1 * r2;
    
    // Track derivative for distance estimation
    float denom = max(dot(z - c1, z - c1), 1e-6);
    scale *= (k / denom);
    
    // Step 1: Base Möbius
    vec2 base = c2 + complexDiv(vec2(k, 0.0), z - c1);
    
    // Step 2+3: Normalize and rotate
    vec2 unit = (base - c2) / r2;
    unit = complexRotate(unit, gen.rotation);
    
    return c2 + r2 * unit;
}

The inverse transformation undoes these steps in reverse order:

vec2 inverseSchottkyMap(vec2 z, Generator gen, inout float scale) {
    vec2 c1 = gen.interior.center;
    float r1 = gen.interior.radius;
    vec2 c2 = gen.exterior.center;
    float r2 = gen.exterior.radius;
    float k = r1 * r2;
    
    float denom = max(dot(z - c2, z - c2), 1e-6);
    scale *= (k / denom);
    
    // Step 1: Undo rotation
    vec2 unit = (z - c2) / r2;
    unit = complexRotate(unit, -gen.rotation);
    vec2 baseOut = c2 + r2 * unit;
    
    // Step 2: Inverse Möbius
    return c1 + complexDiv(vec2(k, 0.0), baseOut - c2);
}

Numerical Iteration and Distance Estimation

To render the limit set, we iterate the IFS defined by the four maps $A, A^{-1}, B, B^{-1}$ :

IterationResult iterateSchottky(vec2 z, Generator genA, Generator genB) {
    vec3 color = vec3(1.0);
    float scale = 1.0;
    
    for (int i = 0; i < MAX_ITERATIONS; i++) {
        bool moved = false;
        
        // Apply whichever transformation has z in its domain
        if (isInsideCircle(z, genA.interior)) {
            z = schottkyMap(z, genA, scale);
            color = mix(color, genA.color, 0.5);
            moved = true;
        }
        else if (isInsideCircle(z, genA.exterior)) {
            z = inverseSchottkyMap(z, genA, scale);
            color = mix(color, genA.color, 0.5);
            moved = true;
        }
        // ... similar for genB ...
        
        if (!moved) break;
    }
    
    return IterationResult(color, z, scale);
}

The scale parameter accumulates the magnitude of the derivative along the orbit. Since Möbius transformations have derivative $|f'(z)| = \frac{r_1 r_2}{|z - c_1|^2}$ , we multiply this at each step.

For distance estimation, we compute the raw distance from the final position to the nearest circle, then correct by the accumulated derivative:

float correctedDistance(IterationResult result, Generator genA, Generator genB) {
    float d = distanceToSchottkyCircles(result.finalPosition, genA, genB);
    float positionCorrection = 1.0 / (10.0 + length(result.finalPosition));
    return d / (result.scale * positionCorrection);
}

The position correction prevents artifacts when orbits escape to infinity—in the limit Schottky case with tangent circles, this rarely happens as orbits tend to spiral toward cusps instead.

Animation and Variation

By varying the radii $r_A(t)$ and $r_B(t)$ independently, we maintain tangency while morphing the shape of the limit set:

float rA = 0.5 + 0.12 * sin(0.8 * iTime);
float rB = 0.5 + 0.10 * cos(0.6 * iTime);

Because computeRotationForTangency automatically adjusts the rotation angles to maintain the parabolic condition, the limit set remains a closed curve even as the circles breathe.

The group stays “on the edge” between classical Schottky (disjoint circles, Cantor set limit set) and degenerate Schottky (overlapping circles, no longer a Schottky group at all). This edge case produces the most visually striking fractals.

Code Notes

A few implementation details worth noting:

Complex arithmetic: We represent complex numbers as vec2 and implement division, rotation, etc. manually. GLSL has no native complex type.
Numerical stability: The max(denom, 1e-6) guards prevent division by zero when we hit circle centers exactly. In practice this almost never happens due to floating-point imprecision.
Distance estimation: The derivative tracking could be made more sophisticated using automatic differentiation, but the simple chain rule multiplication works well for these smooth maps.
Color mixing: The mix(color, gen.color, 0.5) accumulates generator colors along the orbit. This gives a visual indication of which words in the group were used to reach each point.