The Medusa Algorithm for Polynomial Matings

The Medusa algorithm takes as input two postcritically finite quadratic polynomials and outputs the quadratic rational map which is the mating of the two polynomials (if it exists). Specifically, the output is a sequence of approximations for the parameters of the rational map, as well as an image of its Julia set. Whether these approximations converge is answered using Thurston's topological characterization of rational maps. This algorithm was designed by John Hamal Hubbard, and implemented in 1998 by Christian Henriksen and REU students David Farris, and Kuon Ju Liu. In this paper we describe the algorithm and its implementation, discuss some output from the program (including many pictures) and related questions. Specifically, we include images and a discussion for some shared matings, Lattes examples, and tuning sequences of matings.


Introduction
The study of the dynamics of rational maps of the Rieman sphere is greatly facilitated by the fact that a wide variety of dynamical phenomena can be illustrated using only the quadratic family P c (z) = z 2 + c. Of course most general theorems about rational maps have examples in the quadratic family, but further, in some cases the dynamics of a quadratic polynomial appear within a rational map. The most basic example of this phenomena is through polynomial-like behavior. In addition, there are several ways to combine two (or more) quadratic polynomials to produce rational maps whose dynamics can be described via a combination of the quadratic polynomial dynamics. Probably the first such example was a polynomial mating discovered by Adrien Douady [Dou83].
In order to define matings, first we must step back to quadratic polynomials. It is simple to write a computer program which, given a c, will compute (approximately) the orbit of any given point under the quadratic polynomial P c . To illustrate the overall behavior one draws the filled Julia set, K c , the set of points whose orbit under P c does not tend to ∞. This also illustrates the Julia set, J c , the topological boundary of K. (See §2, Figure 1 for a sample J c .) We may examine experimentally the dynamics of one map at a time with such a program.
The next natural step is to understand how the dynamics changes with a change in the parameter, c. We organize the parameter space by defining M , the Mandelbrot set, as the set of all c in C for which the Julia set J c is connected (see §2, Figure 3). By Fatou's fundamental dichotomy theorem, this is equivalent to the set of all c such that the orbit of the critical point 0 under P c lies in K c . Thus it is also a simple matter to generate a picture of M , and a program which will draw the Julia set J c when a parameter c in M is selected. After a brief investigation with such a program, one sees intriguing patterns, and a relationship between M and the Julia sets of its children, the quadratic polynomials.
In addition to the definition of M , many basic results in the theory of the iteration of rational functions support the premise that the behavior of the critical orbit is crucial for describing the dynamics. The dynamics are most amenable to analysis when the polynomial P c is postcritically finite (PCF), i.e., the orbit of the critical point 0 is finite. A key technique in giving a mathematical description of the patterns of quadratic polynomials turns out to be combinatorics. For a postcritically finite quadratic polynomial, we can build a labelled graph, called a spider, which gives a combinatorial description of the dynamics of the polynomial. This is described in §2.2.
The reverse problem, of starting with a combinatorial spider and producing a quadratic polynomial P c (i.e., producing a parameter c) whose dynamics are given by that model, is solved by the spider algorithm. The spider algorithm is an iterative procedure, based on Thurston's topological characterization of rational maps [DH93], and is described fully in [HS94].
The main subject of this paper is the Medusa algorithm, which takes two combinatorial spiders, glues them together in a certain manner (hence the name Medusa), then runs a sort of double spider algorithm which, if it converges, produces a rational map which is the mating of the two quadratic polynomials associated with the originally inputted spiders, see Theorem 3.9.
John Hamal Hubbard designed the Medusa algorithm, based on Thurston's theory ( [DH93]) and the foundational theory of polynomial matings developed by Douady, Hubbard, Shishikura, Rees, Tan Lei and others ( [Dou83,Ree92,Lei92,Shi00], see §2.3). The computer program implementing the algorithm was written under Hubbard's direction by David Farris, Christian Henriksen and Kuon Ju Liu, in a 1998 summer research experience for undergraduates program. The full source code for Medusa is available for download at [Dyn]. Some progress has been made in the study of polynomial matings since 1998, however there are still many intriguing questions. The goals of experimental software like Medusa are to help form conjectural answers to existing questions, as well as inspire new questions. After explaining the algorithm and implementation, in the final section of this paper we provide several examples of images we created using Medusa, which serve to illustrate and examine several of the phenomena of matings. Specifically, we include images and a discussion for some Lattès examples, shared matings, and tuning sequences of matings. We hope this paper will energize future researchers to study polynomial matings, and we expect Medusa is of service in advancing the field.
Organization of sections. In §2 we provide needed prerequisite material on the dynamics of quadratic polynomials and polynomial matings. In §3, we describe the Medusa algorithm and its implementation, proving Theorem 3.9. The final section, §4, contains examples of output from the program related to a few areas of interest in the study of matings.
Acknowledgements. The authors thank Dierk Schleicher, Adam Epstien and Tan Lei for inspiring discussions and advice on how to write this paper. All images of Julia sets of quadratic polynomials were generated with the Otis fractal program [Kaw].
2. Background 2.1. Notation. We writeĈ = C ∪ {∞} for the Riemann sphere, i.e., the one point compactification of the complex plane, endowed with the complex structure with respect to which the identity restricted to C is a chart, and z → 1/z a conformal isomorphism. We write S 2 forĈ viewed as a topological manifold, i.e., not equipped with a canonical complex structure.
2.2. Quadratic polynomials and combinatorics. If K c is connected, then there is a unique conformal isomorphism such that ψ c (∞) = 1. This map conjugates w → w d to P c . The curve R t (c) = R t = {ψ(re 2πit ) : r > 1} is the external ray of angle t. For a postcritically finite polynomial the filled in Julia set K c is locally connected, and then ψ c extends continuously to the boundary. If we parameterize the circle by R/Z, then the map ψ c on the boundary becomes and γ c is a semiconjugacy of multiplication by two to f , i.e., γ c (2t) = P c (γ(t)). Then γ c (t) is called the landing point of R t (c). Call γ c the Carathéodory map of P c . See Figure 1 for a picture of a Julia set and some external rays.
Given a postcritically finite quadratic polynomial, P c , choose θ c ∈ R/Z so that R θc is the external ray associated with the critical value, c. That is, R θc lands at c, if c ∈ J c . Otherwise the critical point is periodic. If the critical point is fixed, take θ c = 0. If the critical point is periodic of period n > 1, the critical value is contained in the immediate basin U of a superattracting cycle and there exists a pair of rays landing at the root of U whose closure seperates the critical value from the other points in the critical orbit. Take θ c to be one of the two angles corresponding to this pair of rays.
Given a rational number θ ∈ R/Z, following Hubbard and Schleicher ([HS94]) we define the standard θ-spider, S θ ⊂Ĉ by: See the image on the left in Figure 2 for an example, it shows the spider for one of the Julia sets of Figure 1. One may view this as a spider, with legs the rays emanating from the unit circle which are in the orbit of θ under angle doubling, and body the point at infinity.
Since γ c semi-conjugates P c to angle doubling, γ c maps S θc to the union of R θc and its images under P c , plus the point at infinity. Note if θ is rational, then it has finite orbit under angle doubling, so the spider has a finite number of legs. Similarly, if P c is postcritically finite, then θ c will be rational. We denote the endpoints on the unit circle of the spider legs by z j = e 2iπ2 j−1 θ .
The spider illustrates the critical orbit. Using this diagram we can also create a sequence called the kneading sequence of θ which records information about the order of the critical orbit in this diagram. Take the plane containing the spider S θ , and cut along the line composed by the rays of angle θ/2 and (θ + 1)/2. Label by A the open half of the plane containing θ, label the other open half B. See the right hand image of Figure 2. Label the ray of angle θ by * a , and the ray of angle (θ + 1)/2 by * b . For any angle t, its θ-itinerary is the infinite sequence of labels from (A, B, * a , * b ) corresponding to the position in the labelled plane of the points in the forward orbit of t under angle-doubling. The kneading sequence of θ, denoted k(θ), is the θ-itinerary of the angle θ. Note a symbol * n appears in this sequence if and only if θ is periodic under angle doubling. Left: the spider for θ = 1/6. The critical orbit is (1/6 → 1/3 → 2/3 → 1/3). Right: the kneading sequence for this spider is K(1/6) = A AB. This spider models f (z) = z 2 + i, whose Julia set is shown in Figure 1.
In this paper, we are interested in combining and comparing quadratic polynomials. In order to keep track of the dynamics of various maps we are studying, we use the discovery of Douady and Hubbard (see [DH82]) on how θ c relates to the position of c in the Mandelbrot set. They show the Mandelbrot set, M , is connected, with simply connected complement inĈ, hence there is a unique conformal isomorphism Ψ M :Ĉ − M →Ĉ − D which fixes ∞ and such that Ψ M (∞) = 1. Then Ψ M defines external rays outside of M , by images of straight rays outside of the disk. It happens that for any rational angle θ = p/q, the map Φ M extends radially to the boundary, to define a landing point c(θ) for the ray of angle θ. Given a postcritically finite polynomial P c to which we associate the angle θ c , then the parameter ray of angle θ c will either land at c (in the preperiodic case) or at the root of the hyperbolic component of M that has c as a center (in the periodic case). For example, for the basilica, f (z) = z 2 − 1, the external rays associated with the critical value −1 is of angle 1/3 and 2/3. The parameter rays of angle 1/3 and 2/3 lands on the Mandelbrot set at the root point of the bulb containing the basillica (the real bulb). Figure 3 shows the Mandelbrot set and some external rays.
2.3. Mating quadratic polynomials. Let f n (z) = z 2 + c n , n = 1, 2 be two quadratic polynomials, with Julia sets J n . Assume each J n is locally connected, and γ n is the Carathéodory map of f n . Define K = K 1 K 2 / to be the quotient space of the disjoint union of K 1 and K 2 in which for each t ∈ R/Z, we identify γ 1 (t) with γ 2 (−t). In other words, we obtain a topological space K by gluing K 1 and K 2 together along their boundaries via γ 1 (t) γ 2 (−t). Consider this definition while viewing Figure 1. In Figure 3. The Mandelbrot set, i.e., the set of all c in C for which the Julia set J c is connected, shown above in black, together with the external rays: 0, 1/511, 1/7, 10/63, 1/6, 3/14, 1/5, 1/4, 169/511, 1/3, 255/511, 1/2, 2/3, 5/6. general one might imagine K as some bizarre balloon animal (possibly with infinitely many body segments), but we will see below that in many cases, K is simply a sphere. On the space K, define the map f 1 f 2 by f n on K n , n = 1, 2. Since γ n semiconjugates f to multiplication by two on J n , this map is well-defined and continuous (no matter how bizarre the space K may be).
If there is a quadratic rational map F which is topologically conjugate on C to f 1 f 2 on K, then F is called a mating of f 1 and f 2 . We denote this relationship by F ∼ = f 1 f 2 , and in this case say the mating of f 1 and f 2 exists. The conjugacy h : K →Ĉ is required to be an orientation preserving homeomorphism which is holomorphic on the interiors of each K n . It is believed that if F exists, it is unique up to Möbius conjugation.
Note that a mating of any quadratic polynomial f 1 with f 2 (z) = z 2 yields F ∼ = f 1 .
Results of Rees, Shishikura, and Tan Lei ([Ree92, Lei92, Shi00]) show that whether the mating of two PCF quadratic polynomials f 1 and f 2 exists can be answered in terms of the location of c 1 and c 2 in parameter space. The fundamental existence theorem is: Theorem 2.1. If f 1 , f 2 are PCF quadratic polynomials, TFAE: • K is homeomorphic to the sphere S 2 ; • there exists a quadratic rational map F which is the mating of f 1 and f 2 ; • c 1 and c 2 do not belong to complex conjugate limbs of the Mandelbrot set, M .
We refer the reader to Milnor's book [Mil99] for detailed background on the dynamics of polynomial maps of C, and his article [Mil04] for a more complete discussion of the definition of mating and its subtleties, a discussion of many foundational results on matings, and a detailed analysis of an interesting example of mating.

