Hide again ...

Background:

In the iteration of quadratic polynomials f_c(z) = z² + c , the parameter c is fixed and the dynamic variable z is evolving when the mapping is applied again and again. Both c and z are complex numbers. The filled Julia set K_c contains those values z that are not going to the point at ∞ under iteration. The dynamics is chaotic on its boundary ∂K_c , which is the Julia set. The mapping f_c(z) = z² + c is conjugate to F(z) = z² in a neighborhood of ∞. The conjugacy is provided by the Boettcher mapping Φ_c(z), which satisfies Φ_c(f_c(z)) = F(Φ_c(z)). We have

(1)  Φ_c(z) = lim_n→∞ (f_cⁿ(z))^1/2n

for a suitable choice of the 2ⁿ-th root. This is made precise by the following formula, which is valid for large |z| with the principal value of the roots:

(2)  Φ_c(z) = z ∏_n=1^∞ (1 + c/(f_c^n-1(z))²) ^1/2n = z ∏_n=1^∞ (f_cⁿ(z)/(f_cⁿ(z)-c)) ^1/2n.

In the dynamics of quadratic polynomials, there is a basic dichotomy:

• When the critical orbit stays bounded, then the Julia set is connected, and Φ_c(z) extends to a conformal mapping from its exterior to the exterior of the unit disk. External rays and equipotential lines in the dynamic plane are defined as preimages of straight rays and circles, respectively. By definition, the parameter c belongs to the Mandelbrot set M.
• When the critical orbit escapes to ∞, the Julia set is totally disconnected, and the parameter c is not in M. Now Φ_c(z) does not extend to the critical point z = 0, e.g., but it is well-defined at the critcal value z = c.

The mapping Φ_M(c) is defined by evaluating the Boettcher mapping at the critcal value z = c, i.e., Φ_M(c) = Φ_c(c). Adrien Douady has shown that it is a conformal mapping from the exterior of the Mandelbrot set M to the exterior of the unit disk, thus M is connected. By means of Φ_M(c), external rays and equipotential lines are defined in the parameter plane as well. See the image above. When the argument θ of an external ray is a rational multiple of 2π, then the ray is landing at a distinguished boundary point: a rational dynamic ray is landing at a periodic or preperiodic point in the Julia set. A rational parameter ray is landing at a Misiurewicz point or a root point in the boundary of the Mandelbrot set. See the following papers by Dierk Schleicher (arXiv:math.DS/9711213) and John Milnor (arXiv:math.DS/9905169).

Computation of the external argument:

External rays are characterized by the fact that the argument arg_c(z) of Φ_c(z), or the argument arg_M(c) of Φ_M(c), equals a given value θ. For large |z|, formula (2) implies

(3)  arg_c(z) = arg(z) + ∑_n=1^∞ 1/2ⁿ arg(1 + c/(f_c^n-1(z))²) = arg(z) + ∑_n=1^∞1/2ⁿ arg(f_cⁿ(z)/(f_cⁿ(z)-c)).

Here arg(z) denotes the argument of the complex number z, i.e., its angle to the positive real axis. This angle is defined only up to multiples of 2π. E.g., a point z on the negative real axis has the arguments ±π, ±3π, ±5π ... The principal value of the argument is the unique angle with -π<arg(z)≤π. This definition is used because a function should be single-valued. However, a discontinuity is introduced artificially when a point z is crossing the negative real axis, say going from i through -1 to -i: its argument should go from π/2 through π to 3π/2, but the principal value goes from π/2 to π, jumps to -π, and goes to -π/2.

Formula (3) means that the orbit z_n=f_cⁿ(z) of z and the arguments arg(z_n/(z_n-c)) are computed. The function arg(z/(z-c)) is discontinuous on the straight line segment from the critical point z = 0 to the critical value z = c, which is drawn in red in the figure: when z belongs to this line segment, then z/(z-c) will be on the negative real axis. When the Julia set is connected, the location of the discontinuity can be moved onto a kind of curve joining 0 and c within the Julia set, such that the argument is continuous in the whole exterior:

