On Benford's Law and the Coefficients of the Riemann Mapping Function for the Exterior of the Mandelbrot Set

We investigate Benford's law in relation to fractal geometry. Basic fractals, such as the Cantor set and Sierpinski triangle are obtained as the limit of iterative sets, and the unique measures of their components follow a geometric distribution, which is Benford in most bases. Building on this intuition, we aim to study this distribution in more complicated fractals. We examine the Laurent coefficients of a Riemann mapping and the Taylor coefficients of its reciprocal function from the exterior of the Mandelbrot set to the complement of the unit disk. These coefficients are 2-adic rational numbers, and through statistical testing, we demonstrate that the numerators and denominators are a good fit for Benford's law. We offer additional conjectures and observations about these coefficients. In particular, we highlight certain arithmetic subsequences related to the coefficients' denominators, provide an estimate for their slope, and describe efficient methods to compute them.


Introduction
The Mandelbrot set M was first introduced and drawn by Brooks and Matelski. By analyzing the family of functions f c (z) = z 2 + c, Douady and Hubbard began the formal mathematical study of the Mandelbrot set as the set of parameters c, for which the orbit of 0 under f c remains bounded. We study Benford's law in relation to the Mandelbrot set to both investigate the distribution's extension to fractal geometry and search for patterns in the Mandelbrot set.
In 1980, Douady and Hubbard were able to prove the connectedness of M by constructing a conformal isomorphism Φ : C \ M −→ C \ D between the complement of the Mandelbrot set and the complement of the closed unit disk [1]. Using the Douady-Hubbard map Φ, we can define related conformal isomorphisms, where M −1 = {1/c : c ∈ M}, by setting Ψ = Φ −1 and Θ(c) = 1/Ψ(1/c). One of the most heavily studied questions in complex dynamics is whether or not M is locally connected (MLC). By a theorem of Caratheodory [2], these two maps can be extended continuously to the unit circle if and only if the Mandelbrot set is locally connected. As such, we focus on studying these maps and their respective Laurent and Taylor expansions, as outlined in [3] and [4], respectively: In Section 3, we outline the methods [3][4][5] we used to compute the a m and b m coefficients. The computation time grows exponentially, so methods of improving computation are explored. Using recursion, we were able to compute the first 10240 coefficients.
In Section 4, we discuss Benford's law along with statistical testing to determine whether the coefficients obey a Benford distribution. Given a base b ≥ 2, a data set is Benford base b if the probability of observing a given leading digit, d, is log b (d + 1)/d (see [6,7]). We can write any positive x as S b (x)b k b (x) , where S b (x) ∈ [1, b) is the significand and k b (x) is an integer. If the probability of observing a significand of s ∈ [1, b) is log b (s), we say the set is strongly Benford (or frequently just Benford). In most cases, a data set is demonstrated to be Benford through statistical testing. There are few straightforward proofs for Benfordness, all of which rely heavily on understanding the structure and properties of the data. Well understood sets such as geometric series and the Fibonacci numbers have explicit proofs for Benfordness, but the structure and properties of coefficients we study are still the subject of active research in complex dynamics. Therefore, we rely on statistical testing for our results. We consider the standard χ 2 distribution and the sequence of the data's logarithms modulo 1 for our statistical testing, and a standard goodness of fit test demonstrates that the numerators and denominators are a good fit for Benford's law, while the decimal representations are not.
Section 5 deals with conjecture, observations, and theorems related to the coefficients. Sections 5.1 and 5.2 are meant to tie together the most important of these that we have found from various authors for the a m and b m coefficients, respectively. In Section 5.3, we present new results and conjectures on the a m and b m coefficients from our work. Theorem 7 [3,8] states that they are 2-adic rational numbers; in other words, they are of the form p/2 −ν , where p is an odd integer. The integer ν is, by definition, the 2-adic valuation ν 2 of a m or b m . Therefore, we focus on the denominator's exponents −ν(a m ), −ν(b m ). Setting m = 2 n m 0 , with m 0 odd and n fixed; the subsequences {−ν(a m )}, {−ν(b m )} appear to be near-arithmetic progressions. We present the results observed in the following conjecture.

