## Exploding the Dark Heart of Chaos

Chris King March 2009 Last Revised Jul 2016

PDF

Mathematics Department

University of Auckland

### Contents

Fig 1: Manifestation of chaos in population boom and bust. Centre: the bifurcation diagram with the bulbs of the

corresponding Mandelbrot set, showing identical periodicities. Top row: Corresponding Julia sets as boundary
between inner and outer basins highlighted by their level sets. Right Period 2 Julia set by modified inverse iteration.

Bottom row. Dynamics by iterating y = f(x) and x = y. At 4.5 the domain fragments into a fractal Cantor set.

### 1. Real ÔBoom and BustÕ Origins of the Dark Heart

To understand the basis of the black heart in complex numbers it is useful to start with a straightforward process in real numbers, which is also at the heart of both chaos theory and world futures – population boom and bust. Suppose we have a population that breeds each spring according to exponential natural growth, so that the next seasonÕs population would be . If we also limit the growth to be proportional to the remaining food area , we end up with the so-called logistic iteration  [1.1].

Sound (.wav) file sweep of the dynamics of fig 1 from period 2 at 3.36 to deep chaos at 4

We can follow this process graphically, by first calculating  and then letting , as shown in the various small graphs in fig 1(b). As the parameter r varies upwards from 1, a whole series of changes occur in the dynamics at discrete values called bifurcation points. Firstly the population reaches equilibrium – a fixed point where , then at  we suddenly find a good year and a lean year – the birth of a period 2 attractor, followed by an infinite sequence of period doublings 4, 8, 16, 32 etc. until in the limit at , we get our first sight of classical chaos – the iteration, although deterministic, appears to wander over the domain in a seemingly random manner. If we plot the limit points the iteration ends up in, for increasing values of r, we get the Feigenbaum bifurcation diagram of attractor form, with the repeated pitchforks and dappled chaotic fan in fig 1(a). The third prong of the pitchfork is the old period 1 fixed point attractor, which has now become repelling towards the new period 2 attractor, as shown dotted in the first fork.

We can tell we are in chaos, because the process displays the Ôbutterfly effectÕ, coined by Edward Lorenz – that if the weather is chaotic, the subtle disturbance of a butterfly flying in the botanical gardens in Hawaii could become the seed of a tropical cyclone hitting Fiji – otherwise known as sensitivity on initial conditions – that arbitrarily small initial differences can lead to global divergences. We can even plot the average rate of exponential spreading, using the Lyapunov exponent, which is positive in chaos, but plunges negatively towards  as the ordered periodicities move from attracting to super-attracting fig 1(a).

Fig 2: Above: Detail of the bifurcation diagram shows the period 3 window (and every other window) has on the left an

abrupt transition from chaos by intermittency crisis and on the right a period doubling sequence (3, 6, 9 É). The importance of the critical points is shown by plotting successive graphs of iterates  of the critical value with varying r, which correspond to the darkened curves in the chaotic region because orbits are drawn together there (inset). Below: The explosion of the complex Julia set from a closed curve to disconnected Cantor dust at a single point, as we exit the right hand cusp of the cardioid of fig 3 shows intermittency crisis at work in the complex case.

In the midst of chaos there are also a series of windows of other periods. In the centre right is a period 3 window, followed to the left by a sequence of all odd periods 3, 5, 7 etc and then, in a reflected series of period doublings, series of even multiples 6, 10, 14 etc. These fall into a number sequence determined by SarkovskiiÕs theorem, which applies to all continuous real functions: Attracting periods further to the right in the following series of all natural numbers have hidden periodicities of all the values to the left:                                                                                          [1.2].

Each window has an abrupt transition from chaos to order, followed as it closes by an infinite period doubling sequence into chaos again. Finally we see our first sign of fractals fig 1(c). When we go past , the process becomes unstable, because a central interval of points escapes [0,1] in the first step and tends to , and the two smaller remaining intervals each have a sub-interval escaping 2 steps, and so on, so that, in the limit, we are left with a totally disconnected residual fractal, called a Cantor set, which has an endlessly repeating motif of self-similarity, in which any part has the same recursive structure as the whole, in the manner of a snowflake.

The question naturally arises: Why are we getting so many different dynamical answers here? Why so many systems doing apparently different things? Is there no natural unification at all, in this mix of order and chaos? This is where the beauty of complex numbers enters into the picture.

In the top row of fig 1 are the corresponding chaotic Julia sets of the same complex number iteration, with their x-coordinates vertical. For period 1 the Julia set intersects the real line transversely at only two points, one repelling fixed point and the point which reflects into it, leaving all the other real points attracted to the attracting fixed point in an ordered manner. By contrast for r = 4, the Julia set has become an interval lying directly along the x-axis, resulting in real chaos on the number line.

Fig 3: Lightening up the dark heart of the Mandelbrot set of : Click images to enlarge (Centre-right) illuminating the closest approach the iterates of the origin (critical point) make to the origin inside the set, as well as its exterior level sets, surrounded by some of its associated Julia sets. (Centre top-left) Partitioning into which iterate of the critical point makes the closest approach. (Center lower-left) Internal level sets of periodic attractors. Left: Modified inverse iteration portrait of Julia set 2. The other surrounding sets are distance estimator method portraits. Applications to generate these high resolution B&W images are included in the Dark Heart Package..

### 2. Illuminating the Writhing Dark Heart of Complexity

Dark Heart Viewer: Application and Source

If we iterate using the plane of complex numbers, instead of the real number line, we find that every possible r value has a fractal set of points on which it is chaotic, called the Julia set, surrounded by one or more ordered basins tending to fixed and periodic points, as shown in fig 3, for the quadratic function                                                                                                                             [2.1].

The Mandelbrot set forms an atlas, or Ôparameter planeÕ of the Julia sets plotted by the level set method, which we will explain below. Inverse iteration turns the repelling Julia set into an attractor by solving  [2.2] repeatedly, starting from a repelling fixed point, (there is always at least one) to get  points on . It can be further refined by focusing on points that are weakly attracting, using the derivative to eliminate repeated and strongly attracting iterates. This helps get around the fact that despite a Julia set being repelling overall, there are exponential variations in probability density of the iteration, which causes the inverse process to fail to represent weakly attracting areas, no matter how long the iterative method runs.

Now we get to the topological nub of the problem! Some Julia sets are connected as in 1 and 4 in fig 3, and have internal basins tending to a period 1 and 3 attractor respectively shown as dots.  However others are totally disconnected Cantor sets – 2D versions of the one we saw in fig1(c), illustrated by 5 and 10.   Continuity of the Julia sets is the question that led to the discovery of the Mandelbrot set as the atlas set of c values for which the Julia set is connected, as shown in fig 3. At first sight, it is far from obvious how one could possibly characterize the set of points for which a complex 2-dimensional fractal will be topologically connected, and it took from Julia in the 1920s to Mandelbrot in the 1980s for the full fractal nature of the black heart to become apparent, instigated by an ASCII image in *Õs the actual discoverer sent to him.

We can derive the central period 1 cardioid of this set fairly easily.  For obvious reasons to do with net shrinkage of neighbouring points to the attractor, it turns out that a periodic attractor has absolute derivative  [2.3], so if we look for the boundary of where we have period 1 – i.e. , with  [2.4], and solve for c in terms of z, where we get  [2.5], the equation for the cardioid. We can also immediately see that the positions of the bases of the period-q bulbs will correspond to , giving their positions as  [2.5a].

To find the equation of the period 2 bulb, we look for a fixed point of the second iterate: . Solving for c we have:

or

The first solution is the period 1 cardioid since . The second solution gives .

[2.6]

A circle of radius ¼ centered at -1 as per fig 3.

However to get the full picture, we need to figure how to turn a deeply enigmatic topological criterion into a computable algorithm, and here is where the heart of fractal complexity really hit the fan.

If we consider [2.7], we find that numbers outside the unit circle get larger with the square of the radius (modulus) and tend to , while those inside tend to 0.  The points on the unit circle are trapped at modulus 1 and become the Julia set, forming a repelling crest between the two basins of the attractors 0 and . Since the angle (argument) of each point is doubled under [2.7], on the unit circle we have [2.8], and after n steps .

This process is chaotic, satisfying three chaotic properties [2.9]:

(1) Sensitive dependence: Since angle is repeatedly doubling, the process has the butterfly effect – there exist points arbitrarily close which get mapped globally apart after a finite number of steps through.

(2) Topological transitivity: Any small open region (angular interval) is eventually expanded over the whole circle overlapping every other small open region so it has topological mixing.

(3) Dense periodicities: The dynamics is permeated with (repelling) periodic points as close as we like to any other point. Solving  [2.10], we get , or  [2.11], i.e.  evenly spaced points of period n, which are all repelling, since .

These are the three axioms of chaos demonstrated.

Fig 4: Inverse images of a large circle in (a) the connected and (b) disconnected cases (Peitgen et. al. 1988).

Now letÕs go back to  and look at the function on a very large circle, as in fig 4. Here  is much larger than c, so itÕs essentially mapping a circle to an even larger circle. The point at infinity is a super-attracting fixed point and the Julia set is the boundary of its basin of attraction. Now take successive inverse images of this circle. By continuity, these must be continuous closed curves and they have to come to  one of 3 cases: (i) a single topological circle (simple closed curve) with two points rotating around it as the two solutions of   [2.12] as c varies on the original circle, or (ii) two separate circles each with one solution or (iii) a transitional case of a figure 8 (or figure !). We canÕt have any other cases like two circles overlapping at 2 points because this would cause a repeated root at each crossing and hence 4 roots when we only have a degree 2 quadratic with 2 roots, whose repeated root, from [2.12], is 0.

If the inverse images are all topological circles, then the filled in Julia set  must be connected, as in fig 4(a), as the circles and their inverse images form an effective conformal mapping of the plane minus a disc onto the plane minus the filled in Julia set. On the other hand, if we do have a first figure 8, we have a repeated root i.e. , and  [2.13], so c is on the last circle whose inverse image (a figure 8) is connected and cÕs inverse image is 0 on the figure 8, as shown in fig 4(b). In the connected case we never reach c, so it is inside the filled in Julia set and thus 0 must stay confined inside when iterated.  By contrast, in the disconnected case, 0 will iterate to  and eventually escape any fixed large circle. Notice that is the critical point of , as  [2.14]. The critical point is the acid test because it is the last point to escape iterations, since the function is stationary there.

Fig 5: (Inset) The Mandelbrot set can be proved connected by extending the mapping beyond c to the entire complement of M, resulting in a conformal mapping between the complement of M and the complement of a disc, corresponding to radii and angles of the standard potential function around the unit disc. This is illustrated by using a binary decomposition in which the level sets are chequerboarded by whether the escaping n-th iterate is in the upper or lower half plane - i.e. has positive or negative imaginary part. The decomposition external angles being ivervse powers of 2 all lead to eventually fixed end points, except for the period 1 main cusp. The external angles of given periodicities can be followed by moving left and right e.g. period 2 cusp 1/3 folows 01 repeating period 3 cusp 1/7 is 001 repeating 0=L, 1=R.

The quadratic Mandelbrot set can also be shown to be a connected set by extending the above conformal mapping arguments for Julia sets in the region outside c to a conformal mapping over all regions where c is outside M, as shown in fig 5.

We now have a computable method for finding the Mandelbrot set of any one-parameter quadratic family: Iterate [2.15], for each point c in the complex plane. If it remains inside a given large circle after a fixed number of iterations, colour the point black, otherwise colour it with the number n of steps to escape. This is the level set method. It leaves the Mandelbrot set black with its complement a coloured topographical dynamic landscape – hence the Ôdark heartÕ. It is also used to picture the exterior of a Julia set by starting from each point  in the plane for fixed c.  We will use this method to investigate the properties of Julia sets and the Mandelbrot set of a variety of functions. One can also, with some greater difficulty, attempt to shade the interior basins of the Mandelbrot and Julia sets, by counting how many steps the iteration takes to get within a small neighbourhood of any attracting periodic points, as illustrated in the Julia sets of fig 1 and the Mandelbrot set of fig 3.  However this fails if the Julia sets contain neutral points.

One can also colour the complement of Julia and Mandelbrot sets with a continuous potential, by comparing the potential of the n-th iterate escaping a large circle with the standard 2D logarithmic potential  [2.16], using the function

[2.17]

We can also picture the rays corresponding to external angles, , measured in fractions of a revolution, of the dynamic rays of potential function by colouring c black or white, according to whether the escaping  lands in the upper or lower half plane, called a binary decomposition. The  factor comes from the repeated angle doubling of  which is dominant on large circles. The method turns fig 5 into a checkerboard pattern highlighting segments of both rays and level curves.

