**Kleinian and Quasi-Fuchsian Limit Sets: An Open Source Toolbox
Genotype:
1.1.49 PDF**

The Intrigues and Delights of Kleinian and Quasi-Fuchsian Limit Sets

Mathematical Intelligencer 2020 Republished in: The Best Writing on Mathematics 2020 Princeton Univ. Pr.

Chris King

**Download** the complete **Matlab toolbox, C scripts and Mac app** with web resource Toolbox contents page.

**Full live linked pdf of this article**

This page contains an in-depth article with a Mac app, a cross-platform C implementation and a Matlab toolbox that depict limit sets of Kleinian and quasi-Fuchsian groups. In all the depth-first examples I have used the c-script for its speed and efficiency, using the matlab functions to aid analysis, provide trace parameters and run stochastic approximations when dealing with intractable situations.

**Contents**

- Overview of the Computational Package
- Diverse Mathematical Foundations
- The Depth First Search Algorithm
- Rational Maps and Period Counting
- Finding the Elusive Traces
- The Maskit Slice as a Fractal Solution Space
- Endless Varieties among Limits of Limits
- Lost in Four Dimensional Space
- Finding Conjugate Limit Sets
- Rational Periodicity becomes Irrational Degeneracy
- The Total Eclipse of Double Degeneracy
- The Spiral Nemesis
- The Holy Grail at the Edge of Chaos
- Discreteness Enters the Chaotic Milieu
- The Non-Free Case and Automatic Group Tables
- Generating Finite State Automata
- Symmetric Normalisations and Hyperbolic Planes
- The Stochastic Algorithm
- Escape-time and Extra Dimensions
- Program Listings

**Overview of the Computational Package**

The package provides a computational toolbox in an OS-independent form which is freely available to those seeking to understand these systems and explore them, including the key examples in "**Indra's Pearls: The Vision of Felix Klein**" by David Mumford, Caroline Series and David Wright, Cambridge Univ Pr. 2002 ISBN 978-0-521-35253-6. The images seek to fully explore and elucidate the underlying mathematical systems, but anyone seeking to produce fractal art via 3D or ray-traced versons could easily extend the toolbox algorithms in Matlab to do this.

The depth first search algorithm, derived from Indra's Pearls, is fully articulated in the** ****Matlab function**, which is also compiled in a fully functional native gcc **C script**, which is a 'smoking' hundred times faster. The C version is also encapsulated in the Mac OS app **DFS Viewer**. The Matlab toolbox includes supporting functions, has the capacity to find all the coefficients of high degree trace formulae and easily address new research questions. However all features can also be explored without Matlab, using DFS Viewer on Mac and the key native C scripts on any gcc C compatible platform. This makes the algorithm widely and freely available to the general public, as well as maths students and researchers.

The Mac app **DFS Viewer** has a variety of controls to aid successful viewing, most of which are also included as parameters in the **C script** and **Matlab** versions:

The **C script** has been configured to be generic and universal – able on any platform to generate very high resolution colour or B&W images and output them as .ppm files which are easily converted to png or jpg formats, so for all users, this is the most accessible and time efficient way to generate high quality images. It is 93 times faster on the right-hand example in fig 7 (15.2 secs vs 1422 secs for levmax = 1000 and epsilon = 0.0001). It can also be easily adapted to generate movies, as batch files, and has a timing meter, based on the successive tree branches of the first four search levels, which gives a good indication of the progress of a high precision computation.

A list of key fractional and degenerate traces explored in this article is also provided, which can be input into any of the versions. Also included are Matlab and C script versions of **tracepq** and **tracepqr** which can be used to find traces of new examples.

The as well as a variety of suplementary functions aiding further research, the Matlab toolbox also includes a stochastic algorithm which can depict chaotic cases, a colour filling escape-time algorithm, and a 3-D example.

**An Intersection of Diverse Mathematical Areas**

Kleinian and quasi-Fuchsian limit sets form an intersection of diverse branches of mathematics, becoming a marriage of group theory, the topology of handles and pinched manifolds formed by the transformations on the Riemann sphere, hyperbolic geometry, fractal tilings, and the theory of chaos and fractal dimensions

Kleinian and Quasi-Fuschsian limits sets are generated by a pair of Mobius transformations in the complex plane, *a *and *b*, along with their inverses *A = a *^{–1}* *and *B* = *b*^{ –1}. Mobius transformations are the analog of affine Euclidean transformations (translations, scalings, rotations) on the Riemann sphere (1 in fig2), the complex plane with a single point at infinity added, forming the north pole. Mobius maps are fractional linear transformations *m*(*z*) = (*pz*+*q*)/(*rz*+*s*), which compose by matrix left-multiplication:

Mobius transformations, as noted, correspond to complex versions of the rotations, scalings and translations of affine Euclidian maps.They can be decomposed into similar elementary transforms and are classified into elliptic, parabolic, hyperbolic and loxodromic types. Parabolic transforms are like transalations and have only one fixed point on the Riemann sphere, like the point at infinity of a translation. Loxodromic maps are complex dilations. They have a repelling and an attracting fixed point and the map expands in a spiral form from the repeller to the attractor. Elliptic maps are conjgugate to pure rotations. Pairs of Mobius transformations induce a group structure of transformations under composition. The limit sets consist of the asymptotic points the rest of the Riemann sphere is eventually mapped onto under the group transformations. Conjugacy powerfully simplifies the classification, because any map *M *=* C*^{-1}*NC* generates isomorphic behavior, so we can regard the groups as equivalent and reduce the analysis to conjugacy classes, or convert our mappings into convenient conjugate forms. The interplay of the group's generators *a*, *b*, *A*, *B *become defining words in the construction of the limit set algorithm, were the exact action of each generator on the limit set becomes apparent, as in fig 1(right) and fig 4. The traces of parabolic maps, defined by *p* + *s* above (the sum of the diagonal elements) lie on the real line between –2 and 2. Those of elliptic maps lie on the unit circle excepting 1 and –1. Traces of hyperbolic maps lie on the real line outside [–2,2] and loxodromic maps lie anywhere in the complex plane outside these ranges, although hyperbolic are sometimes included in loxodromic.

A brief summary of some of the mathematical ideas underlying these limit sets is illustrated in fig 2. Mobius mappings preserve circles and lines in the complex plane Lines are simply circles going through the north pole at infinity. When two Mobius transformations are combined, the actions of the matrices and their inverses generate a symmetry group, just as the 3 rotations and three reflections of an equilateral triangle form a symmetry group. These groups can take varying forms. If we take two disjoint pairs of circles and map the inside of one onto the outside of the other on the Riemann sphere, we also generate a topological transformation of the surface as shown in (e), generating a 2-handled torus – the pretzel. When we apply these two transformations to the surface, successive images of the pairs of Schottky circles become mapped fractally as shown in (b) forming a fractal dust or Cantor set – the limit set of the group. As the circles come together and kiss or touch as in (c) the limit set can become a circle. If the commutator is parabolic wiith value -2 the double torus becomes pinched as in (f) grandma's recipe can be used to generate a quasi-Fuchsian group consisting of a fractal quasi-circle forming a topologically closed loop or Jordan curve. The set containing the converging circles in fig 2, on which the process makes a discrete tiling is called the **ordinary set** and the complementary asymptotic set the **limit set**. These transformations generate a quasi-Fuchsian limit set, when it is a Jordan curve, which is Fuchsian if the limit set is a geometrical circle as in (c).