Conjecture 1.
Let m be written as 2 n m 0 as above, with n = n fixed. Then, the sequence {−ν(a m )} n=n is asymptotically linear, with slope We also present an efficient way to compute the denominator's exponents for the cases n = 0, 1, 2.
Our work offers a new approach to some classical problems in complex dynamics, and we do this through our statistical testing. These coefficients have been studied extensively, so we compile relevant observations from disparate authors to use as a basis for offering new conjectures and results. The idea of studying Benford's law in complex dynamics and fractal sets is also original; to our knowledge, no other authors have attempted to study Benford's law in this setting. Our approach offers both a new and unique way to study important series in complex dynamics, and it provides motivation for number theorists and statisticians to study Benford's law in the new field of data, namely fractal sets.

Preliminaries in Complex Dynamics
We give a brief introduction to complex dynamics. For more detailed proofs see [9] and [10].
Let f :Ĉ →Ĉ be a rational map. The Julia set J f associated to the map f may be defined as the closure of the set of repelling periodic points of f . For a rational map of degree 2 or higher, the Julia set J f is non-empty. We now restrict our attention to polynomial maps of degree d ≥ 2, which have a superattracting fixed point at infinity. We can thus define the filled Julia set K f as the complement of the basin of the attraction of infinity: Making use of the above, it is possible to redefine J f as the boundary of the filled Julia set.
Lemma 1. Let f be a polynomial of degree d ≥ 2. The filled Julia set K f ⊂ C is compact, with boundary ∂K f = J f equal to the Julia set. The complementĈ \ K f is connected and equal to the basin of attraction A(∞) of the point ∞.
It follows that the Julia set J f of a polynomial f is precisely the boundary of the basin of attraction A f (∞).
We may now characterize the connectedness of J f . This is determined entirely by the activity of the critical points of f .

The Mandelbrot Set
We now focus primarily on the family of quadratic functions of the form { f c (z) = z 2 + c} c∈C . Since f c has a single critical point 0, it follows from Theorem 1 that J f c is connected if and only if the orbit f n c (0) | n ∈ N is bounded. This motivates the following definition of the Mandelbrot set, M. Definition 1. M ⊂ C is the set of all the parameters c ∈ C such that the Julia set J f c is connected. Equivalently, M is the set of all the c such that the orbit of 0 under f c remains bounded:

Remark 1.
It is possible to generalize this definition and most of the following results to the family of unicritical degree d polynomials f c,d (z) = z d + c, where d is an integer d ≥ 2. In this case, M d is called the multibrot set of degree d. For simplicity, we focus only on M = M 2 , which has historically been the object of greatest interest. For more information on M d , see [3,8].
It is possible to demonstrate that the interior of M is nonempty. We utilized an escape time algorithm and computer graphics to obtain the visualization of M presented in Figure 1. When the first computer images of M were generated, Benoit Mandelbrot observed small regions that appeared to be separate from the main cardioid, and conjectured that M was disconnected, which was later disproved.
This result was first proved by Douady and Hubbard [1] by explicitly constructing a conformal isomorphism Φ :Ĉ \ M →Ĉ \ D. Douady and Hubbard's proof is significant not only for the result, but also since it provides an explicit formula for the uniformization of the complement of the Mandelbrot set.
A large amount of research has been devoted to the local connectivity of the Mandelbrot set, which is generally regarded as one of the most important open problems in complex dynamics. We recall that a set A in a topological space X is locally connected at p ∈ A if for every open set V ⊂ X containing p, there is an open subset U with p ∈ U ⊂ V such that U ∩ X is connected. The set A is said to be locally connected if it is locally connected at p for all p ∈ A.
As above, let Φ :Ĉ \ M →Ĉ \ D be the conformal isomorphism constructed by Douady and Hubbard; notice that the map Ψ :Ĉ \ D →Ĉ \ M is the Riemann mapping function of C \ M. We consider it the Laurent expansion at ∞ Another possibility is to consider Θ(z) := 1/Ψ(1/z), which is the Riemann mapping of the bounded domain C \ {1/z : z ∈ M}. We have the corresponding Taylor expansion for Θ at the origin: For the general Multibrot set M d , we refer to the coefficients with the notation b d,m and a d,m . To underline the importance of these maps, we reference a lemma from Caratheodory [2]. Therefore, if the map Ψ or the map f can be extended continuously to the unit circle ∂D, then M is locally connected. To demonstrate this extension, it would be sufficient to prove that one of the two series converges uniformly on D.
There have been numerous attempts to prove this result. For example, Ewing and Schober demonstrated that the inequality 0 < |b m | < 1/m holds for every m < 240, 000. A bound of the type |b m | < K/m (1+ ) would lead to the desired result; however, this would imply that the extension of Ψ is Hölder continuous, which it is not, as expressed in [11]. Proving that |b m | < K/(m log 2 (m)) would prove that the series converges absolutely; however, modern computations suggest that such a bound does not exist [11]. Therefore, the MLC conjecture and its consequences remain an object of active study. Another consequence of MLC is related to a topologically equivalent description of ∂M. In particular, the boundary of the Mandelbrot set can be identified with the unit circle S 1 under a specific relation ∼, known as the abstract Mandelbrot set [12]. More information and other implications of MLC and the Density of Hyperbolicity may be found in [13] and [1].

