An effective algorithm to compute Mandelbrot sets in parameter planes

In McMullen (2000) it was proven that copies of generalized Mandelbrot set are dense in the bifurcation locus for generic families of rational maps. We develop an algorithm to an effective computation of the location and size of these generalized Mandelbrot sets in parameter space. We illustrate the effectiveness of the algorithm by applying it to concrete families of rational and entire maps.

In non linear science two fruitful ideas has been considered over the years to have a better understanding of the dynamics of mathematical models. On the one hand, the decomposition of the phase space of a particular dynamical system into several pieces, so that the restricted dynamics is more tractable; and on the other hand, the exploration of the changes that occur in the phase space under small perturbation of the original dynamical system. Under the general framework of smooth dynamical systems these ideas lead, respectively, to the study of invariant manifolds and bifurcation diagrams.
Accordingly, when working on mathematical models it is of high interest to know in advance, before studding a concrete dynamical space, the parameters for which the model may present undesirable dynamics. In other words, we want to characterize in the parameter space the set of parameters for which in dynamical plane contains open sets of initial conditions incompatible with the dynamics we are interested with.
In holomorphic dynamics, i.e. discrete dynamical systems generated by the iterates of holomorphic maps, the phase space is the complex plane and it is decomposed into two completely invariant sets: the Fatou and the Julia set. Moreover inside the Julia set of many rational and entire maps is very common to observe copies of the Julia set of the quadratic family. In a similar way in the parameter space we distinguish between the hyperbolic and the bifurcation parameters, and is very common to observe copies of the Mandelbrot set in the parameter space of many families of rational and entire maps.
In this sense the building block in holomorphic dynamics is the Quadratic family given by Q c (z) = z 2 +ξ where ξ is a complex parameter. The phase space or dynamical plane is divided into two completely invariant subsets: the Fatou set and the Julia set. The Fatou set, denoted by F(Q ξ ), is defined as the set of points z 0 where the family of iterates {Q n ξ (z 0 ) , n ∈ N} is normal in some neighbourhood of z 0 , and its complement is the Julia set, denoted by J (Q ξ ). The Fatou set is formed by points with tame dynamics while its complement, the Julia set, is fulfilled by points with chaotic dynamics.
For the quadratic family, the central object of the parameter plane, or ξ−plane, is the Mandelbrot set (see Figure 1(a)), denoted by M 2 , defined as the set of ξ−parameters such that the critical orbit {Q n ξ (0)} n≥0 remains bounded as n tends to ∞ (see also the definition of M d below). Easily, one can show that the Mandelbrot set is a compact set. Each one of the black bubbles is called a hyperbolic component of M 2 and corresponds to ξ-parameters for which the dynamical plane has an attracting periodic orbit of the same period. Moreover, in each bubble the Julia set moves continuously with respect to ξ, having a (unique) center for which the corresponding attracting periodic orbit has zero multiplier, i.e., it is superattracting. In contrast, the bifurcation parameters, consists of the set of parameters ξ such that the Julia set does not move continuously respect to ξ. The Misiurewicz parameters, i.e., ξ−parameters such the critical orbit (the orbit of z = 0) lands in a repelling periodic orbit, form a dense subset of ∂M 2 . Proving that the interior of the Mandelbrot set consists only on hyperbolic parameters is still a challenged open question. For an introduction to complex dynamics and the quadratic family we refer to [DH84,DH85a,McM94,Mil06,CG93,Tan00] and references therein.
For more general (uni-parametric) families of holomorphic maps a similar approach has been considered to split the dynamical plane into the Fatou and the Julia sets as well as to define the bifurcation locus. In this situation is very common to observe in the dynamical plane copies of the Julia set of the quadratic family and in the parameter plane copies of the Mandelbrot set. One of the cornerstone tools for this assertion is the theory of polynomial-like mappings, developed by A. Douady and J.H. Hubbard in the landmark paper [DH85b]. They rigorously explained why it is so common to observe, in computer experiments, homeomorphic copies of polynomial Julia set in dynamical planes of rational or transcendental families, and small homeomorphic copies of (generalized) Mandelbrot set in the corresponding parameter planes. Remember that for a given d ≥ 2 the generalized Mandelbrot set is defined as In Figure 1(a)-(b) we show the parameter planes of z d + ξ for d = 2 and d = 3, respectively. The black bubbles, as before, are hyperbolic components and the boundary of M d is the bifurcation locus. This central contribution due to A. Douady and J.H. Hubbard has been, and it is still being, applied widely. In 2000 C. Mcmullen [McM00] push forward the arguments above by showing that the existence of polynomial dynamics in families of rational maps is, in fact, universal. More precisely, he proved that small generalized Mandelbrot sets are dense in the bifurcation locus for generic families of rational maps, and so pieces of the Julia set for those non-polynomial maps resembles the Julia set of polynomials.
This result has important theoretical implications. For instance, let P be a polynomial of degree d ≥ 2. A natural root finding algorithm to solve the equation P (z) = 0 is given by the Newton map It is well known that the zeros of P correspond bijectively to the fixed points of N P which are all (locally) attracting. Roughly speaking, if there are no other stable dynamics for N P different from the attracting basins of the attracting fixed points, we can ensure that N n P (z 0 ) will numerically converge to one of the roots of P for almost all initial conditions z 0 ∈ C. However, if other stable dynamics coexist in the dynamical plane this will no longer be true. Precisely the existence of these bad initial conditions in dynamical plane is in direct correspondence with the existence of generalized Mandelbrot sets in parameter plane. So, The parameter plane of z 2 + ξ.  the natural question is where are they? Where, in parameter plane, are located the universal Mandelbrot sets predicted by Mcmullen? The main goal of this article is to give an effective algorithm to compute hyperbolic components (location and size) of the copies of the Mandelbrot set lying in the parameter plane of a uni-parametric holomorphic family (see section 2). It is worth to be noticed that to implement the algorithm it is not required to run over Mcmullen's paper (see section 1), and, more interesting, can be applied to general holomorphic, not necessary rational, maps. We illustrate, by means of different sort of examples, the use of the algorithm in Section 3.