The topology of the surface can undergo other pinching and puncturing transformations as in (g), resulting in Kleinian groups such as the Apollonian gasket (partial Schottky computation in (d) and full limit set fig 1 (left), which is a limiting case of neighbouring quasi-Fuchsian sets illustrated in fig 1 left, as the traces* *tend from 2.2 down to 2. The groups generate tiling symmetries on the hyperbolic plane as illustrated by Maurice Escher (i) and the symmetries become reflected in the forms of the circle packings generated in many of the limit sets where the second matrix trace is 2. The limit sets are also intimately connected with hyperbolic geometry as examplified by the modular group of symmetries illustrated in the Escher print (i), with topological pinching (h) also echoed in the fractal structures of the circular limit sets were *TB* = 2.

In fig 3 one can see the effects of changing *Ta* and *Tb* from the 2, 2 of the gasket to complex values. While *Ta* alters the inner region to form spiral fronds, as in fig 1 right, *Tb* transforms the outer circle into exterior fronds. When both have complex values, the interactive effects can become unexpected, as in the case of + / + (left), which is not simply a partial reflection of the centre figure and is also about ten times larger. Changing the signs of real and/or imaginary parts yeild four symmetry versions of each figure, consistent with the symmetries of the Maskit plane and the double roots for finding *Tab*.

**The Depth First Search Algorithm – **Program listings in colour**: **(a) C script (b)** **Matlab

The Depth First Search algorithm is the quintessance of depiction methods of quasi-Fuchsian and Kleinian limit sets and is the key instrument for studying the diverse key examples. The DFS algorithm does a full depth search of all the key examples in the book, including the non-free cases, which require an automatic group schedule to take account of degenerate words. It includes intelligent line linking using the attracting fixed points of the commutators and special words to take account of the Apollonian Gasket and other examples with whorls with narrow necks. The essential inputs are the traces (*p*+*s*) of the two generating Mobius transforms *a* = [*p q*;*r s*], sending *z* to (*pz*+*q*)/(*rz*+*s*), plus other parameters to control the accuracy of the plots, noted in the discussion.

The gallery of high resolution 1100 x 1100 examples demonstrate the quality of the algorithm, with examples whch test it to the limits of degeneracy, 'accidental' coincidences of double cusps, and situations in which the group is not free due to relations between the generators.

The DFS algorithm has an ingenious strategy designed to work with transformations with parabolic commutators that produce quasi-Fuchsian groups which can be scanned in a circular manner by exploring the tree of words defined by compositions of the transformations.

Firstly "grandma's algorithm" is used to generate 2x2 matrices for Mobius transforms with parabolic commutators (*aba*^{-1}*b*^{-1} etc) from the two traces, *Ta* and *Tb *that are unique up to complex conjugacy. In the case of the Apollonian gasket, both *Ta* and *Tb *are 2 so *a* and *b* are also parabolic.

One can see how grandma's algorithm is extended by reading the subfunction grandma() in the program listing.

The depth-first search algorithm to depict the limit sets is ingenious and extremely elegant. The aim is to traverse the algebraic space of all word combinations of the generators *a*, *b*, *a *^{-1}, *b*^{ -1} in a way which draws a continuous piecewise linear approximation to the fractal of any desired resolution. The algorithm does a depth-first exploration of the tree of addresses of all possible iterations of the transformations, proceeding around the addresses in a manner that ensures we steadily traverse the maze of address space and in particular, the fractal topological quasi-circle of the limit set, by ignoring cancelling paths such as *aa *^{-1}, storing the matrices of composed transformations going down the depth search, so these don't have to be re-calculated when we back up part way to explore another vacant branch. The cancellations are made by always going down the tree corkscrew-wise, making a right-hand turn after each level descent, and then exiting the level when we reach the inverse of the transform we entered the level by after three left turns to explore the other tunnels. It makes a branch terminal check by applying the composed transformation to the fixed point of the appropriate commutator defined by the last transform e.g. *t*(1)**t*(2)*...**t*(*n*-1)*a operates on the attracting fixed point of *bA*^{}*B*^{}*a*, which effectively defines the infinite path composed of the *n* terms and the infinite repeat of the right-hand commutator. It then compares this "new" point with the last "old" point found the same way (equivalent to using the fixed point of the opposite commutator *BAba*) and if they are less than *ε* apart draws a line between. This approach neatly overcomes a critical flaw in random iteration algorithms, that visit some parts of the limit set exponential rarely and thus effectively never get drawn in a finite time. This means that it doesn't need to explore further down the tree once it has a satisfactory linkage of the scale we are depicting, avoiding exponential runaway. Nevertheless on a dual core intel Mac at a resolution of *ε* = 0.0001, each image takes around an hour to produce, although moderately lower resolution images will display in a matter of seconds or minutes.

Special words are introduced to take account of two anomalies, Firstly in the Apollonian gasket example to close the top and bottom circles and their fractal offspring. Secondly, this process is also invoked to avoid pinching off the whorls of limit sets with large appendages with narrow necks as in fig 4. To do this, in addition to the attracting fixed points of the four sets of commutators we also include those of the original four generators, which, as can be seen from the first image in fig 1 and the image on the right from Indra's pearls, are at the very tips of the whorls. Fig 4 illustrates the problem for the second example in fig 3, where two "new" and "old" candidates at the neck are *ε-*close, but the "tip" of the whorl is further away and joining them would cut the whorl at the neck. This means we can then decide to fill in a line only if the distances from our "new" and "old" points to the local whorl tip we call "mid" are both less than *ε*. In the case of the gasket we fill "new" to "mid" to "old" as the 'tip' is right in the centre of the gap, but for other sets just "new" to "old" provided both *ε*'s are small enough.

Further examples aided by special words to get vey high resolution images efficiently arise when dealing with the cusp examples below, where other elements are chosen to be parabolic, such as *a*^{15}*B* in the 1/15 example. We then proceed by making a table of all the parabolic elements for each generator, including the commutators and joing "new" and "old" points only if all the elements arising in the tree exploration sequence between these two are within *ε*. These examples are not included as they would each apply to only a single case and the algorithm is already giving high reolution images of these in real time.

The many non-free examples in figs 19 – 22 go further and have generators connected by an algebraic relation such as (*aba*^{-1}*b*^{-1})^{n} = I, which couples them and causes many of the words in the search tree to become null, so we have to set up an automaton generating the table of an automatic group to correctly portray the set. This requires a new array "state" which iteratively generates successive word states by recursively entering the lookup table FSA, forbiding entry to tunnels involving any of the null word combinations, so we don't endlessly chase our tail and get nowhere fast.

I have so included a function specialwords() in the Matlab implementation so that anyone who wishes to set up a special words directory for other more complex examples can do so easily, as the plotting routine is already included to process the special word extensions, so all that is needed to be done is to replace the array entries in specialwords() with the ones one has defined from the permutations of the parabolic element. As an example in fig 6, I have included IPs 8.16 designed to show the method in the case of *aaB*. This ends up with the look up structure shown in fig 6 left, while the image on the right illustrates two levels of resolution of the limit set algorithm. The words in the tree maze can take multiple forms due to equivalences and cancellations of inverses. For example the commutator *ABab* is equivalent to *BAba* because its single parabolic fixed point is the same as its inverse. This means that a word address like *aaBaab*[*ABab*] = *aaBaab*[*BAba*] = *aaBa*[*abAB*]. One also has to bear in mind these examples apply only to one example each and there is a significant time cost in checking a large number of extra cases for higher denominator words, and their many permutations all having to be below threshold drives the computation further down the exponentiating levels of the search tree. The example above needs some final tuning in its parameters, indicated by irregulatiries in the low resolution image, but otherwise, the Matlab implementation gives a faithful rendition of a full spread of examples of limit sets in the book.

**Rational Maps and Period Counting**

Group theorists take a special interest in the spectrum of properties of limit sets where one trace factor is 2, both because they provide a space of examples which can be analysed and because many other limit sets can be shown to be conjugate to these, so below are several key examples.

In fig 7 you can see examples of rational cusps 1/15, 2/19 and 7/43, with the central chain of circles numbered to highlight the relationship with the denominator. In the left figure the fixed point of the parabolic element *a*^{15}*B* is in green and circles 1 and 2 are thus mapped to themselves. B carries C1 to C16 and then *a*^{15} carries it back to C1. In a similar way, B carries C2 to C17 and then *a*^{15} carries it back to C2. This means that *a*^{15}*B* maps both the disks C0 and C1 to themselves, the fixed points of *a*^{15}*B* must thus be on both these circles meaning there is one fixed point (parabolic) on the tangent point. b is parabolic with fixed point -*i* and pushes points out from its fixed point along clockwise circular trajectories. a is loxodromic and pushes the disks in the numbered spiral chain from right to left across the shrinking and spiralling into the sink on the left.The chain of 15 circles connecting to the opposite large circle 16 map to one another and the mauve circles are in turn emerging from the action of loxodromic *a* in a spiral from the repeller and proceeding in the orange numbered spiral to the attractor. In the 7/43 case, you can also see the way the chain is running in cycles of seven, corresponding to the numerator, and cycles of 2 in the 2/19 case.

These mappings give a good idea of the discrete action of the underlying group, where the isometries map discs to discs so that the images of a given point under the group transformations form a discrete set and the images of the discs effectively form a hyperbolic tiling of the Riemann sphere. In non-discrete cases, as in fig 25 and fig 26, where the stochastic algorithm is simply displaying the orbits of one or more points on the limit set under random sequences of transformations, one can see that the images of such points will become chaotically mixed, so that the isometries of a given point are no longer discrete sets. Any discrete group introduces a tiling of the plane, or higher dimensional space e.g. with the annular tiles exemplified by the hyperbolic tilings in fig 2(2,3).

More accurate plots of the small circles can be gained by including the parabolic elements such as* a*^{15}*B* for 1/15 or *a*^{3}*Ba*^{3}*Ba*^{2}*Ba*^{3}*Ba*^{2}*B * for 5/13 as special words, along with all their permutations, linking "new" and "old" only if all of the steps in the chain between, involving these parabolic critical points are within *ε*.

To explore a given type of limit set we need to be able to calculate thre traces which will have particular properties and periodicities in terms of the generators. We can work out the correct trace of parabolic words recursively from *T*(*p*/*q*) of lower denominator fractions, as they combine in the manner of mediants and Farey tree.

The traces combine in two key relations: (1) The "grandfather" identity *Tmn* + *Tm ^{-1}n* =

Key to this is the nature of any composition of the generators which is parabolic. As an exercise in finding the trace coefficients of some of these examples lets look at the cases 1/*n* where the parabolic element is *a*^{n}*B*. If we are not trying to use the special words algorithm at this point, we only need to find the correct value of *Ta.* Using the extended Markov and grandfather identities, we can immediately see

*–2 = TabAB* = (*Ta*)^{2} + (*TB*)^{2 }+ (*TaB*)^{2}–* Ta.TB.TaB* – 2, or (*Ta*)^{2} + 4^{ }+ (*TaB*)^{2}–* 2TaTaB *= 0, from which we can deduce that *TaB *= *Ta* + 2*i* . Then using the grandfather identity on *a* and *B *, with *TB* = 2, we get a series of equations *Ta *^{n }*B* = *Ta*.*Ta*^{ n-1}* B* –* Ta*^{ n -2}* B*, making it possible to write a recursive process to define *Ta* as an *n*th degree polynomial. We can do this readily in Matlab and at the same time solve the polynomial numerically by a Newton-type method as follows:

function s=recurscoeffs(stps,opt)

syms ta

a=cell(1,stps+1);

a{1}=2; % This is
TB

a{2}=ta+2i; %This is TaB

for i=3:stps+1

a{i}=expand(ta*a{i-1}-a{i-2}); % Ta^(i)B = Ta.Ta^(i-1)B-Ta^(i-2)B

end

f=symfun(a{stps+1},ta); % Make it a function

if opt==0

s=eval(vpasolve(f==2,ta)); %Solve to get the roots of f=2

else

s=eval(vpasolve(f==-2,ta));

end

This will give a list of *n *roots. Theory confirms that among all possible solutions to the trace equation (for a fixed p/q), exactly one will be discrete, so only one of these will have a credible root one can try. You can directly input the values such as s(7) into depthfirst as parameter, having run this function, giving full precision. Example values are 1/7 1.846276 – 0.09628*i*, 1/15 1.958591 – 0.01128*i*, 1/23 1.98179 – 0.00320*i*, 1/31 1.98987 – 0.00131*i*, and 1/60 1.99726820 – 0.00018236, tending towards *Ta = *2, (or *μ* = 2*i* ) in fig 11 centre.

Going beyond the 1/*n* fractions to *p*/*q*, we can also calculate the trace *Ta* for general fractions, because parabolic expressions like *a*^{3}*Ba*^{3}*Ba*^{2}*Ba*^{3}*Ba*^{2}*B* can also be generated recursively e.g. from *a*^{3}*B* and *a*^{2}*B* by similar relations. To be more precise, examination shows that the word *w*_{(p+r)/(q+s)} = *w*_{p/r}.*w*_{q/s}, in the manner of mediants and the Farey tree, so we can use the grandfather identity as follows:

*T**w*_{(p+r)/(q+s)} = *Tw*_{p/r}.T*w*_{q/s} – *T*(*w*_{p/r}^{-1}.*w*_{q/s}). These relations simplify in rather beautiful ways if we have generated a list of Farey tree predecessors. For example

*Tw*_{(3+5)/(8+13)} = *T**w*_{8/21} = *Tw*_{3/8}.T*w*_{5/13} – *T*(*w*_{3/8}^{-1}.*w*_{5/13}) = *Tw*_{3/8}.T*w*_{5/13} – *T**w*_{2/5}, because *w*_{3/8}^{-1}.*w*_{5/13} = (*a ^{3}Ba^{3}Ba^{2}B*)

The Fibonacci fractions used to illustrate this above are are all linked as immediate neighbours in the Farey tree (right). Now it turns out that any other pair of fractions that are Farey neighbours also have generator word sequences which will correctly cancel in the same way as above. Thus to find the correct trace polynomial for an arbitrary fraction, we simply need to determine the path down the Farey tree from the precursors to get to the fraction. So the function **tracepq()** firstly calculates the sequence of left and right handed moves to get from 0/1 and 1/0, (which represent *Ta* and *TB* = 2 respectively, with 1/1 corresponding to *TaB* = *Ta*+2*i *) to the particular fraction whose trace we are searching for, using the subfunction **fareyseq()**. It then builds up the trace of *p*/*q *iteratively using the grandfather identity. I previously had a much more cumbersome algorithm recursecoeffs() using recursion on arbitrary splits, which is still listed in the toolbox. The trouble with the top-down recursion is that you have to find a split which avoids degenerate fractions and is in Farey neighbours. Consequently some larger fractions failed to solve out correctly and the recursion is liable to exponential runaway.

A list of some fractional traces.

**Coefficient recursion** tracepq() program listing in colour

A sample listing of tracepq() for 2/7 iis shown below. The 7 roots are listed in increasing order of their real part, so that one can easily search the lowest few roots for the one which portrays correctly using depthfirst().

>> s=tracepq(2,7,0,1);

2 -1.53746746-0.21870978i

5 -1.40222578-0.58650847i

4 -0.39303611-0.25387457i

3 -0.38479895-1.27170824i

1 0.79526188-1.15961696i

6 1.25146576-0.06152298i

7 1.67080064-0.44805900i

We can then follow with a portrait using depthfirst()

depthfirst(s(7),2 ,0.0001,1000,1,3,0,1100);

The parameters in depthfirst select s(7) the 7th solution listed, 2 for *Tb,* followed by 0.001 accuracy scale 1000 maximum tree depth, 1/3 image scaling, sav = 3 both prints progress (≥2) and saves the image (mod(2)=1), option 0 is for a standard group without special words for the generators (see listing), and 1100 is the image size. Printing progress means you can see a series of stages of exploring the tree so you know progress is taking place for a detailed plot that may take an hour to complete. The report lists the time and olace of each time it falls back below level 4, and if the process reaches levmax it prints a dot. This can happen either because the limit set has issues or because it is larger than the scaling. Ideally one wants levmax large enough so that a nice curve on screen doesn't reach it for a given epsilon. One can also use sav=-1 as a test mode which will report whether the process reached levmax by setting outfl false. This will happen with a bad limit set and cause the process to strand, so this mode can be used within a function to test for a good limit set. The function scantraces() scans all the trace solutions *i*/*n* for a given *n* and automatically finds the one correct root which will produce a low resolution image and declare outfl=true to provide all the *i*/*n* solutions automatically. These can then be plotted in sequence using batchlimsets().

**The Maskit Slice as a Fractal Solution Space**

In figs 9 and fig10 below are shown the Maskit slice, consisting of all solutions to these polynomial equations for the trace of the commutator being both 2 and -2.

A group in which no simplifications other than cancelling out adjacent transformations and their inverses (*aA* and *bB*) is called free. Groups for which large products stay away from the identity I are called discrete. In the non-discrete, case *m*(*z*) can eventually come very close to *z*, resulting in a chaotic recurrent orbit which prevents the areas outside the limit set being consistently tiled by the transformations, resulting in the Cantor dust limit sets we can see using the stochastic algorithm. We are here looking for the double-cusp groups which hover on the boundary between discrete and non-discrete. The points in the white upper region all represent discretegroups, while 'most' points in the shaded region are not discrete. The upper (and lower) cuspy curve is called the Maskit boundary. *μ*-values on the boundary give groups right on the transition between order and chaos.

Above are shown all the fractional rotation states modulo 29, the Maskit slice up to *q* ≤ 51 generated using the script slicec, which collects all the roots found by tracepq(), saving them in mus.mat The conjugate roots constituting parabolic solutions can then be plotted using the function plotmus() . The "Maskit slice" is the plane of *μ*(*p*/*q*) highlighting all polynomial roots where the generating *p*/*q* word, such as a^{15}B for *μ*(1/15), is parabolic. In the Maskit parametrization, whose limit set corresponding to *Ta* = 1.91 + 0.1*i*, *Tb* = 2 is illustrated in fig 23, we have one parameter *μ* defining a Mobius transformation which is paired with translation by 2, which is parabolic with fixed point ∞. The two generators arise immediately from the parameters: *a*(*z*) = *μ* + 1/*z* and *b*(*z*) = *z* + 2. In the standard parametrization we get the same limit set in conjugate form from *Ta* = –*i**μ*(*p*/*q*), *Tb* = 2. We will thus refer to such *Ta* simply as *μ*(*p*/*q*)).

Theory shows that no value of *μ* with imaginary part between –1 and +1 could give a discrete group, while values with imaginary part more than +2 or less than −2 would always be discrete. The points in the white upper region all represent discrete groups, while 'most' points in the shaded region are not discrete. Some isolated *μ*-values in the shaded central region give discrete groups, but such groups are never free. For example *Ta *= 1 and *Ta* = 2^{1/2} , which emerge as solutions to *μ*(1/6) and *μ*(1/8) give non-free groups with *a *^{6} = I and *a* ^{8} = I respectively and isolated limit sets, fig 24, surrounded by non-discrete values.

**Endless Varieties of Limits of Limits**

In fig 11 are shown a sequence of limit sets of 1/*n* tending in limit to 1/i which should theoretically be the Apollonian gasket. Note that in the right superimposed figure for 1/31 in fig 11, the gasket coincides with the upper and lower circles and their fractal subcircles center. Note as well, the dark circles which have emerged in both 1/23 and 1/37, and are present as cicles interrupted by the central chain of circles in 1/15 (fig 8 left), forming a chain between the fractal loxodromic sources and sinks that lies in the limit along the large Apollonian circles. The limiting limit set results in a four-fold rotational symmetry, rather than the two-fold symmetry of the gasket.

If we examine the matrices of the gasket arising from grandma's recipe we have:

Hence we can see that (*a ^{n}B*)

One can take the limit of 2/(2*n*+1) and derive yet another distinct limit limit set which also has fourfold symmetry but has pairs of middle-sized circles on each diagonal rather than a single one, as shown left for *μ*(2/121). And we could do this for 3/(3*n*+1) and 3/(3*n*+2) and so on so we have an infinite number of different limits of limit sets, all tendi*ng to μ*(0/1).

**Lost in Four Dimensional Space**

When we move away from *Tb *= 2, and search for *Ta*, we no longer have a 2D Maskit slice of the parabolic generator complementing Tb = 2, but essentially a 4D space, in which we can't find all polynomial roots as before because the grandfather identities no longer produce polynomials. Fig 13a and discussion shows a way around this problem.

The two sets in fig 13a were generated by finding *Ta for Tb* = 2, using tracepq() for 1/7 and then letting *Tb* = *Ta *and using tracepqs() for *Tb* ≠ 2 to get a new parabolic *Ta* for 1/5. The generalized equation is no longer a polynomial, so the root solving will give only one root, but you can enter a suggestion as the last parameter to prompt it to search near the original value of *Ta*. As above, you can see hidden aspects of the periodicities in the limit set on the left in the five and seven crossing circles, which was the actual solution generated, and to a lesser extent in the congugated version right.

Chapter 9 of Indra's Pearls shows that limit sets with parabolic traces can be shown to be conjugate to a sister limit set with a parabolic *Ta* and *Tb* = 2, by a transformation of generators to a new pair (*u*, *v*) expressed as words in (*a*, *b*) so that *u*, *v* are both parabolic and (*a*, *b*) can also be expressed as words in (*u*,*v*), causing the new group to be a conjugate of the original. For example if *u* = *w*_{1/2} = *aaB* and *v* = *w*_{2/3} = *aaBaB* then
we have:

*aaB*(*aaBaB*)^{-1}*aaB* = *aaB*(*bAbAA*)*aaB* = *aaBbAbAAaaB* = *a*

(*aaB*)^{-1}*aaBaB*(*aaB*)^{-1}(*aaB*)^{-1}*aaBaB* = (*bAA*)(*aaBaB*)(*bAA*)(*bAA*)(*aaBaB*) = *B*.

If *a* and *B* are parabolic then so are *u* and *v* and finding their traces, which we can calculate by multiplying out the products of *a* and *B* and entering the traces into grandma's algorithm we can get matrices *u'* and *v'* with the same traces as the words* u* and *v*, with the commutator parabolic. Thus given any double cusp limit set defined by *Ta* and *Tb* = 2, we can find a double cusp limit set with *u'* and *v'* words in *a* and *B.* We can explore this numerically with the above example. If we apply tracepq(8,13), we get *Ta* = 1.62217905-1.28170820*i. *If we then put this *Ta *into grandma() with *Tb* = 2, we get two matrices *u' *and *v'. *If we now calculate the following traces for *w*_{1/2} and *w*_{2/3}, which are close to the middle of the Farey journey from 1/0 to 8/13, directly with Matlab, we find *TaaB = Tu'u'V'* = 1.55210536 – 0.91396228*i* and *TaaBaB *= *Tu'u'V'**u'V'* = 1.55210536 + 0.91396228*i,* conjugates to 8 decimal places! We can in turn plug this into either depthfirst algorithm and we get the results shown below in fig 13b.

We can see this process in action in the context of searching for doubly degenerate limit sets in the conjugacy between fig 15 (left) and fig 16 (centre), where the construction is key to undestanding how we can in fact get double degeneracy. That said, demonstrating the images in fig 13 are conjugate to ones with *Tb* = 2 is another matter. All we have proved here is a special case. In fact, if we took the alternative grandma's root to the original we would get a similar limit set to the one centre.

**Rational Periodicity becomes Irrational Degeneracy**

In fig 15 is an illustration of the singly-degenerate limit set defined by Troels' point, named after Troels Jorgensen –*μ*(*γ*)=*μ*((–1 + 5^{1/2})/2), of the golden mean limit of the Fibonnacci fractions representing the lowest point on the Maskit slice. The only noticeable difference between the right hand image in fig 14 and the limit set of Troels' point in fig 15 is right on the boundary of the horizontal circles, where *μ*(89/144) still has infinitessimal hammer heads defining the horizontal circles.

The Troels' point example above is deemed singly degenerate because half the region (in this case inside the bounding circle minus its circular replicates, has been engulfed by the limit set in a space-filling fractal of fractal dimension 2. Such points exist because the double cusp values on the boundary of the Maskit slice are countable and the (continuous) boundary has been proven to be uncountable so irrational vaues clearly exist. In fact, almost all values are irrational so degeneracy is the rule! Notice our manner of approach to the limit alternately moving left and right down the Farey tree 2/3, 3/5, 5/8, 8/13 etc. alternating on either side of the golden mean and tending to it in limit. Looking at the sequence in fig 14, one can see why this limit set is degenerate and englufs the free interior. Each successive Fibonacci fraction runs to a smaller pair of horizontal circles, which in limit engulf each of the horizontal circles and their fractal replicates. The fractally self-similar nature of the limit sets of Mobius maps of fractional linear transformations, means that the scale approaches a geometric space-filling limit.

From the fundamental transformations, one can that under L or R moves, the traces transform as follows:

(*Ta, TaB,* *Ta.TaB* – *TB*) <*– *L* –* (*Ta*, *TB*, *TaB*)* – *R*–> * (*TaB*,* TB,* *TBTaB* – *Ta*)

and these compose, so that for RL we get (*Ta*, *TB*, *TaB*)* – *RL* –> * (*Ta.TaB* – *TB**, TaB,** T ^{ 2}aB*.

The examples we have used of fractional traces *μ*(*p*/*q*) have been based on finding the roots the *q*-th degree polynomial using the polynomial solver vpasolve() in Matlab. This provides a complete list of all possible roots without prompting, but tends to break down as the degrees pass about 400, but a carefully crafted Newton iteration as used in Indra's Pearls, such as the one included as **tracepqr()**, which is provided in both Matlab and a native C version, so it can be used without Matlab on a Mac. When suitably tweaked, and given a suitable priming value, this can capture much higher fractions using an initial trace value such as a Farey-related lower fractional trace. One can input successive Fibonacci fractions as inputs to the next and easily reach *μ*(987/1597), *Ta* = 1.616888739503 -1.294321525461*i*, very close to the limiting value above, which was also generated this way in Indra's Pearls. The spirally degenerate examples in fig 17 have had their traces discovered by the Newton method.