• Suppose that z is crossing the red line from right to left, and moving to the fixed point α_c and beyond. The principal value arg(z/(z-c)) is jumping from π to -π at the red line, goes up to about -0.68π at α_c and stays continuous there.
• The modified function arg(z/(z-c)) is continuous with value π, as z is crossing the red line. It goes up to about 1.32π when z is approaching α_c from the right, and jumps to about -0.68π as z is crossing α_c to the left.

When the modified function arg(z/(z-c)) is used instead of the principal value for computing arg_c(z) according to (3), this formula will be valid not only for large |z| but everywhere outside the Julia set. In practice, the Julia set may be approximated by two straight line segments from the fixed point α_c to 0 and c, respectively: when z_n is in the triangle between 0, c, and α_c, adjust the principal value of arg(z_n/(z_n-c)) by ±2π. The triangle is given by three simple inequalities, cf. mndlbrot::turn(). In the two images below, the different colors indicate values of arg_c(z), and thus, external rays. The computed value of arg_c(z), and thus the color, is jumping on certain curves joining preimages of 0, where some iterate z_n is crossing the line where arg(z/(z-c)) is discontinuous. In the left image, the principal value of arg(z/(z-c)) is used, and in the right image, the argument is adjusted. The erroneous discontinuities of the computed arg_c(z) are moved closer to the Julia set.



The same modification seems to work for the computation of the argument arg_M(c) in the parameter plane as well, although the Julia set is disconnected. The left image below shows the discontinuities arising from the principal value of arg(z/(z-c)) in (3), they are appearing on curves joining centers of different periods. In the right image, the argument is adjusted, and the discontinuities are moved closer to the Mandelbrot set.

Drawing “all” rays:

To get an overview of all rays, compute the argument at every pixel and map it to a range of colors. The images above where obtained with Mandel's algorithm No. 6. The code of mndlbrot::turn() has been ported to Java by Evgeny Demidov. See also the algorithms 7 and 8, where no discontinuity appears, but only the ends of certain rays are visible. Adam Majewski has suggested that boundary scanning (below) can be used to draw several rays simultaneously.

Drawing a single ray:

• In the dynamic plane, external rays can be drawn by backwards iteration: It is most effective for a periodic or preperiodic angle. You must keep track of points on the finite collection of rays with angles θ, 2θ, 4θ ... Say z_l,j corresponds to a radius R^1/2l and the angle 2^jθ. Then f_c(z) maps z_l,j to z_(l-1),(j+1). This point, which was constructed before, has two preimages under f_c(z) . The one that is closer to z_(l-1),j is the correct one. This criterion was proved by Thierry Bousch. The ray will look better when you introduce intermediate points. Backwards iteration is used in Mandel again since version 5.3, both in the plane and on the sphere. The current implementation consists of QmnPlane::backRay().

The following methods apply with small modifications to the dynamic plane and the parameter plane as well:

• Choose a sequence of radii tending to 1, e.g., R_n = R^1/2n with R large. Search corresponding points on the ray, which satisfy Φ_c(z_n) = R_n e^iθ or Φ_M(c_n) = R_n e^iθ, respectively. To this end, the Newton method is employed. Using the approximation (1), the function N(z) = z - P(z)/P'(z) is iterated to find a zero of P(z) = f_cⁿ(z) - R e^i2nθ, or analogously in the parameter plane. Intermediate points should be used to make the curve between the points smooth, and to reduce the danger of the Newton method failing when the previous point is not in the immediate basin of attraction for the new point. Presumably, this method has been developed by John Hamal Hubbard or John Milnor, and it has been implemented by Arnaud Chéritat, Christian Henriksen and Mitsuhiro Shishikura as well. It is used in Mandel since version 5.4, the current implementation consists of QmnPlane::newtonRay() and QmnPlane::rayNewton(). These were ported to Java as well. The paper by Tomoki Kawahira gives a detailed description and an error estimate for the approximation (1).
• 
  There are two approaches to draw the external ray of angle θ by computing the external argument (3): Argument scanning means to check every pixel, which is very slow, but might be improved by a recursive subdivision of the image. Argument tracing means to construct a sequence of triangles along the ray. (I have learned this technique from Chapter 14 in this book, where boundary tracing is used to draw escape lines.) First two pixels at the border of the image are found, such that the external argument is <θ at P₀ and ≥θ at P₁. Then compute the external argument at the third vertex P₂:
  • If it is <θ, then the ray is intersecting the edge from P₁ to P₂. Reflect the triangle along this edge to obtain the new vertex P₂, rename the old vertex P₂ to P₀, and keep P₁.
  • If it is ≥θ, then the ray is intersecting the edge from P₀ to P₂. Reflect the triangle along this edge to obtain the new vertex P₂, rename the old vertex P₂ to P₁, and keep P₀.
  Repeat this process with the new triangle … and draw the ray along the sequence of triangles. This basic idea can be improved in several ways:
  • To make the ray look smooth, draw line segments instead of single pixels. Choose the way of reflection such that the distance to previously drawn pixels is increased.
  • Check this distance also to avoid going in circles around the endpoint. (This is not yet implemented in the current version of Mandel.)
  • When computing the external argument according to (3), adjust the principal value of arg(z_n/(z_n-c)) by ±2π, if z_n=f_cⁿ(z) or z_n=f_cⁿ(c) is in the triangle between 0, c, and α_c as explained above.
  • Save the intermediate arguments arg(z_n/(z_n-c)) for the next pixel to check for jump discontinuities by comparison, and to adjust them. (The previous technique shifts these closer to the boundary, but does not remove them.) A drawback of this approach is that it restricts the number of iterations to the length of the array, which means in particular that the end of the ray may be lost.
  • When the image shows a very small part of the Julia set or Mandelbrot set, maybe the starting point will not be found in spite of the previous techniques. It can be searched for on the boundary of a virtual enlarged image, but this will be impractical when it is necessary to enlarge the image by several orders of magnitude.
  Argument tracing has been used in Mandel since 1997, the current implementation consists of mndlbrot::turnsign() and QmnPlane::traceRay().
• Johannes Riedl has used a method, where the next point is chosen from a finite collection of points at a prescribed distance from the last one.
• The power series expansion of Φ_c^-1 or Φ_M^-1 could be used as well to draw external rays, but I expect it to converge extremely slowly close to the boundary.
• Noting that the direction of an external ray is determined uniquely from the fact that it is orthogonal to the equipotential line, external rays can be drawn by solving a differential equation numerically. I am not sure if this method will be numerically stable, and I have received different claims about it.

The inverse of Φ_c(z) or Φ_M(c) can vary considerably on a small interval: the angles 1/7 (of period 3) and 321685687669320/2251799813685247 (of period 51) differ by about 10^-15, but the landing points of the corresponding parameter rays are about 0.035 apart. On this scale, the rays cannot be drawn well by tracing the argument, unless high-precision arithmetics is used. The following image shows the parameter ray for 1/7 and two rays of period 51 in the slit betwen the main cardioid and the 1/3-component. The argument is not distinguishable (left), but the Newton method has no problem (right). In the dynamic plane, backwards iteration works as well.



Note that most of the methods above can be modified to draw equipotential lines. In this case there is no problem with the ambiguity of the argument. Boundary scanning is slow and the Newton method has three problems: finding a starting point, drawing several lines around a disconnected Julia set, and choosing the number of points depending on the chosen subset of the plane. Therefore I am using boundary tracing with starting points on several lines through the image. See QmnPlane::green(). Still, sometimes not the whole line is drawn, or not all components are found.