Fig 6: (1) Parts more complex than the whole – the fractal black hearts, connected to the main body by dendrites, get warm and fuzzy on blowup, since the dendrites become locally space-filling as the Farey denominators tend to infinity. On descending fractal scales the Mandelbrot set boundary fills 2-D space and thus has fractal dimension 2. (2) Locally varying complexity. The Mandelbrot dynamics in each locality (a)  captures the dynamics of each corresponding Julia set (b). (3) Embedded Julia set (a) in the intermittency crisis region of the main period 3 island on the period 2 dendrite with (b) central small island. (c) Corresponding Julia set with fractal embedded homologs (d) each with a period 1 Julia set in their core (e) corresponding to the period 1 region of the small Mandelbrot island (b).

As an atlas of all the varied complex dynamics in the entire family, the Mandelbrot set is sometimes regarded as the most complicated mathematical object in existence. The local dynamics at any point c on the Mandelbrot atlas closely represent those of the corresponding Julia set  (Tan), as shown in fig 6.  Although it was only proven in 1991 (Shishikura), the Mandelbrot set boundary has fractal dimension 2, the same as the plane it sits in, essentially because the number of dendrites grows to , in such a way that the boundary locally forms a space-filling tree, because the Farey denominators (see next section) of the dendritic islands increase without bound, as we descend down fractal scales (fig 6). It is not known whether this results in a boundary of non-zero area.

A convenient way of completing the complex plane so that the iterations are consistent with a single point at infinity is to close it off into a sphere, called the Riemann sphere, with the single point at infinity at the north pole, the unit circle as the equator, and the origin at the south pole, mapping all points in the plane by stereographic projection from the north pole. This way the point at infinity can be treated like any other point.

If we now turn back to the equation we started with – the logistic iteration  [2.18] – although the global form of its Mandelbrot set is apparently different, both the period multiplying bulbs and the dendritic islands, and the local dynamics and forms of the Julia sets are all homologous with those of the standard function .

Fig 7: Although the parameter plane of  consists of two circles forming a pair of the quadratic period doubling bulbs on the standard Mandelbrot set, its dendritic islands have the same cardioid as .  The inside of the set is shaded by approach to its internal periodic attractors.

Dark Heart Viewer: Application and Source includes

To generate an image of the parameter plane we replace the origin with the correct critical point. Solving  for c in terms of  gives:

[2.19]

This gives two solutions representing the left and right unit circles in fig 7 centered on 0 and 2, respecting the essential symmetry in the function between z and 1 – z, while the roots are 0 and 1 and the critical point is 1/2. Each of these two major regions of the parameter plane is of course a perfect copy of each of the bulbs on the periphery of all versions of the black heart. We can derive the circular equation 2.6 for the period 2 bulb of  (see appendix 1), however, the higher periodicity bulbs, despite being obvious circles, have intractable solutions, because quadratic period n corresponds to solving a polynomial of degree , which is impossible classically for degree greater than 4.  Effectively the parametrization has created a bilinear map, giving a pair of degree 1 circular ÔheartsÕ, which only show their true quadratic nature on fractalization, but give us a quick window into the equations for the circular bulbs.

Intriguingly, the equation [2.5] connecting the unit circle to the cardioid,  is also expressed between the parameters of  and  in the form  [2.20], precisely conjugating the sequence of their bifurcations on the real axis. For example first period doubling  3 maps to -3/4 and the tip of the main dendrite maps 4 to -2.

### 3. The Dark HeartÕs Magic Numerology

In fig 8 are shown three classical routes from order to chaos.  The first â€“ period multiplying â€“ a generalization of period doubling, requires an infinite sequence of bifurcations, arriving in limit at a strange attractor (Schršder) and chaos. The strange attractor for the real mapping of fig 1 is shown as  and the corresponding Julia set above has a central black area highlighting this strange attractor, contrasting with the generally ÔrepulsiveÕ nature of the Julia set under the forward mapping. By contrast with the real number period doubling route, a journey to the dendrites of the dark heart can involve any sequence of periodicities determined by which bulb is entered. This complicates the expression of SarkovskiiÕs theorem, which also doesnÕt actually describe the full order of the periodicities, as we shall see shortly. As we move along any dendrite of any compound period, we will go through a transition similar to the real transition after  in fig 1. The dendrites have the base n periodicity of the local process, but as we move out along them, they are fractally permeated by a tree of dendritic Mandelbrot islands with varying fundamental periods, each of which is capable of generating new components of any periodicity on its bulbs.

Fig 8: (a) The three routes to chaos illustrated on the real axis period 3 dendritic island of the black heart of the level set method, (b) the Farey tree of winding numbers and dendrites of the bulbs,  (c) the DevilsÕ staircase of mode-locking.

The second – intermittency crisis – illustrated in fig 2, involves switching between periodicity and chaos in a single bifurcation, as shown in fig 1 at 3.8282 just to the left of period 3 because part of the graph of  crossing  has become tangent and then separated allowing intermittent slippage of the period 3 into chaos. The same transition is responsible for the exploding transcendental Julia sets we will look at later.

All dendritic processes on the Mandelbrot set go through repeated fractal transitions through black hearts, firstly following reverse intermittency-crisis, and then period-multiplying, as we move out along the fractal dendrites. All parts of the complex structures in the dendrites consist of fractal branching arcs of the periodicity of the bulb from which they protrude, interconnected by black hearts, where the periodicity on the dendrite undergoes further multiplication. In the related Julia set, the black hearts are replaced by structures of the periodicity, or transition dynamic, corresponding to the Julia set, also possessing the basic 180o rotational symmetry of quadratic Julia sets.

Fig 9: Mode-locking bulbs follow an decreasing sequence of fractional rotations by 1/n, with all the remaining rational fractions generated in between as mediants.

In the third process – mode-locking – any irrational local dynamic close enough to a rational fractional rotation undergoing dynamical feedback becomes locked to the periodicity of that rotation, forming a series of intervals of mode-locked states, with only a residual set of points in between retaining their unperturbed irrational motion. These mode locked intervals form a continuous fractal ascending (monotone increasing) function, constant on intervals surrounding rationals, called the DevilÕs staircase as in fig 8(c). Mode-locking is manifest in many processes where periodicities interact, including the non-mode-locked orbits (to Jupiter) of the remaining asteroids, because the mode-locked ones were thrown into chaotic orbits and collided with planets, the ordered mode-locking of the MoonÕs day to the month and MercuryÕs day to 2/3 of its year.

As shown in figs 3, 9 the bulbs on the Mandelbrot set follow the fractions on the Farey tree in fig 8(b), adding fractions as mediants  [3.1], as can be seen by counting the number of their dendrites, which also corresponds to the periodicity of the attractor corresponding to the main bulb.  Mediants correctly order the fractional rotations between 0 and 1 into an ascending sequence providing a way of finding the fraction with smallest denominator between any two other fractions.  A way of seeing why this is so is provided by using a discrete process to represent the periodicities or fractional rotations. For example if we have 2/3 =110 and combine it with 1/2=10 by interleaving, we get 11010, or 3/5.

Numbers like the Golden Mean  [3.2], which by virtue of its defining relation  [3.3], is the limit of the ratios of successive Fibonacci numbers 1, 1, 2, 3, 5, 8, 13, 21 É for which  [3.4] are non-mode-locked. The Farey tree leads to Golden Mean numbers, if we alternate left and right as we go down, following a series of Fibonacci fractions. Numbers g, such as these, avoid becoming mode-locked because their distance from any fraction of a given denominator q exceeds a certain bound:

[3.5].

The golden numbers can most easily be described in terms of continued fractions, which, when truncated represent the closest approximation by rationals:

[3.6].

Golden numbers end in a series of 1Õs thus having slower convergence to fractions of a given denominator than any other numbers.  The Golden Mean itself is  [3.7]. More generally, the Farey Tree has straightforward natural rules of parental and descendent inheritance, using continued fractions and any fraction or quadratic rational has a repeating continued fraction (Schršder).

All the Julia sets on the cardioid boundary are caught on bifurcation out of period 1 and fall into rational and irrational types. The rationally indifferent or mode-locked Julia sets such as 2 and 8 in fig 3 (called parabolic) have neutral points which are in the Julia set (compare the attracting fixed point in 1), because they have n repelling branches, as well as n attracting petals, forming an n-branched saddle point. After bifurcation the n petals form the components of a period n basin, corresponding to the periodicity of the adjacent bulb, as in 4. By contrast the irrationally indifferent Julia sets have a neutral point in their internal basin, around which all the points in the basin are irrationally rotating.

One of the most challenging features of the Mandelbrot set is the bewildering arrangement of the dendritic Mandelbrot ÔislandsÕ which are strung like IndraÕs beads along the dendrites, playing a key role in introducing higher periodicities and causing the fractal dimension of the boundary to be space-filling in the plane.

One way of visualizing the precise locations of all the 2n periodic bulbs and dendritic islands of a given period n is to use the wave function method described in appendix 1 with a radial function coloured by sines and cosines of a Gaussian distribution, viz .

Fig 9b: Locations of the periodic bulbs and Mandelbrot dendritic islands for periods n=2-12, located by Gaussian wave function method. Images of non-prime n include all factor periods. The top-left image shows a superposition of periods 1-101.

Full size Quicktime movie of the periodic loci of the first 101 periods of the Mandelbrot set

Full size image of the superimposed wave Mandelbrot top right in fig 9b
Matlab m-file of: Gaussian Wave Periods Mandelbrot

Douady and Hubbard discovered a number of intriguing relationships depending on the critical point that index the features of the bulbs and their dendrites in terms of the critical points. These depend on the external angles, the angle of the ray from the unit disc that corresponds to each ray in figs 5,10. There are intriguing ways to calculate these angles both for the bulb cusps and key points on the dendrites.

Inside each period-n bulb is a central super-attracting point, where  [3.7] for any point on the period-n cycle. But this means that one of the derivatives  has to be zero, so  the critical point. Hence the critical point is periodic with period n at this c value. Hence there is a point c in each bulb where , resulting in an equation of degree . For example for period 3, we get   [3.8] , giving the period 1 heart and the three locations of the small period 3 dendritic Mandelbrot island on the main dendrite and the two period 3 bulbs above and below the main heart.

Once we exit the bulbs and enter the dendrites, the fundamental periods become repelling and the critical point becomes eventually periodic to a repelling periodic orbit of period k after n steps, at Misiurewicz points , the tips and branching points of all dendrites. For example for  we have  [3.9], giving c=0 (where the critical point is fixed) and c=-2 the tip of the main dendrite.

Fig 10: A collection of dynamic rays with their external angles (Pastor et. al., Romera et. al.).

It is possible to calculate the external angles of a variety of key points on the boundary of the Mandelbrot set by an ingenious method which depends on the squaring nature of the quadratic to form exterior angles over powers of 2 using the upper and lower half plane as a discriminating target as noted in the binary decomposition.

If we consider a critical point in either a periodic orbit or an eventually periodic orbit, we can encode the dynamics as a binary ÔdecimalÕ setting a 0 for each iterate in the upper half of the local traverse of the 0/1 and 1/2 rays, which are asymptotic to  the positive and negative real axis. When we expand this decimal as a power series it gives the external angles of the cusps of the bulbs and dendritic Mandelbrot islands and the Misiurewicz dendrite tips and branching points.

Fig 11: Calculation of the external angles of several key period-3 points. We start with the first iterate (1) and assign 0 or 1 if a given iterate is above or below the traversal between the connections between the 0/1 ray and the 1/2 ray: (a) The at the tip of the largest period 3 dendrite. (b) The  at the tip of the second period 3 dendrite.  (c) the main branching point of the period 3 dendrite. (d) The super-attracting period of the period 3 bulb gives the exterior angles of the two cusps, confirmed by cusp Julia set rays above.  (e) The exterior angles of the cardioid main cusp on the period 3 Mandelbrot island on the period 2 dendrite. The ray 0/1 is the preimage of edge 01. Ray 1/2  is the preimage of 0/1. In (b) 03 is the preimage of 14 so the preimage of 01 must extend from 03 at an angle. In (c) the ambiguity in the orbit [001***] provides the third exterior angle.

Rays emerging from cusps and dendritic islands all have external angles fractions of the form  [3.10] because of the repeated binary decimal associated with their periodicity. E.g.   [3.11].  Furthermore these lead to every odd denominator fraction as a result of EulerÕs generalization of FermatÕs little theorem  where a, n are coprime and  is the number of integers less than n coprime to n. For example for 9 we have 1, 2, 4, 5, 7, 8 give , so  and  is actually an external angle of the period 6 bulb.

By contrast, those from Misiurewicz points have even denominators because the non-periodic initial iterations cause an irreducible power of 2 in the denominator.

There are also some intriguing irrational angles. The sequence of period doubling points have sequences  [3.12], tending to the Morse-Thule fractal sequence in, which we invert each binary decimal and append it to the right. There are a corresponding sequence of rays emerging from the Misiurewicz points ,  [3.13] also tending to the same limit equivalent to  in fig 1. Further codings can be found in Pastor (2002).