**The Total Eclipse of Double Degeneracy**

One can take the degeneracy issue up a notch to doubly (or totally) degenerate if we can find a pair of traces engulfing both the complements e.g the outside and the inside sets. The trace parameters for a close approximation can be been derived using generator transformations, as described above.

We start with a high order Fibonacci fraction close to the Troels' point golden mean, generated from 8 left-right passes of Farey mediant addition from 0/1 and 1/0 down the Fibonacci fractions, for example using the function fareyadd() below, use trace recursion, to get the limit set of –*iμ*(987/1597) and apply a generator transformation to (*w*_{21/34}, *w*_{13/21}) the middle pair of terms in the left-right sweep, using grandma's algorithm as noted above, we end up with the two nearly conjugate traces: T (*w*_{21/34}) = 1.501347474086 – 0.865385203320*i* , Tr (*w*_{13/21}) = 1.501347474169 + 0.865385203218*i*. However this limit set (fig16 centre) is still singly degenerate, and has to be, since it is conjugate to the Troels' point limit set (fig 15). In actual fact the set looks closer to being conjugate to one of the right hand sets in fig 15, (which it is, since it is conjugate to 987/1597).

function [p0,q0,p1,q1]=fareyadd(s)

% e.g.[p0,q0,p1,q1]=fareyadd([0 1 0 1])

