Lower bounds on the Hausdorff dimension of some Julia sets

We present an algorithm for a rigorous computation of lower bounds on the Hausdorff dimensions of Julia sets for a wide class of holomorphic maps. We apply this algorithm to obtain lower bounds on the Hausdorff dimension of the Julia sets of some infinitely renormalizable real quadratic polynomials, including the Feigenbaum polynomial $p_{\,\mathrm{Feig}}(z)=z^2+c_{\,\mathrm{Feig}}$. In addition to that, we construct a piecewise constant function on $[-2,2]$ that provides rigorous lower bounds for the Hausdorff dimension of the Julia sets of all quadratic polynomials $p_c(z) = z^2+c$ with $c \in [-2,2]$. Finally, we verify the conjecture of Ludwik Jaksztas and Michel Zinsmeister that the Hausdorff dimension of the Julia set of a quadratic polynomial $p_c(z)=z^2+c$, is a $C^1$-smooth function of the real parameter $c$ on the interval $c\in(c_{\,\mathrm{Feig}},-3/4)$.


Introduction
The Julia set of a rational map on the Riemann sphere can be defined as the set of all its Lyapunov unstable points. Informally speaking, this is the set of points near which the system behaves chaotically, i.e. has unpredictable long term behavior. A natural question is how big and complex this set can be. One of the classical measures of geometric complexity of a Julia set (or any other fractal set) is its Hausdorff dimension.
The study of the Hausdorff dimension of Julia sets of rational maps in general, and of quadratic polynomials in particular, has a long history. It was proven that for large classes of rational maps, including the hyperbolic and parabolic maps, the Julia set has Hausdorff dimension strictly less than two (e.g., see [Urb94,McM00,GS98,GS09]). At the same time there exist rational maps whose Julia sets are proper subsets of the Riemann sphere and have Hausdorff dimension equal to 2. For instance, M. Shishikura showed that such maps are generic in the family of all quadratic polynomials p c (z) = z 2 + c with the parameter c on the boundary of the Mandelbrot set [Shi98].
In what follows, for every c ∈ C, we let p c : C → C denote the quadratic polynomial p c (z) = z 2 + c.
The Julia set of p c will be denoted by J c , and dim H (J c ) will stand for its Hausdorff dimension.
The first examples of quadratic polynomials p c with the Julia set of positive Lebesgue measure (and hence, also with dim H (J c ) = 2) were constructed by X. Buff and A. Chéritat [BC12]. Other interesting examples of quadratic polynomials were constructed by A. Avila and M. Lyubich in the sequence of papers [AL06, AL08,AL22]. Namely, they showed that there exist Feigenbaum maps (i.e., infinitely renormalizable quadratic polynomials p c with bounded combinatorics and a priori bounds) of the following three types: (i) dim H (J c ) < 2; (ii) dim H (J c ) = 2 and the Lebesgue measure of J c is zero; (iii) the Lebesgue measure of J c is strictly positive.
We note that even though, according to M. Shishikura [Shi98], there are plenty of complex parameters c ∈ C, for which the Julia set J c has the highest possible Hausdorff dimension (i.e., dim H (J c ) = 2), it is not known, whether any of these parameters can be real. In fact, not much is known about how large the dimension dim H (J c ) can be when the parameter c is restricted just to the real line: it was shown in [LZ13] that there exist c ∈ R, for which dim H (J c ) > 4/3, and this seems to be the highest lower bound that has been rigorously confirmed. Our first result improves this lower bound.
We recall that for a rational map f :Ĉ →Ĉ, the hyperbolic dimension dim hyp (J(f )) of the Julia set J(f ) is defined as the supremum of the Hausdorff dimensions of all forward invariant hyperbolic subsets of J(f ). It follows immediately from this definition that dim hyp (J(f )) ≤ dim H (J(f )).
Let c Feig ≈ −1.401155189 be the Feigenbaum parameter, i.e., the limit of the period doubling bifurcations along the real axis, starting from the main cardioid of the Mandelbrot set. We prove the following: Theorem 1.1. The hyperbolic dimension of the Julia set for the Feigenbaum polynomial p c Feig satisfies dim hyp (J c Feig ) > 1.49781.
We note that the best known upper bound on the Hausdorff dimension of the Julia set J c Feig is dim H (J c Feig ) < 2 (see [DS20]). Together with Theorem C from [AL08] this implies that dim hyp (J c Feig ) = dim H (J c Feig ).
The proof of Theorem 1.1 involves computer assistance. All computations are done with explicit estimates of the accumulated errors, so the obtained lower bound is rigorous.
Using the exact same machinery as in the proof of Theorem 1.1, we obtain lower bounds on dim hyp (J c ) for other parameters c. We illustrate this fact in Figure 1(a), where we plot lower bounds for 1000 evenly spaced parameters −2 ≤ c ≤ 2. Some parameters, corresponding to infinitely renormalizable real quadratic polynomials with stationary combinatorics, are of particular dynamical relevance; we single them out in the following theorem and compute the corresponding lower bounds with higher accuracy. Notice that any such parameters c is uniquely determined by a unimodal permutation s on the set {0, 1, . . . , n − 1}, where n is the smallest number, for which  The quantitative results of Theorem 1.1 and Theorem 1.2 are illustrated as red dots in Figure 1(b). Conjecturally, all Julia sets from Theorem 1.2 have Hausdorff dimension strictly less than 2 (see [Dud20]), therefore, their Hausdorff and hyperbolic dimensions are expected to coincide.
Finally, we apply our computational machinery to verify that for any real parabolic parameter c ∈ (c Feig , −3/4) we have dim H (J c ) > 4/3. Combining this with a very recent result of Ludwik Jaksztas and Michel Zinsmeister (see [JZ20]), we obtain the following theorem, conjectured in [JZ20]: Refining and extending the computations required to prove Theorem 1.3, we can construct a piecewise constant function whose graph forms a lower bound on the hyperbolic dimension over the entire parameter domain [−2, 2]. An example of this is illustrated in Figure 2(a), which constitutes a graphic representation of over 21000 small theorems (one lower bound for each subinterval). Note that the graph generally appears less smooth as the parameter c decreases. Indeed, the Hausdorff dimension is expected to behave in a rather irregular manner for parameters c ∈ R outside of hyperbolic components, as indicated in Figure 1 and Figure 2(b). (See also [DGM20] for the study of the Hausdorff dimension dim H (J c ) when c ∈ R is close to −2.)  The paper is structured as follows: in the first part of the paper (Sections 2 and 3) we present a general method for computation of lower bounds on the hyperbolic dimension of the Julia sets (see Theorem 2.6). Our method is a rather straightforward modification of C. McMullen's eigenvalue algorithm for expanding maps [McM98]. However, we show that our technique is applicable to a wider class of holomorphic maps, including non-expanding ones. We also expect that the same method can be used for effective computation of lower bounds on the dimension of the limit sets of some Kleinian groups.
In the second part of this paper (Sections 4 and 5) we apply our techniques to obtain rigorous lower bounds on the hyperbolic (and hence, also Hausdorff) dimension of the Julia sets that appear in Theorems 1.1, 1.2 and 1.3, as well as in the discussion above. The involved computations are done using computer-assisted means with explicit estimates of the accumulated errors.
Last but not least, let us mention that another algorithm for estimating the Hausdorff dimension of the Julia sets was constructed by O. Jenkinson and M. Pollicott [JP02]. However, this algorithm, as well as the original algorithm of C. McMullen [McM98], is designed for expanding maps, hence, is not directly applicable for establishing the results of the current paper. Estimates on the Hausdorff dimension of various other dynamically defined objects were obtained by many other authors (see for example [PV22], and the references therein).
2. Lower bounds on the Hausdorff dimension 2.1. Partitioned holomorphic dynamical systems. We will consider a rather wide class of holomorphic dynamical systems, possessing certain kinds of Markov and transitivity properties.
Definition 2.1. A partitioned holomorphic dynamical system F onĈ is a finite collection of holomorphic maps (1) f j : U j →Ĉ, where j = 1, . . . , m, such that (1) the domains U j are pairwise disjoint open sets inĈ, and each f j is a proper branched (or unbranched) covering map of U j onto its image; (2) (Markov property) for every pair of indices i, j, the image f i (U i ) either contains U j , or is disjoint from it; (3) (transitivity) for each domain U j , there exists a finite sequence of maps and maps U surjectively onto the union ∪ m i=1 U i . The domains U j of the maps f j will be called the tiles of F.
We note that, according to this definition, every rational or polynomial-like map is a partitioned holomorphic dynamical system that consists of a single map f . Furthermore, the restriction of a rational map on a Yoccoz puzzle of some fixed level can also be viewed as a partitioned holomorphic dynamical system with the number of maps being equal to the number of puzzle pieces.
Our definition of a partitioned holomorphic dynamical system can also be viewed as a generalization of a Markov partition for a conformal dynamical system from [McM98]. There is also a strong resemblance with the complex box mappings, introduced in [KvS09].
For a partitioned holomorphic dynamical system F that consists of the maps f 1 , . . . , f m as in (1), let denote the union of the domains U j of the maps f j and respectively the union of the images f j (U j ). We say that U and F(U ) are the domain and respectively, the range of F. Note that due to the transitivity property in Definition 2.1, we have the inclusion U ⊂ F(U ).
Since the domains U 1 , . . . , U m are pairwise disjoint, we can view the partitioned holomorphic dynamical system F as a single map Let C(F) be the set of all critical values of F. That is, Every point z ∈ U has a finite or infinite forward orbit O(z) = {z, F(z), F •2 (z), . . . } under the dynamics of F. The postcritical set of F, is the closure of the union of all critical orbits of F.
Remark 2.2. Assumption on F: In the remaining part of the paper we will always assume that F : U → F(U ) is a partitioned holomorphic dynamical system consisting of maps (1) defined on the tiles U 1 , . . . , U m , and such that F(U ) \ P C(F) is a nonempty open set. Due to the transitivity property of partitioned holomorphic dynamical systems, it follows that U j \ P C(F) is also a nonempty open set, for each j = 1, . . . , m. Furthermore, for simplicity of exposition, we will also assume that the number of tiles U j is at least two.
2.2. Geometric tree pressure and Hausdorff dimension. For every positive integer k ∈ N, a point z ∈ U \ P C(F) and a real number t > 0, consider the sum where the summation is taken over all w ∈ U , such that F •k (w) = z. The derivatives here are considered with respect to the spherical metric onĈ. The geometric tree pressure P (z, t, F) is then defined as For rational maps of the Riemann sphere, the geometric tree pressure was studied in [Prz99]. Various equivalent definitions were also considered in [PRLS04] and [Prz21]. It is not difficult to show that if for every j = 1, . . . , m, the set U j \ P C(F) is path connected, then the pressure P (z, t, F) does not depend on the point z ∈ U \ P C(F). Indeed, under these conditions, by the Koebe Distortion Theorem (see Lemmas A.1 and A.2), the pressure P (z, t, F) does not depend on the choice of z within the same domain U j . Then, the transitivity property of Definition 2.1 implies that P (z, t, F) is the same for any choice of z in U \ P C(F).
In case if F is a rational map on the Riemann sphere (without any assumptions on the geometry of the postcritical set), it was shown in [Prz99] that P (z, t, F) is independent of the choice of z ∈ C \ P C(F). We do not attempt here to prove a similar statement for all partitioned holomorphic dynamical systems F. Instead, it is sufficient for us to consider the pressure function and define the Poincaré exponent t F as follows: Lemma 2.3. The inequality P (2, F) ≤ 0 holds. Hence, the set {t > 0 : P (t, F) ≤ 0} is nonempty and the Poincaré exponent t F is well defined.
Proof. Let z ∈ U \ P C(F) be an arbitrary point, and choose a sufficiently small spherical disk D ⊂ U \P C(F), centered at z. Then, according to the Koebe Distortion Theorem (see Lemmas A.1 and A.2), for any k ∈ N and any w ∈ F −k (z), the term |(F •k ) ′ (w)| −2 is commensurable with the spherical area of the connected component of F −k (D), containing w. Summing up over all such connected components, we conclude that Θ k (z, 2, F) is uniformly bounded in k. Thus, P (z, 2, F) ≤ 0.
Remark 2.4. If ∞ ∈ U and U is a bounded subset of C, then all computations in this, as well as in the following sections can be done in the Euclidean metric instead of the spherical one. Since these two metrics are equivalent on U , this does not change the pressure function P (t, F). All further proofs repeat verbatim.
Let us mention an important relation between the Poincaré exponent t F and the dimension of the Julia set of F, in case if F is a rational map. More precisely, if f :Ĉ →Ĉ is a rational map of degree d ≥ 2, let J(f ) ⊂Ĉ denote its Julia set. We recall that the hyperbolic dimension dim hyp (J(f )) of the Julia set J(f ) is defined as the supremum of the Hausdorff dimensions of all forward invariant hyperbolic subsets of J(f ). It follows immediately from this definition that The following fundamental result is a generalized version of Bowen's formula (see [Bow79,Prz99]).
Theorem 2.5. Let F : U → F(U ) be a partitioned holomorphic dynamical system that is a restriction of a rational map f :Ĉ →Ĉ of degree d ≥ 2. Assume that F viewed as a single map has the same topological degree d as the rational map f . Then We note that even though hyperbolic and Hausdorff dimensions of a Julia set are not necessarily equal in general, they are known to be equal for certain classes of rational maps f including hyperbolic ones, Collet-Eckmann rational maps [Prz99], and Feigenbaum maps whose Julia sets have zero Lebesgue measure [AL08].
2.3. Construction of a lower bound on t F . A key difficulty one can face when working with a Poincaré exponent t F is that they are often difficult to compute. In this subsection we provide a quantity δ F that is usually much more tangible for direct rigorous computations than the Poincaré exponent itself. In Section 3 we prove that δ F t F . Afterwards, in Section 4 we discuss application of the latter result to the real quadratic family. Finally, in Section 5 we provide details on computer assisted proofs of Theorems 1.1, 1.2 and 1.3.
Assume, F is a partitioned holomorphic dynamical system, consisting of the maps (1). For each index j = 1, . . . , m, let d j be the topological degree of the (branched) covering map f j . Now, for any t ≥ 0, we introduce the so called, McMullen's matrix M (t) = M (F, t) whose (i, j)-th element m ij is defined as where the supremum above is taken over all z ∈ U i , such that f i (z) ∈ U j , and the derivatives are computed in the spherical metric (see also Remark 2.4).
For a square matrix A, let ̺(A) denote the spectral radius of A. That is, ̺(A) is the maximum of the absolute values of the eigenvalues of A. For a partitioned holomorphic dynamical system F, we define Finally, we can state the lower bound on the Poincaré exponent t F : Theorem 2.6. For any partitioned holomorphic dynamical system F, satisfying the assumptions from Remark 2.2, the number δ F is well defined, and δ F ≤ t F .
We finish this subsection with several practical remarks.
Remark 2.7. First, we note that any McMullen's matrix M always has nonnegative entries. Furthermore, due to the transitivity condition in Definition 2.1, any McMullen's matrix M is primitive. That is, there exists k ∈ N, such that the k-th power, M k has only positive entries. The spectral radius of a primitive matrix can be rigorously estimated both from above and from below by a version of the Perron-Frobenius Theorem, obtained by Collatz in [Col42] (see Theorem 5.6).
Remark 2.8. The following general fact turns out to be helpful in obtaining tight rigorous bounds on δ F : Then the spectral radius ̺(M (δ)) is a concave up function on R.
The proof of Lemma 2.9 is provided in the Appendix.
Remark 2.10. Finally, the lower bound δ F on t F can be improved by considering a refinement of F.
Definition 2.11. Let F be a partitioned holomorphic dynamical system, consisting of the maps (1). Another partitioned holomorphic dynamical systemF :Û →F (Û ), consisting of the mapŝ is a refinement of F, if each tileÛ j is contained in some tile U i , the corresponding mapf j is the restriction of f i toÛ j , and for each point z ∈F (Û ), the setsF −1 (z) and F −1 (z) coincide.
It is obvious from this definition that ifF is a refinement of F, then t F = tF . On the other hand, we show in Remark 3.2 that δ F ≤ δF . Hence, one can hope to get a better lower bound on the Poincaré exponent t F by applying Theorem 2.6 to a refinement of F, rather than to the system F itself. We implement this approach for real quadratic polynomials in Section 5.