Algorithms and Complexity
There are algorithms to compute both the a m and b m coeffecients. While these algorithms work for a generic degree d ≥ 2, we focus on d = 2, which is historically the most interesting case, since it is the one associated with M. For simplicity, we denote b m = b 2,m and a m = a 2,m . The behavior for other values is similar. Derivations for the explicit form of the b m and a m coefficients may be found in [3,4] respectively. There is also a formula to switch between the coefficients, outlined in [14].
where the sum is over all non-negative indices j 1 , . . . , j n such that and C j (a) is the generalized binomial coefficient The a m coefficients can be obtained from the b m using the formula or they can be directly calculated as in the following theorem.
where the sum is over all non-negative indices j 1 , . . . , j n such that While the above theorems give the explicit forms of the coefficients, the following theorem provides a recursive method to find b m , which is more suitable for computers. Once we find b m , we can apply the relationship between a m and b m outlined in theorem 4 to find a m . More details can be found in [5]. where the following holds true.
For example, to calculate the first several b m coefficients, we can use Theorem 6 to obtain:

Direct Computation
Our initial algorithm to generate these coefficients was to directly compute them, and we include our original methodology for reference as we cross checked our results. Most of the data was generated through the recursive algorithm, and details are provided in Section 3.2.
We wrote a program in Python to compute b 2,m and a 2,m based on the formulas given in Theorems 4 and 5, and we obtained the first 1024 exact values of both coefficients' sequences.
Our methodology for computing the m th coefficient was to first generate the solutions j 1 , . . . , j n to the Diophantine Equations (4) and (6). Then, we plug them into the exact formula of a m = a 2,m and b m = b 2,m to find the m th coefficient.
We improved the method to solve the Diophantine equations by first setting an upper limit on the degree for which to obtain coefficients, generated the solutions for the highest order coefficient, and created data structures for dynamic storage. We stored each individual solution as a tuple of length n, where the k th entry denoted the value of the coefficient j k . Every solution for the upper bound was then given a reference in a linked list, which we can use to find the highest order coefficient. The solution stored for the upper bound can then be modified through decrementing the value for j n for each tuple and then deleting the reference to the tuple and deallocating the memory in the linked list when the value for j n reaches zero.
To deal with the time sink in generating the binomial coefficients, we utilized multi-core parallel computing. Each coefficient can be computed independently after we obtain the solutions to the Diophantine equations. We structure our code for concurrent computation and use generator expressions so that we can use multiple cores where our code is executing simultaneously. In a multi-core setting, each core deals with one coefficient at a time. When one coefficient calculation is conducted, the core takes the next awaiting task that is not being taken by other cores. We chose a high-performance server machine and ran our code in a parallel environment with 72 valid cores. We obtained the first 1024 coefficients with a CPU time of 166 hours and a total run time of 7 days.

