# Introduction

I am a physicist by trade, but I've always felt myself to be part geometer. I cut my graduate-school teeth on string theory, and my doctoral work was concerned with supergravity and topological field theory. During the last few years of my teaching career, I became very interested in hyperbolic manifolds. My first exposure was via Perelman's proof of Thurston's Geometrization Conjecture. It became clear after a couple of pages that I was in way over my head.

What followed was several years of struggling through papers and books by, among many others, Thurston, Ratcliffe and Marden. Finally, I found the book "Indra's Pearls", by Mumford, Series and Wright, wrote a program, and began reproducing their plots. Along the way, I became curious about a number of things: why does it take so long for groups just "south" of the Maskit boundary to display their non-discreteness? Is there a way to find non-free groups in the Riley slice? What's it like outside those slices of Teichmuller space?

For the most part, I use the phrase "Teichmuller space" to refer to the Teichmuller space associated with groups of two generators. All are, of course, subgroups of the Moebius group on E2, the complex plane plus the point at ∞, and will be treated as matrices in PSL(2,C). I will refer to the generators as "a" and "b", and to their inverses as "A" and "B", respectively.

This work is an attempt to summarize what I've learned while pondering those questions. What I have learned comes mostly in the form of conjectures, much as experimental physicists might conjecture explanations for their observations. I am no mathematician, though I seem to be playing at it more and more these days. When I quote results, I'll note the references, most of which I am still far short of fully understanding. If you see me make any errors, or have any comments about what I've written here, please e-mail me at kenneth.koehler@uc.edu.

# The experimental apparatus

All of the limit sets I'll display here were generated in the same fashion:

1. The generators are constructed either using their actual components, or from one of the following methods (I did a fair amount of experimentation to get to this point!):

1. using the traces tr(a), tr(b) and tr(ab) (Mumford, Series and Wright):
P = tr(a)2 + tr(b)2 + tr(ab)2 - tr(a)tr(b)tr(ab)
Q = ± √(4-P)
R = ± √P
a = {{tr(a)/2, (R+tr(b)) (iQ*tr(a)+2tr(ab)-tr(a)tr(b)) / (2(tr(b)2-P))}, {(tr(a)2-4) (tr(b)2-P) / (2(tr(b)+R) (iQ*tr(a)+2tr(ab)-tr(a)tr(b))), tr(a)/2}}
b = {{(tr(b)-iQ)/2, (tr(b)+R)/2}, {tr(b)2+Q2-4) / (2(tr(b)+R)), (tr(b)+iQ)/2}}

(the sign on R is chosen to avoid singular matrices; the other ± sign is chosen to make all fixed points lie on or in the unit circle, if that is possible)

I worked this out myself following Mumford et. al., and my result is slightly different but is still correct. A careful eye will notice that my limit sets are conjugate reflections of theirs.

Note that this method produces singular matrices when either P=tr(b)2, or P=4 and tr(ab)=tr(a)tr(b)/2.
2. using the Maskit parameter μ:
a = {{-iμ/2, i ± μ/2}, {-(μ2+4)/(4i ± 2μ), -iμ/2}}
b = {{1-i, 1}, {1, 1+i}}

(where the ± sign is chosen to make all fixed points lie on or in the unit circle; note that with μ=2i and the negative sign choice, a becomes the identity matrix; this corresponds to the (0,1) "gasket", whose limit set was computed using the components

a = {{1, 0}, {-2i, 1}}
b = {{1, 1-i}, {1, 1+i}})
Again, this is different from Mumford et. al. There are two solutions to the Markov identity
tr(a)2 + tr(b)2 + tr(ab)2 = tr(a) tr(b) tr(ab)
and I have used the other one.
3. using the Riley parameter λ:
a = {{1, 2}, {0, 1}}
b = {{1, 0}, {λ, 1}}
4. using the traces t(b) and t(ab) (assuming a is parabolic):
a = {{1, i}, {0, 1}}
b = {{1, i(tr(b) - 2) / (tr(ab) - tr(b))}, {-i(tr(ab) - tr(b)), tr(b)-1}}
(here tr(ab) ≠ tr(b))
5. using the traces t(a) and t(b) (assuming the product ab is parabolic):
P = √(4 - tr(a)2)
a = {{(tr(a) ± i P)/2, 0}, {±(P (tr(b)-2) ±i (tr(a) tr(b)-4))/2, (tr(a) - ±i P)/2}}
b = {{1, i}, {i (2 - tr(b)), tr(b)-1}}
(the sign choice here is an input)
For method 4, whenever I wanted to plot a limit set for tr(ab)=tr(b), I offset tr(ab) slightly. In order to be able to consistently compare plots, whenever using method 5 for tr(a)=tr(b), I similarly offset tr(a).
2. After fixing a and b, I then computed A, B and the commutators abAB, AbaB, BabA and ABab, and their inverses.
3. If a "special word" (a string of a's and B's which is parabolic) has been input, it and its inverse are computed.
4. Each of the generators, commutators and the special word have two fixed points:
z = (g1,1 - g2,2 ± √(tr(g)2-4)) / (2g2,1)
(gi,j are the matrix components)
These are computed (the fixed points of their inverses are the same).
5. The plot area is a fixed size in pixels; if zmax is the largest absolute value of the fixed points, the coordinates of the corners are set to ±zmax±i zmax.
Outside of the Maskit slice, zmax varies considerably, and will be specified for those limit sets below.

Whenever a method produces multiple conjugates which yield finite values of zmax, if one of the conjugates yields a zmax of 1, that is chosen; otherwise, the conjugate with the smallest zmax is chosen.

6. If zmax is 1, the program randomly iterates left multiplication by a, b, A and B (but avoids the inverse of the last multiplication, to avoid wasting time), and after each iteration it maps the unit circle using the resulting product g:
z = 1/(-g2,2/g2,1)*
center = g.z

(where "*" indicates complex conjugation and "g.z" indicates

(g1,1z + g1,2)/(g2,1z + g2,2)).
Each new circle (disc) with a radius larger than two pixels and lying within the plot area is colored in light gray (the radius is reduced by two pixels to guard against rounding errors in the auto-detection of non-discrete groups; see step 7). After (by default) 100 iterations (corresponding to a word length of 100), the process starts over.

If the program is used interactively, the plots are displayed in real time using OpenGL; if the program is being used in batch mode, to create a png file of the plot, the "circling" stops after a total of one million iterations.

At this point the plot is a collection of discs which are each equivalent to the exterior of the unit disc; for a cusped group, these discs represent an approximation to one half of the regular set of the group.
7. The iteration process described above is now applied to every fixed point computed in step 4. Each fixed point is associated with a different color, and once a point is plotted another cannot be plotted there. The maximum word length here defaults to 200. If the program is generating a png file, the total number of iterations is (by default) eight million.

The resulting plot now includes a dust approximation to that portion of the limit set of the group lying within the plot area:

A sample plot (the (2,15) gasket group), and a blow-up of the central region showing the color-coding of the fixed points:

........

If any limit set point has landed on a previously-mapped circle, the program flags the group as non-discrete.

8. Finally, the program outputs the number of points and circles plotted.

Adjustable parameters include the corner coordinates and the maximum word length l (which defaults to 200). If l is modified, the total number of fixed point iterations is set to l3, but the maximum number of circle iterations is left at one million.

The program uses the number of seconds past January 1, 1970 as a seed for the random number generator. Curious about how much variation occurs in the limit sets due to the stochastic nature of the algorithm, I created limits sets for the (2,5) gasket group (see below) about five days apart. The two limit sets were identical. For these plots, the plot area was 800 pixels2. Wondering whether increasing the plot size would alter the situation, I ran two more with a plot area of 3200 pixels2, a few minutes apart. Again, the results were identical.
While all of the limit set plots were done in this fashion, the Riemann surface plots, the real trace rays and all color contour plots were done using Mathematica 7.0. Gimp was used to pinch curves on the Riemann surfaces, and ImageMagick was used to build the montages. LibreOffice was used for simple data analysis.

