Julia set
In complex dynamics, a Julia set is a fractal subset of the complex plane defined for a holomorphic function f, specifically the closure of the repelling periodic points of f, where the iterates of f exhibit sensitive dependence on initial conditions and chaotic behavior.[1] Named after the French mathematician Gaston Julia, who first described these sets in his seminal 1918 memoir on the iteration of rational functions, the Julia set J(f) complements the Fatou set, partitioning the Riemann sphere into regions of normal and non-normal dynamics under iteration.[2]
For the quadratic family f_c(z) = z^2 + c with c \in \mathbb{C}, the filled Julia set K_c comprises all points z_0 whose forward orbits under iteration remain bounded, i.e., \sup_n |f_c^n(z_0)| < \infty, while the Julia set J_c forms the topological boundary \partial K_c.[3] The geometry of J_c depends critically on c: if c belongs to the Mandelbrot set \mathcal{M}—the set of parameters for which K_c is connected—then J_c is a connected Jordan curve; otherwise, J_c is totally disconnected, resembling a Cantor dust.[4] Small perturbations in c can transform J_c from dendritic forms like the Douady rabbit (for c \approx -0.123 + 0.745i) to basilica-like structures (for c = -1), highlighting the sets' intricate self-similarity and Hausdorff dimension often exceeding 1.[4]
Although Julia's work laid the theoretical foundation, the visual complexity of these sets was not appreciated until the 1970s and 1980s, when Pierre Fatou's earlier studies on iteration were extended and Benoit Mandelbrot used computers to generate images, revealing their fractal nature and linking them to broader phenomena in chaos theory.[3] Julia sets are invariant under f and its inverses, serving as the "skeleton" of dynamical systems where repelling behavior dominates, with applications in studying universality in quadratic dynamics and modeling phenomena like dendritic growth in physics.[1]
Background and History
Complex Dynamics Prerequisites
In complex analysis, a holomorphic function is a function that is complex differentiable at every point in some open domain of the complex plane.[5] These functions satisfy the Cauchy-Riemann equations and possess powerful properties, such as being infinitely differentiable and representable by power series in their domains of holomorphy.[6] To study global behavior in complex dynamics, the complex plane \mathbb{C} is often extended to the Riemann sphere \hat{\mathbb{C}} = \mathbb{C} \cup \{\infty\}, which compactifies the plane by adding a point at infinity and endows it with the topology of a sphere via stereographic projection.[7] Holomorphic functions on the Riemann sphere are rational functions, providing a natural setting for iterating maps without boundary issues at infinity.[8]
The study of complex dynamics centers on iterating a holomorphic map f: \hat{\mathbb{C}} \to \hat{\mathbb{C}}, typically a rational function of degree at least 2, and examining the behavior of orbits \{f^n(z)\}_{n=0}^\infty starting from initial points z \in \hat{\mathbb{C}}.[9] A fixed point of f is a point z such that f(z) = z, while a periodic point of period p \geq 2 satisfies f^p(z) = z but f^k(z) \neq z for $1 \leq k < p.[10] The stability of these points is determined by the multiplier \lambda = (f^p)'(z): if | \lambda | < 1, the point is attracting; if | \lambda | > 1, it is repelling; and if | \lambda | = 1 with \lambda \neq 1, it is neutral (or indifferent).[10] The basin of attraction of an attracting periodic point (or cycle) is the open set of points whose orbits converge to that cycle under iteration.[11] Repelling points, in contrast, have orbits that diverge from them, often leading to chaotic behavior nearby.[12]
Central to understanding iteration behavior is Fatou's theorem on normal families of holomorphic functions, which states that for a rational map f of degree at least 2, the family of its iterates \{f^n\} forms a normal family—meaning every sequence has a subsequence converging uniformly on compact subsets to a holomorphic limit—precisely on the largest open set where the dynamics are "regular."[9] This theorem, proved by Pierre Fatou in 1920, underpins the partition of the Riemann sphere into regions of stable and chaotic dynamics, with Julia sets emerging as the boundaries separating these domains.[9]
Discovery and Development
The study of what would later be known as Julia sets originated in the early 20th century through the pioneering work of French mathematicians Pierre Fatou and Gaston Julia on the iteration of rational functions in the complex plane.[13][14] Fatou initiated this research in 1917 with a memoir applying Paul Montel's concept of normal families to analyze the global behavior of iterated functions, publishing key results in the Comptes Rendus de l'Académie des Sciences between 1919 and 1920 that laid foundational ideas for holomorphic dynamics.[13] Julia, building on similar ideas, submitted his doctoral thesis in 1918, which was published that year in the Journal de Mathématiques Pures et Appliquées as a comprehensive 200-page treatment of iteration theory for rational functions, introducing concepts central to understanding the boundaries of basins of attraction under iteration.[14]
Despite the mathematical depth of these contributions, Fatou and Julia's work on iteration received limited attention during the mid-20th century, largely because the intricate boundary structures they described—now recognized as Julia sets—could not be visualized without computational tools, rendering the abstract results inaccessible to broader mathematical exploration.[15] This neglect persisted until the advent of digital computing in the 1970s, when Benoit Mandelbrot, working at IBM, began generating computer images of these sets, first producing recognizable plots of Julia sets in late 1979.[15]
Mandelbrot's visualizations, culminating in high-quality images by 1980, dramatically revived interest in Julia and Fatou's theories by revealing their fractal nature and self-similar complexity, which Mandelbrot connected to his broader framework of fractal geometry as detailed in his 1982 book The Fractal Geometry of Nature.[15] This rediscovery not only reestablished Julia sets as a cornerstone of complex dynamics but also influenced emerging fields such as chaos theory, where iterative processes model unpredictable behaviors in physical systems, and computer graphics, enabling realistic rendering of natural forms through fractal algorithms.[15]
Core Definitions
In complex dynamics, the Julia set of a holomorphic function f: \hat{\mathbb{C}} \to \hat{\mathbb{C}}, where \hat{\mathbb{C}} denotes the Riemann sphere, is studied through the iteration process defined by z_{n+1} = f(z_n) with initial point z_0 \in \hat{\mathbb{C}}.[4]
For a rational function f of degree at least 2, the Julia set J(f) is the closure of the set of repelling periodic points of f, where a periodic point z of period n satisfies f^n(z) = z and is repelling if |(f^n)'(z)| > 1.[9]
When f is a polynomial, infinity is a superattracting fixed point, and the filled Julia set is defined as K(f) = \{ z \in \mathbb{C} \mid |z_n| \not\to \infty \text{ as } n \to \infty \}, with the Julia set given by J(f) = \partial K(f).[4]
A prominent family of examples is the quadratic polynomials f_c(z) = z^2 + c for c \in \mathbb{C}, where the Julia set J(c) lies in the complex plane and captures the chaotic dynamics on its boundary.[16]
Equivalent Descriptions
One equivalent characterization of the Julia set J(f) for a rational function f: \hat{\mathbb{C}} \to \hat{\mathbb{C}} of degree at least 2 is as the closure of the repelling periodic points of f. A periodic point z of period n is repelling if |(f^n)'(z)| > 1, and the set of all such points is dense in J(f), with their closure precisely equaling J(f). This description highlights the chaotic, expanding nature of the dynamics on the Julia set, as repelling periodic points are fundamental to the unstable behavior under iteration.[17][18]
Another formulation arises from the theory of normal families: the Julia set J(f) consists of those points z \in \hat{\mathbb{C}} such that the family of iterates \{f^n : n \geq 0\} is not normal in any neighborhood of z. Here, normality means the family is locally equicontinuous or, equivalently by Montel's theorem, that there exists a subsequence converging uniformly on compact subsets to a holomorphic function or to infinity. This characterization complements the Fatou set F(f), defined as the largest open set where the iterates form a normal family, making J(f) the complement of F(f).[19][17]
The Julia set can also be described as the boundary of the Fatou set: J(f) = \partial F(f). Since F(f) is open and completely invariant under f, its boundary coincides with the closure of the set where the dynamics exhibit sensitivity to initial conditions, separating stable regions from chaotic ones. This boundary property underscores the topological role of J(f) in partitioning the Riemann sphere.[17][19]
For polynomial maps f: \mathbb{C} \to \mathbb{C} of degree at least 2, the Julia set is the boundary of the basin of attraction of infinity, denoted A_f(\infty) = \{ z \in \mathbb{C} : f^n(z) \to \infty \text{ as } n \to \infty \}, so J(f) = \partial A_f(\infty). Infinity acts as a superattracting fixed point, and this basin forms a connected component of the Fatou set, with the filled Julia set K(f) = \mathbb{C} \setminus A_f(\infty) having boundary J(f).[17][19]
Dynamical Properties
Properties of the Julia Set
The Julia set J(f) of a rational map f of degree at least 2 is a non-empty, compact, and perfect subset of the Riemann sphere, meaning it contains no isolated points and is closed under taking limits. This structure arises from the closure of the repelling periodic points, ensuring the set's completeness and topological richness.
A fundamental topological dichotomy holds: for polynomials, particularly quadratics, the Julia set is either connected or totally disconnected, resembling a Cantor set in the latter case.[20] This classification depends on the boundedness of critical orbits, with connectivity preserved only if no critical point escapes to infinity.[20]
Julia sets are not always locally connected, even when connected; counterexamples include certain infinitely renormalizable quadratic maps constructed via iterated tuning, where pinch points prevent path-connectedness between nearby components. Such non-local connectivity complicates boundary descriptions and visualization.
The Hausdorff dimension of hyperbolic Julia sets equals their hyperbolic dimension, the supremum of dimensions of repelling Cantor subsets, and reaches 2 in many cases, such as those with positive Lebesgue measure, indicating space-filling behavior.[21] For instance, sequences of quadratic maps near certain Mandelbrot set parameters yield Julia sets of dimension arbitrarily close to 2.[21]
For hyperbolic rational maps, the dynamics on the Julia set are ergodic with respect to the unique conformal measure of maximal entropy, implying strong mixing properties and uniform expansion away from critical points.[22] This ergodicity underpins the statistical uniformity of orbits on J(f).[22]
Small perturbations in the parameter c for quadratic maps z^2 + c can induce drastic topological transitions in the Julia set, such as from connected to totally disconnected or shifts in connectivity loci.[20]
Julia Set versus Fatou Set
The Fatou set F(f) of a rational function f: \hat{\mathbb{C}} \to \hat{\mathbb{C}} is defined as the largest open subset of the Riemann sphere \hat{\mathbb{C}} on which the family of iterates \{f^n : n \geq 0\} is normal, meaning every sequence of iterates has a subsequence that converges uniformly on compact subsets to a holomorphic function or to infinity.[23] This normality is equivalent to the family being equicontinuous on F(f), due to Montel's theorem applied to families of meromorphic functions omitting at most three values.[24] In contrast, the Julia set J(f) = \hat{\mathbb{C}} \setminus F(f) consists of points where the iterates exhibit chaotic behavior, lacking such uniform convergence properties.[9]
The connected components of the Fatou set, known as Fatou components, are classified into periodic, preperiodic, and wandering types, though the latter do not occur for rational maps. Periodic Fatou components are those invariant under some iterate f^k, and by the works of Fatou, Julia, Siegel, Branner-Hubbard, and others, they fall into four main categories: attracting or superattracting basins, where iterates converge to an attracting periodic cycle; parabolic basins, associated with rationally indifferent fixed points; Siegel disks, where dynamics are conjugate to irrational rotations on the disk; and Herman rings, annular regions with similar rotational dynamics.[25] Preperiodic components map eventually to periodic ones under iteration. Sullivan's no wandering domains theorem proves that for any rational f, every Fatou component is eventually periodic, resolving a long-standing conjecture and ensuring the Fatou set has no wandering components—open sets whose orbits under f remain disjoint forever.
The Julia set serves as the boundary of every Fatou component, enclosing the stable regions of F(f) where dynamics are tame, while J(f) captures the chaotic frontier. Since rational maps are open by the open mapping theorem, f maps open sets to open sets, implying that f(F(f)) \subseteq F(f) and f(J(f)) \subseteq J(f); moreover, the sets are completely invariant, as preimages preserve the normality or chaos.[23] For polynomial maps, the basin of infinity—consisting of points whose orbits escape to \infty—forms a distinguished unbounded Fatou component, often simply connected and containing all escaping points.[26]
Examples and Families
Basic Examples
Julia sets arise prominently in the study of the quadratic family of maps f_c(z) = z^2 + c, where c is a complex parameter; the filled Julia set K_c consists of points whose orbits under iteration of f_c remain bounded, while the Julia set J_c is the boundary of K_c, forming the boundary between escaping and non-escaping regions.[27]
A basic connected example occurs for c = -2, where J_{-2} is the closed interval [-2, 2] on the real axis, a simple dendritic structure extended along the real line without branches or interior points.[28] In this case, iteration of points within [-2, 2] produces bounded orbits that oscillate within the interval, while points outside escape to infinity. The fixed points, solutions to z = z^2 - 2, are z = 2 and z = -1, with multipliers \lambda = 2z = 4 and \lambda = -2, respectively; both have |\lambda| > 1, rendering them repelling.[27]
For a totally disconnected example, consider c = 0.25 + 0.5i, which lies on the boundary of the parameter region yielding connected sets; however, slight perturbations outside produce dust-like structures, but strictly for values like c = 0.5 (real), J_{0.5} is a Cantor set, consisting of uncountably many isolated points with no connected components larger than a singleton, resembling fractal dust scattered along the real line.[29] Here, all orbits on the Julia set are bounded but dense in small neighborhoods, while nearby points escape rapidly to infinity. The fixed points solve z^2 - z + 0.5 = 0, yielding z = \frac{1 \pm \sqrt{-1}}{2} = 0.5 \pm 0.5i, with multipliers \lambda = 2z = 1 \pm i; both satisfy |\lambda| = \sqrt{2} > 1, so they are repelling.[27]
For real parameters c \in (-2, 0.25), the Julia sets J_c are connected and exhibit varied qualitative shapes: near c = -2, they resemble thin rods or line-like dendrites; around c = 0, they form rounded, cauliflower-like boundaries enclosing a disk-like filled set; and approaching c = 0.25, they become pinched and quasi-circular with increasing intricacy at the parabolic cusp.[28] These forms illustrate how small changes in c transform the boundary dynamics while maintaining connectivity.
Quadratic Julia Sets
The quadratic family of maps is defined by f_c(z) = z^2 + c, where c is a complex parameter varying over the plane \mathbb{C}. The associated Julia set J_c is the boundary of the filled Julia set K_c = \{ z \in \mathbb{C} : \sup_{n \geq 0} |f_c^n(z)| < \infty \}, and its connectivity depends critically on c. Specifically, K_c (and thus J_c) is connected if and only if c lies in the Mandelbrot set M = \{ c \in \mathbb{C} : K_c is connected \}; otherwise, K_c is totally disconnected (a Cantor set) and coincides with J_c.[20][4]
The dynamics of f_c are governed by the orbit of its sole finite critical point at z = 0, as the behavior of all other orbits is subordinate to this one under the quadratic structure. The orbit of 0 remains bounded if and only if c \in M, which ensures the connectivity of J_c; unbounded escape of this orbit leads to a disconnected Julia set.[4][20] This critical orbit criterion provides a practical means to determine membership in M and classify the topology of J_c.
The Mandelbrot set M possesses a distinctive bulb structure, comprising the main cardioid—a central region where f_c admits an attracting fixed point—and an infinite collection of period-n bulbs attached to it, each corresponding to parameters with an attracting cycle of minimal period n > 1. These bulbs adorn the cardioid in a self-similar fashion, with the largest being the period-2 bulb to the left and smaller ones spiraling outward for higher periods. For example, the case c = -2 yields a connected dendritic Julia set without attracting cycles.[20][30]
Hyperbolic components form the open interior of M, consisting of all c for which f_c is hyperbolic (i.e., possesses an attracting cycle to which the critical point converges). Each such component is bounded by a curve where the multiplier of the attracting cycle has modulus exactly 1, marking the transition to non-hyperbolic dynamics; the component's root (where the multiplier is 1) and center (where it is 0) define key structural points.[20]
In the Douady-Hubbard theory, the combinatorial structure of M and its Julia sets is illuminated through external rays—curves from infinity parameterized by angles \theta \in \mathbb{R}/\mathbb{Z}—which land at specific boundary points. Rational-argument external rays of M land precisely at the roots of hyperbolic components or at Misiurewicz points (where the critical orbit is preperiodic); at these roots, the rays coincide with the landing points of internal rays from the component's dynamics, establishing a bijection between parameter and dynamical planes.[20]
For a fixed point z^* of f_c, satisfying z^* = (z^*)^2 + c, the multiplier is given by
\lambda = f_c'(z^*) = 2 z^*.
Attracting fixed points occur when |\lambda| < 1, which parameterizes the main cardioid via the relation c = \frac{\lambda}{2} (1 - \frac{\lambda}{2}). This multiplier governs local stability and bifurcations at the cardioid's boundary.[20]
Generalizations and Extensions
To Higher-Degree Polynomials
Julia sets can be generalized to monic polynomials of degree d > 2, such as f(z) = z^d + c where c \in \mathbb{C}. For this unicritical family, there is a critical point at z = 0 with multiplicity d-1, as the derivative is f'(z) = d z^{d-1}, vanishing to order d-1 at the origin. In this family, the dynamics are determined by the orbit of this single critical point. While the unicritical family z^d + c captures an analog via the Multibrot set, general monic degree d polynomials have d-1 finite critical points, leading to a higher-dimensional parameter space.
The parameter space for these unicritical families is known as the Multibrot set M_d, the analog of the Mandelbrot set for degree d, consisting of parameters c for which the orbit of the critical point 0 remains bounded under iteration of f_c(z) = z^d + c. For d > 2, the Multibrot set exhibits greater structural complexity than the quadratic Mandelbrot set, featuring more hyperbolic components and higher-period bifurcation loci due to the higher degree and multiplicity.[31]
A key dynamical property carries over from quadratics: the Julia set J(f) of a polynomial f of degree d \geq 2 is connected if and only if the orbits of all finite critical points remain bounded. If any critical orbit escapes to infinity, J(f) becomes a Cantor set of uncountably many components. This criterion generalizes the quadratic connectivity condition, where boundedness of the single critical orbit suffices.[32]
For monic polynomials with connected Julia sets, the basin of infinity \Omega(f, \infty) is simply connected. This follows from the existence of a conformal isomorphism (via an adjusted Böttcher function) mapping \Omega(f, \infty) onto the open unit disk, possible precisely because no finite critical points lie in the basin. In contrast to the quadratic case, the presence of multiple critical points in general polynomials allows for richer combinatorial structures, such as captures of one critical orbit by another, contributing to the increased complexity of hyperbolic components and periodic cycles in the parameter space.[33]
A representative example is the cubic family f(z) = z^3 + c, which has a single critical point at 0 with multiplicity 2. The Multibrot set M_3 is the connectivity locus where the orbit of 0 remains bounded. In contrast, the general cubic family, such as f(z) = z^3 + a z + b, has two distinct critical points (solutions to $3z^2 + a = 0), requiring bounded orbits for both to ensure connectivity, and features a two-dimensional parameter space with more intricate structure than the quadratic Mandelbrot set.
To Non-Polynomial Maps
The generalization of Julia sets to non-polynomial holomorphic maps encompasses rational functions of degree d \geq 2 and transcendental entire functions, introducing distinct dynamical behaviors on the Riemann sphere or extended complex plane. For a rational map f(z) = p(z)/q(z), where p and q are coprime polynomials of degree d, the map has exactly $2d-2 critical points, counted with multiplicity, as the derivative f'(z) is a rational function of degree $2d-2. The Julia set J(f) is the closure of the repelling periodic points of f, and it coincides with the closure of all repelling periodic points under iteration.[34][9]
Fatou components for rational maps include attracting and parabolic basins, as well as rotation domains such as Siegel disks (where f is conjugate to an irrational rotation on a disk) and Herman rings (annular rotation domains, possible only for degree at least 3). Unlike polynomial maps, rational functions act on the compact Riemann sphere \hat{\mathbb{C}}, so there is no dedicated superattracting basin of infinity; instead, infinity behaves as an ordinary point in the dynamics. By Sullivan's no-wandering-domains theorem, every Fatou component under a rational map is eventually periodic, precluding wandering domains where forward images remain disjoint.[9]
For transcendental entire functions, such as exponentials or trigonometric maps, the Julia set J(f) is the set where the family of iterates \{f^n\} fails to be normal, often exhibiting more complex structures due to the non-compact domain. These functions possess infinitely many critical points and asymptotic values (points omitted in the range, including \infty), with singular values defined as the critical values and asymptotic values; the postsingular set, the closure of the forward orbits of singular values, governs much of the dynamics. A canonical example is the exponential map f(z) = \lambda e^z with $0 < \lambda < 1/e, where \lambda is the sole finite singular value (as 0 is an asymptotic value), and the escaping set consists of dynamic rays landing at repelling points or the singular value's orbit.[35][35]
Transcendental entire functions admit phenomena absent in rational maps, including Baker domains—Fatou components U where f^n(z) \to \infty for z \in U but f|_U is not onto an attracting cycle—and wandering domains, whose forward orbits consist of infinitely many distinct components. Baker domains were first identified in the context of normality domains for entire functions, and examples include certain quadratic entire maps like f(z) = z^2 + c \sin(z) for specific c. Wandering domains, conjectured by Fatou, are ruled out for rational maps by Sullivan but exist for transcendentals, as in f(z) = z^2 e^{z + \lambda} for small \lambda \neq 0. If the Fatou set of a transcendental entire function is nonempty, then J(f) has empty interior in \mathbb{C}.[36][36]
Trigonometric entire functions provide illustrative examples; for f(z) = \sin(z), the Julia set is the entire complex plane \mathbb{C}, as the iterates exhibit no normality anywhere due to unbounded critical orbits and dense escaping points. Similarly, the sine family \lambda \sin(z) for \lambda > 0 yields Julia sets of positive area, often Cantor bouquets of curves, highlighting the rich topology possible in transcendental dynamics.[35][35]
Computation and Visualization
Plotting Algorithms
The most straightforward method for plotting Julia sets involves forward iteration, also known as the escape-time algorithm. In this approach, a grid of points is defined in the complex plane, typically within a bounded region such as a square centered at the origin. For each grid point z_0, the iteration z_{n+1} = f(z_n) is performed starting from z_0, where f is the generating function (e.g., a quadratic polynomial f(z) = z^2 + c). The process continues until either |z_n| exceeds a predefined escape radius (often 2 or larger, depending on f) or a maximum number of iterations (e.g., 100–1000) is reached. Points that escape are colored according to the iteration count at which escape occurs, creating a smooth gradient that highlights regions near the Julia set boundary; points that remain bounded after the maximum iterations are typically colored black or a solid hue to indicate membership in the filled Julia set K(f), the set of points with bounded orbits under f. This method effectively visualizes repelling Julia sets but can be computationally intensive for fine grids or high iteration counts.[37]
The inverse iteration method (IIM) addresses limitations of forward iteration, particularly for connected Julia sets, by working backwards from the exterior of the set. It begins with a set of known exterior points, such as a large disk or circle far from the origin where |z| is sufficiently large to ensure escape under forward iteration. For each such point z, the equation f(w) = z is solved for preimages w, yielding multiple solutions (e.g., two for quadratic f). These preimages are iteratively computed over many steps (often hundreds), with random or systematic selection among branches to densely sample the plane. As iterations proceed, the preimages converge to the Julia set boundary, which acts as an attracting set in the inverse dynamics, allowing efficient filling of interior regions without exhaustively testing every grid point. This technique excels at rendering detailed boundaries for polynomials where the Julia set is connected, though it requires careful branch selection to avoid biases toward certain angles.[38][39]
For enhanced boundary smoothness and antialiasing, the distance estimation method (DEM) approximates the distance from each point to the Julia set boundary, enabling subpixel refinement. Starting from a point z_0 outside the set, forward iteration computes z_n = f^{\circ n}(z_0) and the accumulated derivative z'_n = \prod_{k=0}^{n-1} f'(z_k), with z'_0 = 1. The estimated distance is then d \approx \frac{|z_n| \cdot \ln |z_n|}{|z'_n|}, derived from the Green potential function G(z_0) \approx \frac{\ln |f^{\circ n}(z_0)|}{n} and |G'(z_0)| \approx \frac{|z'_n|}{|z_n|}, so d = \frac{G(z_0)}{|G'(z_0)|}. This provides a lower bound on the true distance with controlled error for large n and smooth boundaries. Coloring or shading based on this distance yields antialiased images, particularly useful for quadratic Julia sets where boundaries are C^1-smooth away from critical points. DEM reduces artifacts in zoomed views but demands higher iteration counts near the boundary for accuracy.[40]
To plot the filled Julia set K(f), the boundedness test is central: during forward iteration, points whose orbits remain within the escape radius after the maximum iterations are classified as interior to K(f), while escapers are exterior. For quadratic maps f(z) = z^2 + c, an escape criterion guarantees that if |z_n| > \max(|c|, 2), the orbit diverges to infinity, allowing reliable termination; this test extends to higher-degree polynomials with adjusted radii based on coefficients. Challenges in rendering include aliasing from discrete grids, which DEM partially mitigates via distance-based interpolation, and precision loss in deep zooms, where standard double-precision floating-point arithmetic fails due to accumulated rounding errors in iterations, necessitating arbitrary-precision libraries (e.g., GMP) that increase computation time exponentially with zoom depth.[41][42]
Pseudocode Implementations
The standard method for generating images of quadratic Julia sets employs forward iteration of the map z_{n+1} = z_n^2 + c, where c is a fixed complex parameter and initial points z_0 are sampled on a grid in the complex plane. For each grid point, the orbit is iterated until |z_n| exceeds a threshold (typically 2, as escape to infinity occurs for |z| > 2 under this map), or a maximum iteration count is reached; the escape iteration count determines the coloring, with non-escaping points colored black to indicate bounded orbits.[37]
The following pseudocode illustrates this escape-time algorithm for a filled Julia set, where the filled set consists of points with bounded orbits (colored black if no escape within max iterations), distinguishing it from the normal Julia set, which is solely the boundary of escaping and non-escaping regions.
function generate_filled_julia(width, height, c_real, c_imag, max_iter, escape_radius=2.0):
image = array of size (width, height) initialized to 0 # 0 for black (bounded)
for x in 0 to width-1:
for y in 0 to height-1:
z_real = scale(x, width, -2.0, 2.0) # Map grid to complex plane, e.g., [-2,2] x [-2,2]
z_imag = scale(y, height, -2.0, 2.0)
c = complex(c_real, c_imag)
z = complex(z_real, z_imag)
iter = 0
while |z| <= escape_radius and iter < max_iter:
z = z^2 + c
iter += 1
if iter == max_iter:
image[x][y] = 0 # Bounded orbit: filled set interior
else:
image[x][y] = iter # Escape time for coloring
return image
function generate_filled_julia(width, height, c_real, c_imag, max_iter, escape_radius=2.0):
image = array of size (width, height) initialized to 0 # 0 for black (bounded)
for x in 0 to width-1:
for y in 0 to height-1:
z_real = scale(x, width, -2.0, 2.0) # Map grid to complex plane, e.g., [-2,2] x [-2,2]
z_imag = scale(y, height, -2.0, 2.0)
c = complex(c_real, c_imag)
z = complex(z_real, z_imag)
iter = 0
while |z| <= escape_radius and iter < max_iter:
z = z^2 + c
iter += 1
if iter == max_iter:
image[x][y] = 0 # Bounded orbit: filled set interior
else:
image[x][y] = iter # Escape time for coloring
return image
This approach produces the filled Julia set K_c = \{ z \in \mathbb{C} : |z_n| \leq 2 \ \forall n \}, whose boundary is the Julia set J_c; for the normal Julia set, only boundary points could be highlighted via additional refinement, but the pseudocode above focuses on the filled version for practical visualization.[37]
Optimizations enhance efficiency, such as early escape checks where iteration halts immediately if |z_n| > 2 (leveraging the fact that orbits diverge thereafter) and periodicity detection to identify repeating cycles in bounded orbits, avoiding further computation. Periodicity can be detected by storing visited points in a set and checking for repeats, though floating-point precision limits its reliability in practice.[43]
For multi-Julia sets, which explore variations across multiple parameters c (e.g., perturbed maps or slices through the parameter space), the algorithm extends by nesting loops over a grid of c values, generating a composite image or animation from individual Julia sets. Perturbed maps might involve adding small noise to c per iteration or point, but the core remains forward iteration per c.
function generate_multi_julia(c_width, c_height, ...): # Similar args as above
multi_image = array of size (c_width * width, c_height * height)
for cx in 0 to c_width-1:
for cy in 0 to c_height-1:
c_real = scale(cx, c_width, c_min_real, c_max_real) # e.g., c in [-0.8, 0.2] + i[-0.2, 0.2]
c_imag = scale(cy, c_height, c_min_imag, c_max_imag)
sub_image = generate_filled_julia(width, height, c_real, c_imag, max_iter)
copy sub_image to multi_image at offset (cx * width, cy * height)
return multi_image
function generate_multi_julia(c_width, c_height, ...): # Similar args as above
multi_image = array of size (c_width * width, c_height * height)
for cx in 0 to c_width-1:
for cy in 0 to c_height-1:
c_real = scale(cx, c_width, c_min_real, c_max_real) # e.g., c in [-0.8, 0.2] + i[-0.2, 0.2]
c_imag = scale(cy, c_height, c_min_imag, c_max_imag)
sub_image = generate_filled_julia(width, height, c_real, c_imag, max_iter)
copy sub_image to multi_image at offset (cx * width, cy * height)
return multi_image
This produces a parameter-space visualization akin to a Mandelbrot slice, with each sub-image a Julia set for that c.[43]
Inverse iteration reverses the process, starting from points near infinity or repelling fixed points and applying branching inverses to fill preimage sets approximating the Julia set; for f(z) = z^2 + c, the two inverses are w = \pm \sqrt{z - c}, chosen randomly or exhaustively up to depth N to generate level sets. This method is particularly effective for connected Julia sets, as preimages accumulate densely on the boundary.
function generate_inverse_julia(c_real, c_imag, max_depth, num_samples, grid_size):
c = complex(c_real, c_imag)
preimages = set() # Store points on Julia set approximation
# Start from repelling fixed point or large |z| (e.g., infinity approximation)
start_points = sample_random_points_outside(grid_size, radius=10.0) # e.g., 1000 points
for start_z in start_points:
current = start_z
for depth in 1 to max_depth:
# Choose branch randomly (0 or 1 for + or - sqrt)
branch = random(0, 1)
inv = sqrt(current - c) * (1 if branch == 0 else -1)
add inv to preimages if within grid bounds
current = inv # Backtrack
image = mark_grid_points(preimages, grid_size) # Density or presence on grid
return image
function generate_inverse_julia(c_real, c_imag, max_depth, num_samples, grid_size):
c = complex(c_real, c_imag)
preimages = set() # Store points on Julia set approximation
# Start from repelling fixed point or large |z| (e.g., infinity approximation)
start_points = sample_random_points_outside(grid_size, radius=10.0) # e.g., 1000 points
for start_z in start_points:
current = start_z
for depth in 1 to max_depth:
# Choose branch randomly (0 or 1 for + or - sqrt)
branch = random(0, 1)
inv = sqrt(current - c) * (1 if branch == 0 else -1)
add inv to preimages if within grid bounds
current = inv # Backtrack
image = mark_grid_points(preimages, grid_size) # Density or presence on grid
return image
The approximation improves with higher N and more samples, converging to J_c as the limit of preimage boundaries.[44]
The distance estimation method (DEM) provides smooth distance-to-boundary values for enhanced coloring or ray marching, iterating both the map and its derivative: z_{n+1} = z_n^2 + c, dz_{n+1} = 2 z_n \, dz_n (with dz_0 = 1), then estimating distance as |z| \log |z| / |dz| upon escape.
function julia_dem_distance(z_real, z_imag, c_real, c_imag, max_iter, escape_radius=2.0):
c = [complex](/page/Complex)(c_real, c_imag)
z = [complex](/page/Complex)(z_real, z_imag)
dz = [complex](/page/Complex)(1.0, 0.0) # Initial [derivative](/page/Derivative)
iter = 0
while |z| <= escape_radius and iter < max_iter:
dz_new = 2 * z * dz
z = z^2 + c
dz = dz_new
iter += 1
if iter == max_iter:
return infinity # Or large value for interior
else:
return |z| * log(|z|) / |dz| # [Distance](/page/Distance) estimate to boundary
function julia_dem_distance(z_real, z_imag, c_real, c_imag, max_iter, escape_radius=2.0):
c = [complex](/page/Complex)(c_real, c_imag)
z = [complex](/page/Complex)(z_real, z_imag)
dz = [complex](/page/Complex)(1.0, 0.0) # Initial [derivative](/page/Derivative)
iter = 0
while |z| <= escape_radius and iter < max_iter:
dz_new = 2 * z * dz
z = z^2 + c
dz = dz_new
iter += 1
if iter == max_iter:
return infinity # Or large value for interior
else:
return |z| * log(|z|) / |dz| # [Distance](/page/Distance) estimate to boundary
This yields sub-pixel accuracy for exterior points, useful for anti-aliased rendering of Julia set boundaries.[40]
Potential Function
In the context of polynomial Julia sets, the potential function \phi(z) for a monic polynomial f of degree d \geq 2 is defined on the basin of infinity A_\infty(f) as
\phi(z) = \lim_{n \to \infty} \frac{1}{d^n} \log^+ |f^{\circ n}(z)|,
where f^{\circ n} denotes the n-th iterate of f and \log^+ x = \max(\log x, 0).[45] This limit exists and defines a continuous, non-negative function that vanishes identically on the filled Julia set K(f) (and hence on the Julia set J(f)) but is strictly positive throughout A_\infty(f).[9]
The function \phi is harmonic on A_\infty(f), satisfying Laplace's equation \Delta \phi = 0 there, and extends continuously to zero on the boundary \partial K(f). Near infinity, it admits the asymptotic expansion \phi(z) \sim \frac{1}{d} \log |z| as |z| \to \infty.[9] Moreover, \phi obeys the functional equation \phi(f(z)) = d \cdot \phi(z) for all z \in A_\infty(f), reflecting the scaling induced by iteration.[46]
This potential function is intimately connected to the Böttcher coordinate, a conformal map \psi: A_\infty(f) \to \{w \in \mathbb{C} : |w| > 1\} defined near infinity that conjugates f to the model map w \mapsto w^d, with \psi(z) \sim z as |z| \to \infty. Specifically, \phi(z) = \log |\psi(z)| outside K(f), linking the potential to external coordinates in the dynamical plane.[9]
The inverse of the potential provides a notion of the real iteration number: for t > 0, the level set \phi^{-1}(t) = \{z \in A_\infty(f) : \phi(z) = t\} consists of equipotential curves corresponding to rays at potential level t, which parameterize the continuous escape rate to infinity.[46] These equipotential curves are particularly useful in applications, such as rendering smooth color gradients in visualizations of Julia sets by interpolating between discrete iteration counts.[9]
Field Lines and External Rays
In the study of Julia sets for holomorphic maps, field lines refer to the gradient lines of the Green's function, or potential function, which describe the direction of steepest ascent in the external potential field surrounding the filled Julia set K(f). These lines are orthogonal to the equipotential curves, where the potential is constant, forming a network analogous to electric field lines around a charged conductor.[47]
External rays are a specific family of these field lines, parameterized by an angle \theta \in [0,1), that extend from infinity toward the Julia set J(f) while maintaining a constant argument \theta under the conformal mapping \phi_f that linearizes the dynamics at infinity. Each ray R_\theta is defined as the curve \phi_f^{-1}(r e^{2\pi i \theta}) for r > 1, where the parameter r = e^t decreases from \infty to 1 as t, the potential, decreases to 0, so R_\theta(t) = \phi_f^{-1}(e^t e^{2\pi i \theta}).[20][47]
For quadratic polynomials f_c(z) = z^2 + c with connected K(f_c), the external rays are conjugate via the Böttcher function to the angle-doubling map on the unit circle, meaning the dynamics on the rays mirrors the doubling of angles \theta \mapsto 2\theta \mod 1. In this setting, Douady and Hubbard established that external rays land at repelling periodic points on J(f_c), with rays of rational angle \theta = p/q in lowest terms landing at periodic points of period dividing q.[20][48]
The landing behavior of rays provides a combinatorial model for the Julia set known as the pinching model or lamination: the boundary J(f) can be constructed by starting with the circle of rays at infinity and pinching together pairs of points with angles whose rays land at the same point on J(f), yielding a lamination whose quotient is homeomorphic to J(f) when it is locally connected.[20]
However, for Julia sets that are not locally connected, such as those containing Cremer points, some external rays fail to land on J(f), instead accumulating on the boundary without terminating at a single point; this non-landing is a direct consequence of the failure of local connectivity.[47]