Recursive Computation
The direct computation runs in exponential time, and it is generally impractical for generating large degree coefficients. Therefore, we switched to a recursive method to generate these coefficients.
The method is described in [5,11,14], and we outline the formula for the computation in Theorem 6. This method is efficient because it is able to reuse information from the previous coefficients to compute the next one.
We wrote a Python program to implement the recursion to find b m and then use Equation (5) to find the corresponding a m . We were able to obtain 10240 coefficients for both series within 82 hours with a single core. We are also able to cross-check our computation results with the direct computation method before starting the statistical analysis. All codes and results can be found at https://github.com/DannyStoll1/polymath-fractal-geometry (accessed on 19 September 2022). Detailed instructions can be found in the README file.

Benford's Law
Frank Benford's 1938 paper, The Law of Anomalous Numbers [6], illustrated a profound result, in which the first digits of numbers in a given data set are not uniformly distributed in general. Benford applied statistical analysis to a variety of well-behaved but uncorrelated data sets, such as the areas of rivers, financial returns, and lists of physical constants; in an overwhelming amount of the data, 1 appeared as the leading digit around 30% of the time, and each higher digit was successively less likely [6,7]. He then outlined the derivation of a statistical distribution which maintained that the probability of observing a leading digit, d, for a given base, b, is log b (d + 1)/d for such data sets [6,7]. This logarithmic relation is referred to as Benford's law, and its resultant probability measure for base 10 is outlined in Figure 2. Benford's law has been the subject of intensive research over the past several decades, arising in numerous fields; see [7] for an introduction to the general theory and numerous applications.  Benford's law appears throughout purely mathematical constructions such as geometric series, recurrence relations, and geometric Brownian motion. Its ubiquity makes it one of the most interesting objects in modern mathematics, as it arises in many disciplines. Therefore, it is worthwhile to consider non-traditional data, such as fractals, where the distribution may appear. Basic fractals such as the Cantor set and Sierpinski triangle are obtained as the limits of iterations on sets, and their component measures (the lengths in the Cantor set and the areas in the Sierpinski triangle) follow a geometric distribution, which is Benford in most bases. Building on these results, it is plausible that more complicated fractals obey this distribution as well. We studied the Riemann mapping of the exterior of the Mandelbrot set to the complement of the unit disk, along with its reciprocal function to determine their fit to Benford's law. These mappings are given by a Taylor and Laurent series, respectively. The coefficients are of interest as their asymptotic convergence is intimately related to the conjectured local connectivity of the Mandelbrot set, which is an important open problem in complex dynamics.