Tree pressure over pseudo-orbits
In this section we give a proof of Theorem 2.6. Assume, k ∈ N and z ∈ U are such that z i = F i (z) is defined for all i = 0, 1, . . . , k and z k ∈ U . Let j 0 , j 1 , . . . , j k ∈ [1, m] be the sequence of indices, such that for each i = 0, . . . , k, we have z i ∈ U j i . We will call this sequence, the itinerary of z of length k + 1. Then for each i = 0, . . . , k − 1, the inclusion z i ∈ U j i ∩ F −1 (U j i+1 ) holds, and one can define the expression It follows immediately from the Chain Rule that the following inequality holds: For every z ∈ U \ P C(F) and t > 0, consider the expression where the summation is taken over all w ∈ U , such that F •k (w) = z. It follows immediately from (4) that the inequality holds for every z ∈ U \ P C(F), k ≥ 1 and t > 0. One can also observe that the expressions θ k (z, t, F) are independent from z within the same tile U j . Indeed, according to the Markov property of partitioned holomorphic dynamical systems, if z,z ∈ U j \ P C(F), then for all k ∈ N there is a bijection between the finite sets F −k (z) and F −k (z), where the bijection is taking each point w of F −k (z) into a pointw of F −k (z) with the same itinerary of length k + 1. The latter implies that Q(w, k) = Q(w, k), and the claim follows. For any z ∈ U \ P C(F) and t > 0, consider the following version of the pressure function: In particular, (5) implies that The next lemma states the relation between the values of θ k (z, t, F), p(z, t, F) and McMullen's matrices. For each index j = 1, . . . , m, let e j ∈ R m be the column vector whose j-th coordinate is 1 and all other coordinates are 0. Let also 1 ∈ R m be the column vector, whose coordinates are all equal to 1.
Lemma 3.1. Let t > 0 be a real number and let M = M (F, t) be the corresponding McMullen's matrix. Then, for any pair of integers k ∈ N, j ∈ [1, m] and any z ∈ U j \ P C(F), the following holds: θ k (z, t, F) = 1 T M k e j . As a consequence, we have p(z, t, F) = log ̺(M (F, t)). In particular, p(z, t, F) is independent of z ∈ U \ P C(F).