p0=0; q0=1; p1=1; q1=0;

for i=1:length(s)

if s(i)==1

p0=p0+p1;

q0=q0+q1;

else

p1=p1+p0;

q1=q1+q0;

end

end

To actually get a doubly degenerate group we aim to find a seqence of generator fractions *p _{n }*/

This limit set has intriguing properties linked to hyperbolic spaces and knot theory. In addition to the symmetries of *a*, *b* , there is an additional induced parabolic symmetry *c* , which arises because the LR move which transforms (*a*, *B*) into (*a*^{2}*B, aB*) must itself arise from a Mobius map which conjugates the limit set to itself, has the same fixed point as the parabolic commutator, and commutes with it, although performing a distinct 'translation' flipping the sets bounding the jagged line connecting the two principal centres (IPs p 339, 388). Robert Riley has shown that this system of transformations results in a gluing of 3D hyperbolic space via the discrete tiling to become the topological 3-sphere with the figure 8 knot (right) removed (doi:10.1016/j.exmath.2013.01.003). You can discover more of this knotty area by searching for "Not Knot" on YouTube.

Just as we did alternating left and right passes of successive fractions to get the Troels' point trace value, we can likewise look for sequences of moves which induce a spiral degeneracy. Starting at the pair of fractions 0/1, 1/0, we take ten left steps by means of Farey addition to arrive at 0/1, 1/10. Then a right step leads us to the pair 1/11, 1/10. Applying the same 10+1 steps again leads to 12/131, 11/120 and then to 143/1561, 131/1430. A Matlab and C script version of **fareyadd** are each provided to turn LR sequences into such pairs of associated fractions. Continuing this way several more times, we end up at –*iμ* = 1.9264340533 – 0.0273817919*i*. You can locate the sequence of fractional cusps on the boundary of the Maskit slice as at left and verify that they ascend near the highest cusp and then successive cusps rotate in exploded blowups of the boundary verifying their spiralling limit sequence. We can see why the engulfing takes place in terms of the sequence of Farey fractions, each stage of which leads to a smaller circle, which is broken by one further stage along the sequence seeking a smaller circle by breaking symmetry.This symmetry-breaking results directly from each higher fraction representing a deviation from the previous one so that the previous level horns do not unite but form a new engufling set of spirals.Two stages in construction of the spirally degenerate limit set are illustrated centre in fig 17. At each stage the limit set engulfs a step furtherinto the two large horizontal circles, with a spiral split again just beginning in the right hand image, confirming the process would continue in the same manner as for the golden mean Troels' point. Just as the LR pattern leading to the golden ratio is a limit of continued fractions approaching a quadratic irrational , so the value of the L^{10} R spiral degeneracy approaches a quadratic irrational giving *μ*(0.0916079) as the limit value of *Ta*:

**The Holy Grail at the Edge of Chaos**

We can also use the above techniques to define doubly spirally-degenerate limit sets based on the same development of double degeneracy applied to the spirally-degenerate example above. To see how this works, let's consider instead of alternating LRLR..., as in the doublydegenerate case above, the sequence LLRLLR..., which is a step to the L^{10} R in the above spiral case. The sequence LLR takes us as follows (*a,B*) > L > (*a*,*aB*) > L > (*a*,*aaB*) > R > (*aaaB*,*aaB*), so we need *TaaaB* = *Ta.TaaB* −* TaB* = *Ta, **TaaB* = *Ta.TaB* −* TB* = *TB *and *T*^{ 2}*a* + *T ^{ 2}B* +

Going another step further with L^{3}R... we have *Ta* = *TaTB* – *TaB*, *TB* = *TaTa *^{2}*B* – *TaB*, *Ta*^{ 2}*B* = *TaTaB* – *TB* plus the Markov identity, leading to *TB* = (*T *^{3}*a* – *Ta*)/(*T*^{ 3}*a* – 2*Ta* – 1), and *T*^{ 8}*a* - 4*T*^{ 6}*a* - 3*T*^{ 5}*a* + 4*T*^{ 4}*a* + 7*T*^{ 3}*a* + 3*T*^{ 2}*a* = 0, giving *Ta *= 1.727154868638 – 0.228442123457*i*, *TB * = 1.572495224561 + 1.575566700023*i*. The limit set fig 7a right shows the emergence of spiral degeneracy with both an external and an internal set of spiral fronds.

