Smooth colouring is the key to the Mandelbrot set

At Cambridge Geek Night 6 on Monday I gave a brief talk about what has been occupying some of my spare time in recent weeks. Since Benoît Mandelbrot died last month I have been playing around with making pretty pictures of Mandelbrot sets, culminating (so far) in a movie of a deep zoom targeted near z = -0.15 + 1.04i. The software and materials for my talk are at https://dotat.at/prog/mandelbrot/.

One of the nicest renderings I did when I was a teenager was a 300dpi laser-printed binary decomposition of the level sets of the exterior of the Mandelbrot set. It took a whole afternoon to generate using BBC BASIC on my Archimedes, which was the best tool for numeric code I had available. I wanted to be able to produce smoothly coloured renderings like the ones in the fractal art books by Heinz-Otto Peitgen and his collaborators, but a few things prevented me. Firstly the lack of a computer that could display true colour images :-) and secondly the difficulty of the mathematics. If I remember correctly, I just about grasped enough of the distance estimation method to understand that it requires significantly more work per iteration, which also made it unattractive on the slow computers I had to hand.

Revisiting the Mandelbrot set now, I have the luxury of computers that are thousands of times faster and thousands of times more colourful and millions of times more accurate with their arithmetic. The Internet also provides more accessible resources on the mathematics of the Mandelbrot set. [A page by Linas Vepstas[https://linas.org/art-gallery/escape/escape.html) was particularly helpful with the smooth colouring problem. Even better, it turns out that his simple recipe arises from the mathematics that gives us our deepest understanding of the Mandelbrot set.

First some basic notation to make the formulae less cluttered.

\begin{matrix} f_{c} (z) & = & z^{2} + c \\ z_{n} & = & f_{c}^{n} (z) & = & f_{c} (\dots f_{c} (z)) \\ r_{n} & = & | z_{n} | & = & | f_{c}^{n} (z) | \end{matrix}

To generate a Julia set, c is a fixed parameter and z varies across the complex plane. For the Mandelbrot set, c varies across the plane and z starts its iteration at 0. In both cases we iterate until z escapes to infinity, in which case it is outside the set, or remains bounded, in which case it is inside. Because we are most interested in the behaviour of z as it is iterated, we usually leave the parameter c implicit.

The simplest Julia set, when c = 0, is just the unit circle. In this case we can write down a direct formula for z_n:

\begin{matrix} f_{0} (z) & = & z^{2} \\ z_{n} & = & f_{c}^{0} (z) & = & z_{n - 1}^{2} & = & z_{n - 2}^{4} & = & z_{n - 3}^{8} & = & z^{2^{n}} \end{matrix}

There is something called a Böttcher map which transforms the exterior of Julia set c = 0 into the exterior of other connected Julia sets and back again. It’s like a Riemann map, but it operates on the complement of the set.

\begin{matrix} f_{c} (z) & = & 𝚽_{c}^{- 1} (f_{0} (𝚽_{c} (z))) \end{matrix}

The definition of the map is as follows. Observe that Φ₀(z) is the identity function.

\begin{matrix} 𝚽_{c} (z) & = & lim_{n \to ∞} \sqrt[2^{n}]{z_{n}} \end{matrix}

Adrien Douady and John Hubbard showed that the Böttcher map Φ_c(0) is a conformal mapping from the exterior of the Mandelbrot set M to the exterior of the unit disk, and therefore the Mandelbrot set is connected.

The Böttcher map allows us to treat the exterior of the set as an electrostatic field. Imagine an infinitely long crinkly cylinder whose cross-section is the set. If you electrically charge the cylinder it creates a field in the surrounding space. The Böttcher map tells you how the simple circular equipotential curves and radial field lines of the Julia set c = 0 correspond to those surrounding the Mandelbrot set and the other connected Julia sets.

This is a powerful tool for understanding the Mandelbrot set. The binary decomposition makes a lot of the field structure visible, since the boundaries between the black and white areas are formed from segments of field lines and equipotential curves. If you express the angle of a field line as a binary fraction of a whole turn, you get the sequence of black and white segments that the field line crosses in the binary decomposition. Field lines with rational angles correspond to pinch points on the Mandelbrot set.

For our purposes we are less interested in field lines and more interested in the Douady-Hubbard potential function. We use the Böttcher map to warp the space around the set into its simplest form. This allows us to use the standard formula for the electrostatic potential around an infinite wire, which is the log of the distance from the wire. The base of the logarithm is related to how highly charged the wire is; we can choose it for our convenience.

\begin{matrix} 𝜙_{c} (z) & = & log | 𝚽_{c} (z) | \\ = & lim_{n \to ∞} log | z_{n}^{1 / 2^{n}} | \\ = & lim_{n \to ∞} log r_{n}^{1 / 2^{n}} \\ = & lim_{n \to ∞} \frac{1}{2^{n}} log r_{n} \end{matrix}

Since we don’t have time to count to infinity, in practice we define a bail-out radius R and stop iterating points that get further than that from the origin. The iteration count n is used for the flat blocky level-set colouring scheme.

\begin{matrix} r_{n - 1} & \leq & R & < & r_{n} \end{matrix}

The way to get a smooth colouring is to ask what fractional iteration value ν would make the point land exactly on the bail-out radius.

\begin{matrix} n - 1 & \leq & 𝜈 & < & n \end{matrix}

We can use the potential function to formalize this idea. If R is large, varying it will not make much difference to the potential we calculate for a point. So to derive 𝜈 we assume that we get the same result from the two expressions for the potential in the first line below. (I think this is what causes the small errors and discontinuities discussed by Linas.)

\begin{matrix} \frac{log R}{2^{𝜈}} & ≃ & \frac{log r_{n}}{2^{n}} \\ 2^{n - 𝜈} & = & \frac{log r_{n}}{log R} \\ 𝜈 & = & n - {log}_{2} (\frac{log r_{n}}{log R}) \\ = & n - {log}_{2} (log r_{n}) + {log}_{2} (log R) \end{matrix}

The final term in that equation is a constant, since R is a fixed parameter of our implementation. We might as well throw it away since it does not add any useful information.

\begin{matrix} n ʹ & = & n - {log}_{2} (log r_{n}) \end{matrix}

And thus we get the equation to use for smoothly colouring of the exterior of the Mandelbrot set.