All computations were done in double precision. Because this work is based on numerical computation, it is important to know what hardware was used. These calculations were done on two Intel Core Duo CPUs, an E8400 stepping level 10, microcode level A07, and an E7500 stepping level 10, microcode level A0B. These processors use 64-bit doubles, with 52 bits for the significand, giving more than 15 digits of precision; in both, computations in the CPU use 80-bit doubles with 64 bits for the significand, giving more than 19 digits of precision. Both OS are Slackware 14.2, running kernels 4.14.11 or later (kernel and library upgrades occurred on the E8500 during the course of this work). The program was compiled with gcc -Wall -Wconversion -g -Og, which is supposed to optimize while facilitating debugging (some variables were optimized out).

While working on this project I made an interesting observation, based on plots like these:

Non-discrete region for tr(a) = 2, tr(ab) = 4+5i, point density and zmax for same:

(tr(b) varies from i at bottom left to 8+9i at top right)

• The first plot shows which groups were nondiscrete, based on simple visual analysis of over 6,000 limit plots.
• The second plot encodes the Log10 of the number of points displayed in each limit plot I created for those traces, with Nmin = 93 encoded as red (outer edges), and Nmax = 948,747 encoded as magenta (center).
• The third plot shows Log10(zmax) for those traces, with zmax ranging from 1.131 (red) to 1768.21 (magenta).

The point density closely mirrors the discreteness of the corresponding groups, although it hides what I suspect is a fractal boundary for the non-discrete groups (based on what we'll see below in the Maskit and Riley slices). The zmax plot clearly indicates the center of the non-discrete region, but it is clear that zmax drops off quickly as distance from the value of tr(ab). I've used the point density visualization technique in all of the sectors of Teichmuller space that I have investigated. The plots were made using Mathematica. To standardize the plots, all of the limit sets used to do this had fixed word length=100, iteration limit=106 and zmax=1.

Teichmuller space is a manifold with boundaries, so locally it is homeomorphic to Rn (n is discussed below). Sampling is therefore required in any study of individual limit sets. In order to make the first plot of the three above, I had to plot 6,561 limit sets. These were sampled at intervals of 0.1 and 0.1i. But in general it is unclear whether this, or any particular, level of detail is sufficient to display the nature of any finite region of the space. In order to get a feel for this, I looked at the same region, but with the more course intervals of 0.5 and 0.5i. These two plots have been scaled to the same image resolution for comparison:

........

Modulo isolated instances, the limit sets appear to be "continuous functions" of their coordinates in Teichmuller space. Because of the large number of limit sets involved, in most cases I have opted to sample at intervals of 1 and i.

Since we will be plotting hundreds of thousands of limit sets, some discussion of performance is warranted. When running the program interactively, the user is responsible for deciding when to terminate the algorithm. Whether mapping circles or fixed points, it is usually clear that after a short while the time between "new" points (circles) becomes very long. When running the program in batch mode, the program as described arbitrarily terminates after l3 iterations, where l is the word length. The following plot describes the results of plotting the (2,5) gasket for word lengths varying from 100 to 1000:

The time required to complete l3 iterations increases approximately quadratically with l. For all this effort there is essentially no improvement in the number of circle maps, and little improvement in the number of fixed points mapped after l=400.

I then ran a number of tests, with the word length varying between ten and 3,200, and maxit varying from 100 to 100 million. The first tests were for the (2,5) gasket (left), and the gaskets at μ=1+2i (middle) and μ=1+√3i (right):

Horizontal axis is log10maxit; lower plot is log10runtime(s), upper plot is log10points; color varies with word length according to the following key:

With the original termination criteria and a word length of 200 (which implied an iteration limit of 8 million), the (2,5) gasket took 21 seconds and generated 50,538 points. With the new criteria, a word length of 200, and maxit set to 1,000, the plot completed in 2 seconds and contained 53,179 points. While there was little improvement in the point tally for word lengths above 200, it is worth noting that a word length of 3,200 and a maxit value of 10,000 generated 65,985 points in 10 seconds.

I also ran tests for the quasifuchsian group for tr(a)=-2-i, tr(b)=-3-i (right; word length = 3,200, maxit = 10,000):

........

Looking at the data, it seems that there is little need to consider word lengths above 50, or maxit values above 10,000. While larger word lengths and maxit values seem to yield higher point counts, the difference is hardly worth the processing time.

But that changed when I tested the Riley slice: run time still increased with maxit, but point totals did not always increase with word length. I ran tests for λ=1+i,...,5+5i. For the first three values of λ, the point counts for word lengths greater than 10 were

1. 3313 ± 6, σ = 2.9114
2. 426 ± 2, σ = 0.7559
3. 230 ± 1, σ = 0.5345.

For the last two values of λ (and maxit=108), the point counts as a function of word length were:

(λ=4+4i in red, λ=5+5i in blue)

The peak point counts were for a word length of 100 for λ=4+4i and 50 for λ=5+5i, which took approximately 20 minutes each. And for those two limit sets, where I expected discrete groups I saw instead the first random signs of non-discreteness, here circled in red (click on the image for a 5-image sequence showing λ=1+i,...,5+5i):

(all images in the sequence have word length 50 and maxit value 108)

Lest this be the result of numerical instability, I checked the group at λ=i for a maxit value of 108 and word lengths of both 50 and 10,000; neither plot displayed any randomness. It seems that large point densities are a sufficient condition for non-discreteness, but they are not necessary.

# The structure of Teichmuller space for genus 2

Before concentrating on genus 2, let me review how I got here from hyperbolic manifolds.

Any complete hyperbolic 3-manifold M is the quotient of H3 by a Kleinian group G, with fundamental group π1(M) isomorphic to G. Its limit set Λ(G) (of all fixed points of G) lies on the S2 at ∞ (which is taken to be the "boundary" of H3, and is isomorphic to E2). The complement of Λ(G) in E2 is Ω(G), the regular set.

A Kleinian group is a discrete torsion-free group of isometries of H3 (and is therefore a discrete subgroup of Moeb(E2)), which acts properly discontinuously (the inverse image of any compact set is also compact). The latter implies that a discrete group has a countable number of elements. Kleinian groups are usually taken to be non-elementary - they leave an infinite number of fixed points on the sphere at ∞ (dim(Λ) > 0). The components of Λ are either topologically circles, fractal sets, or totally disconnected (every component is a point).
(Marden, sections 2.2-2.5)
Specifying that G is torsion-free means G has no elliptic elements (if G has elliptic elements, M is an orbifold, with conical punctures).
If G has n generators and no parabolic elements, Λ(G) is totally disconnected, and the boundary ∂M(G) = Ω(G)/G is a Riemann surface R of genus n. Such a surface has 3n-3 disjoint, non-parallel simple closed geodesics (which divide ∂M(G) into 2n-2 pairs of pants). The surface can be cut along any such loop and reconnected with a twist. The Teichmuller space of M, which describes the set of hyperbolic structures on M, can be parameterized by the lengths and twist angles of each loop. It therefore has real dimension 6n-6; since it possesses a complex structure, we often say it has complex dimension 3n-3.
(Thurston Notes, section 5.3)
If G contains parabolic elements, it can have at most 3n-3, and for each parabolic element, the corresponding loop is pinched to zero length.

Teichmuller space is the quotient of the space of metrics admissible to R by the action of diffeomorphisms homotopic to the identity, by a homotopy which takes the boundary into itself at all times. (Thurston, section 4.6) There is therefore not a 1:1 correspondence between points in Teichmuller space and the space of all hyperbolic manifolds; conjugate groups g and g' (where g'=hgH for some h in Moeb(E2)) describe the same manifold. When the program chooses a ± sign, it is actually choosing between conjugate groups. For the hyperbolic, quasifuchsian and cusped groups below, there are two limit set plots, one for each sign choice. For the other limit sets, I include one which has been conjugated by