From Thurston's Algorithm to the Medusa Algorithm
Thurston's algorithm is a proof that given a branched covering g of the sphere there exists a rational map F that is Thurston equivalent to g unless there exists a Thurston obstruction. The proof can be made into an iterative procedure computing a sequence of complex structures and rational maps F n which, when properly normalized, converges to F . In this section we see that we can take g to be a model of the mating of two quadratic rational maps, and extract finite dimensional but crucial information about the complex structures produced by Thurston's Algorithm so that the sequence F n can be recovered. This is the heart of the Medusa Algorithm. Because of the finite dimensional information needed to run the algorithm, it lends itself to actual computation.
Normalizing matings. Assume f 1 , f 2 are postcritically finite quadratic polynomials and F ∼ = f 1 f 2 . Each f n has one critical point 0, which lies in K n . Thus F has two distinct critical points. By conjugating F with a mobius transformation we can arrange that the critical point coming from f 1 is at the origin, the other critical point at infinity and the two glued-together beta fixed points are at 1. Therefore we know that any such mating belongs to the following family of maps.
Note that every rational map of degree two is conjugate to (at least one) member of F.
The following innocent lemma, which is trivial to prove, is of fundamental importance to why there is such a thing as the Medusa Algorithm.
The lemma shows that there is some magic to quadratic rational maps. Normalized in the way described, we just need the position of the two critical values (and which correspond to which critical point) to uniquely determine the map. We don't need any extra combinatorial information.
Proof. We prove the lemma in the case where u, v are different from infinity. The case where either u or v equals infinity is just as easy and left to the reader. First notice that F : has the desired properties, so we need to show that this is the only such map in F. Since the origin and infinity are critical points, we can write has rank 3. It follows that every solution to the three equations can be and therefore F is uniquely determined.
In the following we will write F u,v for the map given by the lemma.
The Standard Medusa. We now build a model for the mating F = f 1 f 2 of the two postcritically finite quadratic maps f 1 , f 2 . We start by defining the standard Medusa.
Definition 3.3. Let θ 1 , θ 2 ∈ Z be the two rational numbers we associate to f 1 and f 2 , as in §2.2. Define the (θ 1 , θ 2 ) standard Medusa M(θ 1 , θ 2 ) ⊂ S 2 to be the union of the unit circle S 1 , the interior legs and the exterior legs Defined in this way we have that z → 1/z maps M(θ 2 , θ 1 ) bijectively to M(θ 1 , θ 2 ). The endpoints of the interior legs we denote by x j , and the endpoints of the exterior legs we denote by y j , hence x j = 2 exp(2iπ2 j θ 1 ), j = 1, 2, . . . , and We can think of the standard Medusa as a coupling of two standard spiders S θ 1 , S θ 2 , where the bodies have been cut away, then the two are glued along the cut. See Figure 4 for a schematic diagram of this process. Thurston Matings. Recall that two postcritically finite branched coverings F : S 2 → S 2 and g : S 2 → S 2 with postcritical sets P F and P g are called Thurston equivalent if there exists orientation preserving homeomorphisms φ and ψ such that φ restricted to P F maps bijectively onto P g and ψ −1 • φ is isotopic to the identity on S 2 rel. P F .
We proceed to define a branched covering g of S 2 by itself that in nondegenerate cases is Thurston equivalent to the mating F = f 1 f 2 . Let g| M(θ 1 ,θ 2 ) be the angle doubling map r exp(iφ) → r exp(2iφ). Extend g smoothly to a degree two branched covering of the sphere so that: (1) g : D → D is a degree two branched coveing with critical value at x 1 , and (2) g : S 2 \D → S 2 \D is a degree two branched covering with the critical value at y 1 . Denote by ω 1 the critical point of g in D and by ω 2 the critical point of g in S 2 \ D. Notice that ω i coincides with an endpoint of a leg if and only if θ i is periodic under angle doubling, θ → 2θ mod 1.
Notice that if we redefine g outside the unit circle to by setting it equal to z → z 2 here, we obtain a map that is Thurston equivalent to f 1 . Similarly, if we instead redefine g inside the unit circle so it restricts to z → z 2 here, we obtain a mapping that is Thurston equivalent to f 2 . Hence it is reasonable to view g as our branched covering model of the mating F. Shishikura [Shi00] guarantees convergence in the nondegenerate case: Definition 3.4. Let f 1 , f 2 be PCF quadratic polynomials not in complex conjugate limbs of M . If the two critical orbits of F ∼ = f 1 f 2 are disjoint, then f 1 and f 2 are called strongly mateable.
Thurston's algorithm is an iterative process that will give us a sequence of rational maps converging to F when F and g are Thurston equivalent. Using g as our model map, it works as follows. Let σ 0 : S 2 →Ĉ be an orientation preserving homeomorphism mapping ω 1 to 0, ω 2 to ∞ and fixing 1. Recursively define σ n and F n as follows for n = 1, 2, . . . Interpret σ n−1 as a global chart defining a complex structure on S 2 . This complex structure can be pulled back by g. Indeed, since g is a local homeomorphism everywhere except at ω i , i = 1, 2 we can just compose restrictions of g with σ n−1 . The complex structure defined in this way can be uniquely extended to the missing points ω 1 , ω 2 . By the uniformization theorem S 2 equipped with the pullback complex structure is conformally equivalent toĈ. So let σ n : S 2 →Ĉ be the conformal isomorphisms and normalize it so ω 1 is mapped to 0, ω 2 to ∞ and 1 is fixed. By construction F n defined by the composition σ n • g • σ −1 n is holomorphic. The sequence of maps constructed can be illustrated by the commutative diagram shown in Figure 5.
In principle Thurston's algorithm solves our problem, the sequence of generated maps rational maps should converge to our mating. However the set of possible complex structures on S 2 is beyond actual computations, so we need to adapt the algorithm to allow for this. This is exactly what Hubbard's Medusa Algorithm does for us.
The Medusa Algorithm. Notice that each map F n in Thurston's algorithm (in the strongly mateable described) is a degree two rational map fixing 1 and having the origin and infinity as critical points. In other words, F n ∈ F. By Lemma 3.2 we just need to know to where 0 and ∞ are mapped to identify F n . Hence we don't need all the information contained in the sequence of complex structures to find F n , it is enough knowing σ n−1 restricted to the standard Medusa M(θ 1 , θ 2 ). Motivated by this we make the following definition. Notice there is a natural projection π from the complex structures on S 2 onto M. Given a complex structure Σ we know by the uniformization theorem that there exists a conformal isomorphism σ : (S 2 , Σ) →Ĉ which we can normalize so that ω 1 maps to 0, ω 2 to infinity and 1 is fixed. We let π(Σ) equal the equivalence class of σ| M(θ 1 ,θ 2 ) in M(θ 1 , θ 2 ).
One can show that there is a natural bijection between M(θ 1 , θ 2 ) and the Teichmüller space of S 2 \ {x 1 , x 2 , . . . , y 1 , y 2 , . . . , 1}, so Medusa space is a finite dimensional complex manifold in a natural way.
Mappings in Medusa space can be lifted. More precisely we have the following lemma.
Lemma 3.7. Let s n−1 ∈ M 0 (θ 1 , θ 2 ) be given. Set u n = s n−1 (x 1 ), v n = s n−1 (y 1 ) and let F un,vn ∈ F be the unique mapping as in Lemma 3.2. Then there is a unique mapping s n ⊂ M 0 (θ 1 , θ 2 ) such that the following diagram commutes. Proof. Since the simple closed curve γ = σ n−1 (S 1 ) seperates one critical point 0 and its image u n = F un,vn (0) from the other critical point ∞ and its image v, the preimage γ of γ by F un,vn is a simple closed curve and F un,vn : γ → γ is a two to one covering map. Identify the fundamental group on S 1 with Z so that a curve having index 1 with respect to 0 correspond to +1 ⊂ Z. Do similarly for γ and γ . Then the induced map g * : Z → Z is multiplication by two. Since s n−1 extends to a homeomorphism that maps ω 1 to 0 (s n−1 ) * : Z → Z is the identity. Finally, F un,vn maps the bounded component ofĈ \ γ onto the bounded component ofĈ \ γ which implies that (F un,vn ) * : Z → Z is multiplication by +2. Hence (s n−1 • g) * : π 1 (S 1 ) → π 1 (γ ) has the same image as (F un,vn) ) * : π 1 (γ) → π 1 (γ ). It follows by a fundamental theorem of algebraic topology that there exists a covering map s n : S 1 → γ so that g •s n−1 = s n−1 •F un,vn on S 1 , and this lift is unique when we require that s n (1) = 1. We can extend s n to M(θ 1 , θ 2 ) by lifting each leg seperately, in the way that agree with how s n is defined on the circle. In this way we have obtained a homeomorphism s n mapping M(θ 1 , θ 2 ) to its image, and we must show that s n ∈ M 0 (θ 1 , θ 2 ). However, since F un,vn maps the bounded (unbounded) part ofĈ \ γ to the bounded (unbounded)Ĉ \ γ the image of an interior (exterior) leg is interior (exterior), so we can extend s n to a orientation preserving homeomorphism of the sphere as required.
We still need to show uniqueness of s n . For s n to be an element of M 0 (θ 1 , θ 2 ) we must have s n (1) = 1 and that uniquely determines s n on S 1 . Knowing s n on the unit circle means we know to where the base point of the legs must lift and therefore there is only on extension to M(θ 1 , θ 2 ) such that g • s n−1 = s n−1 • F un,vn . Finally suppose that s n−1 and s n−1 represent the same element in M(θ 1 , θ 2 ) and let s n , s n ∈ M 0 (θ 1 , θ 2 ) be the two unique lifts. By assumption there exists an istopy connecting s n−1 to s n−1 , through maps in M 0 (θ 1 , θ 2 ). This isotopy can be lifted to an isotopy connecting s n and s n . Each map in the isotopy maps 1 to 1 so as before we can prove that it is an element of M 0 (θ 1 , θ 2 ). Let a starting point S 0 ∈ M(θ 1 , θ 2 ), be given. The Medusa algorithm consists of repeatedly applying Lemma 3.7 to get a sequence S n ∈ M(θ 1 , θ 2 ) and rational maps F un,vn ∈ F for n = 1, 2, . . . . The beauty of the algorithm is that we produce the same sequence of rational maps that Thurston's algorithm produces.
In practice, the algorithm seems to converge without assuming the maps are strongly mateable. Thus we expect that a stronger theorem holds; namely, it should be the case that anytime f 1 and f 2 are PCF quadratic polynomials in complex conjugate limbs of M , the Medusa algorithm should converge to the mating. The case not covered by Thurston's theorem is when two polynomials that are not in complex conjugate limbs have a mating with only one critical orbit. In this case naively running the Medusa algorithm produces a sequence of Medusas which does not converge (rather tends to the boundary of the Teichmuller space), but the obstruction points (the critical orbits becoming identified) are all pushed together upon iteration of the algorithm, hence the sequence of rational maps seems to converge to the mating. To prove this stronger result one could investigate how the maps in the Medusa algorithm are converging as the boundary of the Medusa space is approached. We expect the techniques of Nikita Selinger's PhD thesis [Sel10] on convergence at the boundary of Teichmuller space could be adapted to solve this question, and leave this future result to the interested reader.
3.2. The Implementation. The point of the Medusa algorithm is that it lends itself to implementation as a computer program. The implementation is an adoption of the implementation of the spider algorithm to the more general setting of quadratic rational maps.
To initiate the program, the user inputs two rational angles θ 1 , θ 2 . The implementation defines an initial Medusa s 0 : M(θ 1 , θ 2 ) →Ĉ, say close to the identity.
To describe our matings, we define a chart on F by letting R a,b : z → In this way we parametrize all the maps in F. Supposing that F ∈ F maps 0 to u and ∞ to v, we let a = v(u−1) We represent a mapping s : M(θ 1 , θ 2 ) → C by several lists of points inĈ. One list represent the image of the unit circle, and the other lists represent the images of the legs. Also we always let the list of points representing the image of the unit circle start with the point 1.
We adopt the convention that two consecutive points in the image of the unit circle or in a leg is connected by an arc of circle. For the points on the image of the circle or on the interior legs the circle chosen is that through s(y 1 ), and the arc of circle chosen is the one connecting the two points and omitting s(y 1 ). For consecutive points on the exterior legs adopt the convention that they are connected by the arc of the circle through the points and s(x 1 ). The arc is the one that connects the two points and omits s(x 1 ).
Clearly, with the information contained in the lists of points and the convention just mentioned we can reconstruct, not s, but the isotopy class of s.
An iteration consists of finding the class of the pullback of s n−1 (as in Lemma 3.7). As in the implementation of the spider algorithm we break the process down into three steps: a pullback step, a rectifying step and a pruning step.
Pullback. Given s n−1 as lists of points as described we first find F un,vn = R an,bn . This corresponds to solving (2) 1 − a 1 − b = u n = s n−1 (x 1 ) and a b = v n = s n−1 (y 1 ).
In other words Notice that R an,bn is the composition of a Mobius transformation with z → z 2 . Hence, pulling back a point consists of first pulling it back by a Mobius transformation M n and then by the square. The question that needs to resolved is, what branch of the squareroot do we need to choose. First we pullback the points corresponding to the image of the unit circle. Suppose that we have pulled back a point z k and obtained the point w k and want to pullback the next point in the list z k+1 . Pulling back first by the Mobius transformation we get that the circle through z k , z k+1 and v n becomes a circle through M −1 n (z k ), M −1 n (z k+1 ) and ∞ i.e. a line. Since the arc of circle connecting the two points was chosen to be the one that did not contain v n the pullback of the arc of circle by the Mobius transformation becomes simply a line segment between M −1 n (z k ), M −1 n (z k+1 ). The preimage of a line by the square is a hyperbola, the two branches of which are contained in opposite quarter planes. Hence knowing one preimage w k , we need to choose the square root so that w k and w k + 1 lies in the same halfplane. So the pullback the points corresponding to the circle we construct to lists, A, B. The first element of A is 1 and the first element of B is −1, i.e. the two preimages of 1 by R an,bn . This was the first step. Next we iterate through the remaing points in the list. The k'th step consists in finding the two preimages of z k , call them w k and w k . If the last inserted point in the list A lies in the same quarter plane as w k then we insert w k in A, and w k in B. Otherwise we insert w k in B and w k in the A. It is easy to verify that the points in the list A are images of the points on the unit circle of angles in the interval 0 ≤ θ < π, whereas the points in B correspond to angles θ with π ≤ θ < 2π. Having pulled back all the points we can concanate the two list so we get one list (starting with the point 1) representing the image of the circle by s n . Notice that this list contains twice the points of the one we have just pulled back.
Next we pull back the interior legs. The leg corresponding to angle θ is the preimage of the leg corresponding to angle 2θ. If 0 ≤ θ < π the point in the list A list that is the preimage of the anchor point of the leg of angle 2θ will be the anchor point of the new θ leg, otherwise it will be the corresponding point in the list B. Hence we have already computed (and can localize) the pull-back of the first point in the leg. Hence, as before we can pull back the rest of the leg, we need to chose the square root so the consecutive points lies in the same halfplane.
Pulling back the outer legs is essentially the same, except that now pulling back by M n two consecutive defines a arc of circle, where the circle goes through 0. However since z → z 2 commutes with z → 1/z we can write the squre as the composition of 1/z, z 2 and then 1/z again. Hence pulling back by M n and then making the change of coordinates w = 1/z we are back in the same situation as the one we were facing when pulling back the interior legs.
In this way we obtain a list of point representing the map s n . However, the points are now connected by arcs of hyperbolas and not arcs of circles. The next step, rectifying, remedy this situation.
Rectifying. Perhaps a better word for the second part of an iteration would be circlyfying. We want to bring us back to the starting position where consecutive points in the lists are connected by arcs of circles. This is the most delicate part of the implementation. What we want to do is replace the arcs of hyperbolas with arcs of appropiate circles without changes the isotopy class of the corresponding element in M(θ 1 , θ 2 ). So given two consecutive points z 1 , z 2 we want to see if there is a homotopy from an arc of hyperbola to an arc of circle so that the intermediate curves does not cross any of the distinguished points s n (x 1 ), s n (x 2 ), . . . , s n (y 1 ), s n (y 2 ), . . .}. It is rather tedious so we will only outline how it is done. The circle and the hyperbola are two (real) quadratic curves and we first find their intersection. This can be boiled down to finding the roots of a degree 4 equation in one real variable. However, since we know that z 1 and z 2 lies on both curves, we can do a division of polynomial and the remaing points (if any) can be found by solving a quadratic equation. The most difficult case when the branch of hyperbola containg z 1 and z 2 intersect the circle in four points. Then the union of the circle and the branch of hyperbola cuts the plane into six parts. By elementary geometric reasoning, one can find exactly to which of the six parts a given point belongs, and this knowledge is enough to decide if the homotopy exists.
If the homotopy exists then we can move on, but if it doesn't we need to do something. What we do is to subdivide the arc of hyperbola in two halves, z 1 , ζ and ζ, z 2 and recursively rectify each half. In case we are not dealing with a leg terminating at a distinguished point, then by compactness the distinguished points are a definite distance away from the arc of hyperbola between z 1 , z 2 . Given any > 0, any fine enough subdivision of the arc of hyperbola, z 1 , ζ 1 , ζ 2 , . . . , ζ k , z 2 will satisfy, that if we replace the parts of hyperbolas with arcs of circles we will stay with a spherical neighborhood of the original arc of hyperbola. Hence, we are able to rectify after adding only a finite number of points. In the case that the arc of hyperbola terminates in a distinguished point z 2 then we are dealing with the image of a leg. It is not difficult to see that we do not change the isotopy class of s n by allowing the homotopy to cross z 2 . In practice, this means that when rectifying a leg, we do not consider the endpoint of the leg a distinguished point, and we are sure that we can rectify adding only a finite number of points.
Pruning. After pulling back and rectifying, we have new lists of point representing s n , but the number of points representing the image of the unit circle has at least doubled. This means that unless we do something we will run out of memory in a finite number of iterations.
What we do is pruning which amounts to checking each point z 2 that is not the attachment point or terminal point of the leg whether it can be removed without changing the isotopy class of the represented map. In practice this means checking whether two arcs of circles, one through z 1 and z 2 the other through z 2 and z 3 , can be replaced by an arc of circle going from z 1 to z 3 without changing isotopy class. Using a Mobius transformation to change coordinates the question becomes whether a line segment (w 1 , w 2 ) and a line segment (w 2 , w 3 ) can be homotopied to a line segment (w 1 , w 3 ) without crossing distinguished points, a question that can be easily answered.
Drawing the Julia set. In addition to producing a sequence of maps R an,bn converging to the mating, the Medusa algorithm can be used to draw successive approximations to the Julia set of the mating. At the beginning of the program, a "painted" sphere K 0 is created, with each point in the upper hemisphere painted black, and each point in the lower hemisphere painted white (or clear). At each iteration of the algorithm, given parameters a m , b m Figure 6. Each of the three columns above shows Maple output of the actual Medusas used in the iteration of the Medusa algorithm for the mating of 1/7 with 1/3 (rabbit mate basilica). In each column, the top figure is the Medusa on the sphere, the lower figure is the Medusa projected onto the plane. Leftmost is the initial Medusa, central is after 2 steps, rightmost is after 20 steps. and a painted sphere K m−1 (i.e., a sphere with each point marked one of black or white), the program computes the pull back of K m−1 by R −1 am,bm , to create K m .
When the sequences (a m , b m ) converge, then R am,bm converges to R a,b ∼ = f 1 f 2 , and K m converges to K, with white or clear marking the Julia set of f 1 , and black the Julia set of f 2 .
For example, let c 1/4 be the parameter which is the landing point in the Mandelbrot set of the external ray of angle 1/4 (c 1/4 ≈ −0.228 + 1.115i). This is a tip point on the rabbit bulb. The mating of z 2 + c 1/4 with itself exists, and is studied in detail in [Mil04]. In this case the Julia set of the mating is the entire sphere, so the approximations K n drawn by Medusa are particularly interesting. Figure 7 shows approximations K 6 , K 10 , and K 14 for this mating. Also see §4.3 for other similar examples.
The full source code for Medusa is available for download at [Dyn]. There are still a few bugs, most notably: when mating with a p/q where q is even, the algorithm will converge properly for a few steps, then start diverging.