Statistical Testing for Benford's Law
A common practice for evaluating whether a data set is distributed according to Benford's law is to utilize the standard χ 2 goodness of fit test. As we are investigating Benford's law in base 10, we utilize 8 degrees of freedom for our χ 2 testing. (There are nine possible first digits, but once we know the proportion that is digits 1 through 8 the percentage that starts with a 9 is forced, and thus we lose one degree of freedom). If there are N observations, letting p d = log 10 (d + 1)/d be the Benford probability of having a first digit of d, we expect p d N values to have a first digit of d. If we let O d be the observed number with a first digit of d, the χ 2 value is If the data are Benford, with 8 degrees of freedom, then 95% of the time the χ 2 test will produce a value of at most 15.5073; this corresponds to a significance level of α = 0.05. We perform multiple testing by creating a distribution of χ 2 values as a function of sample size up to the m th coefficient. This is standard practice for studying Benford sequences, and this is done to account for periodicity in the χ 2 values, which is typical for certain Benford data sets, such as the integer powers of 2. To account for this multiplicity, we also incorporate the standard Bonferroni correction. The overall testing is conducted at the level of significance of α = 0.05, while giving equal weight in terms of significance to each individual test by conducting them at a significance level of α/m. The rationale is to keep the significance of the overall test constant with respect to the number of tests performed. As we increase the total number of tests performed, we wish to increase our correction accordingly. In total, we perform 10,045 tests for the a m and 10,046 tests for the b m coefficients as we compute the χ 2 statistic each time we add a new non-zero coefficient to our data set. This brings our corrected threshold value to 38.9706 and 38.9708, respectively. This corresponds to α = 0.05/10,045 = 4.978 × 10 −6 for each individual hypothesis in the a m dataset and α = 0.05/10,046 = 4.977×10 −6 for each individual hypothesis in the b m dataset. Each data point is not independent, as the χ 2 values are computed by using a running total of the data, and as such, each point is built on the previous one. This results in a high correlation between the data, and the Bonferroni correction likely overcompensates for the increase in type I error. Still, it is one of the most plausible methods of dealing with the increase in multiplicity, since it is one of the simplest and most conservative estimates, and a value above the Bonferroni correction provides strong evidence that the data are not Benford. Using independent increments to compute each χ 2 statistic for Benford's law would fix the issue of independence, but is not recommended since periodic behavior can be missed if the increments are chosen poorly.
We considered the distribution of the χ 2 values to account for random fluctuations and periodic effects. In addition, we provide the provide the p-value of our computed χ 2 statistic to provide the type I error rate for our conclusions, and we compute the powers of the χ 2 tests relative to our null hypothesis that the data are Benford by using the noncentral chi-squared distribution [15]. We also conducted simulations to estimate the sampling error relative to our null hypothesis. We wish to see how likely it is that a random sample falls in the rejection region for our testing; based on our significance level of α = 0.05, we expect the sampling error to be roughly 5% if the data are Benford. We randomly sample 1000 coefficients from our data sets with replacement and calculate the χ 2 value for this sample data. We repeat this simulation 1000 times, and take the ratio the values in the rejection region to the total number of sample statistics calculated, to estimate the sampling error.
An equivalent test is to consider the distribution of the base 10 logarithms for the absolute value of the data set modulo 1; a necessary and sufficient condition for a Benford distribution is that this sequence converges to a uniform distribution [16]. To quantify the uniformity of this distribution, we again consider the standard χ 2 test. We only perform a single test for the total data set, so we do not need to account for multiplicity. Specifically, we split the interval [0, 1] into 10 equal bins. If the data are uniform, we expect that each bin obtains 1/10 of the total data. Therefore, for each bin in the a m data we expect a value of 10, 045/10 = 1004.5 and for each bin in the b m data we expect a value of 10, 046/10 = 1004.6. Since there are 10 possible observations, we have 9 degrees of freedom for the data. If the percentage of the data in the first nine bins is known, then the percentage of data in the last bin is forced, and we lose one degree of freedom. We again take α = 0.05 and for nine degrees of freedom, this corresponds to a χ 2 value of 16.919. These χ 2 results are generated by cells 16, 17, and 18 by the Jupyter notebooks amLogData.ipynb and bmLogData.ipynb, respectively, which may be found under the Data Analysis folder. We also provide the associated p-values and powers of the χ 2 statistic relative to the null hypothesis that the data are uniform.
The coefficients we studied are 2-adic rational numbers, so we considered the distributions of the numerators, denominators, and decimal expansions separately. We considered only the non-zero coefficients, since zero is not defined for our probability measure, and certain theorems and conjectures outlined by Shiamuchi in [4] already describe the distribution of the zeroes in the coefficients. Our goal is to identify which components of these coefficients are a good fit for Benford's law through statistical testing. Table 1 provides examples of the coefficients computed. When the coefficient is 0, numerators are set to 0 and denominators to -for readability. We then use them to compute the exact values in decimal expansion for a m and b m .