h = {{i, -1}, {1, -i}}/√2
which moves the points at ±∞ to (0,i) and the origin to (0,-i).

It should now be clear why I have restricted myself to the genus 2 case: the corresponding Teichmuller space has complex dimension 3. This makes it impossible for creatures like myself, inhabiting a world with 3 real spatial dimensions, to visualize any Teichmuller space. As we will see, at best one can look at cross-sections, and try to imagine how they fit together.

I'm also very interested in Moebius subgroups which do not generate hyperbolic manifolds. These groups are non-discrete; an infinite sequence of group operations converges to the identity, and the limit set is the entire complex plane. If you can characterize the transition from discreteness to non-discreteness, you can identify a boundary in Teichmuller space.

Another interesting type of Moebius subgroup is the non-free group. Here, a finite sequence of group operations equals the identity (which of course means there the group contains elliptic elements). In this case, the limit set is a set of isolated points (dust), and the quotient Ω(G)/G is an orbifold. Non-free groups seem to be isolated points in the larger space in which Teichmuller space can be embedded.

I'll now specialize to genus 2. The two-holed torus can be pinched in four ways:

Since there can be at most three pinches for genus 2, it is clear that in general, there are more ways to pinch a torus than can be "used" to create hyperbolic manifolds.
The surface R (and hence the corresponding hyperbolic manifold) can be parameterized by three complex numbers. I will take them to be tr(a), tr(b) and either tr(ab) or tr(abAB). Which ones are specified determine which "sector" of Teichmuller space they describe:

1. If all the elements of G are hyperbolic, R is a two-holed torus, and the Λ is a totally disconnected infinite set of points (dust):

(method 1, tr(a) = tr(b) = 2.8, tr(ab) = 3.65+1.1i ;
zmax = 1.609 (left), 1.219 (right))

2. Now let's start pinching curves. If word "a" becomes parabolic, R is a twice-punctured torus (from a topological point of view, the two "ends" of the pinch are separated and a point is removed at each end):

(method 4, tr(b) = 3+2i and tr(ab) = 2.01+3.01i ;
zmax = 2.641 (left), 4.308 (right))

Note that we do not consider "b" parabolic as a separate case; exchanging generators results in a conjugation.
Considering what we just mentioned about conjugation, it would seem that while exchanging a and b would not produce a different hyperbolic manifold, it would describe a different region of Teichmuller space. However, a and b are simply labels; one does not expect Teichmuller space to care about generator labels.
Note the recurring patterns in the limit sets, indicative of translational symmetry caused by a fixed point at ∞. All limit sets, with any pinching, may be conjugated to a similar form.

3. If word "ab" becomes parabolic, R is again a twice-punctured torus:

(method 5, tr(a) = -1.01-1.01i and tr(b) = -1+2i ;
zmax = 1.219 (left), 1.129 (right))

These and the previous pair of limit sets appear to include closed circles. This is an artifact arising from both the stochastic nature of the algorithm and the resolution of the images. Quasicircles only appear in limit sets when R separates into two components ("connected" by a pinch). In genus 2, this only occurs when abAB is parabolic.
4. If word "abAB" becomes parabolic, R is a pair of once-punctured tori and Λ is a quasicircle separating the regular set into two components:

(method 1, tr(a) = tr(b) = 2.2, tr(ab) = 2.42+1.955403i ;
zmax = 1 (left), 4.122 (right))

G is called a quasifuchsian group.
In each case where a single curve is pinched, we have reduced the problem of "seeing" Teichmuller space to that of a 4-dimensional space; still difficult, but perhaps approachable.

5. If words "a" and "b" become parabolic, we have a cross-section of Teichmuller space called the Riley slice. R is a sphere with four punctures:

(method 3 (λ = 0.05+0.93i), tr(ab) = 2.1+1.86i;
zmax = 1.913 (left), 2.81 (right), both here multiplied by 12 to get a larger perspective)

6. If words "a" and "ab" become parabolic, R is again a sphere with four punctures:

(method 5, tr(b) = 3+3i ;
zmax = 1.161 (left), 5.515 (right))

7. When words "a" and "abAB" become parabolic, G is a single-cusp group; R is the union of a triply-punctured sphere and a singly-punctured torus:

(method 1, tr(b) = 3, tr(ab) = 3+2i;
zmax = 1 (left), 2.236 (right))

With single-cusp groups, one half of the regular set has closed up into a set of tangent discs, all of which are identified; the other half is simply connected.
8. If words "a", "b" and "ab" all become parabolic, we have the modular group, which is rigid (up to conjugation, it is unique); R is a pair of thrice-punctured spheres:

(method 3 (λ = -2);
zmax = 2.618 (left), 1 (right))

9. Finally, when words "b", "abAB" and a "special word" all become parabolic, R is again two thrice-punctured spheres and we have the groups of the Maskit boundary (which is really a special case of sector 6):

(method 2; μ = 0.766588417465459-1.642138768653476i);
this image corresponds to the (2,5) group, where in addition to b and abAB, the special word a3Ba2B is also parabolic;
this means the group is doubly-cusped; the "a" puncture represents a pinch with five complete twists around the "a" meridian, and two full twists around the "b" meridian;
zmax = 1 (left), 3.882 (right))

The purely-hyperbolic groups constitute the vast majority of Teichmuller space; the other cases are, after all, special cases. Let's examine the special cases in more detail, starting with the most constrained:

Using the Maskit parameter μ and
method 2 above,
• tr(a) = -i μ
• tr(b) = 2
• tr(ab) = i (2-μ)
• tr(abAB) = -2
Each group G(μ) is uniquely determined by μ, but G(μ) is conjugate to
• G(μ + 2k), k∈Z
• G(μ*)
• G(-μ*) (Mumford et. al.), and therefore (as we focus on the region 0 ≤ re(μ) ≤ 2),
• G(2-μ)
So we only need focus on im(μ) ≥ 0. For im(μ) ≥ 2, all G(μ) are free and singly-cusped, so our primary region of interest is 0 ≤ re(μ) ≤ 2 and 0 ≤ im(μ) ≤ 2. In the montage below, we have sampled this region in increments of 0.1 and 0.1i; each sample limit set is centered on its coordinate:

Samples from the Maskit slice (click for larger image):

(zmax = 1 for all Maskit slice limit sets)

The red points trace the Maskit boundary; all groups on the boundary are gaskets (doubly-cusped, with a special word of trace 2) or degenerate (the black point marks a degenerate group), and all groups above the boundary are discrete and singly-cusped. All groups below the boundary are either non-discrete or non-free; examples of the latter are at μ=i, and of course μ=2+i. The samples with im(μ)=0 (at the bottom), whose limit sets are a circle, appear to have the two sides of the regular set identified (μ=0 is a non-free group; tr(a)=0, so a2 is the identity matrix). This also seems to be the case for the (1,4) gasket group at μ=1+i.

Since the only discrete groups with a geometric meaning are above or on the boundary, groups below the boundary are not part of the Teichmuller space for the genus 2 torus.