Small Mandelbrot sets in the bifurcation locus
To provide a solid theoretical content to the algorithm we will describe in next section, we firstly discuss the key results in [McM00]. Even though we will apply the algorithm to transcendental maps we use a rational framework, as it is done in [McM00], for the theoretical discussion. At the end of the section we will briefly argue how to deal with the general scenario.
We consider a holomorphic family of rational maps over the unit disc D as a map requiring that for each parameter t ∈ D, the rational map f t :Ĉ →Ĉ satisfies deg(f t ) ≥ 2. The Bifurcation locus, denoted in what follows by B(f ) ⊂ D, it is defined as a set of t-parameters such that the Julia set J (f t ) does not vary continuously (in the Hausdorff topology) over any neighbourhood of t. Without lost of generality we assume that 0 ∈ B(f ). Indeed, we can always arrange the base point using the change of parameters t := t 0 + t. See Subsection 3.2. A local bifurcation is a holomorphic family of rational maps f t over the unit disc such that t = 0 ∈ B(f ), and for simplicity we say a local bifurcation f . A point c t ∈ C such that The following definitions and proposition play a crucial role in the description of B(f ) (see details in [McM00]).
Definition 1. Let f be a local bifurcation. A marked critical point c is a holomorphic map c : D →Ĉ such that f t (c t ) = 0. To fix the notation we will write (f, c) as the holomorphic family of rational maps with a marked critical point. We say that the marked critical point c of f is active if the sequence of iterates {f n t (c t )} n≥0 fails to be a normal family in any neighbourhood of t = 0.
Definition 2. Let f be a local bifurcation. We say that (f, c) is a (degree d) Misiurewicz bifurcation if the following conditions hold. He shows that any Misiurewicz bifurcation produces a cascade of generalized Mandelbrot sets (corresponding to values of t n tending to t = 0). The precise statement is given below. We use diam(A) to denote the diameter of a set A ⊂ C and a n = O(b n ) if |a n | < C|b n | for some constant C > 0. Theorem 1.2 (Theorem 4.1 in [McM00]). Let (f, c) be a (degree d) Misiurewicz bifurcation. Then, the parameter space D contains quasiconformal copies M n d of the degree d Mandelbrot set M d , converging to the origin, with ∂M n d contained in the bifurcation locus B(f ). More precisely, for all n 0 there are homeomorphisms Although details in the proof of the above result are far from being straightforward, the idea of the construction is quite simple and useful to define the algorithm in the next section.
For small values of t we have two crucial dynamical facts. On the one hand, the fixed point p 0 := f 0 (c 0 ) moves continuously, that is, there exists p t such that p 0 := f 0 (c 0 ), f t (p t ) = p t and |λ t | := |f t (p t )| > 1. On the other hand, because of condition M3, we can assume that f t (c t ) is not anymore a fixed point but, by continuity, f t (c t ) is close to the repelling fixed point p t . See Figure 3(b). We claim that for n large enough there exists a small t such that f n t (c t ) = c t and so we have created a periodic point with multiplier 0, so a parameter corresponding to the center of the main cardioid of a generalized Mandelbrot set (degree depending on the multiplicity of the critical point c t ). To see the claim we observe that c 0 is unramified (infinitely many backward images) and p 0 belongs to the Julia set, so there must be a point b 0 near p 0 such that f n 0 (b 0 ) = c 0 for n large. The value of the parameter r in Theorem 1.2 (b) is equal to the multiplicity, in terms of the variable t, of the map f t (c t ) at t = 0; or in other words f t (c t ) = f 0 (c 0 ) + O(t r ). For a complete proof se [McM00].
Remark 1. Since all arguments used in [McM00] are local we can extend those results to other holomorphic families of non-rational maps, like transcendental entire.