Proof. It follows immediately by induction on
where the summation is taken over all w ∈ U i , such that F •k (w) = z. (If there are no such points w, then the sum is considered to be equal to zero.) Then θ k (z, t, F) is equal to 1 T M k e j , which is the sum of all elements in the j-th column of M k . Since the matrix M = M (F, t) is primitive (see Remark 2.7), the second part of the lemma follows from the version of the Perron-Frobenius Theorem, which states that M k /(̺(M )) k converges to a matrix, whose columns are positive multiples of the leading eigenvector.
Since the function t → log ̺(M (F, t)) is continuous, the above inequalities imply that this function attains a zero on the interval [0, 2], hence, the set {t ≥ 0 : ̺(M (F, t)) = 1} from (2) is nonempty and the number δ F is well defined. The inequality δ F ≤ t F follows from (7).
We conclude this section with several observations.
Remark 3.2. IfF :Û →F(Û ) is a refinement of F, then for any z ∈Û \ P C(F) and t > 0, we have since the suprema in (3) for a refinementF are taken over smaller domains compared to the case of F. Thus, ̺(M (F, t)) ≤ ̺(M (F , t)) and in particular, Remark 3.3. Let F k be a sequence of partitioned holomorphic dynamical systems such that F k+1 is a refinement of F k for every k ∈ N. Then for any z ∈Û \ P C(F) and t > 0 the sequence p(z, t, F k ) is increasing and therefore converges. Assume that F 1 (and hence, all other F k ) is a restriction of a rational function f and the topological degree of each F k : U k → F k (U k ) coincides with deg f . Assume that diam(F k ), defined as the maximum of the diameters of all tiles of F k , converges to zero when k → ∞. In [Prz21] F. Przytycki showed that under the above conditions, lim p(z, t, F k ) = P (t, f ), where P (t, f ) is the geometric tree pressure of a rational function f . Taking into account that the sequence p(z, t, F k ) is increasing in k, we obtain that δ F k → t f = dim hyp (J(f )), as k → ∞.