........
1. Note that p and q must be relatively prime. Also, due to step 7 below, groups with p>q are conjugate to groups with p<q.
2. First find a special word, consisting of q "a"s and p "B"s, which consistently wraps the torus. To do this, you can use the following arithmetic prescription: construct the set starting with {1}, adding q or subtracting p so that the elements are always positive and never greater than p+q, until you return to 1. The group element is the product of a's and B's created from writing an a for each subtraction and a B for each addition (in reverse order). For instance, the (2,5) set is
{1, 6, 4, 2, 7, 5, 3, 1}
and the corresponding group word is aaaBaaB.
3. Using the canonical form of the Maskit groups (for simpler computations; the trace is invariant under conjugation):
a = {{-iμ, -i}, {-i, 0}}
b = {{1, 2}, {0, 1}}
find the trace of the special word. For our example, the (2,5) group, the trace is
-i5 - 4μ4 + 9μ3 - 12μ2 + 9μ - 4)
4. Set this trace (which will be a polynomial of order max(p,q)) to 2, and find all of the solutions which satisfy 1≤im(μ)≤2 (no group for μ<1 can be discrete). For the (2,5) group, these solutions are
0.375189 + 1.30024i, 1.06548 + 1.2824i and 0.766588 + 1.64214i.
I used Mathematica for this.
5. In general, the one with the largest imaginary part will be a discrete gasket group, lying on the Maskit boundary.
For any (p,q), only one solution of the trace polynomial will be discrete and lie on the boundary. (Mumford et. al.)
6. For p or q < 0, begin with (-p,1), for which the parabolic word is abp, and (0,1) (parabolic word a). For any (p,q) find relatively prime integers (r,s) for which ps-rq=±1. A new parabolic word will be
w(p+q),(r+s) = w(p,q)w(r,s).
So, for example, w(-1,2) = w(-1,1)w(0,1) = aba.
(This construction works for positive p and q as well.)
(p,q) and (-p,-q) will produce the same word.
7. I have not proved this, but for every example I have tried (0≤p≤21 and 1≤q≤32, p<q), (p,q) and (p+kq,q) (k∈Z) are conjugate. In addition, from a previous exploration for |p| and |q| both ≤ 5, the groups (-p,q) and (p,q) are conjugate. These results appear to arise from the facts that

1. G(μ) and G(μ+2k) (k∈Z) are conjugate, and
2. G(μ) and G(-μ) are conjugate.
8. (0,2) is conjugate to (0,1), and all (0,k>2) are non-free; their parabolic word = -I (here I is the identity matrix).
The majority of the groups shown in the sample plot above are non-discrete; in principle, they should all be completely colored in. The stochastic algorithm I used to produce all of these limit set plots will only hit all of the points in the complex plane if it is run for an infinite amount of time. The point density plot for the Maskit slice is:

(Nmin = 2,590, Nmax = 550,550)

The slice was sampled at intervals of 0.01 and 0.01i. All groups sampled below the black line were auto-detected as non-discrete. The "light spots" indicate non-free groups; the red line at the bottom corresponds to the "identified" groups on the real axis described above. The plot is markedly non-symmetric about the line Re(μ)=1 because of the tendency for groups with higher p/q to "fill" more slowly with the stochastic algorithm.

To understand the transition from discrete single-cusp groups to gasket group to non-discrete groups, I plotted an 11-image sequence for Re(μ)=1 and Im(μ)=4, 3.5, 3, 2.5, 2, 1.9, 1.8, 1.732050807568877 (the double-cusp), 1.7, 1.6 and 1.5, where the last three are non-discrete (μ=1+1.7i is the second plot; click on this image and the sequence will appear in a new tab):

You can see that the non-discreteness displays itself as an increasing "overlap" between the upper and lower halves of the limit set. Some of the groups just below the Maskit boundary are very slow to show this overlap. This is probably a function of the limited resolution of the plots. Of the 3,146 groups examined (step 5 above),

• 2,465 groups with (p,q) from (1,5) to (21,31) were identified as non-discrete with a word limit of 200;
• 46 groups with (p,q) from (2,27) to (20,23) were identified as non-discrete with a word limit of 400;
• 53 groups with (p,q) from (2,27) to (17,31) were identified as non-discrete with a word limit of 600;
• 14 groups with (p,q) from (8,25) to (21,26) were identified as non-discrete with a word limit of 800, and
• 26 groups with (p,q) from (4,31) to (21,29) were identified as non-discrete with a word limit of 1000.
Of the 288 limit sets which were not clearly identifiable by the program with word length 1000, 180 were non-discrete (having stared at over three thousand limit sets by this point, the last few hundred were easily separated by a trained eye). The (9,20) gasket, which we will look at in more detail below, has a non-discrete solution to its trace polynomial at μ=0.972815685123106+1.64763216820249i. It is one of the groups which required a word limit of 400 to be auto-detected as a non-discrete group by the program. In those runs, the plot resolution was 800 pixels2. When the resolution was increased to 3200 pixels2, the same group was auto-detected with a word length of 200. (At a resolution of 1000 pixels2, the (6,31) group at μ=0.178840219946299+1.74888845124066i failed to auto-detect non-discreteness at a word length of 102,400 and a maxit value of 108.)
This "overlapping" appears to be ubiquitous; the following sequence for tr(a)=2, tr(ab)=4+5i shows the transition from tr(b)=4.5+3.5i to 4.5+4i (click for a 6-image sequence):

As part of these computations I have of course constructed trace polynomials P(p,q)(μ) for those 3,146 groups. For any such polynomial, the solution of
Im(P(p,q)(μ)) = 0
is a real trace ray: along that curve the trace of the polynomial is real. Keen and Series showed that traveling along those rays from above the Maskit boundary to the point P(p,q)(μ)=2 defining the (p,q) gasket group corresponds to shrinking the length of the "pinching curve" to zero. Here I have plotted real trace rays for the first few free gasket groups after (0,1):

( (2,3) is conjugate to (1,3), (3,4) is conjugate to (1,4) and (4,5) is conjugate to (1,5) )

Gasket groups on the boundary are marked by white dots. The trace polynomials for (1,q≤5) are:

• P(1,2)(μ) = -μ2+2*μ-2
• P(1,3)(μ) = i3-2μ2+3μ-2)
• P(1,4)(μ) = μ4-2μ3+4μ2-4μ+2
• P(1,5)(μ) = -i5-2μ4+5μ3-6μ2+5μ-2)
While the other groups are conjugates, their trace polynomials are different:
• P(2,3)(μ) = i3-4μ2+7μ-4)
• P(3,4)(μ) = μ4-6μ3+16μ2-20μ+10
• P(4,5)(μ) = -i5-8μ4+29μ3-56μ2+57μ-24)
Note that the real trace rays for the (1,4) and (its conjugate) (3,4) meet at μ=1+i, the location of the (1,4) group discussed below. The "landscape" around the (1,3) gasket group is shown below:

(White point at right is the location of conjugate (2,3) gasket group.)

The colors represent Re(P(p,q)(μ)) for values between ±4. Keen and Series were not interested in the analytic continuation of the real trace rays to values < 2 as they appear to have no geometric significance. For a more complicated group, the landscape is much more interesting:

The landscape of the (9,20) gasket group:

(White points indicate solutions of P(9,20)(μ) = ±2; color contours are as above, with additional -2 contours indicated.)

Keen and Series showed that the real trace rays above the Maskit boundary (for the groups on the boundary) foliate the space above the boundary. The above plot indicates that this will be far from true for the analytic continuations, although the significance of these is unclear.

........

(μ = i)

Of the groups I've investigated, (p odd and ≥ 5, q=2p+2) all seem to be non-free, although (up to word length 12,800 and maxit value 108) they do not manifest the circular dust patterns:

(from left, (5,12), (7,16), (9,20), (11,24), (13,28) and (15,32); for all, re(μ) = 1;
from left, im(μ) = √2, 1.553773974030037, 1.618033988749895, 1.652891650281069, 1.673898962244985 and 1.687530463435423)

(p/q ranges from 0 (red) to 21/22 (magenta) for the groups investigated)

But the boundary is continuous, and is indeed fractal. The "irrational" points on the boundary label degenerate groups, in which one or both sides of the regular set are empty (for example, the black point above). That territory is replaced by points in the limit set:

A singly-degenerate Maskit group (left) and a doubly-degenerate group (right), which is not in the Maskit slice:

........

(μ = 0.7056734968-1.6168866453i, and tr(a) = 1.5-√(3)/2, tr(b) = 1.5+√(3)/2, tr(abAB) = 2 (zmax = 3.732);
in both plots, with infinite time, the white regions will be filled with points in the limit sets)

(μ = 1+i)