Examples
In this section we discuss several types of matings with different properties. For simplicity, we will refer to a PCF quadratic polynomial simply by its rational angle θ c = p/q, or sometimes f p/q .

4.1.
Simple examples. We explain our first example of an image of a mating produced by the Medusa algorithm in detail. We will mate the two quadratic polynomials shown in Figure 8: f 1 will be the rabbit, 1/7, and f 2 will be the basilica, 1/3.
Let F = f 1 f 2 = 1/7 1/3. The rightmost sphere in Figure 8 illustrates the Julia set of the mating F . Due to our normalization (Equation ??), the critical point 0 of f 1 is always at z = 0 in the sphere, shown as the south pole, and the critical point 0 of f 2 is sent to z = ∞ in the sphere, shown as the north pole. The portion of the filled Julia set of the mating F Figure 8. From left to right: The Julia set of the rabbit, critical angle 1/7, then the Julia set of the basilica, critical angle 1/3, both shown with both sets of critical orbit rays (1/7, 2/7, 4/7, 1/3, 2/3) for comparison; finally, the mating 1/7 mate 1/3 on the sphere, with 1/3 in black, and 1/7 clear. which corresponds to J(f 1 ) (the rabbit) is shown in clear, and "centered" about the north pole. The portion corresponding to J(f 2 ) (the basilica) is shown in black on the front half of the sphere, and grey on the back half (to indicate that to see this, you are looking through J(f 1 )). However, due to the symmetry of the Julia sets of quadratic polynomials, this image is invariant under 180 degree rotation about the vertical axis, hence the grey image in the back does not convey new information. Also, the fixed point z = 1 (corresponding to the β-fixed points of f 1 , f 2 ), is in the dead center of the image, in the front. Note reversing the order of mating, drawing the image of 1/3 1/7, would have the effect of a 180 degree rotation about the central horizontal axis (from z = 1 to z = −1), and flipping the colors.
Self-mating. The limb of the mandelbrot set enclosed by rays of angle 1/3, 2/3 (see Figure 3) is the only limb which is its own complex conjugate. As such, any PCF quadratic polynomial which is not in that limb can be mated with itself. Such a mating clearly has extra symmetries. The leftmost image in Figure 9 is the rabbit 1/7 mated with itself. We discuss self matings more in §4.4.
Tuning. One simple way to make a mating more complicated is by tuning one of the quadratic polynomials. The result shows up as you would expect. In figure 9, compare the rabbit mate rabbit on the left with the right figure, in which the clear rabbit has been tuned with a basilica. We explore further expectations (and surprises) concerning tunings in §4.5.