For L^{4}R ^{∞}, we have *TB* = (*T *^{4}*a* – 2*T ^{2}a*) / (

One can see from the repeating grandfather identity how this process extends readily to the doubly spirally-degenerate solution of L^{10}R^{ ∞}. The corresponding system of equations iteratively solves *Ta ^{n}B* back to

We can solve out the higher terms unsing the grandfather identity in terms of *Ta*, *TB* and* TaB* as follows:

syms a B aB

a2B=a*aB-B;

a3B=a*a2B-aB;

…

a10B=a*a9B-a8B;

a11B=a*a10B-a9B;

a11B=expand(symfun(a11B,[a B aB]));

a10B=expand(symfun(a10B,[a B aB]));

a11B(a, B, aB) =

a10B(a, B, aB) =

We can then assert *Ta* = *Ta*^{11}*B* and *TB* = *Ta*^{10}*B* , by setting the left-hand sides to *a* and *B* respectively. Since these equations are linear in *B *and *aB*, we can make *aB *the subject of both and cross multiply the two fractions to give the form below in which *TB *can again be made the subject of a fractional equation and we can then assert *Ta* = *Ta*.*TB* – *Tab, *or *TaB *= *Ta*(*TB –* 1) to get TaB as well.

We can again use Matlab to advantage to expand the cross product of the fraction above, to get* TB* as a fraction of powers of *Ta*:

Hence:

We can then substitute these back into the Markov identity and solve. In fact Matlab can solve the entire system outright for any L^{n }R, as shown right in the solvable Maskit 'slices' for
*Ta*, *TB*,
*n *≤ 53, coloured from red through to blue, with again the roots running along the top curve of *Ta *being the actual trace solutions. Notice how differently they behave!

syms a B aB

s=[0 0 0 0 0 0 0 0 0 0 1];

t0=a;

t1=aB;

t2=B;

for i=1:length(s)

if s(i)==0

h=expand(t1*t0-t2);

t2=t1;

t1=h;

else

h=expand(t1*t2-t0);

t0=t1;

t1=h;

end

end

aB1=solve(t0-a,aB);

aB2=solve(t2-B,aB);

B=solve(aB1-aB2,B);

aB=a*(B-1);

f=a^2+B^2+aB^2-a*B*aB;

g=expand(f);

h=simplify(g);

e=eval(vpasolve(g,a));

Setting this to zero and evaluating gives the somewhat surprising solution with large real values: *Ta* = 1.936872603726 + 0.021074052017*i*, Tb = 1.878060670271 – 1.955723310188*i*, but this is clearly the correct set of trace values.

We note however that the Matlab solution breaks down for 'mixed' sequences like LRL^{2} R because we end up with mixed powers of *a*, *B* and *aB*.

**The Non-Free Case and Automatic Group Tables**

The final examples in this section are non-free groups, in which the generators *a*, *b* are chosen so that *C*^{2 }= (*aba*^{-1}*b*^{-1})^{2}= I and hence these relations between the generators become trivial identities, causing the search tree to have new dead ends, which have to be eliminated by generating what is called an automatic group. The generators also no longer have parabolic commutators so we have to find the appropriate value of *Tab* from *Ta*, *Tb* by using the four alarm version of "grandma's algorithm" by solving for the trace *TC *= *T*(*aba*^{-1}*b*^{-1}) = 0 arising from the above relation, for the given values of *Ta*, *Tb*.* C*^{2} = I implies *TC* = 0 since the off-diagonal elements of *C*^{2} i.e. (*p*+*s*)*q* = (*p*+*s*)*r* = 0, the off-diagonal elements of I, so *TC* = *p*+*s* = 0 unless *r* =* q *= 0. If you examine line of the code with the comment %Grandma's three parameter double alarm
*TC* = 0 gives a quadratic in *Tab* and you can try out each solution. There is a short function below to do this. You can likewise find the p/q cusps for the non-free case. The commutator trace has moved from -2 to 0 so the only change needed to tracepq() to calculate these is to the initial solution *TaB* = *Ta* + sqrt(2)*i, *giving the centre figure above for *μ*(1/10)*.* The central circles now no longer touch the outer ones, emphasized by their colouring in yellow, or grey in the coloured images.

We then have to set up an automaton(right) in the form of a table generating the dead ends in every part of the maze by recursively looking up successive words and terminating at neutralizing words such as *b*^{-1}*a*^{-1}*ba* and *aba*^{-1}*b*^{-1}*a*. You can see this right and encoded as FSA , where the 19 rows are given by I, *a*, *b*, *A=a*^{-1}, *B=b*^{-1}, *ab*, *AB, ba, bA, Ba, BA, abA, ABa, baB, bAB, Bab, BAb, abAB*, and *ABab*, with x counting as a death state. The table is used sequentially to take us from any state on the left by right multiplication by one of the four generators. The states are stored sequentially as we move down a level of the maze so they can be looked up and the next recalculated. The end result, centre, has two disconnected sets of circles, coloured white and yellow in the figure. It is also a double cusp in which *b* and *a*^{10}*b* are parabolic. Two further examples from Indra's pearls are shown left and right, each coloured by the first iterate, and three more below to illustrate the diversity, with the disconnected central region in yellow.

Here is a short version of the function nonfree() to extract an appropriate *Tab* for the non-free case. Option selects which root (usually 0):

function tab=nonfree(A,B,opt)

%ta=1.924781-0.047529i;

%tb=2;

ta=A;

tb=B;

p=-ta*tb;

q=ta^2+tb^2-2;

if opt

tab=(-p+sqrt(p^2-4*q))/2;

else

tab=(-p-sqrt(p^2-4*q))/2;

end

You can use tracepq() in the non-free case to generate any fractional twist with *Tb* = 2, followed by non-free on the result to get the required *Tab *and then depthfirst to get the limit set. Below is a session for 3/7. The option is now 2 for non-free() and prt is 1 to print the identities. Non-free has option zero to pick the first root. The parameters in depthfirst() choose s(2) and now have option 2 for non-free and choose t for *Tab*.

>> s=tracepq(3,7,2,1);

3 -1.65805636-0.71965515i

1 -1.61195573-0.46319322i

4 -0.42695172-0.16466462i

5 -0.25623249-1.15299930i

6 1.03890744-0.78635572i

7 1.08561081-0.34873804i

2 1.82867804-0.60703463i

>> t=nonfree(s(2),2,0);

>> depthfirst(s(2),2,0.001,1000,1/3,3,2,1100,t);

The limit sets show the differences in the non-free case where the commuator squared is the identity. The image on the right is the analogue of the right hand image of fig 4. The centre image displays periodic circles in the central region while the exterior is fronded, as in the left image of fig 19. You can also use tracepq() and then tacepqs() with option 2, 3 and nonfree() to produce limit sets in the same manner as fig 18 right.

These figures show the same set of Fibonacci fractions as in fig 15 trending from periodicity to non-free degeneracy. The circles are still evident in the centre image at the crossings.

To generate the non-free cases using the DFS algorithm requires developing a tailored automaton for each case. Fig 26b thus shows some of the usual examples in this case, portrayed using the DFS algorithm. In this case I realized that T*abAB* = 1 can have the solution (*abAB*)^{6} = I, just as *Ta* = 1 did previously. This presents another type of generator solution than T*abAB* = 0, because, here the commutator has *r* = *q* = 0, rather than *p* + *s* = 0. The traces for increasing powers of *abAB* thus follow the sequence T*abAB* = 1 -> -1 -> -2 -> -1 -> 1 -> 2 = TI, without arriving at a trace of 0. However, because (*abAB**)*^{2} = I also satisfies (*abAB*)^{6} = I, the T*abAB* = 0 automaton works well enough in practice for this case to depict the limit sets reasonably well, although convergence is limited at high resolution.

**Generating Finite State Automata** However for other values, we need to generate custom FSAs. The packages KBMAG or MAF can be used generate exact automata for a much wider variety of limit sets. All Kleinian groups have been proven to be automatic, although this isn't the case for groups in general. MAF is available for Windows and as a Mac binary which runs flawlessly from Mac Terminal, or from any unix command line. A variety of these are included in the Matlab FSA folder.

The lower series in fig 26b have had finite state automata (FSAs) generated using MAF where either the generators or commutator or both are elliptic: (1) for *a*^{3} = I and (*abAB*)^{3} = I , (2,3) for *b*^{3} = I and (*abAB*)^{3 }= I and (4-6) for (*abAB*)^{3 }= I. A version of the DFS algorithm depthfirstcom is included, which allows such FSAs to be input into the function, as well as the commutator trace. Sample files with several of these FSAs are included with this package.

Both the matlab depthfirst function depthfirst.m and the C-script have been extended to accomodate input including an external FSA file. The FSA folder has a series of MAF generated files.

The C script depthfirstFSA.c contains the FSAs required to generate all the examples above.

For example for the first non-free case above with T(abAB)=0 i.e. (abAB)^2=I, one can run the following script called kleinian in the bin folder of MAF using Mac Terminal:

*_RWS := rec
(
isRWS := true,
ordering := "shortlex",
generatorOrder := [a,b,A,B],
inverses := [A,B,a,b],
equations :=
[
[(a*b*A*B)^2,IdWord]
]
);*

by using the command **./automata kleinian**. This will produce the following FSA in kleinian.wa suitable for using in the C-script, once 1 has been subtracted from each entry and the entries are edited in TextEdit to correct the layout. The C FSAs are created by subtracting 1 from each entry, in the MAF (and Matlab) versions, as C arrays go from 0 to *n*-1. The order above is [a,b,A,B] not Alun's preferred order [a,A,b,B], so as to be compatible with our depthfirst algorithm.

*int FSA[25][4]={{1,2,3,4},{1,5,-1,4},{6,2,7,-1},{-1,8,3,4},{9,-1,10,4},{6,2,11,-1},{1,5,-1,12},{-1,8,3,13},{14,2,7,-1},{1,15,-1,4},{-1,16,3,4},{-1,8,3,17},{9,-1,-1,4},{18,-1,10,4},{19,5,-1,-1},{6,5,-1,-1},{-1,2,7,-1},{-1,-1,10,4},{20,-1,-1,4},{1,5,-1,21},{1,22,-1,4},{9,-1,23,4},{6,2,24,-1},{-1,-1,3,4},{-1,8,7,-1}};*

This is slightly longer than the table above right for this example, which lists as:

int FSA[19][4]={{1,2,3,4},{1,5,-1,4},{7,2,8,-1},{-1,2,3,6},{9,-1,10,4},{7,2,11,-1},{12,-1,10,4},{1,5,-1,13},{-1,2,3,14},{1,15,-1,4},{-1,16,3,4},{-1,2,3,17},{1,18,-1,4},{9,-1,-1,4},{-1,-1,10,4},{7,2,-1,-1},{-1,2,8,-1},{-1,-1,10,4},{7,2,-1,-1}};

but it processes in depthfirstFSA.c to give essentially the same limit set.

**Symmetric Normalisations and Hyperbolic Planes**

Several of the limit sets in fig 26c show both close relationship with the hyperbolic plane and inversion of the exterior space of the limit set in the centre of every whorl. The system fig 26c(second from right) with *b ^{3}* = I and

The complementary question to the use of FSA automata is how to transform the images from those using grandma's normalization into ones which reflect the rotational symmetry of the relations when the commutator and/or the generator(s) are elliptic and have powers equalling the identity. In the symmetric inverse examples of fig 26c , grandma's recipe is transformed. In the first step one can sometimes use a change of generators from *a* and *b* to* e* = *a*^{-1}*b*^{-1 }= *AB *and* a ^{-1}*. This is invertible since from

As can be seen from the images in figs 26c-g, this displays the limit sets in the *n*-fold group symmetry of their commutator, resulting in 3, 4, 6 and 7-fold symmetries, further ornamented by the symmetries of the generators or generator words, also reflecting the multiple group symmetries of von Dyck triangle groups.

To understand the symmetries produced, we can start with the parabolic case where T*abAB* = -2. We can use the Grandfather identity to show (a) T(*c*^{2}) + T*c *^{-1}*c* = T*c*T*c *^{-1}= (T*c*)^{2} so T(*c ^{2}*) = (T

There is a counterintuitive relationship between the commutator trace value and odd and even orders of an elliptic matrix. While the only non-trivial solution of *a*^{3} = I has trace -1, the use of -1 for the commutator results in 6-fold symmetry in the symmetric inverse normalization (fig 26h), while the value 1, which should be of order 6, has 3-fold symmetry (fig 26c2).

Grandma's normalization is a way of transforming limit sets, especially those with parabolic commutators, which are the analogue of translations, from the linear friezes of Jorgensen and Maskit normalizations to a centered symmetrical image. By contrast, symmetric inverse normalizations seek to picture limit sets of groups with elliptic commutators in terms of the fractional rotations their elliptic commutators invoke.

Alun Williams, who developed MAF, notes that the usual generators of the group don't have a nice automatic structure. Changing to using three elliptic generators [* p*,* q *and* r* ] gives a more elegant structure. The generators can in turn be translated back to our original *a*, *b*, *A*,*B* using the gpxlatwa utility in MAF. The FSA array can be easily generated by executing the highlighted commands in Example 4 in order with Terminal, generating a .xwa file containing a 623 x 4 entry FSA array that also plots correctly.

It is possible, given *a* and *b*, to find a suitable *x* as follows. The Lie bracket *L* = (*ab*-*ba*) turns out to have the property that *L*^{2} =* k* I. One can thus find a suitable *x* by normalising *L*^{2} so that it has determinant 1. Although *x* has order 4 as a matrix, as a Moebius transformation it is an involution (self-inverse). This means that *x* has trace 0 and determinant 1 and furthermore *xx* = -I and *xax* = -*a*^{-1} and *xbx* = -*b*^{-1}. This means that the matrix *aba*^{-1}*b*^{-1} = *abxaxxbx* = -(*abx*)^{2}. Hence *p*=*x*, *q=bx*, and *r=ax* are each involutions as Mobius maps and *abx* is the effective square root of the commutator.

In the examples in fig 26c, we have taken the elementary solutions of (*abAB*)^{7} = ± I using MAF to generate a series of FSAs and applied them to the two trace solutions *TabAB *= -1.8019 = -2cos(π/7), and *TabAB *= -1.2470 = -2cos(2π/7). I have also replicated the *pqr* system successfully in the successive images.

I have found the FSA generating process to be completely idiosyncratic. After a two year pause where all my FSAs produced only approximate images at low depth search values, I accidentally hit on an illogical combination of abAB and and p,q,r generator order and modified symmetric inverse normalizations, which suddenly produced the fully accurate limit sets below, although some such as fig 26d still obstinately refused to do so completely.

I have uploaded a fast C script depthfirstFSA.c and a Matlab function depthfirstcom.m, both of which use the commutator value, symmetric inverse normalizations and FSA input to depict key elliptic and non-free commutator examples. I have used the c-script throughout for depiction due to its speed. This plots using the variable *ep* (epsilon) to draw a line between points, but *levmax* (depth level maximum) does not, contrary to the case in the original script depthfirst.c. Thus by setting *ep* to the minimum pixel width for the image size, elliptic limit sets with disconnected fractal "islands", which are thus not quasi circlular fractals, can still be accurately drawn. If the generated FSA is not perfect, these may have to have their levmax set to 32 or less to compute, but the illustrated examples below, all of which except for 26d, have excellent commutator square root FSAs, level maxima as deep as 4000 will depict rapidly in real time.

The generators and the commutator in these examples need to be carefully chosen to reflect compatible rotation and triangle group compatibilities.

We now examine a system where both *ab* and *abAB* are elliptic of order 7, setting up a dual order 7 symmtery in the symmetric inverse normalization. *Tb* can here be easily calculated for a given *Ta *using the extended Markov identity:

To get a sequential development of the fine details of using MAF with the square root of the modulator, I next took the above system and producred first a simple FSA to deal with (abAB)^7=I as utilized above using the command * ../bin/automata ../21a/kleinian* on the file

*_RWS := rec
(
isRWS := true,
ordering := "shortlex",
generatorOrder := [a,b,A,B],
inverses := [A,B,a,b],
equations :=
[
[(a*b*A*B)^7,IdWord]
]
);*

In each case this command produces a *kleinian.wa* file containing the FSA, which can be viewed with Text Edit and edited for the C-script by reducing all the values by 1.

I then expanded it to include the fact that (ab)^7=I using the same command on *kleinian2* producing a 127 step FSA:

*_RWS := rec
(
isRWS := true,
ordering := "shortlex",
generatorOrder := [a,b,A,B],
inverses := [A,B,a,b],
equations :=
[
[(a*b)^7,IdWord],
[(a*b*A*B)^7,IdWord]
]
);*

However neither of these produced the required accuracy.

I then applied the x, a, b, c version on the following file *kleinian3*:

*_RWS := rec
(
isRWS := true,
ordering := "shortlex",
generatorOrder := [a,b,A,B],
inverses := [A,B,a,b],
equations :=
[
[(a*b)^7,IdWord],
[x*a*x,A],
[x*b*x,B],
[(a*b*x)^7,IdWord]
]
);*

With the added file *kleinian3.sub*

*_RWS_Sub := rec
(
subGenerators := [a,b,A,B],
subGeneratorNames := [a_,b_,A_,B_],
subGeneratorInverseNames := [A_,B_,a_,b_]
);*

Issuing the three commands:

**../bin/automata ../kleinian3 -cos -nokb
../bin/automata ../kleinian3 -nokb
../bin/gpxlatwa ../kleinian3**

then produces *kleinian.sub.xwa*, with the final version 389 step FSA.

Although this didn't produce the accelerated result for this system, the same process applied to the corresponding systems with (*ab*)^3=I and (*abAB*)^7=I in fig 26d2, the tricorn gasket above fig 26d2b and the system in fig 26c, where *b*^3=I and (*abAB*)^7=I, immediately gave FSAs providing rapid accurate depth first search at essentially unimited depth. The addition of the *x*-involution thus gives an extremely rapid and accurate limit set for these systems, which now run at potentially unlimited depths limit of 4000+, rather than the previous 20-32 maximum, confirming the acuity of the method.

Since the equations for determining cusp values are more difficult when the reference value *Tb *= 2 is changed to an elliptic Møbius map, so that the solutions to the trace equations are no longer polynomial, if one has a fast accurate FSA, one can pick a loxodromic* Ta*, modify its real part towards degeneracy and then modify both the real and imaginary parts to navigate the associated Maskit slice towards a local cusp point. Fig 26d2c shows this process in operation.

An extension of the stochastic algorithm, qFstochasticom, is also provided to freely explore arbitrary limit sets without needing to determine a viable FSA.

**All Three Traces Elliptic**

Fig 26e shows a series of examples where we are stipulating three elliptic traces for *ab*, aB and *abAB*, resulting in just a single possible case for each since all the traces are now determined by the periodicities. Alun Williams' example, centre-left displays a 245 group 5-fold symmetry. The centre right limit set is an equivalent portrait with the same trace values in the symmetric inverse normalization reflecting the commutator square root periodicity of four. It has a whorl-centered normalization, while the left hand limit set is centered on an image complementing the circular-boundary. Choosing the alternate root of *Tab* in the stochastic portrait far right stochastic sketch also displays a planar Euclidean 244 symmetry, which appears identical to the one far-left from Alun. This implies the differences between the two centre images are purely a result of a difference in the normalisations. In the centre left figure, the symmetric inverse is order 5, which is the order of *ab* not the commutator, so it has to have been designed for this example, by mapping the limit set around one of the four centres of valency 5. The equations to derive the values of T*a* and T*b* are shown in fig 26f.

One can also derive the limit sets for cases where each of T*a*, T*b* and T*c* are directly elliptic, leading to a further series, which extend beyond the constraints of the triangle groups to number combinations such as 345, with each of the three underlying 3, 4 and 5 symmetries highlighted in yellow, green and orange in fig 26f.

Examples of Alun Williams' elliptic limit sets can be viewed here, here and here.

**Discreteness Enters the Chaotic Milieu**

When the DFS algorithm faces a limit set which is non-discrete and fails to tile the plane in an ordered fractal the search tree can fail to find closely connected points and so will appear to hang, runnning to the end of levelmax each time before going onto the next branch. Many of these limit sets also have interesting properties, and can be explored using two different approaches. The DFS Viewer app has a timer function and an escape routine both of which can interrupt the process if too much time has elapsed or the iteration has become unbounded. This can throw up a surprising diversity of images. However, because the strategy of the DFS algorithm is to draw lines beween closely identifies points in the limit set, its depictions give unrealistic portraits of these examples because a large number of lines are drawn higgledy-piggledy, only the end points of which are potential members of the limit set. As shown in fig 17d the FSA version of the DFS search algorithm does not draw line when the distance measure epsilon is set to the pixel size so fractal pontillistic chaos ensues.

An alternative method to portray these processes is the stochastic algorithm below, which is simply remapping a set of limit points already found on the limit set. This has another drawback, that some regions of the limit set are visited exponentially less often, but the examples in the enxt section show a variety of non-discrete limit sets with chaotic iterative behavior, where the orbits are intermingled.

The stochastic algorithm **qFstochastic()** is inferior at depicting the cases above where we know that the limit set is a fractal quasi-circle of a quasi-Fuchsian group or related examples which converge by the DFS algorithm, but it comes into its own when trying to explore the boundaries at the edge of degeneracy or dive into the chaotic regime where the depth first algorithm breaks down.

The function can depict limit sets in various parametrizations including the Jorgensen and Maskit versions.

The algorithm works by firstly finding attracting fixed points of the transformations, then choosing a partial orbit of the commutator and then repeatedly iterating these points in a random sequence of transfrmations *a, b, a*^{-1} and *b*^{-1}, collecting all these points and then plotting them in the same way as the above program. Somewhat paradoxically for those used to chaotic iterations of Julia sets, the random application of the four transforms preserves the limit set, but it fails to distribute evenly, visiting some points exponentially rarely. Nevertheless it does give an accurate depiction of a couldy version of the limit set and works irrespective of whether the example is quasi-Fuchsian, degenerate or chaotic, so one can explore a wider class of sets and use the algorithm to find new examples to apply to an image via the depth search algorithm.

The Jorgensen and Maskit parametrization provide different normalizations to grandma's recipe. In Jorgensen we define the matrices directly from the traces as follows: *a* = [*Ta*-*Tb*/*Tab* *Ta*/*Tab ^{2}*;

The utility of the stochastic algorithm is emphasized by its ability to depict atypical limit sets such as the non-discrete examples below in fig 25 and those arising from elliptic trace values *such as Ta *= 1 and *Ta* = 2^{1/2}, with *a *^{6} = I and *a* ^{8} = I respectively even under variations of the complementary trace, including doubly non-free examples where *(aba*^{-1}*b*^{-1})^{2}^{ }= I as well.

There is an interesting pattern to the non-free trace values giving *a ^{n}* = I. The solutions to this have a close relationship with those for

n = 2 factor(a{3}-2) = t - 2, t + 2

n = 3 factor(a{4}-2) = t - 2, **(t + 1)^2**

n = 4 factor(a{5}-2) = t - 2, t + 2,** t^2**

n = 5 factor(a{6}-2) = t - 2, **(t^2 + t - 1)^2**

n = 6 factor(a{7}-2) = t - 2, t + 2, **(t - 1)^2, (t + 1)^2**

n = 7 factor(a{8}-2) = t - 2, **(t^3 + t^2 - 2*t - 1)^2**

n = 8 factor(a{9}-2) = t - 2, t + 2, **t^2, (t^2 - 2)^2**

n = 9 factor(a{10}-2) = t - 2, **(t + 1)^2, (t^3 - 3*t + 1)^2**

All the *n* values include the trivial solution *t* = 2 for which *a*^{1}= I. The even *n* values have double solutions *t *= ± 2 and odd values just 2 corresponding to *a*^{1} = I. The remaining roots give pared solutions to *a ^{n}* = I and its factors of

Because of the close relationship with *μ*(1/*n*) recursion, the same set of solutions appears in the *μ*(1/*n*) polynomial factorizations, with one of the root pairs unchanged and the other modified by complex coefficients, derived from defining *Tab *in terms of *TabAB *via Markov:

n = 2 factor(a{3}-2) = t^2 + ki*t - 4

n = 3
factor(a{4}-2) = **t + 1**, t^2 - (1- ki)*t - (ki + 2)

n = 4 factor(a{5}-2) =** t**, t^3 + ki*t^2 - 4*t - 2*ki

n = 5 factor(a{6}-2) = **t^2 + t - 1**, t^3 - (1-ki)*t^2 +(3-ki)*t + (2-ki)

n = 6 factor(a{7}-2) =** t - 1, t + 1**, t^4 + ki*t^3- 5*t^2 - 3ki*t + 4

n = 7 factor(a{8}-2) = **t^3 + t^2 - 2*t - 1**, t^4 - (1-ki)*t^3 - (4+ki)*t^2 + (3-2ki)*t + (2+ki)

n = 8 factor(a{9}-2) =** t, t^2 - 2**, t^5 + ki*t^4 - 6*t^3 - 4ki*t^2 + 8*t + 2ik

n = 9 factor(a{10}-2) =** t + 1, t^3 - 3*t + 1**, t^5 - (1-ki)*t^4 - (5+ki)*t^3 + (4-3ki)*t^2 + (5+2ki)*t - (2-ki)

These polynomial solutions also manifest fractional rotations, so that the solution for *n* = 3 is -1 = -2cos(π/3), the 2 solutions for *n* = 5 are -1.618 = -2cos(π/5) and 0.618 = 2cos(2π/5). The 3 solutions for *n* = 7 are -1.8019 = -2cos(π/7), 1.2470 = 2cos(2π/7) and -0.4550 = -2cos(3π/7) and the 5 solutions for *n* are (-1)^{k}2cos(*k*π/*n*) *k* = 1 - *n-1*.

Consequently the most negative valid real roots form a sequence of traces -2cos(2*pi/*n*) tending to -2 = -2cos(0), with even numbers giving factored values: (*n, Tn*) = (3,-1), (8, -2^{1/2 }= -1.414), (5, -(1+5^{1/2})/2) = -1.618, (12, -1.732), (7, -1.802), (16, -1.847), (9, -1.879), (20, -1.902), (11,-1.919 ), (24, -1.932), (13, -1.942), ..., (61, -1.99734), with the number of central circles between the top and bottom increasing by 1 at each stage. However the values of *n *remain faithful to their non-free relation, so while *Ta* = *T*(1/16) lies between *Tb* = *T*(1/7) and *Tc* = *T*(1/9), and *b*^{7} = I and *c*^{9} = I, *a*^{8} = -I so we get *a*^{16} = I. Where the positive root (opt=4) is chosen for the negative trace values or we can alternatively use the ususal negative root (opt=0) with the minus of the real value. As noted in fig 24, these give elliptical variants of the limit sets, when a parabolic generator trace of 2 is replaced by one of these, and the same values apply to both the standard case where the commutator has trace -2 and non-free cases e.g. where it has trace 0 because the solutions above don't depend on the value of *k*.

Below are examples where the parameters have been chosen to take the set outside the regime of discrete groups, so that the limit set becomes unstable, chaotic, or fractal cantor dust, which can include arbitrary commutator traces, or non-free situations requiring automatic groups for the depth first search algorithm.

**Life at the Edge of Chaos:** The series of images below show how the limit sets vary with the trace of the commutator, when both *Ta* and *Tb *are* 2.* The stochastic algorithm can freely explore the commutator parameter space generated by Grandma's extended recipe as normalizer. As can be seen from fig 26, the edge of non-discrete chaos according to the stochastic algorithm appears to be at T*abAB* = 1.

Higher rotational periodicities become less and less stable and more prone to chaotic non-discrete limit sets as the fraction denominator increases. For example if we look at the fractional rotations for *n*/11 which are also the solutions of a fifth degree polynomial solving (*abAB*)^{11} = ± I, we find these values which should provide discrete limit sets in the same manner as fig 26c now have instabilities. Fig 26g shows the limit sets for these fractional rotations in grandma's normalization for *Ta *=* Tb *= 2. As can been seen the positive values are all non-discrete and most of the negative values also show instabilities that disrupt their intrinsic symmetries. On the far right, symmetric inverse normalizations of two of these values using similar parameters to fig 26c, show that mode locking with lower periods is disrupting the limit set causing it to both be non-discrete and display asymmetric versions of 5-fold and 7-fold rotational symmetry. In the above case the commutator trace of -1.3097 is close to that of the 7-fold system with trace -1.2470 . In the lower case the 11-fold symmetry has been coupled to 5-fold with a perturbation to roughly 5.5 = 11/2-fold chaotic rotation.

These examples display a deep parallel to the quasi-periodic and chaotic regimes of a double pendulum, which is a renowned instance of a conservative dynamical system involving chaotic regimes. In both cases (simple harmonic motion and a Mobius transformation) we have two completely ordered transformations which when coupled in pairs are capable of dynamic chaos due to mode-locking interactions.

This demonstrates the value of the stochastic algorithm to give at least a sketch of the actual dynamics in the non-discrete case. One also needs to note that both the DFS and particularly the IFS algorithm strive to eliminate the onset of chaos, in the case of the DFS by approximating by a piecewise linear graph and in the IFS case by making an initial search of the word tree and selecting only words which further the ordered process, eliminating those which would result in instability and thus strand the algorithm.

**Escape-time algorithm** **m-file**

The orbital algorithm from Jos Leys web site works only in Maskit parametrization and only with tractable examples with Tb=2. The method scans all points in the field and iterates. The point is iterated by moving it via translations by 2 in the parallellograms bounded by the sloping line defined by the real and imaginary parts of Ta, until it ends in the centre section, then iterating it by a if it lies above the elliptical curve chosen below to separate the regions, or a if it is below, again translating if necessary and iterating again. If the point goes out the top, we colour it brown, if below, we colour it green and of it ends in a 2-cycle we colour it blue, separating the regions as shown in the diagram. It thus doesn't draw the limit set but the points on either side, highlighting the limit set by default as in level set methods for Julia sets. This works only for a few easily separated examples and require intial curve fitting of the separating sigmoid curve's slope and exponent.

**quasi3D** is a very simple example of a 3D transformation forming a limit set taken from Matlab central as the begining seed of the stochastic algorithm.

To view a much more thrilling 3D example, see Jos Ley's fly through.