Other (1,q even > 2) groups, as well as (n,n+1) for n odd, appear to have similar solutions, as do (9,20), (11,20) (which is conjugate to (9,20)), (13,28) and (15,28) (which is conjugate to (13,28)):

(μ = 0.773301174241798+1.467711508710224i, 1.226698825758202+1.467711508710224i, 0.89512338225965+1.55249182006188i and 1.10487661774035+1.55249182006188i, respectively)

Here is a plot showing the locations of all of the "identified" groups I have found, with the Maskit boundary for reference:

Omitting the (1,4) group, they are fit reasonably well by the equation

4.82572 x8 - 38.6567 x7 + 124.113 x6 - 203.045 x5 + 178.5 x4 - 83.0921 x3 + 20.855 x2 - 3.88996 x + 1.97609
although I suspect that with more data points, this too might be a fractal curve.

# The Riley slice

Using the Riley parameter λ and
method 3 above,
• tr(a) = 2
• tr(b) = 2
• tr(ab) = 2+2λ
• tr(abAB) = 2+4λ2
Any two groups in the Riley slice with the same values of the invariant
tr(ab) - 2
up to sign, are conjugate (Lyndon and Ullman). This implies that the slice is symmetrical about both the real and imaginary axes. In the montage below, I have sampled the region between the origin and 2+i in increments of 0.1 and 0.1i; as above, each sample is centered on its coordinate (zmax drops from 21.05 near the origin to √2, approximately as one over the square root of the distance from the origin):

Samples from the Riley slice:

(Nmin = 424, Nmax = 482,529; color contours range from red = Log10Nmin to violet = Log10Nmax)

The origin is the modular group; the remaining limit sets on the real axis seem to be conjugate to the modular group, at least for rational λ. Most of the region near the origin contain non-discrete groups, with the exception of what appears to be a non-free group at λ=i/2:

(zmax = 3.864)

The group at λ=i is a conjugate of the Apollonian gasket; this is the only group relation between the Maskit and Riley slices:

(zmax = √2)

In the point density plot below, I have explicitly shown the axial symmetries of the Riley slice. The plot covers the region -2 ≤ Re(λ) ≤ 2 and -1 ≤ Im(λ) ≤ 1. The dotted line outlines the boundary of the region containing non-discrete groups. It is not known that the boundary is fractal, but more detailed plots (Gilman) would suggest that this is true. As we have seen above, there appear to be non-discrete groups outside the plotted region.

Lyndon and Ullman showed that groups of the Riley slice are free for all groups outside the solid lines. The white points appear to be non-free groups:

Non-free group at λ = 0.6163849+0.3962686i:

(zmax = 3.836)

# The slice where both a and ab are parabolic

Using method 5 above, tr(b) is the free parameter, tr(a) and tr(ab) are of course 2, and using
tr(abAB) = tr(a)2 + tr(b)2 + tr(ab)2 - tr(a) tr(b) tr(ab) - 2
(Mumford et. al.), we have
tr(abAB) = tr(b)2 - 4 tr(b) + 6
While this sector is algebraically disjoint from the Maskit slice, the groups with tr(b)=2±2i are conjugate to the Apollonian gasket. The group at tr(b)=2 is conjugate to the modular group at λ=0 in the Riley slice.

Here are samples for this region for -1.5 ≤ Re(tr(b)) ≤ 4 and -2 ≤ Im(tr(b)) ≤ 2 (the group at tr(b)=0 is of course non-free, so I plotted the limit set for trb(b)=0.01 instead, highlighted in pink below):

(zmax ranges from 1.138 to 4.562, with an average of 1.718 and a standard deviation of 0.677;
zmax was ∞ for tr(b) = 0.5)

Though non-discrete, the group at tr(b)=1.5 is very interesting: the non-discreteness bifurcates into disjoint maps of the two limits points of b on either side of the imaginary axis, although the axis itself has maps of other group elements (click on the image for a larger version):

Here are the point counts for tr(b) from 1.00 to 1.99:

The group at tr(b)=1.5 is circled in red. All of the groups from tr(b)=1.01 to 1.99 look essentially the same; only the "exposure" changes. zmax varies from 2.637 at tr(b)=1.01 to 4.562 at tr(b)=1.5 to 200.501 at tr(b)=1.99, rising rapidly as the point count drops after tr(b)=1.9 (where zmax=20.512). We see that the lower point counts are associated with larger zmax, as we would expect for plots with a fixed resolution.

The groups for tr(b)∈{-1.5,-1,-0.5,0,1} show a similar bifurcation but appear to be non-free. The limit set of the group at tr(b)=0.5 is a single point at the origin.

Some of the groups in the sample plots are distinctly similar to Riley groups, and the groups at tr(b)=2±i are conjugate to the Riley group at λ=i/2. While the groups at tr(b)=2±i appear to be the same, the one below the real axis was slower to "fill out". The samples appear to be symmetric with respect to reflection in the real axis (the symmetry appears to be "exact" for the discrete groups, although we have to bear in mind the stochastic nature of the algorithm, and the limited resolution, when we contemplate such a statement).

And here are the point density samples for the region -1.5 ≤ Re(tr(b)) ≤ 4 and -2 ≤ Im(tr(b)) ≤ 2:

(Nmin = 0, Nmax = 951,606)

The isolated light points on the real axis represent tr(b)∈{0,0.5,1}, and the red point represents tr(b)=2, which is conjugate to the modular group. There appears to be a series of non-free groups along the line Re(tr(b))=2 approaching the conjugate to the Apollonian gasket.

# The quasifuchsian "slice"