Shared Matings.
One of the intriguing observations in the study of matings is that it can happen that two distinct pairs of PCF quadratic polynomials give rise to the same mating F . If f 1 f 2 ∼ = F ∼ = f 3 f 4 , and f 1 = f 3 or f 2 = f 4 , then we call F a shared mating. The simplest kind of shared mating is when f 1 f 2 ∼ = f 2 f 1 . For example, the left side of Figure 10 illustrates such a shared mating of the rabbit (1/7) and aeroplane (3/7). Of course, taking a shared mating and performing the same tuning on each quadratic polynomial will produce another shared mating, for example as on the right side of Figure 10.
It is not known whether there is a bound on the number of ways in which a quadratic rational map can be realized as a mating. The quadratic polynomials involved above are: f 1/6 (z) = z 2 + i, a tip point on the rabbit limb; f 5/6 (z) = z 2 −i, the complex conjugate of f 1/6 ; f 5/14 , a tip point of the bulb on the basilica bulb corresponding to the rabbit; and f 1/2 (z) = z 2 − 2, the real tip point of the basilica limb (the leftmost point in the mandelbrot set). The Julia set for each of 1/6, 5/14, 3/14, 1/2 is a dendrite, hence has empty interior. For example, the Julia set of f 1/4 is a dendrite, shown in Figure 7. Below is a characterization of when this occurs.
Fact 4.1. Suppose P c is a PCF quadratic polynomial. Let θ c = p/q be a reduced fraction. TFAE: (1) K c has empty interior; (2) q is even; (3) θ c is strictly pre-periodic under angle doubling.
Thus the mating of any two quadratic polynomials satisfying Fact 4.1 (including the shared mating above) has Julia set the entire Rieman sphere. You can visualize such a mating as a space-filling curve on the sphere (each Figure 10. Upper left: the rabbit, 1/7; Upper right: the aeroplane, 3/7, both shown with both sets of critical orbit rays (1/7, 2/7, 4/7, 3/7, 6/7, 5/7). Lower left: the shared mating the rabbit mate the aeroplane, 1/7 mate 3/7, equivalently, the aeroplane mate the rabbit. Lower right: basilicas in the rabbit mate basilicas in the aeroplane, 10/63 mate 28/63. of the empty interior Julia sets is a curve which is pulled into becoming a space-filling curve). Further, since the Julia set of f 1/2 is a line segment, any mating of the form p/q 1/2 where q is even will create a space-filling Peano curve.
Since the Julia set is the entire Riemann sphere, we cannot very well study such matings by drawing their Julia sets. The harmonic measure supported on the Julia set is an object which deserves further study. One could hope to learn something by examining the approximations to the Julia set drawn by the program Medusa in the steps of the algorithm converging to the mating. See Figure 11. 4.4. Self Matings. Carston Peterson has observed that if f is any PCF quadratic polynomial which is not in the 1/2-limb of the Mandelbrot set (i.e., not in the unique limb which is its own complex conjugate), then the following two rational maps are topologically conjugate: Figure 11. Each of the four images above illustrates a Medusa approximation K 12 to the same shared Lattés mating. Upper left: 1/6 mate 5/14. Upper right: 3/14 mate 3/14. Lower left: 1/2 mate 3/14. Lower right: 1/2 mate 5/6. (Note the two lower figures are mated in reverse order from the shared mating. Just rotate the picture 180 degrees and exchange the colors to see the correct image).
(1) start with f f , then mod out by the obvious symmetry, and (2) f f 1/2 , where f 1/2 (z) = z 2 − 2. This is because for f 1/2 , the Julia set is a line segment, [−2, 2], and every external ray of angle θ has the same landing point as the ray of angle 1 − θ (the ray 0 is horizontal and lands at 2, the ray of angle 1/2 is horizontal and lands at −2).
For example, shown in Figure 12 is the Julia set of f 1/5 , together with the Julia sets of both the self mating of f 1/5 , and the mating of 1/5 with 1/2. Since the Julia set of 1/2 is simply a line segment, note in the figure how this simple segment is twisted to fill up all of the black. 4.5. Sequences of matings, and their limits. One question about matings which has yielded an interested study is: If f 1 and f 2 are quadratic polynomials not in complex conjugate limbs, which are not PCF, when does a mating exist (assuming connected Julia sets)? If f 1 and f 2 are hyperbolic, thus stable perturbations of hyperbolic PCF polynomials g 1 , g 2 , each The Julia set of f 1/5 , which is the center of largest baby Mandelbrot set off of the rabbit bulb, shown with critical orbit rays 1/5, 2/5, 4/5. Upper Right: 1/5 mate 1/5. Lower Right: 1/2 mate 1/5, i.e., mod out the upper figure by the obvious symmetry. Lower Left: an approximation K 16 , to 1/2 mate 1/5. The black is 1/2, so shows the simple line twisting to fill up the alloted space.
with a super attracting periodic cycle, the mating exists as a deformation of the mating of g 1 , g 2 . Several papers have appeared constructing matings between particular non-hyperbolic polynomials (see Haïssinksy and Tan Lei [HL04], Luo [Luo95], Yampolsky, Zakeri [YZ01].) However, Epstein [Eps] has shown that mating does not extend continuously to the boundary of the hyperbolic component (in fact, the set of points in ∂M × ∂M where there is no continouous extension is dense). Epstein's theorem is that an obstruction to continuously extending this map to a mating between the two root points of the hyperbolic components occurs whenever in the mating g 1 g 2 , the immediate basins of the superattracting cycles of g 1 , g 2 touch along a distinguished repelling cycle (excluding g i (z) = z 2 ). For example, this occurs in the mating of the rabbit and the aeroplane, Figure 10. That this is a shared mating is an additional coincidence, not needed for Epstein's theorem. We can use Medusa to see a different type of example of why mating as a map from M × M to the space of quadratic rational maps is not continuous. We examine a few convergent sequences of quadratic polynomials, θ m , ω m → θ, ω, as m → ∞, such that the mating θ m ω m exists for every m, but θ ω either does not exist, or is not the limit of θ m ω m .
Below are some simple examples of sequences with no limit, or the wrong limits.
(1) First consider θ m = ω m = 1 2 m −1 , so θ = ω = 0. Note 0 corresponds to z → z 2 , so θ ω = 0 0 is just z → z 2 , with Julia set the circle. However, Medusa output suggests that the Julia set of θ m ω m is much more complicated than the unit disk. The leftmost image in Figure 13 shows the Julia set of 1/255 1/255 (recall Figure 9 shows the first element of the sequence, 1/7 1/7).