From SarkovskiiÕs theorem, we gained an idea of the way the periodic windows form odd and period-doubling series, which from fig 1 correspond to dendritic Mandelbrot islands in the complex case. However SarkovskiiÕs theorem doesnÕt actually tell us which period windows occur and in which order and there are clearly more than the sequence indicates, for example a further period 5 window to the right of 3. In fact the distribution of these is encoded in the eventual (k,n)-periodicities of the Misiurewicz points, in turn encoded by rational external angles. However we need to have a more direct idea of where these periodicities occur.

Fig 12: Several processes in interplay on the period doubling dendrite determine the order of the periodicities. Yellow links show the first members of the odd Sarkovskii series 3, 5, 7, 9 and in mauve the first two of the period doubled odds 6, 10. Green, linked by increasing steps of 1 from 2 to 10 is the main period winding sequence. In orange is a mirror odd sequence 3, 5 leading to a further fractal winding sequence. Blue links are mediant-generated periodicities as in mode-locking. The winding sequences are generated through period 2 dividing the periods into even (0 mod 2) and odd (1 mod 2) creating two ascending series of Mandelbrot dendritic islands with base periodicities covering all the natural numbers above 2 and their fractally entwined recursions repeating the sequence from higher levels.

To make some sense of this while writing the paper, I investigated the base periods of the collection of dendritic Mandelbrot islands illustrated in figs 12 and 13. The sequence of periodicities in fig 12 standing alone are enigmatic and appear both random and repetitious, however once we elucidate the patterns in periods 3 and 5, we can begin to see both a generalization and several interpenetrating processes determining the overall structure.

Firstly, instead of period 2Õs odd series and even doublings, we have n classes corresponding to the residues mod n, corresponding to the n fractal dendrites, each of which has an ascending sequence in steps of n with the k-th residue remaining constant on a given arm, forming  [3.8]. The base arm consists of the multiples of the base periodicity. This means that an island exists with every base period above the base bulbÕs periodicity, winding in to the hub of the main set of n spokes. A second reflected series occurs running out along each dendrite from the dominant island.  In addition mediant periodicity addition occurs, resulting in fractal extension of the process, as illustrated in fig 13.

Fig 13: Period 3 and 5 dendrites (1/3 and 2/5 a revolution) show the fractal effects of mediant formation and period winding with periods of all orders, strung in ascending series following the base period n of the bulb, in such a manner that the k-th of the n arms are composed of periodicities , with the base periodicity and first period doubling bulb forming members of the sequence. (Left) Period 3 showing in yellow mod(3,1) and mod(3,2) branches 4, 7, 10 and 5, 8, 11 and the periods 6, 9, and 12 on the main dendrite forming a period doubled set. In orange a reflected n step series running outward. In cyan, additive period mediants, 11 of 4 and 7 and 13 of 5 and 8, are located between the associated triple branching and the next Mandelbrot dendritic island, with the n-1 other associated periods forming a fractal version of the main winding sequence. Far left the corresponding sequences of the period 9 sub-bulb follow the same pattern. (Right) Fractal winding sequences on period 5 dendrite, with 2/5 of a revolution, each rotate by two vertices, forming a pentagram spiral, with fractal recursions on each branch (cyan and green). Right: The same period 3 sequence at left depicted directly by period colour using Dark Heart Package.

Given the ray angles forming a continuous circle, there are clearly points on the Mandelbrot set corresponding to every possible binary ÔdecimalÕ angle. So we can characterize all points on the boundary in terms of their sequences, periodic, eventually periodic or irrational. Clearly although the rational numbers are dense the irrational numbers constitute real rather than countable cardinality and thus the countable collections of bulbs, islands, and M-points we are examining only form the occasional stars in a mysterious irrational sky of boundary points, with a few additional irrational limit points thrown in.

If we use the binary sequence corresponding to a given ray, we can immediately see sequences which correspond to Mandelbrot islands, decoding their ordering and periodicity.

For example beginning with the period 2 bulb and its dendrite, we have the sequence of periods 2, 3, 4,É, n, stepping from the period 2 bulb to the main period 3 island and then out, converging to the eventually periodic tip point : . Each of these is a period n island because it is periodic and its external angle places it on the period 2 dendrite.

There are then period doubling versions of this sequence of periods 4, 6, 10, É , 2n converging to the  point: . Although this point has eventual periodicity 1, the map is flipping half-planes here, so it is the period 2 Ôband mergingÕ point with sequence  corresponding to the main M-point junction of n dendrites of a period n bulb. This period doubling sequence of sequences continues, in steps of , with M-point limits converging to the Morse-Thule sequence of the period doubling  point, as noted above.

We can also see periodicities 3, 5, 7, É , 2n+1 forming a sequence of islands to the right of the period 3 island converging to the same Ôband mergingÕ point: , and again period doubled versions of these, consistent with the bifurcation diagram of fig 1.

If we put this together with the even sequence , which also converges to the same Ôband mergingÕ point, we get the main winding sequence of periods 1, 2, 3, 4, É, n of fig 12: .

Fig 13a: Left Periods 1, 2, 3, 4  (black, red, green and blue rays) on the circle map diagram, with the sectors clamed by periods 2 and 3 bulbs shaded. Looking at period 4 clockwise from 0, the first two rays form the cusps of the main period 4 bulb, the next 2 from the period 4 island on the main period 3 dendrite, the next is period 2, leaving two ray pairs in the period 2 sector, the period 4 period doubling rays and the rays from the one period 4 island at the tip of the dendrite.

We can now begin to think about finding the position of a given ray sequence on the boundary of the Mandelbrot set. Now to see where these actually might be we need to consider the way the sequence of lower periodicity bulbs and islands carve out ray parameter space in a manner similar to the devilÕs staircase. The easiest way to see this is to look at the disc from which the rays emerge and to consider the periodicities in terms of the angle doubling circle map corresponding to . Recall from [2.11] that this has a dense set of periodic points of period n, with angles  coinciding with our external angles. For example for n=2 there are two period 2 points  and , mapping to one another, and 1 which is period 1, together forming an equilateral triangle about the origin. Moreover the eventually periodic points to any period are also dense. For example there are  points eventually fixed to 1, since .

In fig 13a, period 1 is in black, period 2 in red, 3 in green and 4 in blue. We can now analyse the positions of all the period 4 bulbs and islands. Counting round clockwise from 0 as in fig 13a, we find firstly the arcs from the main cardoid period 4 bulb (the only one on the cardioid of this period), next because the two arcs ar in the period 3 sector this must be the period 4 island on the main dendrite, next is a period 2 ray, followed by two ray pairs centered on ½, forming firstly the arcs above and below the first period doubling period 4 bulb and second those emerging above and below the main cusp of the one period 4 island on the period 2 dendrite.

More generally, we can decode any period-n sequence by examining whether the ray (or ray pair) falls within any of the sectors of the main bulbs of periods 1, 2, É, n-1. If not it is a main bulb. If it does fall in a lower period sector, we can analyse itÕs address further, using key landmarks on the bulb concerned.  For example, we know the period 4 island is on the main dendrite because, looking at the landmark sequences for the period 3 bulb with cusp rays  and we see the main branching M-point has lateral rays  and .  The first period doubling rays are likewise  and .  Further landmarks noted include the tips and the other main M-point ray. We can immediately see that the period 4 rays  and  lie on the largest dendrite beyond the principal M-point.

The M-points can be treated slightly differently. They are all on dendrites, but can be on dendrites of any period, so we need to decode their position sequentially, using the first n-digits to look for ranges in a recursive set of sectors corresponding to the generalized  points, which is made more problematic by the fact that these points have irrational external angles. On the other hand, as noted in fig 13a the addresses of the principal M-points do encode precisely where they are on the tree, so one can alternatively analyse the structure of the eventual periodicity of the M-point to decode its physical address. Looking at these on the period 3 bulb, we see that the principal tip is  which is precisely what we get if we reduce the periodicity of the cusp rays to  eventual periodicity 1, the secondary tip is  which is the eventual period 2 reduction for the lower cusp ray , and the main branching M-point was analysed above.

Fig 13b: Julia sets showing (magenta) the varieties of rational and irrational iteration of the critical point. (a) Critical point eventually periodic 0 > -2 > 2 (fixed). (b) Chaotic interval subset of dendritic Julia set (b2) Chaotic dynamics on strange attractor of (b). (c) Neutral points on bulb cusp point, with critical orbit stratified in petals (d,e) Irrational recurrent cycles (Siegel disc) (f) Super-attracting critical point on period 3 (periodic) (g,h) Asymptotically attracting to periods 3 and 1. (i) Seahorse near dendrite with escaping orbit (inset). (j) Fractal strange attractor (cyan) at limit of period doubling on pd 3 bulb. For an image of the collective iteration of all critical points see fig 35. For preimages, and Julia dynamics fig 40, 41.

Quicktime movie of modified inverse iteration working

The variety of dynamics of irrational orbits can better be understood by iterating the critical point and Julia set together. In fig 13b are a variety of periodic and eventually periodic orbits with rational external angles as well as boundary points with recurrent cyclic orbits which are not periodic but densely fill a closed curve, as well as chaotic intervals and fractal strange attractors at the period doubling limit point. The critical point can thus be caught in attracting internal basins, on petals at bifurcation on bulb cusps, be caught in the Julia set itself in chaotic or eventually periodic orbits, and finally can escape to infinity after a sometimes long transient finite orbit near the Julia set. A further dynamic way of viewing the active dynamics in these Julia sets is the wave function method in appendix 1. For a collective image of the iterations of all critical points see fig 35.

### 4. A Heart is only Half a Hamburger: Degenerate Critical Points

Dark Heart Viewer: Application and Source includes

A favorite of the standard Mandelbrot function to higher polynomials is to consider simple extensions of the standard function to higher degree i.e.  [4.1]. This certainly gives a one-parameter family of degree d and the results are interesting and beautiful to look at.  It is also a trivial extension of deriving the standard cardioid to find the higher dimensional equivalent simply by solving  together to get  [4.2], which gives a series of dark heart generalizations, to a hamburger, a quartic trefoil of cusps and on to a (d–1)-petalled ÔpolygonalÕ cusp curve, as shown in fig 14(left). The case for degree 1, as noted in the logistic case, and in all the quadratic bulbs, would be a circle, however the above equation has trivial solution 0 for d = 1.

Fig 14: (Right) ÔHamburger MandelbrotÕ of  determines whether the Julia sets are connected or Cantor

sets, but these are simply ÔtriangularÕ versions of the quadratic sets. Although the main body

and all dendritic islands are hamburgers, the bulbs have now all become quadratic ÔheartsÕ.

(Left) The series of multi-cusp polygonal curves for higher powers of z.

We can see immediately on examining the cubic case in fig 14 that we have a hamburger main body, which is repeated in all dendritic islands, but the bulbs have all now become recursive quadratic hearts. The parameter plane is a defining measure of whether the associated Julia set is connected or a Cantor set, but the Julia sets are simply ÔtriangularÕ versions of the same series of Julia sets we saw in the quadratic case.  The reason for this is can be seen if we run the process backward by inverse iteration. We then find  [4.3] is producing 3 complex cube roots of each point symmetrically arranged around the origin in an equilateral triangle where the quadratic inverse iteration was producing 2.

Fig 14b: First three rows; Function, Mandelbrot and Julia views of ,  k = 2-5 and show the principal power atoms forming the stable critical points of a diversity of functions with main body boundaries following well-defined exponential curves.  Fourth row: Mandelbrot views of second order atoms which arise in a region in which two or more dark hearts of critical points overlap giving full continuity (next section). Left to right: , , both of which have a hidden infinite dark heart due to a second critical point a z = 0 mapping to 0, and  k = 2-4 (black inner bodies with full continuity) consist of a fractal curve with quadratic cardioids emerging from the cusps and fractal symmetries of the order k, rather than the k-dimensional bulbs and well-defined body boundary curves of the upper sequence. Notice the last three have outer bodies in the k +1 form but bulbs in the k-form due to the unusual c parametrization.

These form the root possible single critical point kernels which appear in all more complicated examples to follow and they represent the only kinds of solutions we can have, except for an exponential plateau and the inversions of these solutions in the presence of singularities, such as those caused by negative powers of z. In fig 14b are shown each of these along with the period 6 Julia sets of a period 3 sub-bulb of a period 2 bulb on each set to show their dynamics are essentially similar except for the k-fold geometry. All of these have well-defined boundary curves of their central region, so Dark Heart has a script which can generate interesting movies of the changes in the Julia sets as we move around the boundary, or out of the cusps.

.

Transitions between parameter planes of c(z+k)z(z-1).

