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.
Fig 0: Left: Click to see full resolution – Moving set of images of all the Kleinian limit sets of Ta = –iμ(n/61), Tb = 2 generated with the C-script depthfirst61.c This script is nearly 100 times faster than the extended feature Matlab version and will run on any platform. See explanation below. Right: The Mac DFS Viewer app in operation, displaying in colour the 26-th iterate of a quasi-Fuchsian limit set. The fractal fragmentation of the iterates has only reached part way out from the centre and the tips are still single coloured.
Contents
Overview of the Computational PackageOverview 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.
Fig 0b: Convergence towards space-filling degeneracy as the real part of loxodromic Ta = 0.95+0.05i is reduced towards the free discrete group at 1.886+0.05i , each coloured by the 30th iterate.
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:
Ta and Tb set the real and imaginary parts of the two matrix traces, levmax sets an upper limit in the depth search tree, epsilon gives the fineness of the limit set resolution relative to the overall scale used to expand or shrink the image. The timer interrupts the process after a fixed time in seconds the next time it hits levmax, to exit the process if the trace values are super-critical. If only a part of the limit set draws, you need to set the timer for a longer period. Set timer to 0 if you want a high resolution image you know has good traces. Non-free chooses the non-free case discussed below. Necks performs the additional midpoint checks to deal with large narrow-necked fronds as in fig 4, Escape turns on an escape routine if the process begins to go unstable or chaotic. If super-critical values are input the app will then attempt to detect runaway and exit displaying the chaotic iteration up to this point. but can also disrupt depiction of large limit sets. Iterate, if set to a non-zero value n, colours each pixel by the n-th iterate of the four generators – a, b, a-1, b-1. Set to 0 it gives a black and white image, and set to -1 it gives the final level iterate for each pixel. Alt chooses the alternative square root for Tab, which sometimes differs between the Matlab and C versions. If you find you are getting different effectsfor a given set of trace parameters, try choosing alt root, and/or reversing the signs of the real or imaginary parts of the traces.
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.
On a Mac, open depthfirst.c using Text Edit and edit it to change the parameters to suit your example. Then open Terminal and navigate to the folder containing depthfirst.c. You can do this by typing cd< folder name> until you get to the folder, or simply by using cd<space> and dragging the folder icon rfrom the Finder into the Terminal window. Then copy and paste the gcc command at the top of the C script into Terminal and the script will compile and run. After this make a copy of the .ppm file, call it "test.ppm" and copy and paste the sips command into Terminal and voila! – you have a full resolution highly compact colour image – test.png. Alternatively a Matlab user can execute replayC() on the original filename without .ppm. to get te same png image. If you don't have Apple's XCode installed you will need to download it, free from the Apple Store and then install the Terminal command line tools by executing xcode-select --install in Terminal to activate gcc. On older systems, a suitable version of XCode can be downloaded by signing on for a free Apple Developer account. The script can also be edited, compiled and run in XCode as a command line tool.
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:
Fig 1: The limit sets for real joint traces greater than 3 are a circle (Fuchsian). As the traces Ta, Tb reduce 3 through quasi-Fuchsian 2.2 (right) to 2 (left), the fractal circle becomes a Kleinian Apollonian Gasket. The subset of circles shaded light grey are mapped from the largest left circle by the transformations a and BAb and never mapped over the white circles, forming a modular subgroup of the gasket group. The limit set (centre ) with traces Ta = 1.92 + 0.03i, Tb = 2 is colour coded by the first iterate. The right-hand figure shows locations of some key elementary repeated elements.
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-1NC 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.
Fig 2: Diverse mathematics underlies the forms of the limit sets. (b) shows Schottky circles of non-kissing circles as in (2), defined by loxodromic maps on the Riemann sphere (a), taking the inside of one circle over the outside of the paired one, resulting in fractal dust Cantor sets as the limit sets. (c) shows how this becomes a circle when pairs of circes touch. (e) shows the Riemann sphere with the two circle interiors removed ready for gluing in ways which can generate new surfaces. (h) is the construction for the modular group represented by the Escher tiling consisting of the group of matrices with determinant 1. (f) shows the arrangement for a quasi-Fuchsian group, where the touching circles are further drawn together, so that the surface becomes two punctured tori joined at their cusps and the limit set becomes a fractal tpological circle, or Jordan curve. This is manifest algebraically in choosing mappings in which the trace of the commutator is parabolic, with TabAB = –2 causing the pinching. As we move the circles, more can kiss, so that when all do, as in the gasket (d), two more transformations have become 'pinched' so that both a and b are parabolic and we end up with (g). The four tangent circles generating the Kleinian group of the Apollonian gasket partially portrayed in (d) are highlighted in black. One of the four circles is the line y = 0, with centre the point at infinity. The gasket group is called doubly cusped because we have pinched two extra handles because, a and b are also parabolic. It is also sometimes called maximally cusped, because, after all this squeezing, there are no more curves left to pinch. In many of the cases below with Tb = 2, only one generator b is parabolic, although in double cusp limit sets both Tb and a generator word defining the fractional motion are, and in some such as fig 3, neither are, although each quasi-Fuchsian limit set is still a fractal topological circle arising from the parabolic commutator..
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.
Fig 3: Left and centre: Quasi-Fuchsian limit sets Ta = 1.91 ± 0.05i, Tb = 1.91 ± 0.05i If we independently exchange the signs of the real and imaginary parts to find complementary points on the Maskit slice, which is symmetric under changes of sign of both x and y, we will find four symmetry variants of each of these sets.
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-1b-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 bABa, 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.
Fig 4: Special words and the words defining a frond of a limit set.
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 a15B 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-1b-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.
Fig 5: Lacy quasi-Fuchsian limit sets where Tb = 3 Left: Ta = 1.91 + 0.05i , Tb = 3 IPs 8.9 238 Right: Ta = 1.9656+0.99536i, Tb = 3 crop of spirals IPs 8.23 264
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.
Fig 6: Left: Special word arrays for the parabolic element aaB. Right" Special words algorithm working on Ta = 2 + i, Tb = 2 8.16 247.
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 a15B is in green and circles 1 and 2 are thus mapped to themselves. B carries C1 to C16 and then a15 carries it back to C1. In a similar way, B carries C2 to C17 and then a15 carries it back to C2. This means that a15B maps both the disks C0 and C1 to themselves, the fixed points of a15B 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 a15B for 1/15 or a3Ba3Ba2Ba3Ba2B 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 ε.
Fig 7: Double cusp 1/15 Ta = 1.958591030 + 0.011278560i, Tb = 2 both b and a15B are parabolic IPs 9.1 269. In general solving for a given denominator q involves iteratively solving a q-th degree polynomial. Centre: 2/19 cusp Ta=1.9038 + 0.03958i, Tb=2 IPs 9.15 293 a10Ba9B and b are parabolic. Right: 7/43 cusp Ta=1.80785524 + 0.136998688i, Tb = 2 IPs 9.16 294 − 295 plotted to higher resolution using the associated C-program. You can see the denominators clearly in the periodic mappings of the circles, and the numeratos both in the stepping periods along the chain and in the chains between the horizontal circle pair and the vertical pair.
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-1n = TmTn and (2) the extended Markov identity TabAB = (Ta)2 + (Tb)2 + (Tab)2 – Ta.Tb.Tab – 2. These form the founding identities to enable the decoding of the traces of the examples. We can then find Ta by solving a q-th degree polynomial e.g. by Newton's method, after solving for TaB and iteratively asserting constraints such as TaaB = ± 2.
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 anB. 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 + 2i . 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 nth 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.09628i, 1/15 1.958591 – 0.01128i, 1/23 1.98179 – 0.00320i, 1/31 1.98987 – 0.00131i, and 1/60 1.99726820 – 0.00018236, tending towards Ta = 2, (or μ = 2i ) 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 a3Ba3Ba2Ba3Ba2B can also be generated recursively e.g. from a3B and a2B by similar relations. To be more precise, examination shows that the word w(p+r)/(q+s) = wp/r.wq/s, in the manner of mediants and the Farey tree, so we can use the grandfather identity as follows:
Tw(p+r)/(q+s) = Twp/r.Twq/s – T(wp/r-1.wq/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) = Tw8/21 = Tw3/8.Tw5/13 – T(w3/8-1.w5/13) = Tw3/8.Tw5/13 – Tw2/5, because w3/8-1.w5/13 = (a3Ba3Ba2B)-1.a3Ba3Ba2Ba3Ba2B = a3Ba2B = w2/5.=w(5–3)/(13–8).
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+2i ) 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
Fig 8: Four higher fractions with rotation rates close to 1/2, 1/3, 1/4 and 1/5 : 10/19 1.70408769 – 1.01161324i 13/40 1.68667404 – 0.57810830i 10/39 1.71314918 – 0.35283835
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.
Fig 9: In the centre is shown a section of the "Maskit slice" - the plane of μ(p/q) where the generating p/q word, such as a15B is parabolic from Indra's Pearls (the Matlab version in in fig 10 below). A Matlab generated version is shown in fig 10. In the Maskit parametrization we have one parameter μ defining a Mobius transformation which is paired with translation by 2. In the standard parametrization we get the same limit set from Ta = -i μ(p/q), Tb = 2. Left and right are shown the values of -i(μ(p/q), 2 for two known values 1/2 and 3/10 with values 1 + sqrt(3)i and (1 + sqrt(11)i)/2 coloured by the final iterate. Notice the overlapping fractal dark circles in the right hand image. For an explanation of μ in relation to the trace values, see fig 23.
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.
Fig 10: Left: Ta= –iμ(i /30) Tb = 2, i = 0-29 showing all the rotation states mod 30. This can be contrasted both with the version for prime 29 (800x800 moving gif for n/29) and fig 0 where prime 61 is shown. Center: "Maskit slice" as in fig 9 centre on a rectangle between -2 and 2 by generating the symmetrized plot of all roots up to q = 51 coloured from 1 red to 51 blue using tracepq() in the script slicec.m and plotting with pltmus(). Right: Conjugate roots for (i /29) and their superposition. The sequences for n/29 and n/61, where all the fractions are irreducible because they are prime can be contrasted with the 800x800 moving gif for n/30, where there are factors of 2, 3 and 5 in the iterating steps. Watch n/29 and n/30 move together.
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 a15B for μ(1/15), is parabolic. In the Maskit parametrization, whose limit set corresponding to Ta = 1.91 + 0.1i, 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 = 21/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
Fig 11: Limit sets of μ(1/7), μ(1/15), μ(1/31) and μ(1/60) show they are not tending to the limit of μ(1/n) = μ(0/1), the Apollonian gasket, as 1/n tends to 1/∞ = 0/1. The limit set of μ(1/31) is pictured in blue with the gasket superimposed in black to highlight the relationship. The Apollonian gasket has trace Ta = 2 which is the limit of the traces of these limit sets and its matrix a is also the limit of the μ(1/n) matrices. But the limit sets themselves tend in limit to another limit set altogether whose group contains an additional factor as well as that of the Apollonian gasket.You can see this with the gasket superimposed in black on the limit set of μ(1/31) in blue (click to enlarge).
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 (anB)μ(1/n) introduces in limit a new symmtery involving rotation by π / 2, not existant in the gasket group.
Fig 12: μ(2/(2n+1)) for n = 60, or μ(2/121), demonstrating a third limit of limit sets tending to μ(0/1) different from both the above sequence in fig 11 exemplified in limit by μ(1/60) and the Apollonian gasket μ(0/1). Note there are now pairs of circles separating the principal ones rather than one as in μ(1/60). This is the second in an infinite sequence of distinct limits of limits.
One can take the limit of 2/(2n+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/(3n+1) and 3/(3n+2) and so on so we have an infinite number of different limits of limit sets, all tending 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.
Fig 13a: Left: Two sets with Ta = 1.8385 – 0.0964i, Tb = 1.7667 – 0.2201i and Ta = 1.8385 – 0.0964i, Tb = 1.7667 + 0.2201i.
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 = w1/2 = aaB and v = w2/3 = aaBaB then we have:
aaB(aaBaB)-1aaB = aaB(bAbAA)aaB = aaBbAbAAaaB = a
(aaB)-1aaBaB(aaB)-1(aaB)-1aaBaB = (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.28170820i. 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 w1/2 and w2/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.91396228i and TaaBaB = Tu'u'V'u'V' = 1.55210536 + 0.91396228i, 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.
Fig 13b: Left: μ(8/13),2. Centre: (w1/2,w2/3), which is clearly conjugate to the left figure, as confirmed by the alternate root of the same traces illustrated right, giving an identical limit set to the left one except horizontally flipped, and rotated through 45 degrees.
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
Fig 14: Periodicity limiting in degeneracy: Fibonacci fractions μ(3/5), μ(8/13), μ(21/34), μ(55/89) and μ(144/233) tending to the Golden Mean γ = (–1+sqrt(5))/2. The limit sets tend in limit to the degenerate golden mean Troels' point illustrated in Fig 16.
In fig 15 is an illustration of the singly-degenerate limit set defined by Troels' point, named after Troels Jorgensen –μ(γ)=μ((–1 + 51/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.
Fig 15: Singly degenerate Troels' point Jorgensen's group Ta = 1.6168866453 − 1.29432650i, Tb = 2 coloured by first iterate to highlight the boundary lines between the complementary discs left and right being consumed by the degeneracy IPs 10.4 318-319. This value corresponds to the minimum imaginary value on the boundary of the Maskit slice fig 11 centre and corresponds closely to the μ value of the golden mean. It is the mirror image of Ta = 1.6168866453 − 0.7056734968i. See fig 14 for sequence convergence images to the golden mean. Right singly degenerate Ta = 1.7847+0.73852i, Tb=3 Ips 10.1 311
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 2aB.Ta – TaB.TB – Ta). Hence these relations can also be used to decode other systems like L10 R. In these progressions, including the one above, the geometric scalings become determined by the eigenvalues of the transformation matrix. This leads to the discovery that the boundary of the Maskit slice is asymptotically self similar. You can see the nature of this phenomenon in fig 17 lower left inset, where successive sections of the corresponding spiral degeneracy are fractally scaled and rotated versions of the original. In our case above, one of the eigenvalues turns out to be (5 + 211/2)/2 ~ 4.791287848, so the traces of each successive pair of Fibonacci fractions are about five times closer. For example
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.294321525461i, 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.
Fig 16: Left: A singly degenerate conjugate of the Troels' point map above, forming an approximation to the doubly degenerate set right. As in fig 13b, the alternate grandma's root gives a circular limit set, again conjugate to the original even though Ta and Tb are not exact conjugates, and despite the mid-resolution iterates of the DFS algorithm (inset) proceding on an entirely different basis. Right: Doubly degenerate limit set LR∞ coloured by the third iterate drawn on the same scale as left of the space-filling in the plane Ta,Tb = (3±31/2i)/2 = 1.5 ± 0.8660254038i IPs10.11 335. The fixed point o the commutator lies at the cusp between one of the principal centre pairs and the additional transformation, c, which is also parabolic, has the same fixed point and commutes with the commutator, although distinct from it.
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 (w21/34, w13/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 (w21/34) = 1.501347474086 – 0.865385203320i , Tr (w13/21) = 1.501347474169 + 0.865385203218i. 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 pn / qn and –qn / pn so that in the irrational limit r and –1/r , they are degenerate both on the inside and outside. For this relationship to happen, the words get extremely long and there is little variation near the center, which suggests the traces of the original (a, B) should be invariant in the sense that they are the same as those of the next in the list (aaB, aB), since the sequence LR takes us as follows (a,B) > L > (a,aB) > R > (aaB,aB). This system has three defining equations: TaaB = Ta.TaB − TB = Ta from the grandfather identity, TB = TaB, plus the Markov equation T 2a + T 2B + T 2aB = Ta.TB.TaB, which provides three equations to solve for Ta, TB and TaB. Hence T 2a + T 2B + T 2B = Ta.TB.TB, or T 2a = (Ta.– 2)T 2B and Ta.TaB − TB = Ta, or (Ta.– 1)TB = Ta. Solving for Ta gives T 2a – 3Ta + 3 = 0, or T
a=(3± 31/2i)/2, TB =Ta / (Ta-1) = (3-/+ 31/2i)/2, with conjugate trace values of 1.5 ± 0.866025i closely neighbouring the above values in the previous paragraph. We have now induced a doubly degenerate limit set LR∞engulfing the entire complex plane. The limit set displayed in fig 16 right clearly has convergence to degeneracy on both sides of the remaining white islands, confirming its double nature.
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 (a2B, 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.
Fig 17: Spirally degenerate Ta = 1.9264340533 – 0.0273817919i, Tb = 2, IPs 10.6 322-323.The left image shows partial engulfment and the right another stage further. It took 9.7 hours to run in the C script using a levmax of 10000 and an epsilon of 0.0000025 and would have taken 36.9 days in Matlab because of the very slow engulfing convergence for this limit set. One can see that this limit set is pervaded by two different types of double spirals, the main set consisting of the main central, four peripheral and so on as in the centre image and all other similar examples, such as figure 11, and a second set unique to this system, arising from the spiral degeneracy. Inset below left: Self-similarity of the Maskit slice on approaching the values leading to sprial degeneracy. In the asymptotic limit the scaling factor approaches 86.214 + 164.198i, accounting for both the scaling and rotation.
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.0273817919i. 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 L10 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 L10 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 2a + T 2B + T 2aB = Ta.TB.TaB. This system solves out via (1) Tab = Ta(TB – 1) and (2) TB =T 2a/(T 2a – 2), to give T 4a – 5T 2a + 8 = 0 or Ta = ±1.63224188 ± 0.40523272i. (2) then gives TB = ±1.5 -/+ 1.32287565Ii, giving non-conjugate traces, which grandma confirms satisfy TaaaB = Ta and TaaB = TB, producing the limit set, fig 17a left, which has similar characteristics of double degeneracy to the previous example in fig 16 which can be compared with the limit set of (LLR)5 = μ(153/418).
Fig 17a: Left μ(153/418),2 - LLR5 compared with the doubly degenrate LLR∞ Right: μ(91/345),2 - LLLR4 compared with the doubly degenrate LLLR∞. Notice in the expanded voews, both the inner and outer whorls of the doubly degenerate cases reflect the limiting behavior approaching the horizontal 'circles' in the singly degenerate examples.
Going another step further with L3R... we have Ta = TaTB – TaB, TB = TaTa 2B – TaB, Ta 2B = TaTaB – TB plus the Markov identity, leading to TB = (T 3a – Ta)/(T 3a – 2Ta – 1), and T 8a - 4T 6a - 3T 5a + 4T 4a + 7T 3a + 3T 2a = 0, giving Ta = 1.727154868638 – 0.228442123457i, TB = 1.572495224561 + 1.575566700023i. The limit set fig 7a right shows the emergence of spiral degeneracy with both an external and an internal set of spiral fronds.
Fig 17b: Left: The corresponding motifs in single (L4R)3 35/169 and double L4R∞ clearly show their relationship.
Right: Double degenerate L6R∞ and single (L6R)3 63/433 generated using tracepqr() from the trace of 8/55 as seed.
For L4R ∞, we have TB = (T 4a – 2T2a) / (T 4a – 3T2a) and TaB = Ta(TB-1) = T 3a / (T 4a – 3T2a), giving T 10a - 6T 8a + 8T6a + 4T 4a = 0 andTa = 1.79210587 - 0.141970345i, Tb = 1.64779887126 + 1.721433237i, which again compares with (L4R)3 35/169 Ta = 1.75300569-0.22184750i as shown above. This leads to a systematic way to resolve any Ln R ∞.
Fig 17c: The Holy Grail: Left: Doublespiral degeneracy – genuine L10R∞ – deduced as outlined below, compared centre with the previous image of spiral degeneracy.On finer resolutionwith a higher levmax, both the external and internal fronds to either side erupt further in space-filling spirals as shown in the moving gif, but the convergence is very slow. The last image took 28.9 hours running on C and would have taken 112 days on Matlab. Right: The limit set with Ta and TB swapped to show it has equivalent fractal structure, when the 'outside' is exchanged with the 'inside'.
One can see from the repeating grandfather identity how this process extends readily to the doubly spirally-degenerate solution of L10R ∞. The corresponding system of equations iteratively solves TanB back to TaB and TB, giving three equations in Ta, TB, and TaB. Because the first two equations are linear in TB and TaB, one can eliminate these, to give TB and TaB as above, to get a higher degree polynomial in Ta.
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 = Ta11B and TB = Ta10B , 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:
o = expand( (B*a^9 - 8*B*a^7 + 21*B*a^5 - 20*B*a^3 + 5*B*a + a)*(a^9 - 8*a^7 + 21*a^5 - 20*a^3 + 5*a )
- (a^10 - 9*a^8 + 28*a^6 - 35*a^4 + 15*a^2 - 1 )*B*( a^8 - 7*a^6 + 15*a^4 - 10*a^2 + 2))
o = 2*B - 15*B*a^2 + 35*B*a^4 - 28*B*a^6 + 9*B*a^8 - B*a^10 + 5*a^2 - 20*a^4 + 21*a^6 - 8*a^8 + a^10
Hence:
We can then substitute these back into the Markov identity and solve. In fact Matlab can solve the entire system outright for any Ln R, as shown right in the solvable Maskit 'slices' for
Ta, TB,
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));
h = (a^2*(a^12 - 12*a^10 + 54*a^8 - 112*a^6 + 105*a^4 - 37*a^2 + 8))/(a^6 - 6*a^4 + 9*a^2 - 2)^2
Setting this to zero and evaluating gives the somewhat surprising solution with large real values: Ta = 1.936872603726 + 0.021074052017i, Tb = 1.878060670271 – 1.955723310188i, but this is clearly the correct set of trace values.
We note however that the Matlab solution breaks down for 'mixed' sequences like LRL2 R because we end up with mixed powers of a, B and aB.
The Non-Free Case and Automatic Group Tables
Fig 18: Left: Ta=1.934141–0.03668126i Tb=2 coloured by the second iterate. Centre: Sequence of non-free group limit sets for fractions n/29. See both the n/29 rotations together. Right Non-free Maskit slice.
The final examples in this section are non-free groups, in which the generators a, b are chosen so that C2 = (aba-1b-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-1b-1) = 0 arising from the above relation, for the given values of Ta, Tb. C2 = I implies TC = 0 since the off-diagonal elements of C2 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
Fig 19: Left: Ta=1.9247306-0.0449408i, Tb=1.91 + 0.2i, Tab=1.85692801 – 1.24257380i 11.3 363 Centre: Ta = 1.92478160 – 0.04752913i, Tb = 2 and Tab = 1.92480000 – 1.46174256i 11.1 354 μ(1/10) Right: Ta =Tb = 2 Tab = 2.0000 – 1.41421356i 11.3 363
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-1a-1ba and aba-1b-1a. 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 a10b 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
Fig 20: Upper row: 1.906-0.06739i, 2 Near degenerate high period, 1.865-0.70567i, 2 Close to 1/3 of a revolution μ(1/60) Ta = 1.9973 – 0.00025756i, showing the same limiting behaviour as in the standard version does not converge to the non-free gasket above. Spirally degenerate L10R 143/1561 1.934141199416-0.036681260008i.Lower row: Convergence towards space-filling degeneracy as the real part of loxodromic Ta = 0.95+0.05i is reduced towards 1.919855, each coloured by the 30th iterate.
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);
Fig 21: Left: 1.92+0.05i, 1.89+0.12i, Centre Left: 1.94+0.05i, 1.94+0.05i with similar non-free coefficients to the left hand images of fig 4. Centre: A limit set generated from μ(3/5),2 with tracepqs() and then coupling this as Tb with μ(1/7) using tacepqs() and nonfree() to get Tab coloured yellow to emphasize the periodic circles in the central region, Centre right:
Limit set with Ta, Tb having values determined by μ(3/5), μ(8/13) showing two fractional parabolic traces can also be combined. Right: 1.91 – 0.1i, 3
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.
Fig 22: Transition from periodicity to degeneracy: Non-free versions of μ(3/5), μ(8/13), μ(21/34), μ(89/144) and μ(144/233) tending to the Golden Mean. For values see table above.
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 TabAB = 1 can have the solution (abAB)6 = I, just as Ta = 1 did previously. This presents another type of generator solution than TabAB = 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 TabAB = 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 TabAB = 0 automaton works well enough in practice for this case to depict the limit sets reasonably well, although convergence is limited at high resolution.
Fig 26b Above: A series of examples for TabAB=1 -- (abAB)6 = I including periodic 1/2, 1/7, 1/19, close to degenerate 144/233 and a quasi-Fuchsian dragon. Below: A series for TabAB = -1 -- (abAB)3 = I including Ta = 1 a3 = I Tb=2, Ta=2 Tb = -(1+51/2)/2 b5 = I , Ta = 1.93-0.05i Tb = -(1+51/2)/2, gasket, Ta = μ(1/13) Tb = 2 and near degenerate Ta = μ(144/233) Tb = 2.
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 a3 = I and (abAB)3 = I , (2,3) for b3 = 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 b3 = I and abAB7 = I, the subject of help example 4 in the MAF package, also included in the m-files folder, the coloured example in fig 26c can be nicely depicted by a variant of the depth first search. This displays the symmetries reflective of the particular symmetry group in hyperbolic space. As noted in MAF, this limit set contains as a subgroup the von Dyck (2,3,7) triangle group (right on the Poincare disc)), the tiling of the hyperbolic plane by triangles: π/2, π/3, π/7. Its importance stems from its connection to Riemann surfaces of genus g with the largest possible order of its automorphism group, 84(g − 1). This triangle is the smallest figure that can tile H2.
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-1b-1 = AB and a-1. This is invertible since from I = aABb = aeb, we also have a = b-1e-1. Secondly, a congjugacy c = [1,1; -1,1] (or equivalently c = [1+ i, -1- i; -1+ i, -1+ i]/2) is applied to transform our new a and b into cac-1 and cbc-1.
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 TabAB = -2. We can use the Grandfather identity to show (a) T(c2) + Tc -1c = TcTc -1= (Tc)2 so T(c2) = (Tc)2–TI = (Tc)2– 2, so if Tc = -2, T(c2) = (-2)2 -2 = 2 = TI. Hence the parabolic case commutator squares to the identity. If Tc = 0 we get T(c2) = (0)2 – 2 = -2 so T(c4) = TI. Examining T(c3)+Tc -1c2 = Tc2Tc -1, we get T(c3) = Tc2Tc -1–Tc = ((Tc)2– 2)Tc – Tc, or T(c3) = (Tc)3– 3Tc. So when Tc = -1, T(c3) = -1 – 3(-1) = 2 = TI. When Tc = 1, T(c3) = 1 – 3 = -2, so T(c6) = TI.
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 a3 = 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).
Fig 26c2: The limit sets for Ta,Tb =2 TabAB = 1 and 0 in symmetric inverse normalization, showing the limit set adopts the rotational symmetry of the commutator square root (compare with fig 26 in standard grandma portrait). As the elliptic order of abAB increases, the central circle and annular limit set become larger. For parabolic TabAB = -2, it becomes linear, like the Jorgensen normalization.
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 L2 = k I. One can thus find a suitable x by normalising L2 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-1b-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.
Fig 26c: Upper row: Examples where b3 = I and (abAB)7 = ±I using the elementary MAF FSA for b3 = I and (abAB)7 = I in depthfirstFSA.c. First from left: The limit set of Ta = 1.91+0.05i, Tb = 1 with TabAB = -2cos(π/7) in grandma's normalisation. Second: The same set in symmetric inverse normalisation showing 14-fold symmetry. Third: Ta =1.93+0.05i, Tb = 1 with TabAB = -2cos(2π/7) in grandma's normalization, Fourth: The same limit set portrayed in symmetric inverse normalization showing the 7-foldsymmetry of the commutator using an additional rotated superimposed image to enhance the boundary. Fifth: Stochastic blow up showing the exterior domain of the set is inverted inside the centre of each whorl except here the inverted symmetry is three-fold. Sixth: The same parameters but using the alternate square root for Tab. The 7 and 14-fold symmetries result from the alternating signs of the solutions to the cubic equation for (abAB)7 = I, as noted in the solutions using recursion above. Lower row: Foreclosure of the interior by engufling towards degeneracy as the real part of Ta moves from 1.93 down to 1.87, coloured by the 30th interate.
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:
Fig 26d: Ta = 1.95+0.05i, ab is elliptic of order 7 i.e. (ab)7 = I and the square root of the commutator (rqp)7 = I, in symmetric inverse (center) and grandma normalizations (left and centre right for alternate roots) coloured by the 15th iterate, using depthfirstFSA.c. It is notable how different switching from symmetric normalisation to grandma is in figs 26c and 26d. Lower row: A sample of parameter space. The corresponding limit sets with Ta = 2, μ(1/2) and μ(144/233) with Tb adjusted so (ab)7 = I. Because Tb = 1, tracepq does not solve as a polynomial, so the usual μ relatioships don't apply.
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 kleinian1 producing a 79 step FSA:
_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.
Fig 26d2: Left and cantre: Ta= 2 and 1.5899+0.5156i, each with Tb chosen so that (ab)^3=I and (abAB)^7=I rendered using the p, q, r FSA. The the right image shows near degeneracy. Coloured by the 30th iterate. Right: Corresponding system Ta=1.9+0.05i , Tb chosen so that b^4=I and (abAB)^7=I rendered using the p, q, r FSA, coloured by the 25th iterate.
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]
]
);
Fig 26d2b: The tricorn Apollonian gasket generated by using Ta = 1+ i, Tb = 1 - i, TabAB = 1.
Both the commutator square root and ab are elliptic of order 3. Coloured by the 5th iterate.
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.
Fig 26d2c: Left: The third example lower in fig 26bTa = 1.93-0.05i Tb = -(1+51/2)/2 TabAB = -1 in symmetric inverse normalisation, showing it is a totally-disconnected Cantor set of fractal islands. Centre: Reduction of the real part of Ta to 1.898444-0.05i brings it close to degeneracy. At this point we can see it is also close to a value where a9B is parabolic. An adjustment of the imaginary part to 1.8987-5.0798i carries Ta very close to the μ(1/9) cusp value. Each on the 25th iterate..
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 Ta and Tb are shown in fig 26f.
One can also derive the limit sets for cases where each of Ta, Tb and Tc 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.
Fig 26f: Left: The equations to find the commutators used in fig 26e. Right: The limit set Ta = 2cos(π/3), Tb = 2cos(π/4), TabAB = –2cos(2π/5) displays all three symmetries, with the 3-symmetry in orange, the 4-symmetry in yellow and the 5-symmetry in green. The colouring of the limit set into red, magenta, cyan and blue represent the 15th iterate of each of the generators. he larger chunking of the iterates lower right by contrast with the finer upper left hints at the way how the composition of the depth-first algorithm's search topology and the symmetric inverse normalisation has wrapped around the limit set. This limit set involved a relatively huge MAF-generated automatic group table with 1234 non-short circuting assignments.
Examples of Alun Williams' elliptic limit sets can be viewed here, here and here.
Discreteness Enters the Chaotic Milieu
Fig 26i: (1–3) The output of the DFS algorithm in the Mac app when the real part of Ta in the example of fig 0 right is changed to 1.41891069, 1.51891069 and 1.71891069 respectively. 4–5 The FSA version of the DFS algorithm can be adapted to either attem pt to circumscribe a closed curve as the original DFS, by drawing short line segments whever points are closer than a small epsilon and when the depth limit is crossed. In the chaotic state this results in tangled lines as in (1–3). In the FSA version, a line is drawn only if the points are epsilon close. If epsilon is then set to the image's single pixel size, totally disconnected disconnected elliptic limit sets can be fully depicted. The former tends to fill in a few too many points while the tatter tends to draw a few too few.
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.
Fig 23: Left: Standard, Jorgensen and Maskit representations of Ta = 1.91 + 0.1i, Tb = 2. Research in the Maskit parametrization has focused on the two transforms a(z) = μ + 1/z, b(z) = z + 2. In the standard (grandma) normalization, this corresponds to the traces Ta = -iμ and Tb = 2. Right: Ta = sqrt(2)-2i, Tb = 2 appears to be a circle packing in an unstable island off the Maskit slice, so is an example of a non-free group.
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/Tab2; Ta Tb/Tab] and b = [Tb-Ta/Tab -Tb/Tab2; Tb TaTab]. In Maskit, grandma's recipe is conjugated as follows: Let p=-i/(21/2) and c=[(1-p2)/(ip) -p; p ip]. We then let a=c*a*c-1 b=c*b*c-1.
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 = 21/2, with a 6 = I and a 8 = I respectively even under variations of the complementary trace, including doubly non-free examples where (aba-1b-1)2 = I as well.
Fig 24: A variety oftochastic sketches of non-free limit sets induced by choosing traces within the Maskit slices.These can also be generated more accurately using the DFS, but require the correct automation and sometimes changes to the DFS algorithm to function. The coefficients are as shown with 1/n corresponding to the non-free traces shown below. The non-free example has Ta=21/2, Tb=2 with command qFstochastic(sqrt(2),2,10,10,1100,0.4,2,2,nonfree(sqrt(2),2,0,0));. Intriguingly the real solutions to 1/n for the non-free case T(aba-1b-1) = 0 and the standard case T(aba-1b-1) = -2 are the same, although all the other complex traces calculated by tracepq() differ.
There is an interesting pattern to the non-free trace values giving an = I. The solutions to this have a close relationship with those for μ(1/n). Both apply the recursive relation a{1} = 2, a{2} = Ta+ki, a{i} = a{1} a{i-1} - a{i-2} solving for a{n} = 2 = TI where k = 0 for solving an = I and 2 for the standard μ(1/n) process with commutator trace -2. There turn out to be ceil(n/2)-1 real solutions, Ta, each of which is a trace solution for an = I, the most negative of which gives a discrete non-free limit set. Corresponding positive values arise from the alternate root, thus defining points on the real line in the Maskit slice. We can compare the 2 recursions and see the underlying ceil(n/2)-1 pattern in Maltab, using firstly , recursean.m and then recurscoeffs.m
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 a1= I. The even n values have double solutions t = ± 2 and odd values just 2 corresponding to a1 = I. The remaining roots give pared solutions to an = I and its factors of n. For example n = 9 has the solutions for 1, 3 and 9.
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)k2cos(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, -21/2 = -1.414), (5, -(1+51/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 b7 = I and c9 = I, a8 = -I so we get a16 = 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.
Fig 25: Varied stochastic non-discrete examplesTa=1.78+0.05i, Tb=3, Ta=1.6, Tb=1.6, and Ta=1.85+0.05i, Tb=1.85-0.2i with TabAB=-2
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 TabAB = 1.
Fig 26: The function non-free() can also be used with a fourth parameter to find Tab for any value of the commutator trace TabAB. Above are shown a series of values spanning the real line from -2.5 to 1.99 showing the evolution of the limit sets towards non-discrete chaotic instability. There is a transition into chaos as we cross 1. The last value 1.99 approaches 2 where the extended Grandma recipe becomes singular and returns 0/0 Nans or infinities for a, b. The scale also fluctuates extensively as we move along the sequence. Two neighbouring complex trace values are also shown to illustrate the dynamics when the commutator is loxodromic. One should note that the traces do not form a spectrum as such, but rather traces of -1, 0 and 1 depend on powers of the commutator which equal the identity ( 3, 2 and 6 respectively), however we will see in fig 26g that even when we pick a series of values all of which are elliptic commutator solutions the same situation pertians. Both the DFS and stochastic algorithms are dependent on generators which are defined using Grandma's recipe, which is only one of several normalizations for Mobius transformations.
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.
Fig 26g: Corresponding limit sets to fig 26 for the elliptic commutator values 2cos(nπ/11) n = 1-5 show trends towards chaos as we approach 2, and display mode-locking instability driven by lower 5-fold and 7-fold periodicities as shown far right. Although the upper image for -0.8308 looks reasonably discrete, the variant top right with similar parameters to fig 26c (Ta = 1.98+0.05i, Tb = 1) has clearly become 7-fold mode-locked instead of 11-fold periodic and the lower image 5-fold ~ 11/2-fold mode-locked.
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.
Fig 27: Ta with real and imaginary parts 1.912, 1/30 and 1.74, 0.24
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.
Fig 28: 3D Plots of quasi3D
To view a much more thrilling 3D example, see Jos Ley's fly through.