The algorithm
In this section we present an algorithm to compute the location and size of the small generalized Mandelbrot sets in the bifurcation locus B(f ) coming from the Misiurewicz bifurcation using the estimations obtained in Theorem 1.2. Next we will effectively use the algorithm in some examples.
We assume that (f, c) is a (degree d) Misiurewicz bifurcation. In order to be more precise we denote the holomorphic family by f (t, z) := f t (z) to emphasize the dependence on the parameter t and the variable z. In a similar way the marked critical point c t is given by the holomorphic map c : D →Ĉ, so we write c t = c(t). We define the following auxiliary maps, that is G k , k ≥ 0, is a function depending on t giving the k−th iterate of the critical point c(t), for the corresponding map f t . We have the following recursive expression, for k ≥ 0, .

Some computations show that
The center of the generalized Mandelbrot set M d of the family z d +ξ is the parameter ξ 0 = 0, so the center of M d corresponds to the polynomial z d having a superattracting fixed point at the origin. See Figure 1. Similarly, we define the size of M d measuring the distance between two concrete parameters in M d . In order to do this, we consider the distance, denoted by L, between the center ξ 0 of M d and a parameter ξ 1 such that z d + ξ 1 exhibits a superattracting two cycle; thus L = |ξ 0 − ξ 1 |. See also Figure 1. Easy computations show that there are d − 1 possible choices for the parameter ξ 1 given by d−1 √ −1. However, these d − 1 values are symmetrically distributed around the origin with equal distance L = 1 to the origin. By definition the center and the size of the generalized Mandelbrot set M d is cero and one, respectively.
Since the definition of the center and the size of M d are in dynamical terms, we can easily transfer these two definition to any small copy of M d .
Definition 3. Let M n d be a homeomorphic copy of M d (see Theorem 1.2). We define t n := φ n (0), the center of M n d ; or equivalently, the parameter such that the marked critical point c(t n ) has a superattracting orbit of minimal period n, and s n := φ n ( d−1 √ −1) (we choose one of the roots), a parameter such that the critical marked point c(s n ) exhibits a superattracting orbit of minimal period 2n. Then the size M n d is defined by L n = |t n − s n |. Under this notation the main concern is to numerically compute t n (location) and L n (size) of the small generalized Mandelbrot sets in the parameter plane of a holomorphic family. However we notice that, at least theoretically, the computations of those values using the auxiliary maps G k 's is immediate.
Lemma 2.1. Assume the notation of Definition 3. The values of t n and s n are solutions of the following system of equations Of course L n = |t n − s n |.
Proof. The proof is straightforward since G n (t n ) = f n (t n , c(t n )) and G 2n (s n ) = f 2n (s n , c(s n )).
We will use Newton's method to numerically solve these equations. Nonetheless remember that the small copies of the Mandelbrot set we are looking for are dense in the bifurcation locus which implies that it is crucial to have an arbitrarily good accuracy on the initial condition in order to find the correct parameters. In particular we notice that if t n solves the first equation, then it also solves the second equation. To attain this goal we will use a continuation method.
Computing the location t n of M n d . We will compute t n recursively. Assume that we know with high accuracy the location t n−1 . From Theorem 1.2(b) we get |t n | ≈ |λ 0 | −n/r , |t n−1 | ≈ |λ 0 | −(n−1)/r and we conclude |t n | ≈ |t n−1 | · |λ 0 | −1/r . So, we numerically compute t n as a result of the iterative process (Newton's method applied to the equation G n (t) − c(t) = 0) given by We can numerically compute G n (t) and G n (t) using the recursive expressions obtained in Equations (2) and (3).
Computing the size L n of M n d . In this case we also assume that t n−1 and L n−1 are known values and we will compute the size L n = |t n − s n | of M n d . As before we can obtain s n solving (4) using the Newton's method from an initial condition close to s n .
We firstly find an initial condition close to s n . In this case we identify the diameter of M n d with L n and using Theorem 1.2 estimations (b) and (c), we have that so, we conclude that s n ≈ t n 1 + L n−1 t n−1 · |λ 0 | − 1 d−1 · d−1 √ −1 .
Finally, we use this initial condition to solve numerically the equation G 2n (t) − c(t) = 0. We construct the following sequence of iterates (η k ) k≥0 tending to s n as k tends to ∞, We remark that in the above algorithm to compute t n and L n it is crucial to have an accurate initial condition. Assuming that the rational map f t has degree d, we have that the n−th iterate f n t has degree d n . Hence the number of solutions of (4) grows exponentially. From the numerical computation point of view of the location t n and the size L n we can also visualize the results obtained in Theorem 1.2(b-c). Thus for example simple computations show a linear relation between log (t n ) and log (L n ) with respect to n. More precisely, In the next section we illustrate how to implement the above algorithm to compute the location t n and the size L n of a cascade of small copies of the Mandelbrot set M d for different families covering rational as well as transcendental scenarios.
where d ≥ 2 and t is a complex parameter. To simplify the notation we will assume that d is fixed; so, we erase the dependence on d. In fact the parameter d coincides with the degree of the Newton method applied to the polynomial p t at the unique active critical point.
The expression of the Newton's map applied to p t writes as It is easy to check that the zeros of the polynomial p t are finite fixed points of N t which are always attracting. Moreover the point at ∞ is a repelling fixed point of N t with multiplier λ 0 = (d + 1)/d (independently of t). We remark that given any holomorphic family of polynomials q t of degree d + 1 and their associated family of Newton's method N qt , then the point of ∞ is the unique repelling fixed point of N qt and its multiplier is always equal to λ 0 = (d + 1)/d (independently of q t ). The critical points of N t correspond to the zeros of the rational map So the critical points of N t are the zeros of p t and z = 0 since p t (z) = (d + 1)dz d−1 . The zeros of the polynomial are fixed points of the Newton map, so there are not active critical points. Therefore the origin is the unique active critical point of N t and the degree of N t at the origin is equal to d, so deg (N t , 0) = d. We observe that the marked critical point is c(t) = 0 for all the values of t; and we can interpret the parameter d as the degree of this active critical point. The corresponding critical value, N t (0), is equal to 1/t. Furthermore at the parameter t = 0 we have a Misiurewicz bifurcation of degree d (see Definition 4). In Figure 4 we show the parameter plane of the Newton map for different values of d.
The hyperbolic components in the t-parameter plane are the open subsets of C in which the unique active critical point c(t) = 0 either eventually map to one of the immediate basin of attraction corresponding to one of the roots of p t or it has its own hyperbolic dynamics associated to an attracting periodic point of period greater than one. Of course, the bifurcation locus corresponds to the union of all boundaries of those components and possible accumulating points, see Figure 4. For a detailed study of the dynamical system given by (8) we cite [CGJV14].
We choose this example for several reasons. Firstly, there is only one active critical point, c(t) = 0, and changing the value of the parameter d ≥ 2 the active critical point has local degree d. Thus, changing the value of the parameter d we find out generalized Mandelbrot sets M d densely appearing in the bifurcation locus of of N t (see Figure 4). Moreover, N t has a Misiurewicz bifurcation at t = 0 of degree d (see Figure 2). Secondly, and more important, we claim that the family N t can be considered as a model for the study of Misiurewicz bifurcations for Newton's method applied to arbitrary families of polynomials. To see the claim, suppose that q t is a family of polynomials of degree d + 1 and we denote by N qt the Newton's method associated to q t . Assume, without lost of generality, that t = 0 is a Misiurewicz bifurcation of some degree. As we already explain the unique repelling fixed point of N qt is located at ∞   and its multiplier, λ 0 = (d + 1)/d, coincides with the multiplier of the family N t . But from Theorem 1.2(b-c) this universal multiplier controls the location and the size of the cascade of the generalized Mandelbrot sets tending to t = 0 that appears in the parameter plane. This phenomenon explain why the parameter plane of two Newton's method applied to different families of polynomials of the same degree look like so similar. Going back to N t , remember that it has a Misiurewicz bifurcation at t = 0 of degree d. So we have in the parameter plane of N t a cascade of generalized Mandelbrot sets, denoted by M n d where n denotes the period of the attracting cycle in the main cardioid, tending to the Misiurewicz bifurcation parameter, t = 0. We denote t n and L n the location and size of the baby Mandelbrot set M n d of period n. Applying the algorithm developed in section 2 we can compute numerically the centers and the size of these baby Mandelbrot sets. In Table 3.1 we show these values for n ≤ 25 when d = 2 and d = 3. In Figure 2 we also plot the first terms of this sequence for d = 2.  Table 1. Location t n and size L n of the Mandelbrot sets M n 2 tending to t = 0 for the family N t with d = 2 and d = 3. See Figure 2.
Another way to visualise the output of these numerical computations is to use Theorem 1.2 and (7). From there we have The above expressions can be expressed as a linear regressions between the logarithm of the location and size of M n d and n, respectively. We plot in Figure 5 the results of the numerical computation of t n and L n for different values of the parameter d. In both cases, subfigures (a) and (b), the x-axis correspond to the period n of the small Mandelbrot set M n d , while the y−axis corresponds to log |t n | and log |L n |, respectively.
3.2. Example 2. An entire transcendental family. Due to the existence of an essential singularity at infinity, the global dynamics of transcendental entire maps is quite different to the dynamics of rational maps. For instance new types of Fatou components may exist like Baker or wandering domains (see [Dev94] for a brief introduction to transcendental dynamics).  However, according to Remark 1, since the theory of Misiurewicz bifurcations is of local character, it can be easily extended to families of transcendental entire maps. The main goal is to show that the above algorithm runs with no significant modifications. We consider the family of entire transcendental maps given by where t = 0 is a complex parameter. For all the values of the parameter t we have a superattracting fixed point at the origin, i.e., f t (0) = f t (0) = 0, which is also an asymptotic value (a singularity of the inverse map, see [Dev94], pages 181-182). Computing the derivative we have that f t (z) = tz(z + 2) exp(z) and thus there also exists a unique (active) critical point located at c(t) ≡ −2. The local degree of both critical points is 2 since f t (0) = 0 and f t (−2) = 0.
The Fatou set of f t always contains the basin of attraction of 0, that is the points in C that under iteration tends to 0. Although ∞ is an essential singularity of the map f t and so, the global dynamics is far for being polynomial or rational, the existence of a critical point at −2 of degree 2 makes some similarities between the family f t and the quadratic family. Such similarities are explained when plotting the parameter plane of f t where we can see some copies of small Mandelbrot sets. In this parameter plane we plot three different behaviours, the first one is the set of parameters such that the orbit of the active critical point −2 tends to zero (blue), the second one when the orbit escapes to ∞ (red) and finally when the orbit accumulates in an attracting cycle (black). See Figure 6.
In order to study a Misiurewicz bifurcation we need to find parameter values such that the active critical point lands in a repelling fixed point. Thus, we first compute de critical value v(t) = f t (−2) = 4te −2 , and imposing that v(t) is a fixed point f t (v(t)) = v(t) we arrive to the equation Instead of solve numerically this equation we can take advantage to write it down in term of the function g(z) = z exp(z) and its (multi-valued) function W (z) known as Lambert function.  concluding that g t has a Misiurewicz bifurcation of degree 2 at the parameter t = 0. Applying the algorithm developed in Section 2 we have computed the location t n and size L n of the cascade of Mandelbrot sets M n 2 tending to t = 0. To illustrate these computations we use, as before, the regression between the period n of the small Mandelbrot set M n 2 with log |t n | and log |L n |, in Figure 7(a) or Figure 7(b), respectively.
We end up this section by considering a slight modification of the Misiurewicz bifurcation described above (where, remember, the active critical point c(t) lands in one iterate at the repelling fixed point). Following [McM00], we introduce a new parameter k controlling the number of iterates that the active critical point needs to land at the repelling fixed point We solve numerically the above equation obtaining the parameterτ ≈ −19.63162088. For this concrete parameter value, the critical point −2 lands into two iterates to the repelling fixed point f 2 τ (−2) ≈ −5.37511380 with multiplier λ 0 ≈ 1.9462489. In Figure 8 we show in the parameter plane of the map f t the location of the parameterτ and the sequence of Mandelbrot sets M n 2 tending to this generalized Misiurewicz bifurcation.  We observe that by construction the sequence of copies of Mandelbrot set M n 2 have n ≥ 3, since we need two iterates to be close to the repelling fixed and at least another iterate to return to the critical point. In general, when the critical point takes k iterates to land into a repelling fixed point, the small copies of Mandelbrot set have period n ≥ k + 1. Applying again the algorithm developed in the Section 2 we obtain the locationt n and sizeL n of M n 2 of degree n. Finally, in Figure 9 we show the linear regression between the period and the logarithm of the location and size of M n 2 as expected.  Figure 9. Linear regression between log (t n ) and n (left) and log (L n ) and n (right) for the entire map f t .