Benfordness of the Taylor and Laurent Coefficients
We examine the distribution of the first digits of the a m and b m coefficients. As mentioned earlier, we restrict our discussion to the non-zero coefficients. We conduct the χ 2 test and distribution of the base 10 logarithms modulo 1 to evaluate the data. The notebooks to generate the results are found under the Data Analysis folder. The plots of the χ 2 values are shown in   The denominators stay below the original threshold for significance. This is expected, as they consist of a random sampling of a geometric series, which is known to be Benford in most bases [16]. The distributions for the base 10 logarithms modulo 1 of the numerators are slightly skewed. There is a pattern in the a m coefficients that could account for this; we have observed that when m = 2 n , a m = 1/m. This result seems to generalize for a d,m , such that when m = d n , a d,m = 1/m, which can be observed in the tables provided by Shiamuchi in [4], and we have not found a counterexample in our computations. There is regularity in the b m numerators as discussed by Bielefeld, Fisher, and von Haeseler in [11], and it is likely that similar regularities are present in the a m coefficients as well. The χ 2 values for the a m and b m data are 8.482 and 10.203, respectively. These correspond to p-values of 0.486 and 0.334; the powers relative to the null hypothesis are 0.482 and 0.574. As a result, there is not sufficient evidence to reject the null hyopthesis that the data are uniform. The denominators consist of a sampling of integer powers of 2. Since log 10 (2) is irrational, the sequence x n = 2 n is Benford in base 10, and log 10 (2 n ) (mod 1) converges to a uniform distribution [16]. Since the denominators span many orders of magnitude, it is expected that they will similarly converge in distribution. The χ 2 values for the a m and b m data are 6.334 and 4.416. These correspond to p-values of 0.706 and 0.882; the powers relative to the null hypothesis are 0.358 and 0.248. As a result, there is not sufficient evidence to reject the null hypothesis that the data are uniform. It is worth noting that the distributions of the logarithms modulo 1 for a m and b m decimal expansions are skewed towards different halves of the interval. This asymmetry may be related to how the series represent coefficients of reciprocal functions and how they may be computed from each other. The χ 2 values for decimal expansions are 64.261 and 60.757. These correspond to p-values of 2.008 × 10 −10 and 9.580 × 10 −10 ; the powers relative to the null hypothesis are 0.99998 and 0.99994. There is sufficient evidence to reject the null hyopthesis that the data are uniform.
We may also investigate the magnitude of the data by computing the arithmetic mean and standard deviation of log 10 |x n |. It is typical, but not necessary, for a data set to be Benford if it spans many orders of magnitude (see Chapter 2 of [7,17] for an analysis that a sufficiently large spread is not enough to ensure Benfordness). Our findings for the Taylor and Laurent coefficients are summarized in Table 2. The data are generated by cell 8 in the Jupyter notebooks amLogData.ipynb and bmLogData.ipynb, which may be found under the Data Analysis folder. The numerators and denominators span many orders of magnitude, while the decimals do not. The mean for the decimal expansions being negative indicates that the denominators are larger than the numerators, on average. The ratio between the growth rates of the numerators and denominators likely has some form of regularity as well to account for the small standard deviation, but more analysis would be needed to determine the exact relationship. These observations are consistent with the previously discussed conjecture that 0 < |b m | < 1/m for all m, and it is plausible a similar relation holds for the a m coefficients as well. Ultimately, this provides insight into the growth of the coefficients and the shape of the data. Our testing provides evidence for convergence of the numerators and denominators to a Benford distribution, but there is not sufficient evidence for convergence in the decimal expansions.

On the Taylor and Laurent Coefficients
This Section deals with observations and theorems on the coefficients from various authors. Our goal is to compile and highlight important results from disparate sources. We link to the papers where the original observations may be found and their proofs when applicable. Section 5.1 deals specifically with observations related to the a m coefficients, Section 5.2 deals with the b m coefficients, and Section 5.3 deals specifically with new observations we make.
Theorem 3 highlights the relevance of the study of the Riemann mappings Ψ and Φ. Much effort has been put into the understanding of the behavior of both series. We now refer to a few important results and introduce new conjectures on the behavior of the coefficients. One of the most important theorems relating to both sets of coefficients is the following. This is a combination of Theorems outlined by Shiamuchi in [3,8], and it provided our motivation for studying the numerators, denominators, and decimal expansions, separately. We consider only the case d = 2. The majority of the following results hold also for a general integer d ≥ 2, under simple modifications.
It is unknown whether the converse is true. The proof may be found in [18]. The authors have reported that their computation of 1000 terms of a 2,m has not produced a zero-coefficient besides those indicated by the theorem, which is consistent with our observations. The result may also be expressed as the following corollary: Corollary 1. Let m = m 0 2 n with n ≥ 0, and m 0 odd. If 3 ≤ m 0 ≤ 2 n+1 , then a 2,m = 0.
Making use of the 2-adic valuation, it is possible to obtain the following theorem, as outlined in [4]. Theorem 8. We have −ν 2 (a m ) ≤ ν 2 ((2m − 2)!) for all m, with equality attained exactly when m is odd.
In the following, since our interest is only for the 2-adic valuation, we will make use of the notation ν(x) := ν 2 (x), and we immediately obtain the following remark.