Dynamical partitions for real quadratic polynomials
Starting from this section we restrict our consideration to the real quadratic polynomials p c (z) = z 2 + c, where c ∈ [−2, 2].
Let D r (z) ⊂ C denote the open disc of radius r ≥ 0 centered at z ∈ C, and let U 1 , U 2 ⊂ C be the two connected components (half-disks) of the set D 2 (0)\[−2, 2]. For any fixed c ∈ [−2, 2], consider the maps f 1 : U 1 → C and f 2 : U 2 → C defined as the restrictions of p c to the corresponding domains U 1 and U 2 . It is easy to check that the maps f 1 and f 2 together form a partitioned holomorphic dynamical system F c as in Definition 2.1. (Note that the latter is not true for non-real complex parameters c due to failure of the Markov property in Definition 2.1.) A sequence of refinements of F c can be constructed in a standard way: a refinement F (k) c of level k = 2, 3, . . . consists of 2 k maps, each of them being a restriction of the polynomial p c to a corresponding connected component of the set p 1−k c (U 1 ∪ U 2 ). We refer to these connected components as tiles of level k and denote them by P (k) j (c), where j = 1, 2, . . . , 2 k . The tiles are enumerated so that each tile P (k) j (c) is squeezed between the external rays of angles (j − 1)/2 k and j/2 k . For each k as above, the set of all tiles of level k is denoted by P (k) (c). We will often suppress the parameter c in the notation, writing P (k) j and P (k) in place of P (k) j (c) and P (k) (c) respectively, provided that this does not cause any ambiguity.
It follows from Theorem 2.5 that for any k, the Poincaré exponents t F (k) c coincides with the hyperbolic dimension of the Julia set J c , hence, due to Theorem 2.6, we have for any k ≥ 2. According to Remark 3.2, the difference dim hyp (J c ) − δ F (k) c decreases monotonically as k increases, so our strategy for obtaining lower bounds on dim hyp (J c ) consists of estimating δ F (k) c from below for sufficiently large values of k.
The McMullen matrices, corresponding to the refinements F (k) c will be denoted by c (1). We will often suppress the indices and write M (k) (t) or M (t) instead of M We note that due to local connectivity of the Julia sets for real quadratic polynomials, it follows that diam(F (k) c ) → 0 as k → ∞. Hence, the result of F. Przytycki [Prz21], discussed in Remark 3.3, implies that lim for all c ∈ [−2, 2]. Unfortunately, we don't know a good bound on the rate of convergence of δ F (k) c , when the critical point is recurrent (e.g., when p c is an infinitely renormalizable quadratic polynomial).

