**Chris King March-Dec 2009**

**Mathematics Department**

**University of Auckland**

- Dark Heart Viewer: Application - Source - Flight Manual Latest version 1.3.1 with high degree functions and movie making.
- Riemann Zeta Viewer: Application - Source - Flight Manual
- Wave Function Method Viewer: Application - Source - Flight Manual

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

2.
Illuminating the Writhing Dark Heart of
Complexity

3.
The Dark HeartÕs Magic Numerology

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

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

6.
The Dark Heart of Division: Rational Functions

7.
Transcendent
Chaos out of Quadratic Disorder

8.
Falling into the Heart of Deep Transcendental
Variation

9. Getting to the Heart of the Riemann Zeta Function

10. Appendix 1: Combined Methods of Depicting
Julia Sets

11. Appendix 2: Ray Tracing Hypercomplex and Multi-dimensional Chaotc Iterations

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.

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 : (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 plotted by modified inverse
iteration (see below). (Top left) Partitioning into which iterate of the
critical point makes the closest approach. (Bottom left) Internal level sets of
periodic attractors.

Dark Heart Viewer: Application - Source

XCode Mandelbrot/Julia Explorer of Instructions - Source Code

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.

XCode viewer of Binary Decomposition Source Code

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 .

Dark Heart Viewer: Application - 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.

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 180^{o} 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 2* ^{n}* periodic bulbs and dendritic islands of a given period

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.

Wave Function Method Viewer: Application - Source - Flight Manual

**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).

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, É , 2*n*
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, É , 2*n*+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

Quicktime movies of Julia dynamics by wave function: Parabolic, Siegel disc, Superattracting period 3, Period 3 Island, Seahorse

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.

Dark Heart Viewer: Application - 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.

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).

Dark Heart Viewer: Application - Source includes

Parameter planes of 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.

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.

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 *z*^{5}/5+*cz*^{4}/4-*z ^{2}/2*-

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 *z*^{3}-*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 - 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).

Dark Heart Viewer: Application - 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 180^{o} 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 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.

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

Quicktime movie of cCosz roots by Gaussian wave function

Quicktime Movie of cCosz Julia Set Explosion and Period Doubling

Quicktime Slow Motion of the Explosion

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 .

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.

Dark Heart Viewer: Application - 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 - 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 - 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 8^{th} sector, (2) in the 6^{th} and 7^{th}, (3) in the 4^{th} (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 3^{rd} 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
.

XCode Explorer of
. Source code

Quicktime movie of roots by Gaussian wave function

Qucktime movie of periodicities of the parameter plane by radial wave function

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*)=*cz*log(*z*), *f*(*z*)=exp(-*z*^{2})+*c, f*(*z*)=*czz ^{z }*=

Universality of the dark heart examplified (left) in the parameter plane of *f*(*z*)=*cz ^{z+10 }*and (right) the close relationship

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

Video of journeys around f(z)=cz^{z+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.

Gives excellent high performance rendition of the RZ function and its Mandelbrot and Julia sets and those of related functions.

Quicktime movie of bifurcations in the RZ Julia set

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]:

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 :

The
Lanczos approximation is defined as follows:

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

XCode Collatz viewer - Source code

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.

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.

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).

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.

**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. Source Code

**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.

**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.

XCode Viewer for Pickover stalks with Continuous Potential - Source Code

Fractal Domains Explorer

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

Matlab m-files of: Wave function Julia, Wave function Mandelbrot, Gaussian Wave Periods Mandelbrot

**
**

**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 4*v,* 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.

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 Yvo Smellenbergh 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).

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.
Source XCode of: Quadratic attractor fill, Binary Decomposition, Continuous Potential and Pickover Stalks, Distance Estimator, MIIM, Layered Cubic, 2-parameter Herman, Cosine, Sinz+Sin2z/2, Cosz+Cos2z, Cos(Sqrt(z)), Layered Sinz/z

3.
Pixel Bender Source Code Quaternion Julia, 3D Mandelbulb

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

5.
Matlab m-files of: Herman ring, Siegel Disc, Wave function Julia, Wave function Mandelbrot, Gaussian Wave Periods Mandelbrot, Riemann Zeta, Buddabrot

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

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.

MegaPOV http://megapov.inetart.net/download.html

Fractal Domains __http://www.fractaldomains.com/download.html__

Adobe Pixel Bender http://labs.adobe.com/downloads/pixelbender.html

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