Using method 1 above, tr(a) and tr(b) are free parameters, tr(abAB) is of course -2, and
tr(ab) = ½ (tr(a)tr(b) ± √(tr(a)2 (tr(b)2-4) -4 tr(b)2)
The slice of this region corresponding to tr(b)=2 is the Maskit slice, and the group at tr(a)=tr(b)=2 corresponds to the group λ=i in the Riley slice. The groups at tr(a)=2 and tr(b)=2±i are rotations of the corresponding groups in the slice where a and ab are parabolic.

With two free parameters, this sector of Teichmuller space is 2-complex dimensional, requiring four real coordinates for plots; it is clearly more than a "slice". This is obviously impossible to depict accurately here, so I have elected to show "samples of samples"; each large square represents a point in the tr(a) complex plane, and each small square in it represents a point in the tr(b) complex plane for that value of tr(a). In the montage below, the large squares correspond from tr(a)=-3-3i at bottom left to tr(a)=3+3i at top right; the limit sets in the small squares correspond to tr(b) in the same range (click for a larger image):

(zmax is infinite when tr(b) is zero (and is therefore non-free);
elsewhere it ranges from 1 to 12.815, with an average of 2.272 and a standard deviation of 1.228)

Since we know groups with a generator of trace zero are nonfree, I have offset them by 0.01 and plotted those groups instead. The background colors have the following meanings:

• light red - doubly-cusped groups (gaskets or their conjugates)
• light purple - singly-cusped groups
• light green - non-free groups
• light gray - gasket groups with both halves of the regular set identified
• yellow - regular quasifuchsian groups
• white - nondiscrete groups

There are some interesting "approximate" symmetries, up to conjugation, in the samples:

• the background colors, except at tr(a)=±2, when simultaneously reflected in both the real and imagainery tr(a) axes;
• the limit sets with respect to tr(a)→-tr(a) with a simultaneous 180 degree rotation in the tr(b) plane; and
• with respect to simultaneous reflection in both the real and imaginary tr(b) axes.
The 180 degree rotation is of course equivalent to the pair of reflections tr(b)→-tr(b).
Here are the point density samples for the same regions:

(Nmin = 0, Nmax = 937,227)

The groups in the center of these montages, which were offset to tr(a)=0.01, tr(b)=0.01,±1:

and the group at tr(a)=1, tr(b)=0.01:

are worthy of comment because their limit sets feature concentrations of maps of the fixed points of a (red and green) and the fixed points of b (blue and orange).

# The 2-complex dimensional region where a is parabolic

Using method 4 above, tr(a) is of course 2, tr(b) and tr(ab) are free parameters and, using the equation for tr(abAB) above, we have
tr(abAB) = tr(b)2 + tr(ab)2 - 2 tr(b) tr(ab) + 2
The slice of this region where tr(ab)=tr(b)±2i is conjugate to the Maskit slice, the slice where tr(b)=2 is the Riley slice and the slice where tr(ab)=2 is the same as the slice where both a and ab are parabolic.

As above, we have samples of samples:

(zmax ranges from to 0.762 to 6.423 with an average of 1.872, except for limit sets where tr(ab) ≈ tr(b); for those, it ranges from 141.775 to 412.796 with an average of 252.704)

Here each large square corresponds to tr(ab) varying from -3.01-3.01i at bottom left to 3.01+3.01i at top right, and each small square corresponds to tr(b) in the range -3-3i at bottom left to 3+3i at top right; tr(ab) is displaced to avoid (for this method) the singular case tr(ab)=tr(b).

There seem to be "islands" of non-discreteness located where tr(ab) ≈ tr(b).
The samples are qualitatively symmetric with respect to the simultaneous reflections tr(b)→-tr(b) and tr(ab)→-tr(ab). These symmetries are more clearly seen in the point density samples, here for tr(ab) from -4-4i to 4+4i; each sample contains a 3x3 neighborhood in tr(b) of the value of tr(ab):

(Nmin = 93, Nmax = 948,747)

Looking at the "samples of samples", we see that something very interesting is going on when both tr(b) and tr(ab) are real. Here are samples along the tr(b) real axis for tr(ab) = 0.001 (click for a 12-image sequence showing tr(b) = 0.01, 0.02, 0.03, 0.04, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6 and 0.7):

(blue and orange denote the fixed points of b;
zmax decreases from 222.723 near the origin to 3.446 at 0.7 approximately as 1/Re(tr(b)))

Checking the tr(b) imaginary axis, we again see some interesting limit sets (click for a 12-image sequence showing im(tr(b)) = 0.01, 0.02, 0.03, 0.04, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6 and 0.7):

(zmax decreases from 199.062 near the origin to 3.148 at 0.7i approximately as 1/Im(tr(b)))

Going further along either the real tr(b) axis, or along the imaginary tr(b) axis, the limit sets remain qualitatively similar but slowly "fade into the distance".

Here is the region in the tr(b) plane near the origin for tr(ab) = 0.001, with tr(b) varying from -0.01-0.01i at bottom left to 0.01+0.01i at top right:

(zmax = 2000.5 for the limit set at the tr(b) origin (which, appropriately, is a very small version of the limit sets to either side), and ∈ {135.581, 148.994, 183.229, 199.057, 222.723} for the other limit sets)

Considering what we observed in the quasifuchsian sector, we might be tempted to suggest that such concentrations of limit points for a single group element might be caused by the small values of the traces; for the above plots, tr(ab)=0.001; in a sense, ab is very nearly non-free. But looking at larger values along the tr(ab) real axis between 0 and 1 we see similar phenomena (click for a 9-image sequence):

("teal" and "magenta" bands are actually composites of fixed points of the four commutators;
tr(b) varies from tr(ab)-0.01-0.01i at bottom left to tr(ab)+0.01+0.01i at top right;
zmax = 2000.5 at the limit set at the tr(b) origin for each image, and ∈ {134.905, 148.994, 182.32, 199.057, 222.723} for the other limit sets)

More interesting stuff occurs for tr(ab) = √(3)/2, 1, √2, √3 and √5:

........

........ ........

(zmax = 2001 at the limit set at the tr(b) origin for all but the last image set where it is 2238;
for the remaining limit sets, zmax ∈ {134.905, 141.775, 148.994, 157.761, 159.178, 182.32, 199.057, 200.501, 222.65, 222.723, 223.108, 225.108})

Moving further out the real axis for tr(ab) from 1.1 to 2 (click for a 10-image sequence):

(zmax = 2000.5 at the limit set at the tr(b) origin for each image, and ∈ {134.905, 148.994, 182.32, 199.057, 222.723} for the other limit sets)

Beyond tr(ab)=2, the limit sets are again qualitatively similar but "fade away".

Heading along the negative side of the tr(ab) real axis between 0 and -1 (click for a 9-image sequence):

(zmax = 2000.5 - 1000 * Re(tr(ab)) at the limit set at the tr(b) origin for each image,
and decreases approximately linearly with tr(ab) from the interval (196.121, 321.612) to the interval (142.308, 232.723) for the other limit sets)

On the negative side there seems to be interesting stuff going on for tr(ab) = -1, -√2 and -√3:

........ ........

(zmax = 3000.5, 3413.7 and 3731.6 at the limit set at the tr(b) origin for each image; for the remaining limit sets,
zmax ∈ {202.847, 222.58, 274.138, 298.563, 332.723}, {241.068, 242.484, 339.676, 340.923, 342.923} and {263.543, 264.959, 371.302, 372.706, 374.706}, respectively)

Finally, moving further to the left along the real axis for tr(ab) from -1.1 to -2 (click for a 10-image sequence):

(zmax = 2000.5 - 1000 * Re(tr(ab)) at the limit set at the tr(b) origin for each image,
and decreases approximately linearly with tr(ab) from the interval (270.114, 443.834) to the interval (209.574, 343.834) for the other limit sets)

In all of these "interesting" plots I've been showing you, either tr(ab) is small, or tr(b) is "close" to tr(ab). In the quasifuchsian sector, limit sets with concentrations of specific fixed points occurred when tr(b) was small (i.e., b was "almost" non-free). Lest we think that the concentration of fixed points emerges as zmax gets large, recall that our first such concentration was found in the sector where both a and ab are parabolic, at tr(b)=1.5, where zmax was only 4.562.

# The 2-complex dimensional region where ab is parabolic

Using method 5 above, tr(a) and tr(b) are free parameters, tr(ab) is of course 2, and using the equation for tr(abAB) above, we have
tr(abAB) = tr(a)2 + tr(b)2 - 2 tr(a) tr(b) + 2
This region intersects the Maskit slice when tr(a)=2-2i and tr(b)=2, corresponding to μ=2+2i, which is once again the Apollonian gasket group. The group at tr(a)=tr(b)=2 corresponds to the modular group at λ=0 in the Riley slice. And the slice where tr(a)=2 is the same as the slice where both a and ab are parabolic. The groups where tr(b)=tr(a)±2i are conjugate to the corresponding groups in the quasifuchsian sector. Finally, this and the previous region share a slice where in this region, a is parabolic, and in that region, ab is parabolic.

Again, we have samples of samples:

(zmax = 149.955 for tr(a) ≈ tr(b) = 2; for tr(a) ≈ 0 and tr(b) = 2-2i, it is 100;
for the groups where tr(b) = 1, it ranges from 1.008 to 51.104, with an average of 11.32 and a standard deviation of 9.872; and
for the remaining groups, it ranges from 0.5 to 27.029, with an average of 1.492 and a standard deviation of 1.384)

As before, there seem to be "islands" of non-discreteness located where tr(a) ≈ tr(b). There is also the simultaneous reflection symmetry tr(a)→-tr(a) and tr(b)→-tr(b).
Here are the point density samples for tr(a) from -4-4i to 4+4i; each sample contains a 3x3 neighborhood in tr(b) of the value of tr(a):

(Nmin = 93, Nmax = 948,747, the same as for the previous region, to facilitate comparison)

Looking at the samples of samples, we see once again that something interesting is going on when both tr(a) and tr(b) are real, but now the fixed points of a (denoted as red and green) are as dominant as the fixed points of b. Here are samples along the tr(b) real axis for tr(a) = 0.001 and tr(b) from 0.01 to 0.05 and 0.1 to 0.7 (click for a 12-image sequence):

(zmax increases approximately quadratically from 0.725 to 0.967)

Moving further out the tr(b) real axis, here are samples for tr(a) = 0.001 and tr(b) from 0.95 to 1.05 and 1.1 to 2 (click for a 21-image sequence):

........

(The last plot shows zmax(tr(b)))

Here are samples along the imaginary axis for tr(a) = 0.001 and Im(tr(b)) from 0.01 to 0.05 and 0.1 to 0.7 (click for a 12-image sequence):

(zmax increases approximately quadratically from 0.709 to 0.816 from tr(b) = 0 to tr(b) = 2i)

Here we travel further up the imaginary axis for tr(a) = 0.001 and Im(tr(b)) from 0.95 to 1 and 1.1 to 2 (click for a 16-image sequence):

Here are samples in the region around tr(a) = 0.001 with tr(b) varying from -0.01-0.01i at bottom left to 0.01+0.01i at top right:

(zmax averages 0.7094 with a standard deviation of 0.0017)

Looking at more small tr(b) neighborhoods along the tr(a) real axis between 0 and 1 (click for a 9-image sequence):

(zmax increases approximately quadratically from 0.725 to 0.967;
what appears to be purple is actually intermixed red and blue)

Some interesting stuff happens around tr(a) = 1, √2, √3 and 2:

........

........

(zmax ∈ (1, 1.016), (1.307, 1.339) and (1.932, 2.034) for the first three, respectively;
for tr(a) = 2, zmax = 2032.63 at tr(b) = 0, and ∈ {31.239, 35.844, 36.47, 46.077, 63.02} for the other limit sets)

Moving further out the tr(a) real axis from 1.1 to 1.9 (click for a 9-image sequence):

(zmax increases from 1.055 to 3.581 approximately as a sixth-order polynomial)

Going still further for tr(a) from 2.1 to 2.5 (click for a 5-image sequence):

(zmax decreases from 3.871 to 1.987 approximately as a fourth-order polynomial)

Now moving left of the tr(a) origin along the real axis for tr(a) between 0 and -1 (click for a 9-image sequence):

(zmax increases approximately quadratically from 0.587 to 0.694)

More interesting stuff occurs for tr(a) = -1, -√2, -√3 and -2:

........

........

(zmax lies in the intervals/sets (0.578, 0.583), {263.543, 264.959, 371.302, 372.706, 374.706, 3731.55}, (0.518, 0.525) and (0.5, 0.616), respectively)

Looking further along the negative real axis at tr(a) from -1.1 to -1.9 (click for a 9-image sequence):

(zmax increases approximately quadratically from 0.507 to 0.573)

Finally, looking a little further for tr(a) from -2.1 to -2.5 (click for a 5-image sequence):

(zmax decreases approximately quadratically from 0.665 to 0.574)

Virtually everything we noticed about the previous sector has held true here, with the following exceptions:

• the concentrations of fixed points for specific group elements extends farther along both axes;
• the fixed points of a are here as concentrated as the fixed points of b;
• nothing special appears to happen at tr(ab)=√(3)/2 and √5, but does at tr(ab)=2;
• something interesting happens at tr(ab)=-2.

In all of the "interesting" plots in this sector, either tr(a) is small, or tr(b) is close to tr(a). So far, we have noticed this effect in all of the sectors where one group element is parabolic.

# Trying to put it all together

First, a series of observations:

But having come this far, it seemed to me that I had to try sampling all six dimensions. Using method 1, I drew 117,419 limit sets for tr(a), tr(b) and tr(ab) varying from -3-3i to 3+3i, offsetting each to 0.01 when the trace would have been 0.
The careful reader will note that I should have had 493=117,649 limit sets. 41 of the groups did not have a finite zmax (we'll discuss these later on). And as I noted previously, under certain circumstances method 1 can produce singular matrices. This happened for 189 of the sampled groups, when either

• tr(a) = 2 and tr(b) = tr(ab) (48 instances), or
• tr(a) = -2 and tr(b) = -tr(ab) (48 instances), or
• tr(b) = 2 and tr(a) = tr(ab) (49 instances), or
• tr(b) = -2 and tr(a) = -tr(ab) (48 instances)
(There was obviously some overlap). Groups displaced slightly from these values were densely non-discrete.
While limit sets for those groups could have been drawn using a different method, I elected not to do that so that all of the limit sets I discuss below are consistent.
All of the limit set plots for this section were done using the adaptive termination algorithm.
Here is the point density plot for the entire sampled region, with colors assigned on a logarithmic scale (as before) with Nmin=1 represented by red and Nmax=106 represented by magenta (click on any of the following plots for a larger image or sequence):

The seven major blocks across the top correspond to tr(a) varying from -3+3i to 3+3i;
the seven minor blocks across the top of each major block correspond to tr(b) varying along the same range, and
the points across the top of each minor block correspond to tr(ab) again varying along the same range.

There appears to be an approximate symmetry (again, up to conjugation and the stochastic nature of the algorithm) for the simultaneous reflections tr(a)→-tr(a) and tr(ab)→-tr(ab); on the left is the point density for tr(a)=2+i, and on the right is the point density for tr(a)=-2-i with the tr(ab) reflection:

........

There is also a simultaneous reflection symmetry tr(b)→-tr(b) and tr(ab)→-tr(ab), as we will see in more detail in the next samples.
Aside from the plots for tr(a) purely real or purely imaginary, the plots are all very similar. Here is the "sample of samples" of limits sets for tr(a)=2+i:

The simultaneous reflection symmetry tr(b)→-tr(b) and tr(ab)→-tr(ab) is now more obvious. Also note the limit sets at tr(b)=0.01, tr(ab)=-1+2i and 1-2i:

........

For each value of tr(a), the samples for tr(b)=0.01 are all "dense" except the limit sets at tr(ab)=±(Im(tr(a))-Re(tr(a))i). For the dense limit sets, zmax averages 1.673, while for the "sparse" ones, zmax averages 23.481. This looks like an artifact of fixed resolution, but plotting their conjugates (with zmax values 20.501 and 20.03, respectively) the limit sets are densely non-discrete. We'll look more closely at the small traces values below.

Looking at limit sets with tr(b)≠0.01, zmax for all but 21 were ≤ 100; all but 2,667 were between 1 and 10.
Here is the sample of samples for tr(a)=1:

In addition to the symmetry discussed above, and sparse limit sets at tr(b)=0.01, tr(ab)=±i, we have limit sets along the real axis reminiscent of the "concentrated" limit sets in the sector where a is parabolic. In addition, there are new symmetries:

• with respect to simultaneous reflections in both the real tr(b) and real tr(ab) axes;
• with respect to simultaneous reflections in both the imaginary tr(b) and imaginary tr(ab) axes;
• along the real tr(b) axis, with respect to reflection in the real tr(ab) axis; and
• along the imaginary tr(b) axis, with respect to reflection in the imaginary tr(ab) axis.

All of the samples for real tr(a) were similar.

The black squares above denote limit sets which had no finite zmax. Most of the limit sets with infinite zmax could be plotted when conjugated as discussed above, and most of those were non-discrete, similar to the one on the left below. Several stood out, however: the limit sets for tr(a)=tr(b)=±3i, tr(ab)=2, and tr(a)=-tr(b)=±3i, tr(ab)=-2, while still non-discrete, displayed more symmetry (below, right):

........

(zmax ranges from 33.45 million to 381.1 million for the 34 limit sets whose conjugates had finite zmax)

Here are the samples for tr(a)=i:

Here we have limit sets along the imaginary axis reminiscent of the ones in the sector where ab is parabolic. And in addition to the tr(b)/tr(ab) reflection symmetry, we have new symmetries mimicing the ones for real tr(a):

• with respect to simultaneous reflections in both the real tr(b) and imaginary tr(ab) axes;
• with respect to simultaneous reflections in both the imaginary tr(b) and real tr(ab) axes;
• along the real tr(b) axis, with respect to reflection in the imaginary tr(ab) axis; and
• along the imaginary tr(b) axis, with respect to reflection in the real tr(ab) axis.

All of the samples for purely imaginary tr(a) were similar.

We noticed in the sectors where a and ab are parabolic, that there were "islands" of non-discreteness where the other two traces were equal. We seem to see something similar in the samples for tr(a)=2, and a complementary effect for tr(a)=-2:

........

although at our course sampling it is a little hard to see. So I plotted point densities at 4 times the resolution to clarify the situation. These plots cover the same trace ranges but at increments of 0.25 and 0.25i:

.......

As expected, the first plot is essentially a higher resolution version of the point density plot for the sector where a is parabolic. For the second plot, a is still parabolic, but with trace -2; now the islands of non-discreteness occur when tr(ab)≈-tr(b).

Each of these plots contain 254 = 390,625 points, requiring the computation of one limit set per point. The first plot took approximately 106 cpu-days to complete, while the second took approximately 98 cpu-days to complete. Both plots include null points for 641 limit sets producing singular matrices. For every value of tr(b), the set with tr(ab)=sign(tr(a))tr(b) was singular (625 sets). In addition, there were 16 sets where tr(b) was the same for both plots, but tr(ab) differed by sign.

Curious about how the limit sets behave for small traces, I examined limit sets for tr(a)=tr(b)=tr(ab) = 10-6 to 1, increasing in factors of 10:

(The three limit sets shown above are for 10-6, 10-3 and 10-1; click on the image to see the whole sequence.)

As the traces increase by a factor of 10, zmax decreases by a factor of 10, ranging from 5,463,615.3 for traces equal to 10-6, down to 56.77 for traces equal to 0.1. The last frame in the sequence is for traces equal to 1, where zmax=7.6; there are only 12 points plotted in that limit set, so I have circled them to mark their positions. There is an intriguing self-similarity displayed in these plots: even though the traces change by a factor of 10 for each plot in the sequence, the fixed points are in the "same" places; it's as though increasing the traces increases the "exposure", but not the relative positions. For traces less than 0.01, all of the limit set points are either blue or orange (corresponding to fixed points of b) except the point at the origin, which is red (corresponding to a fixed point of a). For traces of 0.01 and 0.1, fixed points of b again dominate the limit sets, but fixed points of different operators appear, mostly concentrated near the origin (the behavior near the origin is probably a side-effect of the resolution.) For the last image (traces=1), the points (from top to bottom and left to right) correspond to fixed points of

ABab
a a
a
abAB AbaB abAB AbaB
b
a a
ABab

Taking a broader look at the effect of smaller traces, here is a sequence of samples of samples for tr(a) varying from 0.01 to 10-5 (the last is shown here):

(zmax < 100 for all but 22 of the limit sets in all of the samples, but
the maximum value of zmax is in the hundreds of thousands for all of them;
the mean zmax for the ones <100 were between 4.78234 and 4.79333)

Again we see the self-similarity, across four orders of magnitude in tr(a).

So what have we learned?

• When one trace is held fixed, there seems to be a simultaneous reflection symmetry in the other two traces.
• When any group element is "close" to being non-free (i.e., its trace is "small"), there is an "island" of non-discreteness, which of course delineates a boundary of the Teichmuller space.
• When one group element is parabolic, there are additional islands of non-discreteness where the other two traces are "close", modulo the sign of the parabolic element.
• When one of the traces is purely real (or purely imaginary), there are additional reflection symmetries with respect to both the real and imaginary axes for the other two traces.
• Boundaries in Teichmuller space seem to be fractal.
• Teichmuller space near the origin seems to be self-similar with respect to dilation of traces.

So as a first approximation, it might be possible to imagine the whole 6-D Teichmuller space for genus 2 as a subset of C3 with regions removed "near"

• where any complex coordinate is zero;
This would appear to divide the Teichmuller space into disjoint components.
• where one complex coordinate is ±2 and the other two complex coordinates are equal (up to the sign of the parabolic trace).
This would appear to have a similar effect. Consider tr(a)=2; the equation tr(b)=tr(ab) defines a 3-dimensional hyperplane which again divides the space into disjoint components.

The Teichmuller space for genus 3 spans six complex dimensions. Since there are nine ways to pinch the genus 3 torus:

one must pick six of the associated traces to parameterize the space.

Kleinian groups for genus 3 have three independent generators. If we set
a = ((1 + 2i, 4), (1, 1 - 2i)),
b = ((1 + 1.5i, 1), (2.25, 1 - 1.5i)) and
c = ((-13 - 4i, 8i), (-12 + 16i, 11 + 4i)),
the limit set is a gasket:

(This group corresponds to Mumford et. al., figure 11.6; the limit set is conjugate to their figure 11.7.)

In this group, a and b have trace 2, and c, aB, bcBC and cbCB have trace -2, so the torus is maximally pinched and therefore a gasket. The topology is four triply-punctured spheres, and the regular set partitions into four disjoint subsets, each corresponding to one of the triply-punctured spheres (each subset is colored differently in the plot above).

For genus 4, the Teichmuller space has nine complex dimensions. There are now 16 ways to pinch the torus:

In general, for genus g, there are g (outer loops to pinch) + g!/(2(g-2)!) (pinches from one loop to another) + g!/(2(g-2)!) (pinches between "corners") = g2 ways to pinch the torus, so one must choose 3g-3 of g2 traces to parameterize the space. This in turn means that there are

(g2)! / ((3g-3)! (g2-3g+3)!)
ways to parametrize the space.

• For g=2 there were four ways to choose three complex traces (which trace of the four to leave out);
• for g=3 there are 84 ways to choose six complex traces;
• for g=4 there are 11,440 ways to choose nine complex traces;
• for g=5 there are 5,200,300 ways to choose 12 complex traces; and
• for g=6 there are 5,567,902,560 ways to choose 15 complex traces.

Let's stop there.

# Appendix: Mathematics as Art

It is by now well-appreciated that mathematics can be a source of visual (as well as auditory) beauty. In many, perhaps most, cases, symmetry plays a fundamental role; one need only look so far as the Appolonian gasket for a prime example.

But equally interesting are some of the more chaotic limit sets. Some of my favorites are:

"Starry Day"

(tr(a) = 0.00001, tr(b) = -1, tr(ab) = -2-2i, zmax = 1.221283, N = 168,554)

(tr(a) = 0.00001, tr(b) = -2, tr(ab) = -3i, zmax = 1, N = 122,772)

"Twin Twisters"

(tr(a) = 0.00001, tr(b) = -3-3i, tr(ab) = 3-3i, zmax = 1.407854, N = 84,185)

"Behind the Eyes of a Robot"

(tr(a) = 0.00001, tr(b) = -2i, tr(ab) = -1, zmax = 20.429912, N = 233,586)

"Christmas Ornaments"

(tr(a) = 0.00001, tr(b) = i, tr(ab) = -2, zmax = 1.618036, N = 6,624)

"Beyond the Infinite 1" (with thanks to Stanley Kubrick)

(tr(a) = 0.01, tr(b) = 0.01, tr(ab) = 3, zmax = 1.009256, N = 812,247)

"Beyond the Infinite 2" (likewise)

(tr(a) = 0.01, tr(b) = 1, tr(ab) = -1, zmax = 7.522355, N = 87,540)

"Chromosomi" (with thanks to Enrico Rava)

(tr(a) = 1, tr(b) = -3-i, tr(ab) = -3, zmax = 1, N = 7,234)

In some of these images I have changed the background color for artistic effect, but they are otherwise unadulterated limit set plots. In the following image, I duplicated the plot with reflection symmetries to create

"Play Ball!"

(tr(a) = 0.01, tr(b) = 2, tr(ab) = 3i, zmax = 1.005981, N = 655,168)