The computer-assisted proof
In this section we will outline the general strategy employed for the computer-assisted part of the proof. We will also present the key steps of the proof, together with some details that are specific to the problem at hand. We begin by recalling the main results, summarizing our conclusions.
We remark that the Feigenbaum parameter c Feig satisfies the bound |c Feig − c ⋆ | ≤ 10 −10 , so the theorem applies to the Feigenbaum map p Feig (z) = z 2 + c Feig . Thus, Theorem 5.1 immediately implies Theorem 1.1. That is, the Hausdorff and hyperbolic dimensions of the Julia set for the Feigenbaum map satisfy (The equality dim hyp (J c Feig ) = dim H (J c Feig ) follows from [AL08] and [DS20].) By using our computational machinery developed to prove Theorem 5.1, we can study other parameters in the same manner. Interesting candidates are the periodic Feigenbaum parameters of which we explore a few below: Theorem 5.2. Given a parameter c ⋆ from Table 1, for any parameter c ∈ R satisfying |c − c ⋆ | ≤ 10 −10 , the hyperbolic dimension of its Julia set dim hyp (J c ) satisfies the corresponding tabulated lower bound.
This immediately implies Theorem 1.2, where we also have added some information regarding the dynamics of the critical orbits. The computational scheme for rigorous numerical approximations of the parameters c ⋆ from the table is discussed in the Appendix.
Feeding even wider sets of parameters to our code, we can compute a uniform lower bound on the Hausdorff dimension of the associated Julia sets:   From this point of view, k is the length of the bit strings used to encode the 2 k tiles at level k. As a consequence, the matrix M (k) has dimension 2 k × 2 k .
Simply put, each non-zero element m (k) i,j of M (k) relates to the maximal distance each tile of the covering has to the origin of the complex plane. To be more precise, given the collection of tiles P Following Remark 2.4, the derivatives p ′ c (z) in (8) are measured in the Euclidean metric. Note that |p ′ c (z)| = 2|z| is independent of the parameter c. On a technical note, by the nesting properties of the tiles, it follows that we have P for a certain index l = l(i, j). Therefore, when computing M (k) we actually use tiles of the next level k + 1.
Due to the underlying dynamics of real quadratic polynomials, a McMullen matrix has a fourfold symmetry as follows (we suppress k for now): Here the four submatrices all carry exactly the same information. Each of the two right-most matrices (labeled G ) is given by rotating the matrix G 180 degrees. The matrix G has size 2 k−1 × 2 k−1 ; it is sparse and double banded. Out of its 2 2(k−1) elements, only 2 k−1 are non-zero, which corresponds to one row's worth. In fact, using our encoding of tiles, only the elements of G indexed by (i, 2i−1) and (i, 2i) for i = 1, . . . , 2 k−2 carry any information; all other matrix elements are zero. For efficiency, the non-zero entries of G can be stored in a single 2 k−1 -dimensional vector; the elements are computed using (8).
Given a McMullen matrix M , together with a positive number t > 0, we define a new matrix M (t) whose elements are given by In other words, M (t) is the matrix formed by raising each non-zero element of M to the power of t. Thus M (t) has exactly the same sparse, banded structure as M , see Figure 3. which, by Theorem 2.5 and Theorem 2.6, we know is a lower bound for the Hausdorff (in fact, hyperbolic) dimension of the Julia set J c . We will continue suppressing the encoding depth k, writing M (δ) in place of M (k) (δ) and δ in place of δ k . Our aim is now to find good lower bounds for δ. We will use a very simple two-step strategy: First, we make an educated guess δ ⋆ < 2 for a lower bound for δ. Next, we prove that (12) 1 < ̺(M (δ ⋆ )).
Observe that according to Lemma 2.3 combined with (6) and Lemma 3.1, it follows that ̺(M (2)) ≤ 1. Together with concavity of the function t → ̺(M (k) (t)) (see Lemma 2.9) and inequality (12), this implies that indeed, For this strategy to work, there are three challenges to overcome: 1. We need some efficient heuristics for making a good (educated) guess for δ ⋆ . 2. Given t > 0, we need to rigorously enclose the matrix elements of M (t). 3. Given t > 0, we need to rigorously bound (from below) the spectral radius of M (t). For now, we will assume that we can overcome challenge 2. This will be explained in detail in Section 5.4.
In order to compute and bound the spectral radius, we will discuss some important classes of matrices.
Definition 5.4. A matrix A is positive if all its elements are positive; we then write A > 0. Analogously, a matrix A is non-negative if all its elements are non-negative: we then write A ≥ 0.
Using this notation we can write e.g. A ≥ B, meaning that A − B ≥ 0.
Definition 5.5. A square, non-negative matrix A is primitive if there exists a positive integer n such that A n is positive.
A related, but somewhat weaker, property is that of being irreducible. All primitive matrices are irreducible. The latter class is important since the Perron-Frobenius theory for positive matrices generalizes to the class of irreducible matrices.
We point out that, by construction, all McMullen matrices M and M (t) are primitive (and thus irreducible).
By a classic result by Collatz [Col42] we have the following: Theorem 5.6. Let A be a non-negative, irreducible, n × n matrix, and let v (0) be an arbitrary positive n-dimensional vector. Defining we have the following sequence of enclosures This theorem is very useful: it encloses the spectral radius of A, and therefore we can use any point r ∈ [λ (j) , λ (j) ] as an approximation of the spectral radius. Furthermore, we have the bound on the maximal approximation error we are making by using r in place of ̺(A). When the matrix A is primitive, it is known that this error tends to zero as j increases; lim j→∞ (λ (j) − λ (j) ) = 0.
We are now prepared to overcome challenge 1: finding a good candidate for δ ⋆ -the lower bound for δ appearing in (11). We request that δ ⋆ satisfies two bounds: where ε is a small positive number. The smaller we take ε, the better our bound becomes. The leftmost inequality of (15), which is the same as inequality (12), ensures that δ ⋆ is a true lower bound (and not just an approximation of δ). Using Theorem 5.6, we can find a good δ ⋆ via a simple bisection scheme. A straight-forward computation indicates that ̺(M (2)) < 1 < ̺(M ( 1 10 )). Based on this, we start the bisection search on the interval ∆ 0 = [ 1 10 , 2]. After j steps we have an interval ∆ j of width less than 2 1−j , whose left endpoint is a valid lower bound δ ⋆ of δ. In our computations, we end the bisection scheme when we have reached ε ≤ 10 −10 in (15).