In the above video, firstly we move from k ~ 2 to 1, passing a crisis in cz(z-1)2, cz2(z-1) and its rocky coastline and cz(z+1)(z-1). Then we spin twice on a tiny circle around the crisis at k = -1 to see the arms rotate. Then we zoom into a bulb and follow the Julia sets as they cross the rays on the bulb, zooming in again to see the tiny rotational disconnections at the very centre at z = 1 explaining the crisis. Finally when we view the Julia sets in full, the disconnections are too small to be seen.

Fig 14c: The parameter planes of c(z+k)z(z-1), passing (a) cz(z-1)2, (b) cz2(z-1) and (c) cz(z+1)(z-1) display radically different dynamics.

Fig 14c shows three key parameter planes from the video above. While in (c), both critical points give the same parameter plane, resulting in first order bulbs although there are two separate criticals, in (b) there is a second order dark heart rocky coastline parameter plane where both criticals give rise to connected dynamics enclosed in a first order plane with double bulbs. By contrast, in (a) there is a unique crisis around the double zero z = 1 again giving a critical that maps to 0 inducing an infinite parameter plane, but in this case every neighbouring value of k has a finite parameter plane with a rotational connectivity crisis caused by the Julia sets winding around the two closely spaced zeros, resulting in alternating connections and disconnections. In the limit at k = -1, although that of z = 1 is instantaneously infinite, due to the double zero, the combined plane (shown in the video above) has first-order bulbs because all adjacent k-values have a finite plane for the critical z ~ 1.

The reason for this overall profile is that, while the quadratic had a non-degenerate critical point, the single critical point of the higher powers of  are all degenerate, and in the odd powered case is neither a maximum nor a minimum. The single critical point guarantees the parameter plane does determine the connectivity of the Julia sets because we have only one critical point, but we are missing at least one degree of freedom in the solutions. We thus need to look at more general cubics with non-degenerate critical points.

Fig 15: Continuity (a) and Cantor (b) parameter planes of non-degenerate (a, b) have cardioid dendritic islands (3, 8, 11, 14) and circular bulbs (main figure and 10). Julia sets are Cantor outside both (4, 5, 7, 12) and connected inside both (1, 6, 15). In the yellow region between, we have 1 fractal process connected and one disconnected (2, 9, 13). Periodicities in the dark hearts of inner two-critical region a and the bulbs of one-critical region b are illustrated as indicated.

### 5. Hearts Layered as Onions of Complexity: Generalizing to Higher Polynomials

Dark Heart Viewer: Application and Source includes

by Gaussian wave function method

We now explore the universality of the heart construction, along with the Mandelbrot atlas, as we fan out to higher polynomials, rational functions and transcendental functions. We first examine the non-degenerate cubic  [5.1] illustrated in fig 15. The inclusion of –z results in a completely new situation (unlike the quadratic case where  simply gives a translate of the standard Mandelbrot).

has real critical points symmetrically about the origin  [5.2]. If we try to use either critical point to form a pseudo-parameter plane, fig 19(2), we get a distorted image comprising a left half like fig 15(a) and a right half like (b) and vice versa, similar to the artifacts when we iterate from a non-critical point in . The Julia sets also do not properly classify in either set.

Transition from czk to czk - z showing separation of the critical points and its effect on the Mandelbrot ensembles.

.

The solution is to form two nested parameter planes, the first (a) consisting of those c values for which all critical points remain finite  [5.3] and the other in which any remain finite  [5.4]. This gives us the two nested parameter planes a, b in fig 15. These now correctly classify the connectivity of the Julia sets and show that three situations pertain [5.5]: (i) connected  (1, 6, 15), (ii)  Cantor   (4, 5, 7, 12) and (iii) one fractal component is disconnected (2, 9, 13).

Repeated Root Theorem:  To see why, we need to generalize the proof of the inverse circles argument for the standard quadratic. For very large z, as before  is dominated by  and maps large circles to even larger circles.  Now if all the inverse images of circles are again topological circles we have a connected Julia set, whose complement is homeomorphic to the plane minus a disc. However we now have a situation in which we can have two successive figure 8 crossings each corresponding to a repeated root of the inverse iteration  (or even in a degenerate case, such as the previous cubic, a triple point).

Fig 16: Stages in the bifurcation of a large closed curve under cubic inverse iteration

We need to try to solve [5.6] for repeated roots . To avoid having to solve tedious cubics, let us examine whether the critical points represent repeated roots of f , by checking if their image has a repeated root.  Substituting into f we have  [5.7].

Now putting this into the above, we have                   [5.8], or

[5.9].

Fig 17: (Left) Interference between  and . Because of critical point sensitivity, the dynamics may become confused below the level of computability. Nevertheless the layered Mandelbrot method correctly classifies connectivity of the Julia sets, as shown in the disconnection between a and b. Blue and mauve regions are in  (Right) Fractal structures showing sensitivity to the period 1 region of the other critical point are linked to the whole by black hearts through immediate intermittency crisis. The Julia set corresponding to c shown in the small inset has a neck region * which has a quadratic Julia bifurcation kernel undergoing intermittency disconnection. The apparent bleakness of the landscape is caused by the complete inhibition of period multiplying except on dendritic island black hearts.

Fig 17b: Lower right: The order of the base periodicities of the dark hearts connecting the rocky peninsulae to the coastline of the central black region containing both critical points in the parameter plane of f(z) = z3- z + c fig(15) follows the same mediant-based winding sequence as that of the outer bulbs. But the effects are much harder to recognize from the Julia sets. For example the period 5 dark heart (left) has a base period juila set appearing to be of period 1 but is actually just one petal of a disjoint set of 5 iterating as shown. The same for the period 3 example (above right). Inset shows the base periodicity Julia petal (actual period 3) from this set and one from a period 3 bulb (actual period 9). In the bottom row, the journey of the region a dark hearts and the region b bulbs of fig 15 is shown as the function trends from f(z) = z3+ c, as in fig 14 towards f(z) = z3- z + c. In each case, the left lobe folds on itself to form a region a dark heart and the right lobe become the region b bulb. Other dark hearts are generated from left lobes of satellites such as the two visible in the lower left image are drawn down to the boundary of region a to also form dark hearts. In the lower inset two Julia sets in this journey illustrate how the initial period 3 bifurcates into the 3 base period basins shown upper right.

.

How the two halves of the period 3 bulb of z3+c separate into an outer bulb and a dark heart of z3-z+c,
with Julia journeys showing off mode 1's capacity to model periods, Cantor sets of escaping flows and irrational rotations.

Hence each critical point does give rise to a repeated root, and if this occurs outside the Julia set, a figure 8 crossing occurs and its image, the critical number,  will escape to .  By the same token, if there are no figure 8 crossings and thus no repeated roots, both critical points must be inside the intersection of all nested circles and so lie inside the filled in Julia set and remain finite, validating the dual parameter plane scheme. Then we again have a conformal mapping of the complement of a unit disc onto the complement of the filled in Julia set, so it is connected. If either critical point does escape, it must do so by there being a figure 8 for this critical point and thus the resulting fractal structure must have at least one Cantor component and thus be disconnected as a whole. Since a double point corresponds to a quadratic disconnection, in the non-degenerate case, bifurcations of connectivity should correspond to black hearts, and locally quadratic bifurcations within the Julia set.

Fig 17b: Mandelbrot parameter plane of z5/5+cz4/4-z2/2-cz, illustrating that the Julia set is connected if and only if

the c value remains finite when iterated from all the five critical points, regardless of whether these values form a

connected parameter plane themselves. Top right and bottom left Julia sets are connected, the onthers disconnected.

The top right set appears to come from an island but is still connected, because

the black coloring confirms all critical points remain finite here and do not escape to infinity..

The 'repeated root' theorem is a topological result, so it applies generally to the Julia sets of any polynomial and even to the kernels of Julia webs of transcendental functions we shall look at below. As an illustration in 2011 with a new version of the Dark Heart Viewer, which can picture any polynomial up to degree five, the same situation pertains. This brings us to the 'core' of the Dark Heart, which is the transition of a single critical point from convergence to divergence. The heart, its bulbs and the correspondence to Julia sets applies across the board, no matter what kind of function we are dealing with. We can find it in diverse functions, from the elementary quadratic, through rational to transcendental functions, extending to gamma, with its characteristics remaining even with Riemann's zeta function.

Quicktime movie of the Dark Hearts of f(z)=(z+c)(z-c)(z+1)(z-1) separating as c runs from 0 to 2.

Both parameter planes of z3-z+c have universally black hearts as dendritic islands and have circular bulbs, rather than the hearts of the degenerate cubic. However  also has other more complex boundary structures and the two structures interpenetrate and interfere with one another on the boundary of  in surprising and often confusing ways.

Fig 18: A sensitive region of  iterated from one critical point (a) still fragments according to the dynamics of the other (composite dynamics b), showing each critical point is sensitive to the dynamics of the other when in the non-dominant region and the other critical point is disconnecting.  (1, 4)  period 1 and 15 ÔdendritesÕ of a point in a non-dominant region, lying in a period 1 region and a period 15 bulb of the other critical point.  (2) Black heart Ôcut off at the auriclesÕ by the boundary of a bulb of the other critical point in a non-dominant region.  (3) Black heart region in  with reflective symmetry has a degenerate symmetric intermittency crisis at its base..

The structure of and its Julia sets is relatively straightforward. The boundary of  consists entirely of standard mode-locking bulbs, although these are sometimes interpenetrated by its subset.  The Julia sets in  are homeomorphic to , where C is the Cantor set and  is a connected quadratic Julia set. Those outside  although they possess more complex asymmetries, are all Cantor sets, since all Cantor sets are homeomorphic to one another.  is clearly not connected, as can be seen from its many globally remote Cantor islands fig 15(14) and the next cubic example, which have no dendritic connection to the whole, which (3, 8, 11), and all the dendritic islands do.

The structure of  is much more perplexing, appearing to have features of Julia sets linked by black hearts as well as other interrupted regular and irregular structures. We will also see that  is probably not connected for reasons that will become apparent below. To understand the structure of  we need to also consider the idea of critical point dominance and sensitivity.  Because we now have two critical points, whose dynamics differs and different regions of the complex plane have c values which are dominated by one critical point or the other.  For the above cubic, with fixed critical points on the real axis, each critical point has broad regions of dominance the half-plane containing it, however dominance can vary fractally over local regions and in regions where the boundaries of their two parameter planes overlap fractally the dynamics can become highly irregular. For example dominance can be seen to switch critical points crossing the four diagonal clefts.

Fig 19: Confirming classification in the Mariana trench: (1) Regions investigated – a  (fig 17) and b (fig 19). (2) Parameter planes of the positive (P) and negative (N) critical point. (3) Irregular structure is irregular in both P and N, suggesting mutual sensitivity. (4) Tip of irregular region 3a shows fractal variation of dominance and sensitivity (P bulbs right, N bullbs left). (5) Tiny black heart island verifies classification when Julia sets are sampled (6a,b) shown (7a,b) with key region enlarged below. Classification is also verified in irregular region (6c,d) in (8c,d).  b and d in are both connected and a and c in  are both one-component disconnections as individual black islands are a series of connected components and not Cantor sets.

Critical Point Dominance, Sensitivity and Mutuality [5.10]:

(a) For c in its domain of dominance , a critical point is insensitive to the dynamics of the other critical point(s) , and produces the standard bulbs and dendrites characteristic of the black heart.

(b) For c inside , critical point  becomes sensitive to the dynamics of :

(i) If c is inside the parameter plane of the other critical point, instead of forming period multiplying bulbs with black heart linked dendrites, it forms black heart linked fractals representing the current Julia set  of the other critical point  at this c value. Period multiplying on the main body of each fractal is completely inhibited, but both intermittency crisis and period multiplying occur traversing the linking black hearts. The fractals are mostly of period 1 forms, since  is mostly in  period 1, but can be of high periodicity when a region of   lies in high period bulbs of , as in figs 20(4), 23.  Dendritic regions of the sensitive point may also become an extended region of  as in fig 21.

(ii) If c is outside it adopts the dynamics of the other critical point, its own dynamics becoming effectively fragmented by the otherÕs Cantor set.

(c) For c outside both domains, the dynamics of each critical point may show irregular behavior due to mutual sensitivity, which cannot be resolved computationally, because it is manifest in the behaviour of both critical points. Nevertheless, because of the repeated root theorem, will correctly predict connectivity of the corresponding Julia set.

The level set method colours  black, in shades of green and yellow through to blue according to the escaping critical point, and the rest of the plane is coloured from red by the most rapidly escaping point.

To highlight the notion of critical point sensitivity, fig 18(a) iterates one critical point in a region where its dynamics has a cut-off at the boundary of , where it begins to systematically follow the disconnected dynamics of the other critical point, highlighting the fact that sensitivity is part of the dynamics and not an artiact of the dominance algorithm. Conversely in fig 21 there is a black band where a dendrite of one critical point has been rendered into a stable connected region by the connected periodicity of the other.