Remark 4.
In the case that m is odd, we may obtain an efficient algorithm to compute −ν(a m ) through ν((2m − 2)!). Following immediately from Theorem 8 by the properties of the d-adic evaluation outlined in [3], we have, Therefore, if we set a value, N, the denominator's exponent for every odd number m < 2 N is given by We may also summarize Theorem 3.1 and Corollary 3.5 from [8] for the case that d = 2.
Under the same assumptions, the result is also true with

Results for the Laurent Coefficients
Similar results hold for the b m coefficients. We use the notation m = 2 n m 0 , where m 0 is odd. The first result was presented in [19] in 1988.
It is still unknown whether the converse of this theorem is true. In [5], the only coefficients that have been observed to be zero are those mentioned in this theorem. The following result, from [11], underlines that a result similar to the one for the a m holds.
Don Zagier has made several observations and conjectures about the exponents of b m . We shall later extend them to the a m coefficients. The original conjectures are outlined in [11]. These correspond exactly with the values predicted by Corollaries 1 and 5, respectively; when a m = 0, the algorithm gives 0 for the denominator's exponent. In general, for each n there is a partial periodicity with period 2(2 n+1 − 1) in m 0 , and equivalently, 2 n+1 (2 n+1 − 1) in m.
Another direction is to calculate the slope of each of the subsequences, since they seem to grow linearly. Our observations make use of the previous remark and have led to the following conjecture.
The numerators tend to follow similar behavior as the denominators. In particular, the modulus of the numerators in the subsequences tend to organize as follows: {a m } n=0 > {a m } n=1 > · · · {a m } n=N . The possibility of bounding the numerators by making use of its associated denominator is a subject of further study.
We now extend Conjecture 1 to the a m coefficients, as follows: Otherwise, it follows the pattern in Table 3.
This suggests a partial periodicity with period 2 4 (2 n+1 − 1) in m 0 , or of 2 n+4 (2 n+1 − 1) in m. As before, it is possible to write m 0 as 2(2 n+1 − 1)k + l, but it is more difficult to identify a general pattern in this case. Table 3. The distribution of (m 0 ) when m 0 ≡ 12 3 has periodicity 16 × 12= 192. From m 0 = 195, it repeats itself, and will have the same (m 0 ) of m 0 = 3 and following.

Future Work
The most natural extension of our work would be to generate more coefficients, which would allow more thorough statistical testing. Given more data, we could look at the coefficients over a certain zoom or average the coefficients over certain subsets. It would be particularly interesting if certain subsets of the coefficients also converge to the Benford distribution. We may also look at the powers of the denominators to observe whether they follow a Benford distribution.
The algorithms for computing the coefficients of the Mandelbrot can also be easily generalized to obtain other abstract Multibrot sets, which could be analyzed using the same methods. We could look at the data in different bases to observe whether Benfordness holds there. It would be interesting to see the numerators and decimal expansions of the coefficients for the Multibrot set of degree d follow a Benford distribution in base d; the denominators will not since they are sampled from a geometric series with a common ratio d, and they are not Benford in the base of the common ratio.
The results of Section 5 also present interesting extensions for future work. In particular, Remark 6 suggests that dividing the coefficients into subsequences to be bounded separately may be the best approach to study the convergence of the Laurent series of the coefficients. This approach, which has not been followed in the past, to the best of our knowledge, could provide valuable results in the study of the local connectedness of M.