Set-valued computations.
For efficiency reasons, all computations carried out in the bisection scheme are performed using normal floating point arithmetic. As such they are not entirely reliable: rounding errors can accumulate, and we must therefore carefully verify that δ ⋆ indeed is a true lower bound of high accuracy. For the lower bounds reported in Theorem 5.1 and Theorem 5.2, we use an encoding depth of 28, which means that M has 2 29 = 536 870 912 non-zero elements. For Theorem 5.3 we use a mix of depths 17 and 18.
In order to certify the computations described in Section 5.2, we have to resort to set-valued numerics. We will work with interval-valued matrices, and all operations will use interval arithmetic with outward directed rounding. An elementary overview of the theory of interval analysis is given in [Moo66,Tuc11].
For now, we assume that we are given interval-valued McMullen matrices M (k) , whose nonzero elements are intervals: m The width of each such interval reflects the uncertainty in the underlying tile construction; this is described at depth in Section 5.4. Again, we will continue to suppress the encoding depth k for clarity.
Given an interval McMullen matrix M we want to find a good lower bound for the hyperbolic dimension of the Julia set J c : To this end, we will use the following result (see [Nus86]): Lemma 5.7. If A and B are irreducible matrices such that 0 ≤ A ≤ B, then ̺(A) ≤ ̺(B).
Recall that, by construction, each McMullen matrix is primitive (and therefore irreducible). This property carries over to the interval-valued versions M (t) too. Let M (t) = min(M (t)) denote the matrix whose elements are formed by taking the lower endpoint of each non-zero (interval) element of M (t). Then, by Lemma 5.7, it follows that We now validate our (educated) guess δ ⋆ by computing ̺(M (δ ⋆ )) using Theorem 5.6, but with all operations carried out in interval arithmetic. This will produce a rigorous enclosure of ̺(M (δ ⋆ )) which can be used to validate the lower bound δ ⋆ . More specifically, we check that ̺(M (δ ⋆ )) > 1, which implies (12), which in turn, implies that δ ⋆ is a true lower bound on δ.
5.4. Rigorous tile computations. Recall the notation from Section 4: let D r (z) denote the open disc of radius r ≥ 0 centered at z ∈ C. Note that D 2 (0) \ [−2, 2] consists of two half-discs. Given a quadratic map p c (z) = z 2 + c with c ∈ [−2, 2] and an integer k = 2, 3, . . . , we will let P (k) denote the set of connected components of p 1−k c (D 2 (0) \ [−2, 2]). The elements of P (k) are called tiles; there are exactly 2 k of them at level k. When c is real, the tiles have a four-fold symmetry: each quadrant is a reflection of the others. When c is complex, the symmetry is only two-fold (see Figure 4), and the restriction of the map p c to each of the tiles no longer provides a partitioned holomorphic dynamical system as the Markov property fails. Suppressing the level k for the moment, recall that, for each tile P i , we want to compute the maximal distance s i = max{|z| : z ∈ P i }. Since p ′ c (z) = 2z, we have the relation m i,j = 1/(2s i ), whenever P j ⊂ p c (P i ), see (8) of Section 5.1. In order to get rigorous enclosures of the elements of the McMullen matrices, we must be able to enclose the distances s i for all tiles at a fixed level. We do this by covering the boundary of each tile with small discs, starting with the two initial half-discs D 2 (0) \ [−2, 2]. From these initial half-discs, we generate the entire collection of tiles at level k is by repeatedly applying the inverse map p −1 c (w) = ± √ w − c. This of course, requires an inverse map that is disc-valued. We will describe how this is achieved in detail in Sections 5.5 and 5.6. For now, let us assume that we can cover the boundary of each tile P i by a finite collection of small discs d i,j , see Figure 5. Computing upper and lower bounds (s i and s i , respectively) for each maximal distance s i is then a simple matter of traversing the collection of discs, and finding the maxmin/maxmax distance to the origin: We end this section with a practical note. Each tile P (k) i ∈ P (k) is encoded via the binary representation of its index (i) 10 = (b k−1 . . . b 0 ) 2 , where i = 0, 1, . . . , 2 k − 1. The encoding is used to select the correct branch of the inverse image: if b j = 0, then the j:th inverse image belongs to the lower half plane. If b j = 1, we select the inverse image in the upper half plane. This produces a labeling of tiles starting at the right-most position of the fourth quadrant and moving clock-wise through the 3:rd, 2:nd, and 1:st quadrant. By the four-fold symmetry, it suffices to compute the first quarter of tiles. They all belong to the fourth quadrant, and are indexed i = 0, . . . , 2 k−2 − 1. 5.5. Circular complex arithmetic. There are several models for set-valued complex arithmetic: the fundamental choice lies in the the way a set of complex numbers is represented. A survey of different approaches can be found in [PP98]. We will opt for circular sets; this will minimize the overestimation when computing the inverse images of the quadratic map.
A set in the complex plane is thus represented by a disc with midpoint m ∈ C and radius r ≥ 0: Computing the inverse map p −1 c (w) = √ z − c involves two operations: subtraction and taking the square root. Subtracting a disk from another is straight-forward. We have (21) z 1 − z 2 = m 1 , r 1 − m 2 , r 2 = m 1 − m 2 , r 1 + r 2 , which incurs no overestimation: the resulting disc is a sharp enclosure of the difference. The square root, however, is more complicated as discs are not preserved under this operation. Given an input disc z we must find a new disc w that contains the set { √ z : z ∈ z}. We will only focus on one of the two branches of the square root (the second branch is simply the negative of the first). First, we write the midpoint of z = m, r in polar coordinates: m = ρe iθ . Assuming (for now) that 0 / ∈ z, we have 0 ≤ r < ρ. From the two positive numbers α 1 = √ ρ + r and α 2 = √ ρ − r, we formρ = α 1 +α 2 2 andr = α 1 −α 2 2 . We then have the enclosure of the square root which is the sharpest possible. When the ratio r/ρ is small, the overestimation is very small. As the ratio approaches one, however, there is a significant difference between the true image of the square root and its disc enclosure given by (22). This is illustrated in Figure 6. When the disc z = m, r contains the origin, the situation changes drastically. Throughout our computations, however, this will only occur in very special situations. We will handle these special cases by extending our circular complex arithmetic. 5.6. Extended circular complex arithmetic. Considering the computations needed to establish our results, we want to extend the circular arithmetic introduced in Section 5.5. In order to compute the tiles efficiently (see Section 5.4) we would like to be able to represent (and compute with) interval segments of the real and imaginary axes. Indeed, part of the boundary of many tiles are subsets of the real and imaginary axes, see Figure 4. To this end, we will extend the notion of a disc to include such subsets of the complex plane. We introduce the following notation: m, r x = {z ∈ C : |z − m| ≤ r and Im(z) = 0}, m, r y = {z ∈ C : |z − m| ≤ r and Re (z) = 0}.
These collapsed discs contain only real or imaginary numbers, respectively. We can now define precise sets of real parameters for which we will prove our theorems. As an example, for Theorem 5.1 and Theorem 5.2 , we will work on very small subsets of the real axis: c = c ⋆ , 10 −10 x . For Theorem 5.3 we use intervals having much larger radii. In each inverse iteration of the quadratic map, we will be subtracting c from a disc. This operation is straight-forward; given two complex numbers z 1 = x 1 + iy 1 and z 2 = x 2 + iy 2 , we define z 1 , r 1 y − z 2 , r 2 x = −x 2 + iy 1 , r 1 + r 2 , (25) z 1 , r 1 − z 2 , r 2 x = x 1 − x 2 + iy 1 , r 1 + r 2 .
There are other combinations of operands, but these three are the only ones we will use.
We will also need to compute the square root of discs of type m, r x . To simplify notation, leť x be a real number and consider the disc x = x, r x . There are three different situations that we must account for.
If 0 ≤x − r, then the disc x consists of non-negative numbers only. We can use the same formula as (22) with θ = 0, and making sure we stay in the correct class of discs: Ifx + r ≤ 0, then the disc x consists of non-positive numbers only. Once again, we can use the same formula as (22) with θ = π/2. This corresponds to a rotation from θ = π to π/2, which brings us to the second class of extended discs: Finally, if x contains zero in its interior, we split the disc into its non-negative and non-positive counterparts, and treat them according to (27) and (28), respectively. This concludes our derivation of disc-valued arithmetic. In the actual numerical computations (see Section 5.7), all floating point operations are performed with outward directed rounding. This ensures that all results are true enclosures and can be used as rigorous bounds. 5.7. Computer outputs. In this section, we will present the numerical results behind Theorem 5.1, Theorem 5.2 and Theorem 5.3. For the basic interval arithmetic, we use the CAPD library ([KMWZ20]), which is well established in the community of computer-assisted proofs in dynamical systems. The code we have developed, however, can easily be adapted to any other interval library. The main program is based on the ideas presented above, and consists of three main stages: 1. Compute the tiles and the McMullen matrix elements. 2. Compute an approximation of a lower bound for the hyperbolic dimension. 3. Validate the lower bound. In order to speed up the computations, we perform a large part on them in parallel. Indeed, each tile can be computed independently, and all matrix-vector products required in the implementation of Theorem 5.6 benefit from parallelization.
Starting with Theorem 5.1, we launch our program aiming at depth k = 28, with an initial covering of each half-disk made up of 199 full discs and one collapsed disc. The output is presented below. Summarizing, the entire run takes about 20 hours, using 12 threads. In forming the intervalvalued McMullen matrix M , we compute 2 26 = 67, 108, 864 tiles. This part of the run is the most demanding: 80% of the run-time is spent on this stage, despite the fact that this is where we gain most on parallelisation. For each tile, we compute rigorous bounds on the distance to the origin using the disc coverings, and then form the matrix entries as explained in Section 5.4. We note that the widest interval appearing in the McMullen matrix has diameter 0.09313.... While computing the tiles, we split every produced disc of type m, r y until the radius is less than 0.01. This improves the quality of the tile covering by reducing the size of certain discs. This, in turn, prevents the matrix elements from becoming too wide.
The second stage, based on Section 5.2 takes over three hours, and produces a good quality "guess" for a lower bound on the Hausdorff dimension. This is attained by a non-rigorous bisection scheme producing increasingly accurate estimates.
Finally, the third stage (again fully rigorous) validates the guess, and proves that the bound given in Theorem 5.1 is a true lower bound.
Continuing to Theorem 5.2, we do exactly the same as above, but on a set of different parameters. We will not display the output from all ten runs; one typical run will suffice. Here we use the same depth k = 28, but reduce the initial covering of each half-disk: we use 39 full discs and one collapsed disc. 1.33929890631116133 to 1.35502280450309742 (all greater than 4/3). The total computation time was 4 minutes and 15 seconds when run on 12 cores in parallel.
Finally, we redo the computations above with a bit more effort, and over a wider parameter domain. To be more specific: we split the domain [−2, 2] into 21000 subintervals (of varying radii), each on which we compute a rigorous lower bound for dim hyp (J c ). This produces a piecewise constant function acting as a lower bound on the hyperbolic dimension of the Julia sets for the full domain [−2, 2], and is illustrated in Figure 2(a). The underlying computations took three days on ten 12-core machines. Using only a subset of this information, we can (again) prove Theorem 5.3 by showing that the hyperbolic dimension of the Julia set satisfies dim hyp (J c ) > 4/3 over the entire interval [−1.53, −1.23], see Figure 7(a). Since this interval contains all real parabolic parameterŝ c ∈ (c Feig , −3/4), the theorem follows.
We end this section with an illustration of a rigorous cover of the Feigenbaum Julia set, see Figure 7(b). Here we collect several technical statements and their proofs.
Lemma A.1. Let U ⊂Ĉ be a topological disk and assume that f : U →Ĉ is a univalent map. If there exists r > 0, such that each of the setsĈ \ U andĈ \ f (U ) contains a spherical disk of radius r, then for any open subsetŨ ⋐ U , compactly contained in U , the distortion of f onŨ , measured in the spherical metric, is bounded by some constant K ≥ 1. This constant depends on the domainŨ and the real number r.
Proof. After precomposing and postcomposing f with appropriate isometries of the sphere, we may assume that both U and f (U ) are contained in a Euclidean disk of radius ρ(r), centered at zero. Here ρ(r) is some increasing function of r. The latter implies that the restrictions of the spherical and the Euclidean metrics on U ∪ f (U ) are equivalent with the equivalence constant depending on r. Now the statement of the lemma follows from the standard Koebe Distortion Theorem.
Lemma A.2. Let F : U → F(U ) be a partitioned holomorphic dynamical system that consists of at least two maps f j as in (1). Then there exists r > 0, such that for any open spherical disk D ⊂ U \ P C(F) ⊂Ĉ of radius r, a positive integer k ∈ N and any connected component D ′ of F −k (D), the setsĈ \ D andĈ \ D ′ contain disks of radius r.
Proof. Observe that each of the topological disks D and D ′ must be contained in one of the (possibly different) domains U j . If there are at least two such domains, then after choosing a sufficiently small r, the statement becomes obvious.
Proof of Lemma 2.9. Let M = {m i,j } 1 i,j n be a n×n primitive matrix and M (δ) = {m δ i,j } 1 i,j n , δ ∈ R. We need to show that the spectral radius ̺(M (δ)) is a concave up function on R.
Let m k (δ) be the entry in the first column and the first row of (M (δ)) k . From Perron-Frobenius Theorem one can deduce that (m k (δ)) 1 k → ̺(δ) when k → ∞. Observe that for every k one has m k (δ) = r j=1 c δ j , for some r ∈ Z + and c 1 , . . . c r > 0 depending on k. Let us show that (m k (δ)) 1 k is a concave up function for every k.
On has: d dδ (m k (δ)) 1 k is concave up. Taking the limit when k → ∞ we obtain that ̺(M (δ)) is concave up.
Rigorous approximations of real Feigenbaum parameters. Our computations rely on the theory of kneading sequences (see [MT88,CE80]). Given a unimodal map f (z) with the critical point at 0 the kneading sequence is the sequence K(f ) = {k 1 , k 2 , . . .} ⊂ {−1, 0, 1} N , where k i = −1 if f i (0) < 0, k i = 1 if f i (0) > 0 and k i = 0 if f i (0) = 0. Given a Feigenbaum map p c (z) = z 2 +c, c ∈ R, with a stationary combinatorics of period n the map p n c (z) is hybrid equivalent to p c (z) near the origin. In particular, there is a homeomorphism h from a segment I ⊂ R containing the postcritical set of f c to a segment I 1 containing zero such that p n c = h • p c • h −1 . From this we conclude that the kneading sequence K(p c ) has the properties (29) k n l m = (k n /k 1 ) l k m = (−k n ) l k m , k s+n = k s if n does not divide s.
This allows us to compute any number of elements of the kneading sequence K(p c ) given k 1 , . . . , k n . Notice that not all sequences k 1 , . . . , k n of length n are realizable as an initial segment of a kneading sequence. To find all realizable sequences we can start by finding all real superattracting parameters of period n. Given a real superattracting parameter a of period n, it follows from [McM96,Lyu99] that there exists a unique period n real Feigenbaum parameter c = c(a) such that k i (c) = k i (a) for i = 1, . . . , n − 1 and k n (c) = k n (a − ǫ) for a sufficiently small ǫ > 0. The value of k n (c) can also be determined as k n (c) = −sign dp n x (0) dx | x=a . Given a point b ∈ R, b = c, one can find out whether b < c by comparing the kneading sequences of p b and p c . Namely, let i be the smallest number such that k i (b) = k i (c). If sign(k i (b)) = sign dp i x (0) dx | x=b then b < c, otherwise b > c. Now, from the above considerations, one can compute an approximation of c by repeatedly cutting a segment containing c into halves, starting with the segment [−2, 0]. Using intervalvalued numerics one can control the accumulated errors and ensure that the approximations of the Feigenbaum parameters in Theorem 5.2 are obtained with appropriate accuracy.