Fig 20: Variety of Julia sets on the boundary of  show interaction of two interacting fractal processes corresponding to the two critical points. These no longer have the rotational symmetry of the quadratic and degenerate cubic cases, as there are two interacting fractal processes, the roots are non-symmetric and the critical points are arranged non-symmetrically to the roots. Any Julia set with one factor disconnected (top left) is homemorphic to a product of a quadratic Julia set and a Cantor set, but when both processes are connected, new non-symmetric fractal structures arise. Regularities in the overall shape arise from limitations in the relations between the two periodicities on .

A key point concerning complexity of the higher polynomials is that each step of degree introduces a new fractal process, so that degree n effectively has complexity increasing by . To see this, one notes that the inverse iteration of the quadratic has two square roots, being a non-linear version of a two component iterated function system (Barnsley), with a degree n polynomial being an n component system in which each pair of components introduce a fractal process. However because the process is actually determined in terms of topological connectivity by the n-1 critical points, there are effectively n-1 fractal processes involved.

We now examine a cubic that forms a one-parameter family with the roots able to adopt general positions in relation to one another in the complex plane:   [5.11]. After competing this work on cubics, and independently developing the idea of layered Julia sets, I discovered that Epstein and Yampolski (1996) have investigated the topology of an even more general two-parameter cubic representation, called the Ôcubic interconnectedness locusÕ, dating from John MilnorÕs 1995 Stony Brook lectures. We will use such a two parameter investigation into the rational function family forming the Herman ring system.

Fig 21: The main upper field shows and  for  with disconnected islands of above and below right and inset with arrow top left. Insets top centre are and  showing principal regions of sensitivity and dominance (s,d). Large central inset, interfering black heart of  and bulb of  showing sensitivity at (a) where a dendrite of the critical point becomes a connected region because the other critical point remains finite there. Below (a, b) Julia sets for the corresponding regions above, confirming disconnection occurs according to the classification. Lower right expanded view of the above inset showing quadratic hearts on both  and . Note that region (a) is a fractal dendrite that has become fully connected (in ) because the other critical point is period 1 there.

Dark Heart Viewer: Application and Source includes

This cubic has one root at z=1 and two at  [5.12], while the critical points are at  [5.13].  The relationship between these can be seen in fig 22, with roots (1,0), (0,-1), (-1,1) in blue and critical points  in red for c=(1+i).

Fig 22: Relationship between cubic roots and critical points.

This cubic, which is a one-parameter reduction of a general cubic, has a universality lacking in the previous example, because the roots and critical points can assume any position in the complex plane, rather than the critical points being fixed on the real line symmetrically about the origin.

The two critical points now have parameter planesand with asymmetric structures and domains of dominance that are varying with c and can thus be assessed only on the local dynamics of the critical point.   is demonstrably disconnected between a central region overlapping with  and a series of quadratic hearts globally disconnected from one another, giving rise to Julia sets which are topological products of a Cantor set and a quadratic Julia set, as before.

Fig 23: Critical point sensitivity demonstrated for period 17. Above, period 17 bulb of  confirmed by counting the external dendrites (right), causes the ÔdendritesÕ on  inside the bulb iterating from the other critical point (enlarged below) to determine a family of period 17 connected Julia sets varying in form as c varies (see variations on the arms below).

### 6. The Dark Heart of Division: Rational Functions and NewtonÕs Method

Dark Heart Viewer: Application and Source includes Newton cubic

We now apply the previous cubic to analyzing NewtonÕs method, as an example of iterations of rational functions rather than polynomials. In NewtonÕs method we try to iterate to a zero of a function by repeatedly finding , then proceeding down the derivative

[6.1]

to the x-axis determining .  If all goes according to plan we converge towards a nearby zero by iterating:

[6.2].

Fig 24: The Julia set of

However, knowing about chaos, we can see that points which fail to converge to one of the roots of the cubic may become chaotic and form a Julia set.

For the previous cubic we have [6.3], a rational function, with a cubic numerator and a quadratic denominator. Differentiating, we find itÕs critical points are

[6.4]

The critical points are thus 0 and the zeros of f, which the iteration is seeking.

We can thus use the origin as critical point to seek a parameter plane for .  This will not classify whether or not the Julia sets are connected, as the Julia sets are known to form an infinite connected set, which includes  because it is a repelling fixed point, since , but the Newton parameter plane does successfully classify the local properties of bifurcations in the fractal kernels of the Julia sets, which contain components isomorphic with quadratic Julia sets.

There are two special values of c for which the roots are especially simple. For ,  [6.5] with roots -1, 0, 1 and for ,  [6.6] whose roots are the complex cube roots of unity  [6.7] forming an equilateral triangle about the origin.

Fig 25: (a) Newton Mandelbrot of , (b) blow up of the Mandelbrot kernel in the largest lozenge in (a). (c, d) Julia sets of c = 0, 1 (e-h) Quadratic Julia kernels on the Mandelbrot island in (b).

The dynamic has three super-attracting fixed points at the roots of the cubic equation, so the level set method colors each of these basins in distinct color series. A remarkable topological property of the fractal boundaries of each of these basins in the case  is that they consist entirely of triple points bounding all three basins. Consequently the Julia set is the boundary of any one of the basins.

We will now look at an independently derived rational function which has two parameter generality and which is famous for its unusual Julia sets.

The function  [6.8] forms a two-parameter family. All members of the family have super-attracting fixed points at 0 and , so have two competing basins of attraction separated by a Julia ring which can have a variety of additional fractal features, including quadratic kernels.

XCode explorer of using two complex parameters. Source code

Fig 26: (1,2,3) Layered a-parameter planes for . (Inset) heart islands of (2) showing (1,2,3) stages in Julia connection. (4,5,6) Layered c-parameter planes for with positions on (2) indicated. (7) double ring heart with stages in Julia disconnection ((a,b,c).  (8) Bulb from far right of (6) with stages in Julia disconnection (a,b,c)

Seeking for critical points:

[6.9]

The critical point at z = 0 is inaccessible as it is a super-attracting fixed point with f(z) = 0 for all a and c, however we can form a layered pair of parameter planes as with the previous cubic example using the other two critical points.  We will form an effective 4-D layered Mandelbrot by looking in two stages at the layered Mandelbrots of varying a for and, having picked a complex a, deriving the layered Mandelbrot for varying c of this a and finally the Julia set .

Fig 27: (Top 1-7) Julia sets of layered c-parameter plane a = 4 (centre-right atlas). (3) is a Herman ring with an irrational annular flow in each lobe. The exterior has 180o rotational symmetry identical with  but the inner basins attracted to 0 do not. The Julia sets 3 and 4 display complementary Cantorization, (3) basin of   iinto the basin of 0 and (4) basin of 0 into the basin of . (Bottom 1-4) Julia sets of layered c-parameter plane  (centre-left atlas). (1) is similar fig 26 (7c) but neutrally shaded to show basin boundaries.

This forms a good test of the method because the critical points are varying across complex values adjacent to the boundary between the basins of attraction of 0 and  which is universal to all values of both parameters.  As can be seen from fig 26, both the layered a and c parameter planes classify the Julia sets into connected, partial and complete disconnection, even on fractally-reduced scales. The Julia sets show both quadratic kernels (fig 27 upper 6, 7 and lower 2), and features of the cubic Newton (fig 27 lower 3, 4).

For  [6.10] this function gives a Julia set called a Herman ring with an irrational annular flow in each component, fig 27(3), which appears only with functions of degree 3 or higher.

.

Fig 27b: The example below for n=5, d=3 showing the three types of Julia sets: A Cantor set of points, A Sierpinski curve and a Cantor set of simple closed curves. There are 5-1=4 baby Mandelbrots. Each of the intervening 4 concentric domains has 5+1 outer subdomains and 3-1=2 inner subdomains adding to the 5+3=8 fold rotational symmetry fo the Julia sets.

.

Finally an intriguing example displaying geometric symmetries. Although the function has n+d free critical points rotationally symmetric about the origin, these all have the same set of parameter planes so one suffices. Both the function and its Julia sets are n+d symmetric, while the parameter plane has n-1 baby Mandelbrot islands connected by a web of Sierpinski curves. Each concentric domain has d-1 inner subdomains and n+1 outer adding to the n+d symmetry (http://math.bu.edu/people/bob/papers.html).

### 7. Transcendent Chaos out of Quadratic Disorder

Dark Heart Viewer: Application and Source includes cCos(z)

Fig 28: (Above) Parameter plane of consists of quadratic hearts and generates Julia sets entirely with quadratic kernels.
(Below) Julia sets corresponding to the dendritic island black heart (a) have quadratic kernels identical to the quadratic examples in fig 3.

.

According to Devaney (Peitgen and Richter 1986) a fundamental change occurs with transcendental functions like , which causes their Julia sets to become the closure of the points in the plane which escape to  , which is an essential singularity for transcendental functions, rather than the boundary, as in the quadratic case, where the Julia set is the boundary of the basin of attraction of .

One can appreciate the nature of the essential singularity at  if we consider the corresponding singularity of at 0. It is clear that  oscillates an infinite number of times as it approaches  on the real axis, in the same way.  The process being chaotic on the closure of the basin of attraction of infinity, rather than its boundary, is a consequence of the fact that  any neighbourhood of  is mapped over the entire complex plane minus a point. This allows chaos to reign on a much larger subset of the complex plane, and provides the basis for considering explosions of the Julia set, e.g. of , as illustrated in the sequence of images in fig 29, where the Julia set representing the plane except the internal basin of attraction explodes into this basin as the intermittency crisis bifurcation occurs.

The polynomial-like nature of transcendental functions is realized by their series representation

e.g.

with non-degenerate critical points, and identical extrema, which are locally quadratic in form to second order, as is confirmed by the above power series about 0. Despite the fundamental changes in the formal definition of the Julia set, inspection immediately confirms that the transcendental cosine Julia sets contain an infinite set of fractal kernels isomorphic to quadratic Julia sets, and that this explosion at a single bifurcation point corresponds precisely to an intermittency-crisis bifurcation (compare fig 2) of a period 1 kernel into a Cantor set kernel followed later in the path shown in fig 29 by a reverse period doubling sequence, permeated by an infinite sequence of such processes as each of the fractal islands pass in, along the x-axis, through the origin and up the y-axis, as c goes from down to 1.5.

This is hardly surprising because  form an infinite series of non-degenerate periodic maxima, essentially similar, to first order, to a perfect identical chain of quadratic maxima and minima. However those not familiar with complex numbers should also note that these functions also contain imaginary exponential behavior because  [7.1] so  [7.2] and  [7.3]. For this reason it is better to use bounds on the real or imaginary parts of the number for the level set method rather than a large circle.

The close correspondence between Cos(z)+c and the cubic function discussed in fig 15

Although transcendental Julia sets fill major regions of the complex plane, a momentÕs reflection confirms that the origin is a critical point for the cosine and  for the sine, so it is perfectly straightforward to carry out a computation using the parameter plane in the same way as the Mandelbrot set, and when we do, we find that the region M of non-escape corresponds to a principal central bulb set with a fractal series of black hearts, including a series lining the x-axis, and that any c value in M corresponds precisely to Julia set kernels of the corresponding quadratic type, as illustrated in figs 28, 29. It is easy to see why this is the case, because the critical points in each case lie at the centre of the principal kernels, thus setting up a parallel situation to the standard quadratic.

Fig 29: (Above) Parameter plane of  using a more sensitive level set scheme than fig 28. (Below) Classic journey of the exploding Julia sets of  form a sequence proceeding from  running between points 1 and 2 on the parameter plane above, following firstly an intermittency crisis (the implosion) and later a reverse period doubling sequence, punctuated by an infinite sequence of the same entire sequence as structures move in horizontally and out vertically, meeting at the centre in crisis, as illustrated in the first two figures of the bottom row. The sequence is not in linear time steps, but has Ôslow motionÕ around the explosion point.

Having discovered the close correspondence between trigonometric functions and polynomials, we will now examine how the critical point techniques might function, when we have more complex functions composed of both transcendental functions and polynomial and rational functions, which may have several, or even an infinite number of distinct critical points in their domain.

### 8: Falling into the Heart of Deep Transcendental Variation

Dark Heart Viewer: Application and Source includes

We first examine  [8.1], illustrated in fig 30, chosen because it is a superposition which has calculatable roots which are fractions of .  Examining [8.2], we can see by inspection that  [8.3], however the latter value is a point of inflection (see inset in fig 30) and , which is a fixed point, so cannot be used, leaving only .

Fig 30: (1)  and (2) compared. Despite having two superimposed functions, (2) still has only quadratic hearts and classic quadratic Julia kernels and is essentially identical to that for (1)  except for the two XÕs with suppressed hearts, associated with the inflection point.  (3) Mandelbrot dendritic islands located at (a, b)  above are both quadratic hearts, which give rise to standard quadratic Julia kernels connected to extended Julia networks, as shown below. (4) Julia sets corresponding to c and d in (1) and (2) show kernels surrounding shifted critical points.

The parameter plane for this function is almost identical to that for , retaining the three-lobed central region and a consistent set of dendritic island quadratic hearts, each of which generate quadratic Julia kernels as illustrated in (a) and (b).

To see how additional critical points might complicate the picture, letÕs take the derivative of the previous function, namely  [8.4].

Dark Heart Viewer: Application and Source includes

Solving for critical points, we have (as confirmed in the inset in fig 31):

[8.5]

Fig 31:  and  for the function  correctly classify the connectivity of the kernels of the corresponding Julia sets. (1) An island of  has a disconnected period 3 set (1a). Central bulb correctly classifies connectedness and complete and partial disconnection (2a-d). A degenerate quartic black clover in one critical pointÕs iteration, seen also in  and  left, correctly displays fourfold symmetric quartic kernels and three distinct Julia fractals consistent with degree four, having 4-1=3 fractal interactions. Further quartic islands abound (4).

We find that if we examine the critical point parameter planes  and  are identical, however  gives a completely different profile. Combining  and  to form  and , as we did with the non-degenerate cubics in section 4, immediately provides an acid test for the connectivity of the contributing fractal processes, which number 2 and sometimes 3.

The combined analysis covers the behavior of all the critical point types, and  and  give a correct classification of the connectivity of the fractal processes contributing to the kernels.  A notable and unexpected appearance in  is the degenerate quartic clover leaf, fig 31(3) (compare fig 14) which consistently has 2-cleft cubic bulbs. This is partly obscured in  and , because it overlaps both, but is still evident on the superimposed figure on the left.  Consistent with this subset, we find the kernels of the associated Julia sets of the three period doubling bulbs (a, b, c) all contain a component with square quartic symmetry and 3 fractal components, consistent with local quartic behavior, which is also clearly reflected in the form of the inset real graph of f.

The connectivity of the contributing fractal processes is also correctly described for the central period doubling bulb fig 28(2), where (2a, 2c) in  have connected kernels, (2b) in  has one fractal component disconnected and (2d), outside both, has a Cantor kernel. Likewise the period 3 bulb of the quadratic heart island (1, 1a) has product of a Cantor set and a quadratic period 3 Julia set.

Fig 32: Parameter plane of , with (a) principal region of and regions of coloured by which critical number they are generated by, (b) series of embedded Julia sets in the centre region of  (c) A corresponding series in the right heart undergoing period doubling. Individual Jula sets corresponding to locations in (a, c).

Dark Heart Viewer: Application and Source includes

The function  is of interest since, although  is double valued, because , the composite is unique valued. The Julia sets all posses a parabolic symmetry, reflecting the fact that the target set  is equivalent to (Durkin, Brown and Halstead).

Although the critical points  form an infinite sequence, the critical values of the function are , so we can use the two critical values to form a pair of parameter planes forming  as before.

In both this and the previous case explosions can be considered, in which one or both critical values are drawn into a finite periodic attractor, including explosions into the fundamental semi-infinite quadratic region, as illustrated in detail in fig 33.

Fig 33: Intermittency crisis explosions compared under identical conditions at 16384 iterations. Top row main cardioid of . Second row period 3 denritic Mandelbrot island on the period 2 dendrite of . Third row Julia set explosion of  as at the beginning of the series in fig 29. Fourth row explosion of Julia set of

exiting the central black heart in fig 32.

The explosion of  is particularly sensitive to the iteration method because it is a composition of a function of degree ½ with a transcendental function of first order degree 2. Mathematically identical formulations of  give radically different portraits of the explosion and different parameter values for the bifurcation point.  Comparison of the bottom row of fig 33, which depends only on the square root function with fig 41 which uses arctan to calculate the angle halving of shows how different these can be as a result of numerical overflow or underflow effects in the iterative computation, with the latter showing discrete effects leading to increasing complexity approaching the apparent bifurcation point. For a full discussion, see appendix 1.

Fig 34: Layered parameter plane of using the first 8 critical points derived using Newton's method and regions of corresponding Julia sets. (a) A region of high sensitivity to successive critical points shows Julia sets with first disconnection (1) in the 8th sector, (2) in the 6th and 7th, (3) in the 4th (see enlargement) and (4) a set outside , with multiple disconnections caused by the escape of the first critical point, llustrating the evolution of the kernels, each of which corresponds to a particular critical point. (b) An irregular bulb eroded by the third critical point with connected Julia sets (1, 2) and one disconnected at the 3rd sector).

We now extend our discussion to an example where there are an infinite number of distinct critical points, with distinct critical values, each of which could contribute to the fractal dynamics, by examining the rational transcendental function  [8.6] illustrated in figs 34 and 35. Although this is undefined at 0 the singularity is removable by substituting the value 1 at 0 since .

Examining the function for critical points, we find (as confirmed in the real function inset of fig 32):

[8.7]

Following the technique of the previous function in which a dominant maximum and a neighbouring minimum were used successfully, we thus look at layered parameter planes, using NewtonÕs method as above to evaluate the first 8 positive critical points of the function and solve for  and as before.

To highlight the effects of the multiple critical points, in fig 34 we have coloured the parameter plane black if we are in  and shade c value for which the i-th critical point escapes but the previous ones do not, in a different colour. The complement of these points are coloured by the maximum dwell so that the individual  are also highlighted, as well as regions where all critical point escape.

This highlights sensitive regions where successive critical points are continuing to erode the boundary of  and enables one to investigate where disconnections are occurring in the Julia kernels. In this case the functionÕs undulating critical points dominated by the central maximum result in a Julia kernel which near zero contains a finite basin extending and expanding in both real directions. As c moves away from zero the dominant critical point escapes resulting in disconnection of this basin. Further disconnections can happen at any of the subsequent sectors of the Julia kernel each of which corresponds to the behaviour of one of the successive critical points. Fig 34 highlights this correspondence by comparing changes in the parameter planes of successive critical points with changes in the corresponding sectors of the associated Julia set.

Investigation of combined bifurcations of several critical points of

Fig 35: Parameter plane of  with confirmation of quadratic period 3 kernels (a,b)
and a disconnected kernel in the complement (c).

As a demonstration of universality, we have  [8.8], which has a single critical point:

[8.9].

As can be seen from fig 35, this has a global parameter plane with a dendritic island black heart, both of which give rise to connected quadratic period 3 kernels, while a complementary point gives a disconnected kernel.

Further Transcendental Functions and the Gamma Function Explored

f(z)=czlog(z), f(z)=exp(-z2)+c, f(z)=czzz =czexp(zlog(z)), f(z)=Γ(z)

Universality of the dark heart examplified (left) in the parameter plane of f(z)=czz+10 and (right) the close relationship

between the gamma function f(z)=cΓ(z) (above) and f(z)=czz-3/2 (below) reflected in formulae such as Stirling's.

Video of journeys around f(z)=czz+k. Three journeys are made. The first along the real axis from -5 to 5.
Then we make the same journey except circuiting around 0 along the lower unit semicircle.
Finally we move from -1 around the circle to -i, up the y-axis to i and back to -1.

### Mac XCode High performance RZ Viewer and Source Code Gives excellent high performance rendition of the RZ function and its Mandelbrot and Julia sets and those of related functions.

Quicktime movie of wave function method on RZ Julia highlighting the zeros in frame 1.
Quicktime movie of Gaussian wave function highlighting the non-trivial zeros in frame 3.

As an exercise in using these techniques at the frontier of functional complexity, and as a basis to generate open source code, I developed Matlab code to perform rapid calculation of the Riemann Zeta function[9.1]:

[9.1]

on the entire complex plane, except the single pole at z = 1, using the Dirichlet Eta function , and the Lanczos approximation to the gamma function to extend the domain to :

[9.2]

The Lanczos approximation is defined as follows:

[9.3]

with  and  calculated independently by Paul Godfrey as constants from the relation:

where   are coefficients of the Chebyshev polynomial matrix.

A satisfactory formula I have used subsequently is [9.3b] trunctated to the same number of products as the sum in [9.2].

Fig 36: (a) The Riemann Zeta function , showing its pole at 1 the trivial zeros at even negative real values and the non-trivial zeros on the line x=1/2. (b) The Julia set of , highlighting eventually fixed points in the internal basin mapping to , to which the zeros are also mapped. Inset in right overlap of eventually fixed point and non-trivial zero showing their proximity, with the eventually fixed point to  lying on the curve where .

The iteration takes the first 2001 terms of  and halves the next term to average contributions near 0 where the terms remain large. The overall process is computationally extremely intensive and the self-organized criticality of , launching the positive half-plane to the neighbourhood of the pole at 1, is liable to stress numerical calculation to the limit. However the results are intriguing and provide support for the dark heart of the parameter plane even under these extreme circumstances.

As illustrated in fig 36, the Julia set of (Woon) forms the boundary between the basin of attraction of  and the attracting fixed point . The first six non-trivial zeros of the function, from  on, famous for the Riemann hypothesis - that they are on all on the line x=1/2 - lie in the basin of attraction of this fixed point. In fact all the zeros of the function do, including the trivial ones at , as all are mapped to 0, which iterates to .

The overwhelming reason for the zeros to be on the line x = ½ is that there is an internal symmetry in the analytic extension of the function to , which was expressed by Riemann  (1859), in the form of the reflection relation around x = ½:

[9.4].

Once expressed by Riemann in the form [9.5], this internal symmetry is made evident, as shown right, comparing .  It is also evident in the reflective dynamics of the iteration around the Julia set of  in fig 37.

Fig 37: Above: The functions  have the same set of non-trivial zeros but the latter also displays the internal symmetry about the roots. Below: The Julia set is not connected because it contains an infinite fractal series of islands in the V's of the negative real half-plane which are 'adiabatically' disconnected from the main body..

The non-trivial zeros are definitely confined to the critical strip x = [0,1], as x cannot be greater than 1 by the Euler product formula [9.1], whose individual terms immediately result in prime sieving of the series formula, when taken to the left hand side, and x cannot either be less than 0 as a result of the reflection formula (except for the tirivial zeros which come from the sin term). All zeros so far computed lie on the line, an infinite number and around 2/5 of the total have been proven to do so but there the defence rests. The number of primes with imaginary part between 0 and T as a result of the reflection formula and Cauchy's argument principle integrating round zeros and poles, depend on [9.6]. From [9.5] it can be seen that this in turn depends on [9.7], which behaves as a Gaussian random variable growing very slowly according to log(log(T))1/2, whose deviations may not have fully surfaced, even at current large computed numbers. This has led some skeptics to suggesting the hypothesis may actually be false.

Fig 37b: Julia set of the Collatz problem as a compex iteration

Given the relative difficulty of a final proof, the possibility remains that the Riemann hypothesis may be an undecidable proposition, because confirmation ultimately depends of computing the position of all zeros, and thus becomes an effective Turing halting problem, like the unresolved Collatz conjecture - whether the sequence starting from an integer n and either halving it if it is even, or multiplying it by 3 and adding 1 if it is odd, converges (to the cycle 4 > 2 > 1). Like the Riemann hypothesis this iteration contains occasional path records which significantly exceed the overall growth of the orbit lengths. Although all numbers appear to converge to the period 3 cycle, some take an exceptional time, deviating from the generally inverse quadratic relationship between number and orbit length. The undelying complexity in this process can also be realized by extending it to the complex numbers and results in a chaotic discrete iteration corresponding to the associated transcendental function shown in the inset, having a fractal Julia set surrounding the basins of attraction of the integers, with other points possessing diverse orbits, including divergent ones coloured as shown at left.

Although the Julia set appears to be connected, this is unlikely, because all the trivial zeros on the negative real axis also iterate to the attracting fixed point and the dappled lacunae in the divergent basins are in a state of self-organized criticality, where divergent values and values convergent to  appear to densely interpenetrate, as emphasized in fig 4b where diverging points are black and those converging to  are coloured. This is likely to be a result of computational overflow at the pole and indeed testing values with MatlabÕs own standard zeta function (presumably derived from the Maple engine as it is in the symbolic toolbox) froze the computation after three iterations .

Fig 37c: (a) Parameter plane of  iterating from the critical value 1 [-40,6] x [-40,40]. (b) Enlarged region [-20.5,-15.5] showing complex bifurcations. (c,d) A tiny fractal ÔlakeÕ barely visible around c = 0, shown expanded [-0.1,1] x [-1.1,1.1] (c) and larger (d) . This corresponds to c for which the positive half-plane diverges due to the pole at 1.Inset in (c) at half the scale is the corresponding ÔlakeÕ iterating from 0 corresponding to c for which the zeros of  diverge. (e) Sample Julia sets on the real line, with colour chosen to distinguish divergence, convergence to the critical strip and convergence into positive half plane - shades of blue for points tending to , red for points asymptotic to a periodicity in the critical range, green asymptotic to a periodicity with real parts entirely positive, and grey to black for non-periodic points (largely absent). The range from -15 to -0.1 is similar to 0.6 with no apparent bifurcations. (f) An enlargement of the region shown # in (a). The sensitive changes in the Julia set from 0 to 0.0125 correspond to a transition off right of the island in (d).

To test the question of parameter planes as a measure of the bifurcations of a family of Julia sets, we now examine the parameter planes of the function . The most outstanding critical value of  is the value 1, to which all z with positive real values tend as . This however is a case of self-organized criticality, as it projects the entire positive half plane directly on to the neighbourhood of the pole at 1. The Julia set  is thus extremely sensitive to small changes in c, undergoing explosions of the positive half-plane for arbitrarily small changes in c.

In fig 37c we show the parameter plane for the critical value 1, in which the black regions indicate the critical value remains finite and colours indicate divergent iteration to , along with corresponding illustrations of Julia sets, which highlight both complex bifurcations on the real line in the interval [-21,-15], and explosions of the positive half plane on either side of 0, at ~-0.005 and ~0.001. Both these features correspond closely to the parameter plane, which shows both complex fractal structure in the former range and a tiny fractal ÔlakeÕ around 0 connected by a dendritic thread of further ÔlakesÕ winding in the imaginary direction to the divergent region.

In fig 37c2 we demonstrate that the same mediant winding sequences appear on the bulbs of the black region as in the standard Mandelbrot set of polynomial and transcendentals functions. The right hand plane coloured black in fig 37c(a) is iterating to the fixed point 1 + z because = 1, so + z = 1 + z. But this is true for all z with large positive real part, so the iteration is fixed on 1 + z.ue We can thus colour the black region according to how many steps it takes to reach within epsilon of a fixed point or periodicity and colour by the number of steps in blue and add redness for the period. This immiedately shows up the periodicities of the bulbs neighbouring the boundary which can be confirmed to correspond to Julia sets with rotational periodicity the same number, confirming the sequences of periods of the bulbs and the mediant relationship in the fractal progression.

Fig 37c2: Shading the bulbs demonstrates their periodicity, as confirmed by the Julia set portraits, which display the corresponding rotational periodicities in their spirals , confirming an upward set of odd periodicities 3, 5, 7, 9 ... and a downward set of integer periodicities 3, 4, 5, 6 ..., each with mediant fractality viz (3,4)=7, (4,5)=9, (3,5)=8, (5,7)=12.

A second parameter plane iterating from the critical saddle at ~ -15.34 shows how the Julia properties of the zeta function can follow the same rules as our original quadratic and finally demonstrate the universality of the dark heart and its bulb structures as an atlas of Julia dynamics.

Fig 37d: Examination of the region [-16,-15]x[-1,-2] of  , marked * in fig 37c(b), The Mandelbrot set of the critical saddle at ~ -15.34 with critical value 0.5206 displays a classic bulb. Selecting a point in the period 3 subregion generates a beautiful example of period 3 Douady rabbits in the corresponding Julia set, confirming the relationship between the two and the capacity of the zeta function to locally model any function.

We can likewise investigate critical points lying on the critical strip between the non-trivial zeros. These generally show less interesting features because their critical values tend to approach 1 so give results similar to those aove for the critical plateau around 1000, however locally they can give more accurate atlas descriptions of the Julia sets of neighboring points.

Fig 37e: Left Mandelbrot set of using the critical point z = 1.3828 + 42.2909i
provides a period 3 Julia set at the location -0.375467 + 42.981618i starred in the period 3 sub bulb.

Finally the jewel in the crown! When we go to z = -17.3739, the first critical point on the x-axis with negative critical value -3.7436 we find, sitting in the inner bay of the central area, a perfect quadratic black heart, almost exactly as our original quadratic heart of the Mandelbrot set of f(z)=z^2+c was at the beginning of the paper. These hearts are perfectly replicated in the fractal bays of all the dendrites up and down the critical strip.

Fig 37f: Full circle. The quadratic heart at the centre of the zeta cyclone of  at the site

marked by the black star for the critical point at -17.3739, with a Douady rabbit Julia kernel

chosen by selecting a point in the period 3 bulb.

Another dynamical view of the zeta function is provided by taking the function as the values of a vector field, the so-called 'holomorphic flow' Broughan and Barnett (2003), in which successive non-trivial zeros have a complex encoding in terms of sources, sinks and rotational flow.

Fig 37g: Ghostly ÔBuddhabrotÕ image plots the collective orbits of the critical point of by adding 1 to each pixel that the orbit from the critical point of each c value on the plane visits. The orbits of c points outside M are iterated in cyan overlying the orbits of c points in M iterated in shades of ochre and magenta Matlab m-file.

### Conclusion:

The approach gives a complete investigative approach to confirming the universality of the Mandelbrot "dark heart of chaos" across the very broad class of complex functions possessing polynomial-type critical points. The paper has demonstrated that multiple critical point analysis, as developed to deal with cubics and higher polynomials extends readily to composite transcendental functions and thus embraces the comprehensive majority of complex functions central to mathematics. The ultimate conclusion is that chaos, represented by the Julia sets, or their bifurcation kernels will always trap the critical points of the function as the last super-stable points of escape into basins of attraction if these kernels are locally homeomorphic to connected sets.

### Appendix 1: Combined Methods of Depicting Julia Sets and Parameter Planes

Fig 38: Quadratic Julia sets depicted by combined methods. (a) A Siegel disc (). Exterior (blue) by level set, Julia (white) by modified inverse iteration (XCode), irrational flow (green) by the sine of the velocity (Matlab m-file). (b) Super-attracting period 3 (Solution of ). Exterior by continuous potential, Julia by modified inverse iteration, interior basins by level set. (c) Seahorse valley dendrite () by distance estimator, exterior by level set. (d) Parabolic bifurcation period1 to 20 ( adjacent to ). Exterior by level set, Julia by modified inverse iteration (XCode), petals by velocity (Matlab m-file).

### Level Set Method LSM

Begin with the point for a fixed c and iterate . If we are depicting the Mandelbrot parameter plane, we begin instead with the critical point and iterate for each . If the point escapes a circle say of radius 10, , we colour it by the number of iterations.  If it remains bounded after a fixed number of iterations, we assume it cannot escape and colour it black (or white in the above examples).

This method will work for the Julia set of any function provided we can determine basins of attraction of fixed or periodic points, to apply the algorithm to. In the case of  this is  and the Julia set is the boundary of the basin of attraction of . In other functions it may be finite and multiple. In the case of transcendental functions, such as cos and exp,  is an essential singularity and the Julia set is the closure of its basin of attraction and the target set may be  or . The method will also work for the parameter plane (Mandelbrot set) of any function provided we can locate and handle the critical points and establish individual parameter planes for each.

Variation 1 Internal Basins:

If the Julia set has an internal attracting periodic point we can test for this by also finding the periodic cycle by first iterating the critical point until it becomes within  of being periodic and then colour a non-escaping point by the number of iterations to bring  within  of a point on the cycle.

XCode Viewer and Source Code

Variation 2 Continuous Potential:

We can derive a continuous potential for quadratics with an even gradation but less sensitivity at the boundary of

by using the continuous potential formula  for an escaping point.

Fig 39: Herman ring  Irrational annular flow depicted

by sine of the discrete velocity, interior basin level sets red, exterior blue. Matlab m-file.

Variation 3 Distance Estimator (Quadratics only):

For dealing with dendritic Julia sets, or the Mandelbrot set, to demonstrate connectivity, we can make two additional tests to detect proximity to :

(a) Test for overflow of the derivative of any point on the escaping orbit:

(b) Test for distance

We can calculate orbit derivatives iteratively as follows  for the Julia set and  for the Mandelbrot set, since we are differentiating with respect to c rather than z. In practice, O ~ 1600000,  ~ 0.1. The derivative overflow tests for highly repulsive dynamics adjacent to .

XCode Viewer for Distance Estimator the contributions of overflow and distance to complement the non-escaping points (black) are shown in green and blue.

Variation 4 Discrete Velocity of non-attracting Basins and Petals:

Compute, for the points that donÕt escape, the average discrete velocity  on the orbit.

Matlab files Herman Siegel

Variation 5 Binary Decomposition: Split the level sets light or dark depending on whether the n-th iterate is in the upper or lower half plane - i.e. has imaginary part positive or negative. This highlights both equipotential level curves and sections of external rays. See fig 5.

Fig 39b (Clockwise from top left) Pickover stalks for the Mandelbrot set, enlarged, a Julia set from the small bulb and detail of internal structure. The two right images are period 7 Julia sets generated by using an anular trap intersecting the period 7 Mandelbrot bulb and coloring using a mod 7 scheme shaded by lof of the distance from the edge to highlight an infinite series of trap regions modulo the period of the Julia set. Jos Leys has some good art work using these methods, but remember they are not highlighting the Julia set, but rather a fractal collection of orbit trap images lying inside, outside and across the Julia set.

Variation 6: Orbit Traps and Pickover Stalks: Choose a target neighbourhood, e.g. narrow 'stalks' along the x and y axes, or circles or discs about the origin as orbit traps. Iterate a point using the level set method, but if it enters the target region, colour it with the minimum number of steps to enter the target instead, shaded by the log of the distance from the edge.

Fractal Domains Explorer

### Inverse Iteration Method IIM

First find a repelling fixed point , , by solving . For , one of the two fixed points is always a repeller as , unless . Now plot the two inverse images  of this point and repeat to form the   n-th inverse iterates. This method requires a heap, or some equivalent data structure, to keep track of the branching tree of inverse iterates. If memory is exceeded we can randomly plot one or other roots and its pre-images. However this method has the problem that it is computationally intractable because the points are exponentially unevenly distributed over the Julia set, due to multifractality (fractal redistribution of the probabilities), resulting in the inverse mapping being strongly contractive to some features leaving others unrepresented, and thus fails to represent significant features of , even with exhaustive computation times.

Fig 39c MIIM coloured by the log of the absolute derivative (click to enlarge)

Variation 1 Modified Inverse Iteration MIIM (see figs 3, 13b, 39c):

We cut off the sub tree from a given  if the derivative .

This eliminates dominant highly contractive regions of the inverse iteration, which have already been registered. We can calculate successive derivatives iteratively .

XCode viewer and source code for heap processing MIIM using pointers (thanks to Yvo Smellenbergh for help transporting the code to OS X)

Fig 40: Wave function method developed during the production of this work portrays the dynamics leading to formation of Julia sets by iterating successive compositions of f(z) with the complex wave function illustrated by colouring the plane shades of red (r) and green (g) according to a conformal function, such as orwhich highlights recursive pre-images of the critical point.

XCode Wave Function Method Viewer: Application - Source - Flight Manual

Conformal Wave Function Method

This method developed during production of this paper to gives a dynamic picture of the action of iteration of complex functions having two input and two output variables. A suitable function (see fig 40) is used to generate conformal interfering waves of oscillating colour, forming a Cartesian or polar grid.  Successive iterates of f(z) are then applied to form a sequence of images , with numbers which have exceeded computational bounds assigned black by default (Matlab does this automatically). The method generates a discrete dynamical movie of the recursive action of f on the plane, showing how the repeated application of f leads to formation of a Julia set, highlighting the internal dynamics of the basins, including attracting and neutral mode-locked periodic points, irrational flows and Cantor processes. The process is a form of inverse iteration since the n-th frame is drawing the n-th pre-image of the domain of the wave function. If polar coordinates about the critical point are used, this highlights pre-images of the critical point.

Fig 40: Wave function method highlights the dynamic periodicities of all regions of the parameter plane. The first frames immediately highlight periodicities of the critical point as images of the origin, showing both locations and periodicities of period 2 to 7 super-attracting bulbs and dendritic islands. Later iterations highlight the periodicities of the period 2, 3, 4, and 5 bulbs starting with frame 60 having all these bulbs showing images of the origin.

Quicktime movie of Mandelbrot periodicities by wave function method

One can also produce a dynamic movie of the internal periodicities on the Mandelbrot set by repeatedly inserting the original c value into each iterate of f, expressed in terms of the critical value. For  this is just c, but for , with critical number ½ and critical value v=c/4, we need to insert 4v, and use a wave function centered on the v value which will fix the critical point, i.e. v=1/2, or c=2. This method immediately highlights the positions of all the period-n bulbs and dendritic islands, as well as demonstrating the periodicity of each.

Chaotic Processes and Discrete Iteration

How discrete processes such as underflow and Moire patterns can add accidential complexity to investigating discrete chaotic processes.

## Appendix 2: Ray Tracing Hypercomplex and Multi-dimensional Chaotic Iterations

We now investigate how Julia and Mandelbrot sets of the 4D Quaternions, bicomplex numbers and other systems, such as the spherical polar Mandelbulb iteration can be investigated in 3D space. Many thanks to both Tom Beddard and for cooperative help and great code which made the examples in this section possible.

A: Hypercomplex Systems

The quaternions Q are defined as:

Q =

Since i, j, k share a symmetrical relationship, a pure unit quaternion u = (0, b, c, d) behaves as a rotation in 3D space , with each of i, j, k corresponding to rotations of 180 about the axes.

This symmetrical representation can be extended to 1+n dimensions (Dang et. al.). One needs to note however that the only true division rings over the reals where elements can be both multiplied and divided and give rise to a full suite of rational and transcendental functions are the complex numbers and quaternions , with a non-associative extension to the octonions O, by effective quaternion 'complexification', according to the Cayley-Dickson construction , with conjugate , in the same way the quaternions were formed from the complex numbers.

Note that in the quaternions, De Moivre's theorem still holds when we express a given element in terms of a pure quaternion unit vector , but the loss of commutativity means that, for example , so there is no unique way to write a polynomial in terms of powers of z, however, like , we can develop transcendental functions such as the log and exp.

Fig 42b: (Top left) Quarternionic Mandelbrot set in 3D is a solid of revolution of the complex 2D Mandelbrot set about the real axis. (Above) Cross-section of this quaternionic Julia set is a 2-D Douady 'rabbit'. (Right): Quaternionic Julia set [source file] using the free cross-platform package MegaPOV projected into . (Lower left) Symmetric projection of the same Julia set [source file] generated much more quickly using the free cross-platform developer package Adobe Pixel Bender. A given Julia set of can be represented in 3D because it can be rotated by an angle in parameter space to a related Julia set . The coordinates can be used to define a Julia set of an arbitrary via the rotated mapping Movie file of Quaternion Julia sets (mp4)

If we are considering a quadratic mapping in Q, we have:

Hence, given the iteration , we have the following four-assignments:

Expressed in terms of vectors, for a level set, or escape-time algorithm, this becomes:

.

Now an apparent paradox arises, in which the Mandelbrot set appears 'simpler' than the Julia set, in the sense that it can be fully represented in a lower dimensional space.

If we are considering the Mandelbrot iteration then r and c are co-directional. If c is written as then and the iteration preserves the complex plane containing the real axis and the axis of u. The quaternionic Mandelbrot set thus consists of a complex 2D Mandelbrot set rotated about the real axis to form a spherically symmetrical set in Q. This can be represented in as a rotation of the 2D complex Mandelbrot set about the real axis, as shown in fig 1(a).

The same situation pertains for a Julia set .

Now let us consider an iteration of the form . It turns out that the dynamics of this iteration are independent of the angle in .

Let then .

Now consider . If we choose so that , then

.

If f(z) is as above then , so the rotated dynamic is the same as the original.

Now let us consider . Given q there exists . The Julia set is simply a rotation of the dynamic in to .

The effect of the function , is to simply rotate the Julia set of the iteration unchanged within , but resulting in very different Julia sets in the expansion to the rest of . It is thus common to investigate the triple and iterate the above function in , projecting the resulting points to to give 3D projections of a representative spectrum of quaternionic Julia sets. These can either be portrayed by inverse iteration, or by using a distance estimator and 3D ray tracing (see below).

Quaternionic versions of transcendental and other functions can be defined and portrayed, just as they have been in this paper for complex functions, using the above methods, according to the following definitions (Halayka 2009):

Octonionic Julia Sets (Griffin and Joshi 1992) and a generalized Mandelbrot set (Griffin and Joshi 1993) can also be defined and portrayed, and display new features associated with their non-associativity, in which (ij)l=-i(jl)≠i(jl).

B: Ray-tracing Principles

Ray-tracing techniques reflect right rays from a source off a tangent plane determined by a normal from a point estimated iteratively to be arbitrarily close to the Julia set using an iterative process similar to that used to generate the sets themselves (Hart et. al.). Methods need to be sensitive to the fact that fractal sets often have non-differentiable boundaries. For more detailed information see (Dang et. al.).

The standard distance estimator used for a complex or quaternionic quadratic iteration is: applying the chain rule to . This is unchanged for , since .

Fig 42c: A succession of unbounding volumes iteratively determined using the distance function.

The algorithm iteratively checks the distance to the set using successively smaller unbounding volumes not intersecting the region, to step along the ray closer, but not into the region until the distance to it falls below a given threshold. A calculation of the tangent plane is then made, either by differentiating to get a normal, or by approximating a tangent plane using neighbouring points 'on' the surface. The shading of the point projected onto the screen plane is then coloured according to the lighting and the orientation of the tangent plane.

To calculate a normal, a variety of methods can be used. A cross product of vectors to neighbouring points, finding the direction of the point z a fixed distance from on the set whose distance estimation is maximal, finding the direction of maximal attraction, or using a gradient, for example calculated by picking six neighbouring points and the same for y and z, or by directly differentiating the iterative function defining the region. Two good examples of complete platform independent source code files illustrating the technique fully are the Quaternion Julia code example and the 3D Mandelbulb code example illustrated above and below.

Fig 42d: A Julia Island [source code] synthesized using the free cross-platform ray-tracing package MegaPOV.

The discrimination radius used to determine how closely neighbouring points on the screen are computed can be made a power law function of distance from the observation point or 'camera', , so that more distant parts of the image are computed in similar relative precision to nearby parts - a process called 'clarity'.

C: Bicomplex Numbers

Fig 42e: Level sets of the Mandelbrot set (left) correspond (centre, right gray) to product factors of the tetrabrot, for various values of .

Bicomplex numbers (Rochon) are defined a little differently from the quaternions:

T =

Like the quaternions, bicomplex numbers can be expressed in terms of a pair of complex numbers:

One can define a generalized Mandelbrot set in T in the same way as the complex numbers:

The so-called tetrabrot is then defined as the cross-section of this in

This set can then be portrayed by colouring the external level sets on the surface of the tetrabrot level set at distance surrounding the actual fractal. As is reduced, the fractal details of the 3D structure emerge, as shown in fig 3.

T is a commutative unitary ring, but it is not a division ring like the quaternions, so it is not capable of developing a full suite of functions.

D: The Mandelbulb and Spherical Polar Iterations

Fig 42f: Three views [source code] of the degree 8 Mandelbulb using the free package Adobe Pixel Bender (Tom Beddard). The author also has an update for MegaPOV for Mac from Yvo Smellenbergh to enable the Mandelbulb to be generated also using MegaPOV. The patch to recompile MegaPOV or POV-Ray is available here.

A new development in 2009 driving the compilation of this appendix has been the extending the basic rules of the quadratic Mandelbrot set in 2D polar coordinates to 3D spherical polars in the Mandelbulb.

The iteration in polar coordinates takes the form:

This can be readily generalized to the 3D mapping :

Ray tracing the Mandelbulb in the above figures requires using the distance formula to iteratively differentiate the above function: to produce a Jacobian matrix capable of determining the tangent plane.

Although spherical polar coordinates do not form a number system, and indeed there is no division ring over the reals except for the complex numbers and quaternions, the above iteration at least is well-defined for a variety of integer powers, and the 3D set forms a 3D fractal set similar to the 2D complex Mandelbrot set, particularly for larger n in the range of 6 and above.

Just as the complex iteration has odd symmetry for even powers and even symmetry for odd powers, corresponding to (n-1)-fold rotational symmetry, as noted in [4.2] and fig 14, the Mandelbulb displays the same symmetry types in both the and dimensions.

Fig 42g: (Left) Four degree 8 Julia sets of the Mandelbulb iteration. (Right) Mandelbulbs of degree -8, 4, 7 and 16. Even degrees have odd symmetry and odd degrees even (compare the degree 7 and 8 side views in figs 42g and 42f).

### Software and Demonstrations:

The software used to investigate this was an intel native Mac OSX XCode application-generating project originally developed as a Mandelbrot explorer for the standard quadratic by Michael C. Thornburgh and converted to a generalized Mandelbrot/Julia function explorer by the author. A Metrowerks C program suite developed by the author was used to generate movies and Julia sequences at specific loci. The software, source code and movie demonstrations is available from a link at http://www.dhushara.com.

All of the viewers work by dragging a rectangle to kexpand a sub-region. If you click you will switch between the parameter plane and the Julia set of that coordinate. Sometimes checking the random box will fill the interior basins. The Herman viewer toggles from a-Mandelbrot to a c-Mandelbrot with the chosen a and then to the ac-Julia set. Reset takes you back to the beginning. The number of threads accelerates the program at the expense of other threaded programs. Increasing the number of iterations gives a slower, but more accurate image.

1. Compiled Mac OS X viewers of each of the exaples (see article topics).

2.

3. Pixel Bender Source Code Quaternion Julia,

4. MegaPOV/POV-Ray Code , Mandelbulb patch (requires recompilation)

5.

### Bibliography:

1.     Barnsley M. (1988) Fractals Everywhere Academic Press, New York.

2.     Brown D, Halstead M (2007) Super-attracting cycles for the cosine-root family Chaos, Solitons and Fractals 31 1191-1202.

3.     Dang Y, Kauffman L, Sandin D Hypercomplex Iterations, Distance Estimation and Higher Dimensional Fractals http://www.evl.uic.edu/hypercomplex/html/book/book.pdf

4.     Devaney R.L. (1986) An Introduction to Chaotic Dynamical Systems Benjamin/Cummings, Menlo Park.

5.     Douady A, Hubbard J. (1985) On the dynamics of polynomial-like mappings Ann Sci Ecole Norm Sup 18, 287-343.

6.     Durkin M. (1998) Observations on the dynamics of the complex cosine-root family J. Differenc Equat Appl 4 215-28.

7.     Epstein A., Yampolsky, M. (1996) Geometry of the Cubic Connectedness Locus I: Intertwining Surgery arXivMath/9608213v1

8.     Griffin C, Joshi G (1992) Octonionic Julia sets Chaos, Solitons & Fractals 2/1 11-24

9.     Griffin C, Joshi G (1993) Transition Points in Octonionic Julia Sets Chaos, Solitons & Fractals 3/1 67-88

10.    Halayka S (2009) Some visually interesting non-standard quaternion fractal sets Chaos, Solitons and Fractals 41 2842â€“2846

11.    Hart J C, Sandin D, Kauffman L (1989) Ray Tracing Deterministic 3-D Fractals Computer Graphics, 23/ 3 289-296

12.    Milnor, John W. (1999). Dynamics in one complex variable. Vieweg, Wiesbaden, Germany. ISBN 3-528-13130-6.

13.   Pastor G,  Romera M, Alvarez G, Montoya F (2002) Operating with external arguments in the Mandelbrot set antenna Physica D 171  52–71

14.   Peitgen H.O. & Richter P.H., (1986), The Beauty of Fractals Springer-Verlag, Berlin. DC

15.   Peitgen H.O. et.al. (1988) The Science of Fractal Images New York ; Berlin : Springer-Verlag

16.   Peitgen H.O., Jurgens H., Saupe D. (2004) Chaos and Fractals: New Frontiers of Science New York: Springer.

17.   Riemann B (1859) On the Number of Prime Numbers less than a Given Quantity http://www.maths.tcd.ie/pub/HistMath/People/Riemann/Zeta/EZeta.pdf

18.   Rochon, D. (2001) Dynamique bicomplexe et theoreme de Bloch pour fonctions hyperholomorphes, Doctoral Thesis, Universite de Montreal.

19.   Romera M, Pastor G, Alvarez G, Montoya F (2004) External arguments of Douady cauliflowers in the Mandelbrot set Computers & Graphics 28 437-449

20.   Romera M, Pastor G, Alvarez G, Montoya F (2006) External arguments in the multiple-spiral medallions of the Mandelbrot set Computers & Graphics 30 460-469

21.   Schršder M. (1993) Fractals, Chaos and Power Laws ISBN 0-7167-2136-8.

22.   Shishikura M. (1994) The Boundary of the Mandelbrot Set has Hausdorff Dimension Two AstŽrisque, 222/7 389-405.

23.   Tan Lei (1990) Similarity between the Mandelbrot set and Julia Sets, Communications in Math. Phys. 134, 587-617.

24.   Woon S (1998) Fractals of the Julia and Mandelbrot sets of the Riemann Zeta Function arXiv:chao-dyn/9812031v1

25.   Broughan K, Barnett R (2003) The Holomorphic Flow of the Riemann Zeta Function Mathematics Of Computation 73/246 987-1004.

Tom Beddard's Blog http://www.subblue.com/blog

Mandelbulb Skytopia Site http://www.skytopia.com/project/fractal/mandelbulb.html

Hypercomplex Iterations http://www.evl.uic.edu/hypercomplex/

Quarternionic Fractal Explorer http://www.theory.org/software/qfe/

3D Fractals Bicomplex Mandelbrot http://www.3dfractals.com/

Buddhabrot Movies http://www.superliminal.com/fractals/bbrot/bbrot.htm

Mu-Ency Mandelbrot Encyclopedia http://www.mrob.com/pub/muency.html