§ INTRODUCTION

In daily life, the word “chaos” often brings to mind a state of complete disorder and randomness. Edward Lorenz, the father of chaos theory, famously summarized chaos as “when the present determines the future but the approximate present does not approximately determine the future.” These two ideas almost entirely capture the meaning of chaos in dynamical systems theory: chaos is characterized by a sensitive dependence on initial conditions, often leading to behavior that appears random. One of the ways in which chaos arises in a system is through sensitive dependence on initial conditions emerging from its temporal dynamics, commonly referred to as the “butterfly effect.” This is the idea that one small disturbance can result in large implications later on, or in the words of Lorenz, a butterfly flapping its wings in Brazil can set off a tornado in Texas. In dynamical systems theory, this magnification of an initial perturbation can be quantified using Lyapunov exponents, which describe the rate at which nearby trajectories diverge from each other. This type of “temporal” chaos has been well-studied in dynamical systems theory; however, there is another type of chaos that has been less studied in previous research. Commonly referred to as unpredictability, this second type of chaos emerges from a geometrical property or structure of a dynamical system, specifically in systems that exhibit a property known as multistability, which is when a system can achieve multiple different stable states depending on its initial conditions. It follows that this geometrical sensitivity occurs when close initial conditions can end up in completely different final states. Interestingly, both these “temporal” and “geometrical” forms of chaos are related to each other and associated with the presence of fractal geometry. A major focus of this paper is on the types of geometric structures that lead to or result from sensitivity to initial conditions. Specifically, we explore methods of detecting, classifying, and quantifying these structures and their relation to chaos.

Another major focus of this paper is the application of dynamical systems to neurobiology and neuronal modeling. Specifically, we research two systems established by Nikolai F. Rulkov in the early 2000s. These systems are simple models of biological neurons capable of modeling a variety of neuronal behaviors, including silence, spiking, bursts of spikes, and chaotic spiking-bursting. The simplicity of the Rulkov neuron maps allows for experimentation that would not be possible with more computationally intensive models. For example, we explore arranging neurons into complex physical structures and coupling them with a flow of current.

In the past, researchers have mainly focused on the Rulkov maps' biological properties, with minimal research being done on the chaotic dynamics and possibility of multistability and complex geometric structures that could lead to sensitivity to initial conditions in these maps. Therefore, the aims of this project and the main source of original research presented in this paper are as follows:

* to analyze and quantify the chaotic dynamics of uncoupled Rulkov neurons,

* to model the behavior and analyze the dynamics of complex Rulkov neuron systems,

* to explore the possibility of the existence of multistability and fractal geometry in the Rulkov maps,

* and if this multistability and geometry exists, to detect, classify, and quantify the geometrical properties associated with chaos and unpredictability.

We believe that this research will provide a novel perspective on how these geometrical properties emerging from and resulting in sensitivity to initial conditions might appear in mathematical neuron systems and biological neurons in general. Another major aim of this paper was to curate the language and methods to be approachable and accessible for a more general audience than the target audience of most papers in chaos theory. For this reason, much of this paper is written in a conversational manner, and concepts are presented in an intuitive way. A reader with only a background in vector algebra and elementary calculus should be able to understand this paper, with most of the higher-level mathematics being delegated to the appendices.

This paper is organized as follows. In Section <ref>, we cover important background information from the theory of dynamical systems. In Section <ref>, we detail much of the past and current literature, as well as our own additions and modifications, relating to the geometrical properties and structures that emerge in dynamical systems, with a focus on those that result from or lead to chaos and unpredictability. In Section <ref>, we discuss the concept of a slow-fast system and develop slow-fast dynamical systems theory for discrete-time systems. To motivate the Rulkov maps, we also discuss neurons and neuronal modeling from a biological perspective. In Section <ref>, we introduce the Rulkov maps, performing a thorough analysis of their non-chaotic and chaotic dynamics, bifurcations, and regimes of neuronal behavior. In Section <ref>, we model the injection of direct current into the Rulkov maps, using the example of a pulse of electrical current. In Section <ref>, we establish the method of coupling Rulkov neurons with a flow of current and perform an in-depth analysis of the dynamics of multiple coupled neuron systems. In Section <ref>, we combine all of our research, applying the mathematical and computational theory from Section <ref> to the Rulkov map systems established in Sections <ref>, <ref>, and <ref> to analyze three systems for geometrical properties associated with sensitivity to initial conditions.

§ BACKGROUND

A dynamical system can be thought of as being a point that evolves, or moves around, according to a specific set of rules meiss. The location of this point is known as the dynamical system's state, which is dependent on and evolves according to a time variable denoted as $t$ [p. 1]arrowsmith. All of the possible states of a dynamical system live in the system's state space, and each state in state space corresponds to a unique point. A simple example of a dynamical system is a ball rolling across a table, which evolves in accordance with Newton's laws of motion. This system's state at a given time $t$ is the ball's location at that time $t$, and the system's state space is the entire surface of the table, which contains all the possible locations of the ball.

In this paper, we will concern ourselves with deterministic dynamical systems, where every state can be determined uniquely from past states, in contrast with stochastic dynamical systems, where randomness is involved [p. 2]alligood. For this reason, we can calculate all of a deterministic dynamical system's future states if we know its initial state meiss. It is important to note that the initial state of a system is not necessarily the state of the system when it was created, but rather, the state at the beginning of any relevant stretch of time [p. 9]lorenz.

There are two main types of dynamical systems: continuous-time and discrete-time. In continuous-time dynamical systems, the time variable is continuous ($t\in\mathbb{R}$).[The symbol $\in$ means “is an element of.”] Consider a system in real $n$-dimensional state space, or $\mathbb{R}^n$, with the state vector $\mathbf{x}$:

\begin{equation}
\mathbf{x} = \begin{pmatrix}
x\e{1} \\
x\e{2} \\
\vdots \\
\end{pmatrix}
= \Big\langle x\e{1},\,x\e{2},\,\hdots\,,\,x\e{n}\Big\rangle
\label{eq:statevector}
\end{equation}

where $x\e{m}$ is the component of $\mathbf{x}$ in the $m$th dimension of space.[$[m]$ is not an exponent.] We can describe the dynamics of a continuous-time dynamical system using a set of differential equations [p. 1]arrowsmith:

\begin{equation}
\frac{d\mathbf{x}}{dt} = \Dot{\mathbf{x}} = \mathbf{g}(\mathbf{x})
\label{eq:diffeq-general}
\end{equation}

For discrete-time dynamical systems, the time variable is discrete ($t\in\mathbb{N}$), and the dynamics are governed by a mapping or iteration function [p. 9]martelli:

\begin{equation}
\mathbf{x}_{k+1} = \mathbf{f}(\mathbf{x}_k)
\label{eq:iterationcond}
\end{equation}

where $\mathbf{x}_k$ is the state of the system at the time $t=k$. In this paper, we will mainly focus on discrete-time dynamical systems.

The forward orbit or trajectory $O^+(\mathbf{x}_0)$ of a discrete-time dynamical system is the set of all its iterates, generated by iterating an initial state $\mathbf{x}_0$ using the iteration function $\mathbf{f}$ an infinite number of times [p. 411]layek:

\begin{align}
\begin{split}
O^+(\mathbf{x}_0) &= \{\mathbf{f}^t(\mathbf{x}_0)\}_{t=0}^{\infty} \\
&= \{\mathbf{x}_0,\, \mathbf{f}(\mathbf{x}_0),\, \mathbf{f}^2(\mathbf{x}_0),\, \mathbf{f}^3(\mathbf{x}_0),\, \hdots\} \\
&= \{\mathbf{x}_0,\, \mathbf{x}_1,\, \mathbf{x}_2,\, \mathbf{x}_3,\, \hdots\}
\end{split}
\end{align}

\begin{equation}
\mathbf{f}^t(\mathbf{x}) = (\mathbf{f} \circ \mathbf{f} \circ \hdots \circ \mathbf{f})(\mathbf{x}) = \mathbf{f}(\mathbf{f}(\hdots \mathbf{f}(\mathbf{x})\hdots)
\end{equation}

$\mathbf{f}$ being repeated $t$ times. A general orbit $O(\mathbf{x})$ of length $j$ is simply a subset of $O^+(\mathbf{x}_0)$ that starts at $\mathbf{x}$ and ends at $\mathbf{f}^j(\mathbf{x})$:

\begin{equation}
\begin{split}
O(\mathbf{x}) &= \{\mathbf{f}^t(\mathbf{x})\}_{t=0}^j \\
&= \{\mathbf{x},\, \mathbf{f}(\mathbf{x}),\, \mathbf{f}^2(\mathbf{x}),\,\hdots\,,\, \mathbf{f}^j(\mathbf{x})\} \\
\end{split}
\end{equation}

An orbit can then be thought of as a subset of state space connected through time evolution.

To establish some basic definitions, a state $\mathbf{x}_s$ is stationary if it satisfies

\begin{equation}
\mathbf{f}(\mathbf{x}_s) = \mathbf{x}_s
\label{eq:stationary}
\end{equation}

If this holds, the point $\mathbf{x}_s$ is called a fixed point of $\mathbf{f}$ [p. 16]martelli. Similarly, a state $\mathbf{x}_p$ is a periodic point of $\mathbf{f}$ if

\begin{equation}
\mathbf{f}^q(\mathbf{x}_p) = \mathbf{x}_p
\label{eq:period}
\end{equation}

where $q\geq 1$. The periodic orbit or $q$-cycle of $\mathbf{x}_p$ is

\begin{equation}
O^q(\mathbf{x}_p) = \{\mathbf{x}_p,\, \mathbf{f}(\mathbf{x}_p),\, \mathbf{f}^2(\mathbf{x}_p),\, \hdots \,,\, \mathbf{f}^{q-1}(\mathbf{x}_p)\}
\end{equation}

where the smallest $q$ satisfying Equation <ref> is the period of the orbit [p. 5]arrowsmith.

§.§ Stability

In simple terms, a given orbit $O(\mathbf{x})$ is said to be “Lyapunov stable” if nearby orbits will remain in a neighborhood of $O(\mathbf{x})$ holmes.[We define a neighborhood of a given orbit $O(\mathbf{x})$ to be an open set (a set that does not contain its boundary) containing $O(\mathbf{x})$.] To formalize this for a fixed point of a discrete-time dynamical system, we will use the definition outlined by bof: a stationary state $\mathbf{x}_s$ is Lyapunov stable if, for every $\epsilon>0$, there exists a $\delta>0$ such that for all $t>0$,

\begin{equation}
|\mathbf{x}_0-\mathbf{x}_s|<\delta \implies |\mathbf{f}^t(\mathbf{x}_0)-\mathbf{x}_s|<\epsilon
\label{eq:stable}
\end{equation}

\begin{equation}
|\mathbf{x}| = \sqrt{\left(x\e{1}\right)^2+\left(x\e{2}\right)^2+\hdots+\left(x\e{n}\right)^2}
\label{eq:magnitude}
\end{equation}

According to the definition by bof, a Lyapunov stable stationary state $\mathbf{x}_s$ is stable, or asymptotically stable, if $\delta$ can be chosen such that

\begin{equation}
|\mathbf{x}_0-\mathbf{x}_s|<\delta \implies \lim_{t\to\infty}\mathbf{f}^t(\mathbf{x}_0) =\mathbf{x}_s
\end{equation}

In other words, $\mathbf{x}_s$ is stable if and only if there is a neighborhood of $\mathbf{x}_s$ where all initial states will asymptotically approach $\mathbf{x}_s$ [p. 5]arrowsmith. If a stationary state is not stable, then we call it unstable.

We will say that $O^q(\mathbf{x}_p)$, a periodic orbit of period $q$, is Lyapunov stable if each point $\mathbf{x}\in O^q(\mathbf{x}_p)$ is a Lyapunov stable point of $\mathbf{f}^q$ [p. 23]martelli. A simple map with Lyapunov stable periodic orbits is the one-dimensional map with iteration function $f(x)=1-x$. In this map, every point except $x=0.5$ is in a Lyapunov stable 2-cycle. For example, $f(0)=1$ and $f(1)=0$, so $\{0,\,1\}$ is a periodic orbit with $q=2$. It is Lyapunov stable because $f^q(x) = f^2(x) = f(f(x)) = x$, which is obviously Lyapunov stable for all $x$.

§.§ Quantification of Chaos

[thick] (0, 0) circle [radius=1];
(0, 0) circle [radius=1.5pt];
[->, thick] (0,0) -> (0.866,0.5);
at (-0.2, -0.2) $\mathbf{x}_0$;
[-Stealth] (-1.5, 1) -> (0.433, 0.25);
[above left] at (-1.5, 1) $\frac{1}{2}p_1(0) = \frac{1}{2}p_2(0) = |d\mathbf{x}_0|$;
[->, very thick] (1.5, 0) -> (3, 0);
[above] at (2.25, 0) $\mathbf{f}^t(\mathbf{x}_0)$;
[rotate around=45:(4.5,0), thick] (4.5,0) ellipse (2 and 0.5);
(4.5, 0) circle [radius=1.5pt];
at (4.3, -0.2) $\mathbf{x}_t$;
[->, thick] (4.5, 0) -> (5.914, 1.414);
[->, thick] (4.5, 0) -> (4.854, -0.354);
[-Stealth] (6.5, 0.707) -> (5.207, 0.707);
[-Stealth] (6.5, -0.118) -> (4.618, -0.118);
[right] at (6.5, 0.707) $\frac{1}{2}p_1(t) \approx |d\mathbf{x}_0|e^{\lambda_1t}$;
[right] at (6.5, -0.118) $\frac{1}{2}p_2(t) \approx |d\mathbf{x}_0|e^{\lambda_2t}$;

Evolution of an initial infinitesimal one-dimensional sphere of perturbations from $\mathbf{x}_0$ by $d\mathbf{x}_0$ after $t$ iterations, with the semi-principal axes of the one-dimensional ellipsoid approximated by a modification of Equation <ref>

According to [p. 8]lorenz, the father of chaos theory, a dynamical system exhibits chaos if it is “sensitively dependent on initial conditions.” To get a better idea of what this means, we will turn to a quantitative definition of chaos provided by [p. 106]alligood: chaos is defined by “a Lyapunov exponent greater than zero.” The Lyapunov exponent is a quantity that characterizes the separation rate of close trajectories brandon. To define this precisely, let us first consider a simple one-dimensional discrete-time dynamical system with an initial state $x_0$ and an iteration function $x_{k+1}=f(x_k)$. If we perturb this initial state by some small amount $\delta x_0$, we can expect the separation of the initial and the perturbed states after some time $t$, denoted by $\delta x_t$, to be approximated by

\begin{equation}
|\delta x_t| = |f^t(x_0+\delta x_0)-f^t(x_0)| \approx |\delta x_0|e^{\lambda t}
\label{eq:lyapunovapprox}
\end{equation}

where $\lambda$ is the Lyapunov exponent. Here, we can see that if $\lambda>0$, any difference in initial conditions will be magnified as time goes on; even if the initial perturbation is infinitesimal, $\delta x_t$ will become significant as $t$ goes to infinity. This is how we quantify the chaos of a dynamical system: the more positive $\lambda$ is, the faster close trajectories will diverge from each other, so the more chaotic a system is. Similarly, the more negative $\lambda$ is, the faster close trajectories will converge on each other, so the more non-chaotic a system is wolf.

The Lyapunov exponent $\lambda$ is calculated from Equation <ref> as $t$ goes to infinity and $\delta x_0$ approaches 0:

\begin{equation}
\begin{split}
\lambda &= \lim_{t\to\infty}\lim_{\delta x_0\to 0}\frac{1}{t}\ln\left|\frac{\delta x_t}{\delta x_0}\right| \\
&= \lim_{t\to\infty}\lim_{\delta x_0\to 0}\frac{1}{t}\ln\left|\frac{f^t(x_0+\delta x_0)-f^t(x_0)}{\delta x_0}\right|
\label{eq:lyapunovunsimp}
\end{split}
\end{equation}

It is shown by brandon that this can be simplified to[See Appendix <ref> for a full derivation.]

\begin{equation}
\lambda = \lim_{t\to\infty}\frac{1}{t}\sum_{i=0}^{t-1}\ln |f'(x_i)|
\label{eq:1dlyap}
\end{equation}

which is simply the average of the logs of the derivatives taken at all the points visited by $O^+(x_0)$.

Now, let us consider an $n$-dimensional system with an initial state $\mathbf{x}_0$ and an iteration function $\mathbf{x}_{k+1} = \mathbf{f}(\mathbf{x}_k)$. Given some initial perturbation in a general direction

\begin{equation}
\delta \mathbf{x}_0 = \begin{pmatrix}
\delta x\e{1}_0 \\[4pt]
\delta x\e{2}_0 \\[1pt]
\vdots \\[1pt]
\delta x\e{n}_0
\end{pmatrix}
\label{eq:perturbation-gendirec}
\end{equation}

we can rewrite Equation <ref> for our $n$-dimensional system as

\begin{equation}
\begin{split}
\lambda &= \lim_{t\to\infty}\lim_{\delta x\e{1}_0,\, \delta x\e{2}_0,\, \hdots\,,\, \delta x\e{n}_0\to0}\frac{1}{t}\ln\frac{|\delta \mathbf{x}_t|}{|\delta \mathbf{x}_0|} \\
&= \lim_{t\to\infty}\frac{1}{t}\ln\frac{|d\mathbf{x}_t|}{|d\mathbf{x}_0|}
\label{eq:lyap-nd-unsimp}
\end{split}
\end{equation}

where $d\mathbf{x}_0$ is an infinitesimal initial perturbation and $d\mathbf{x}_t$ is the evolution of that perturbation after $t$ steps.

If we write Equation <ref> fully as

\begin{align}
\begin{split}
\begin{pmatrix}
x\e{1}_{k+1} \\[4pt]
x\e{2}_{k+1} \\[1pt]
\vdots \\[1pt]
\end{pmatrix} =
\begin{pmatrix}
f\e{1}(x\e{1}_k,\, x\e{2}_k,\,\hdots\,,\, x\e{n}_k) \\[4pt]
f\e{2}(x\e{1}_k,\, x\e{2}_k,\,\hdots\,,\, x\e{n}_k) \\[1pt]
\vdots \\[1pt]
f\e{n}(x\e{1}_k,\, x\e{2}_k,\,\hdots\,,\, x\e{n}_k)
\end{pmatrix}
\end{split}
\label{eq:expandediteration}
\end{align}

where $f\e{m}$ is the iteration function corresponding to the $m$th dimension of $\mathbf{f}$, then we can write the Jacobian matrix $J(\mathbf{x})$ of the system as

\begin{equation}
J(\mathbf{x}) = \begin{pmatrix}
\frac{\partial f\e{1}}{\partial x\e{1}} & \frac{\partial f\e{1}}{\partial x\e{2}} & \hdots & \frac{\partial f\e{1}}{\partial x\e{n}} \\[4pt]
\frac{\partial f\e{2}}{\partial x\e{1}} & \frac{\partial f\e{2}}{\partial x\e{2}} & \hdots & \frac{\partial f\e{2}}{\partial x\e{n}} \\[1pt]
\vdots & \vdots & \ddots & \vdots \\[1pt]
\frac{\partial f\e{n}}{\partial x\e{1}} & \frac{\partial f\e{n}}{\partial x\e{2}} & \hdots & \frac{\partial f\e{n}}{\partial x\e{n}}
\end{pmatrix}
\label{eq:jacobian}
\end{equation}

It is shown by brandon that an infinitesimal perturbation can be iterated by the Jacobian matrix:[See Appendix <ref> for a full derivation.]

\begin{equation}
d\mathbf{x}_{k+1} = J(\mathbf{x}_k)d\mathbf{x}_k
\label{eq:jacobian-iteration}
\end{equation}

It follows that

\begin{equation}
d\mathbf{x}_{t} = J(\mathbf{x}_{t-1})J(\mathbf{x}_{t-2})\hdots J(\mathbf{x}_0)d\mathbf{x}_0
= J^td\mathbf{x}_0
\label{eq:jacobian-alltheway}
\end{equation}

Rewriting Equation <ref> by substituting Equation <ref>, we get that

\begin{equation}
\begin{split}
\lambda &= \lim_{t\to\infty}\frac{1}{t}\ln\frac{|J^td\mathbf{x}_0|}{|d\mathbf{x}_0|} \\
&= \lim_{t\to\infty}\frac{1}{t}\ln\left|J^t\mathbf{u}_0\right|
\end{split}
\label{eq:lyapunitvector}
\end{equation}

where $\mathbf{u}_0$ is a unit vector in the direction of $d\mathbf{x}_0$.

In an $n$-dimensional system, there is a spectrum of Lyapunov exponents $\lambda = \{\lambda_1,\, \lambda_2,\,\hdots\,,\, \lambda_n\}$ ordered from largest to smallest. In the spectrum, each Lyapunov exponent corresponds to the rate of separation of trajectories for a specific starting direction of $\mathbf{u}_0$ wolf. Geometrically, we can imagine all possible $d\mathbf{x}_0$ forming an $(n-1)$-dimensional sphere[In this paper, we say that a one-dimensional sphere is a circle, a two-dimensional sphere is a standard sphere, and so on. Topologically, this is because a circle “looks like” a one-dimensional line up close, a sphere “looks like” a two-dimensional plane up close, and so on.] with infinitesimal radius centered around $\mathbf{x}_0$. As time evolves, this $(n-1)$-dimensional sphere will be deformed into an $(n-1)$-dimensional ellipsoid with $n$ principal axes (see Figure <ref>). If we say $p_i(t)$ represents the length of the $i$th principal axis[The 1st principal axis is the longest, the 2nd principal is the second longest, and so on.] at time $t$, then the Lyapunov exponent $\lambda_i$ corresponds to the growth of $p_i(t)$. Then, from Equation <ref>,

\begin{equation}
\lambda_i = \lim_{t\to\infty}\frac{1}{t}\ln\frac{p_i(t)}{p_i(0)} = \lim_{t\to\infty}\frac{1}{t}\ln\frac{p_i(t)}{2|d\mathbf{x}_0|}
\label{eq:lyap-prin-axes}
\end{equation}

because the radius of a circle is half the length of its principal axis. For an arbitrarily chosen direction of $\mathbf{u}_0$, the perturbation will grow in magnitude according to the maximal Lyapunov exponent $\lambda_1$ because of the exponential nature of Lyapunov exponents brandon. Geometrically, this is because, after a long time, the 1st principal axis of the ellipse will dwarf all of the other principal axes in comparison (see Figure <ref>). Any $d\mathbf{x}_0$ that has a non-zero component in the direction of the 1st principal axis will eventually be overwhelmed by the growth of that component, resulting in $d\mathbf{x}_t$ approaching the direction of the 1st principal axis and growing according to $\lambda_1$.

While Equation <ref> provides good intuition for what is going on, it isn't practical to use it to calculate the Lyapunov spectrum. Instead, brandon shows that the Lyapunov spectrum can be manipulated from Equation <ref> into the form[See Appendix <ref> for a full derivation.]

\begin{equation}
\lambda_i = \lim_{t\to\infty}\frac{1}{2t}\ln\mu_i
\label{eq:lyapunov-eigenvalue}
\end{equation}

where $\mu_i$ is an eigenvalue[For the unfamiliar reader, if $M\mathbf{v} = m\mathbf{v}$ for some matrix $M$, some vector $\mathbf{v}$, and some scalar $m$, $m$ is an eigenvalue of $M$ and $\mathbf{v}$ is its associated eigenvector. [pp. 268-273]linear provide an excellent in-depth explanation of this concept.] of the matrix $J^{t\intercal}J^t$.[$M^{\intercal}$ is the transpose of $M$.] These eigenvalues $\mu_i$ are labelled with their subscript $i$ so that the Lyapunov exponents satisfy $\lambda_1\geq\lambda_2\geq\hdots\geq\lambda_n$. Mirroring the aforementioned definition from [p. 106]alligood, at least one Lyapunov exponent being greater than 0 indicates chaos wolf.

§ THE GEOMETRY OF DYNAMICAL AND CHAOTIC SYSTEMS

As defined in Section <ref>, a chaotic system is a dynamical system exhibiting sensitive dependence on initial conditions. We have already discussed the type of sensitivity to initial conditions known as “temporal” chaos, which emerges from a system's dynamics and can be quantified using Lyapunov exponents. However, another type of sensitivity to initial conditions, a “geometrical” form of chaos, arises in some systems exhibiting a phenomenon known as multistability.

In general, dynamical systems can have many interesting geometrical properties, including attractors, basins of attraction, and fractal basin boundaries. We are interested in these geometric structures because some of them are strongly connected with these “temporal” and “geometrical” forms of chaos. In this section, we will explore some of the geometries of dynamical and chaotic systems and methods of detecting, classifying, and quantifying them.

§.§ Non-Chaotic Attractors

An attractor, as defined by [p. 332]strogatz, is “a set to which all neighboring trajectories converge.” In more specific terms, an attractor can be defined as a set of points $A$ in state space that satisfies the following three properties:

* Trajectories that start in $A$ will stay in $A$. Specifically, any forward orbit $O^+(\mathbf{x}_0)$ with $\mathbf{x}_0\in A$ will satisfy $\mathbf{x}_t\in A$ for all $t>0$.

* $A$ will attract an open set of initial conditions. Specifically, say that the distance from an arbitrary state $\mathbf{x}$ to $A$ is the smallest $|\mathbf{x} - \mathbf{a}|$ such that $\mathbf{a}\in A$. Then, there is some open set $U$ that $A$ is a subset of, or $A\subseteq U$,[The symbol $\subseteq$ means “is a subset of.”] such that if some $\mathbf{x}_0\in U$, then the distance from $\mathbf{x}_t$ to $A$ approaches 0 as $t\to\infty$. Essentially, all forward orbits that start sufficiently close to $A$ will be attracted to $A$.

* $A$ is as small as possible. Specifically, this means there is no set $A'$ satisfying $A'\subset A$[The symbol $\subset$ means “is a proper subset of.” $X\subset Y$ if and only if $X\subseteq Y$ but $X\neq Y$.] that also satisfies the previous conditions.

Attractors are key geometrical features of dynamical systems because, by definition, they represent the typical behavior of a set of a dynamical system's state space. We will begin our discussion of attractors with non-chaotic attractors, that is, attractors that don't exhibit sensitivity to initial conditions.

[thick, domain=-1:1]plot(,1-()^2);
[fill](0, 1.05) circle [radius=0.05];
[above] at (0, 1.05) $\mathbf{x}_s$;

Ball at the top of a hill, an example of a fixed point repeller

[thick, domain=-1:1]plot(,()^2);
[fill](0, 0.05) circle [radius=0.05];
[above] at (0, 0.05) $\mathbf{x}_s$;

Ball at the bottom of a valley, an example of a fixed point attractor

Two systems with fixed points

[dashed, thick] (0,0) circle [radius=3];
[above right] at (2.5,2.5) $U$;
[-Stealth, semithick](2.5,2.5)–(1.5,1.5);
[fill] (0,1.5) circle [radius=0.05];
[above] at (0, 1.5) $\mathbf{x}_p$;
[->, semithick](0.25,1.25)–(1.25,0.25);
[above right] at (0.75,0.75) $\mathbf{f}$;
[fill] (1.5,0) circle [radius=0.05];
[right] at (1.5,0) $\mathbf{f}(\mathbf{x}_p)$;
[<-, semithick](0.25,-1.25)–(1.25,-0.25);
[below right] at (0.75,-0.75) $\mathbf{f}$;
[fill] (0,-1.5) circle [radius=0.05];
[below] at (0, -1.5) $\mathbf{f}^2(\mathbf{x}_p)$;
[->, semithick](-0.25,-1.25)–(-1.25,-0.25);
[below left] at (-0.75,-0.75) $\mathbf{f}$;
[fill] (-1.5,0) circle [radius=0.05];
[left] at (-1.5,0) $\mathbf{f}^3(\mathbf{x}_p)$;
[->, semithick] (-1.25,0.25)–(-0.25,1.25);
[above left] at (-0.75,0.75) $\mathbf{f}$;

A connected set $U$

[fill] (0,2) circle [radius=0.05];
[above] at (0, 2) $\mathbf{x}_p$;
[dashed, thick] (0,2) circle [radius=1.2];
[-Stealth, semithick](1.25,3.25)–(0.45,2.45);
[above right] at (1.2,3.2) $U_0$;
[->, thin](0,2)–(0.849,1.151);
[above right] at (0.35, 1.55) $\epsilon$;

[fill] (2,0) circle [radius=0.05];
[above] at (2,0) $\mathbf{f}(\mathbf{x}_p)$;
[dashed, thick] (2,0) circle [radius=1.2];
[-Stealth, semithick](3.25,1.25)–(2.45,0.45);
[above right] at (3.2,1.2) $U_1$;
[->, thin] (2,0)–(2.849,-0.849);
[above right] at (2.35, -0.45) $\epsilon$;

[fill] (0,-2) circle [radius=0.05];
[above] at (0, -2) $\mathbf{f}^2(\mathbf{x}_p)$;
[dashed, thick] (0,-2) circle [radius=1.2];
[-Stealth, semithick](1.75, -2)–(0.65,-2);
[right] at (1.75,-2) $U_2$;
[->, thin] (0,-2)–(0.849,-2.849);
[above right] at (0.35, -2.45) $\epsilon$;

[fill] (-2,0) circle [radius=0.05];
[above] at (-2,0) $\mathbf{f}^3(\mathbf{x}_p)$;
[dashed, thick] (-2,0) circle [radius=1.2];
[-Stealth, semithick](-2, 1.75)–(-2,0.65);
[above] at (-2,1.75) $U_3$;
[->, thin] (-2,0)–(-1.151,-0.849);
[above right] at (-1.65,-0.45) $\epsilon$;

A disconnected set $U = U_0\cup U_1\cup U_2\cup U_3$

Two possible open sets of initial conditions $U$ that 4-cycle attractor $A = \{\mathbf{x}_p,\, \mathbf{f}(\mathbf{x}_p),\, \mathbf{f}^2(\mathbf{x}_p),\, \mathbf{f}^3(\mathbf{x}_p)\}$ attracts

§.§.§ Fixed Points

As previously defined in Section <ref> (Equation <ref>), a state $\mathbf{x}_s$ is a fixed point of $\mathbf{f}$ if it satisfies $\mathbf{f}(\mathbf{x}_s) = \mathbf{x}_s$. A fixed point $\mathbf{x}_s$ is also an attractor if it satisfies all of the previously listed properties.

To get an intuitive understanding of what makes a given fixed point an attractor or not, we will consider two example systems. First, let us consider a ball on top of a hill (see Figure <ref>). The point at the top of this hill $\mathbf{x}_s$ is our potential attractor. This is clearly a fixed point by Equation <ref>: as time evolves, the state of the system won't change. For the same reason, Property 1 of attractors is satisfied. Now, let us consider an open set that contains $\mathbf{x}_s$. If a ball starts in a neighborhood of $\mathbf{x}_s$, it will roll down the hill and be repelled from $\mathbf{x}_s$ no matter how small that neighborhood is. This means that there is no open set that satisfies the requirements of Property 2, so $\mathbf{x}_s$ is not an attractor. In fact, because states that start in a neighborhood of $\mathbf{x}_s$ are repelled from it, we call $\mathbf{x}_s$ a repeller [p. 17]strogatz.

Our second system is a ball at the bottom of a valley (see Figure <ref>). Again, this point $\mathbf{x}_s$ is clearly a fixed point and satisfies Property 1. However, if we now consider open sets that contain our fixed point, we can see that a ball that starts in a neighborhood of $\mathbf{x}_s$ will oscillate around $\mathbf{x}_s$ before eventually settling at it. Therefore, $\mathbf{x}_s$ satisfies Property 2: there exists a neighborhood of the point such that initial states in that neighborhood are attracted to it. Finally, notice that our attractor $A$ contains only one point: $A = \{\mathbf{x}_s\}$. There exists no proper subset of $A$ besides the empty set, which is obviously not an attractor, so Property 3 is automatically satisfied. Therefore, $\mathbf{x}_s$ is an attractor.

We prove in Appendix <ref> that there is a simple set of criteria to determine whether a fixed point is an attractor, a repeller, or a saddle point.[Saddle points attract some states in a small neighborhood of it and repel others in that same neighborhood. They only exist in systems higher than one dimension.] For a one-dimensional system $x_{k+1} = f(x_k)$ with a fixed point $x_s$, the criteria are as follows:

* $x_s$ is an attractor if $|f'(x_s)|<1$.

* $x_s$ is a repeller if $|f'(x_s)|>1$.

For an $n$-dimensional system $\mathbf{x}_{k+1} = \mathbf{f}(\mathbf{x}_k)$ with a fixed point $\mathbf{x}_s$,

* $\mathbf{x}_s$ is an attractor if the absolute value of each eigenvalue of the Jacobian matrix at $\mathbf{x}_s$ is less than 1. That is, $|\nu_i|<1$ for $i = 1,\,2,\,\hdots\,,\,n$.

* $\mathbf{x}_s$ is a repeller if $|\nu_i|>1$ for $i = 1,\,2,\,\hdots\,,\,n$.

* $\mathbf{x}_s$ is a saddle point if at least one $|\nu_i|>1$ and at least one $|\nu_i|<1$ for $i = 1,\,2,\,\hdots\,,\,n$.

§.§.§ Periodic Orbits

Consider a periodic orbit of the periodic point $\mathbf{x}_p$ that has a period $q$:

\begin{equation}
O^q(\mathbf{x}_p) = \{\mathbf{x}_p,\, \mathbf{f}(\mathbf{x}_p),\, \mathbf{f}^2(\mathbf{x}_p),\, \hdots \,,\, \mathbf{f}^{q-1}(\mathbf{x}_p)\}
\end{equation}

Let us consider what has to be true for this periodic orbit to be attractor $A = O^q(\mathbf{x}_p)$. Of course, Property 1 is immediately satisfied by the definition of a periodic orbit established in Section <ref>. By Property 2, in order for $A$ to be an attractor, there must be an open set such that any forward orbit $O^+(\mathbf{x}_0)$ that starts in this set contains an orbit $O(\mathbf{x})$ that is arbitrarily close to $O^q(\mathbf{x}_p)$. In other words, a state that starts in the open set will eventually fall into a periodic orbit that approaches $O^q(\mathbf{x}_p)$. A simple possibility for an open set of initial conditions that $A$ attracts is a connected set[A connected set is a set where every point in the set can be reached by a continuous path from every other point in the set without leaving it.] that contains all the points in the periodic orbit (see Figure <ref>). Another possibility is a disconnected set composed of neighborhoods of all the points visited by $O^q(\mathbf{x}_p)$ (see Figure <ref>). For example, if we define these neighborhoods as open sets $U_i$ that are the sets of states $\mathbf{x}$ that satisfy $|\mathbf{x} - \mathbf{f}^i(\mathbf{x}_p)| < \epsilon$, where $i = 0,\, 1,\, 2,\,\hdots\,,\,q-1$, then our open set of initial conditions is $U = U_0\cup U_1\cup U_2\cup\hdots\cup U_{q-1}$.[$S_1\cup S_2$ is the union of $S_1$ and $S_2$, which is the set of all the elements contained in $S_1$, $S_2$, or both $S_1$ and $S_2$.] Finally, because we defined the period $q$ of our periodic orbit $O^q(\mathbf{x}_p)$ in Section <ref> to be the smallest $q$ satisfying $\mathbf{f}^q(\mathbf{x}_p) = \mathbf{x}_p$ (Equation <ref>), any $A'$ will be missing at least one point in $O^q(\mathbf{x}_p)$. A trajectory that starts in $A'$ will leave $A'$ when it reaches one of these missing points, so $A'$ fails Property 1. Therefore, a periodic orbit $O^q(\mathbf{x}_p)$ that satisfies Property 2 is an attractor. We will delve deeper into and see examples of periodic orbit attractors in Section <ref>.

[<->] (0,2)–(0,-2);
[<->] (2,0)–(-2,0);
[thick, domain=-1:1]plot(,-1-()^2);
[dashed, domain=-2:2]plot(,);
[above] at (0,2) $x_{k+1}$;
[right] at (2,0) $x_{k}$;

[<->] (0,2)–(0,-2);
[<->] (2,0)–(-2,0);
[thick, domain=-1.323:1.323]plot(,-1/4-()^2);
[dashed, domain=-2:2]plot(,);
[above] at (0,2) $x_{k+1}$;
[right] at (2,0) $x_{k}$;

[<->] (0,2)–(0,-2);
[< ->] (2,0)–(-2,0);
[thick, domain=-1.581:1.581]plot(,1/2-()^2);
[dashed, domain=-2:2]plot(,);
[above] at (0,2) $x_{k+1}$;
[right] at (2,0) $x_{k}$;

Graphs of the quadratic function $x_{k+1} = r-x_k^2$ for values of $r$ before, at, and after a saddle-node bifurcation, with fixed points at the intersection of the function with the line $x_{k+1} = x_k$

§.§ Bifurcations

One of the interesting properties of some dynamical systems is a dependence on parameter. Namely, say a map $\mathbf{x}_{k+1} = \mathbf{f}(\mathbf{x}_k;\,r)$ is dependent on some parameter $r$, where we use the semicolon to separate changing variables from fixed parameters. For most values of $r$, varying $r$ slightly will cause the system's dynamics to change quantitatively, such as a fixed point changing its location slightly. Sometimes, however, the qualitative geometry of a system can change as $r$ is varied. Specifically, there can be an appearance or disappearance of stable orbits, a change in stability, or the emergence of complex geometric structures like strange attractors crawford.[See Section <ref> for details on strange attractors.] Bifurcations are defined as these qualitative changes in a system's dynamics, and bifurcation points are the values of the system's parameters at which bifurcations occur [p. 45]strogatz. Bifurcation theory is a very rich and complex field,[For the interested reader, crawford provides a far more comprehensive introduction to bifurcation theory than we give in this paper.] and there are many different types of bifurcations, but we will focus on two simple ones in this section: saddle-node bifurcations and period-doubling bifurcations.

§.§.§ Saddle-Node Bifurcations

As an example of a system that exhibits bifurcations, let us explore the one-dimensional $r$-dependent quadratic map with iteration function

\begin{equation}
f(x;\,r) = r-x^2
\label{eq:quadratic-map}
\end{equation}

To analyze the qualitative changes in the dynamics of this system, let us first examine the fixed points of this map and their stability. By Equation <ref>, the fixed points of this map satisfy

\begin{equation}
x_s = f(x_s;\,r) = r-x_s^2
\end{equation}

Rearranging and using the quadratic formula, we get that

\begin{equation}
x_s(r) = \frac{-1\pm\sqrt{1+4r}}{2}
\label{eq:fixedpoint-quadratic-eq}
\end{equation}

where we write $x_s$ as a function of $r$ to emphasize its $r$-dependence. By analysis of the discriminant, we can see that there are no fixed points for $r<-1/4$, one fixed point for $r=-1/4$, and two fixed points for $r>-1/4$. Here, as we continuously increase the parameter $r$, two fixed points appear at the bifurcation point $r=-1/4$ in a region where there were no fixed points before. This type of bifurcation is known as a saddle-node bifurcation, named after the appearance of this bifurcation in higher-dimensional systems [p. 48]strogatz.

We can see how these fixed points appear graphically by putting $x_k$ on the horizontal axis and $x_{k+1}$ on the vertical axis. Then, the intersections between the parabola $x_{k+1} = f(x_k;\,r)$ and the line $x_{k+1} = x_k$ are where the fixed points are located. This is shown in Figure <ref>, where $x_{k+1} = f(x_k;\,r)$ is plotted for three values of $r$. In Figure <ref>, $x_{k+1} = f(x_k;\,r)$ doesn't intersect $x_{k+1} = x_k$ anywhere, so there are no fixed points for $r=-1<-1/4$. In Figure <ref>, $x_{k+1} = f(x_k;\,r)$ intersects $x_{k+1} = x_k$ at one point, the point of tangency, and we can see from Equation <ref> that this fixed point is at

\begin{equation}
x_s\left(-\frac{1}{4}\right) = \frac{-1\pm\sqrt{1+4\left(-\frac{1}{4}\right)}}{2} = -\frac{1}{2}
\label{eq:quadratic-saddlenode-bifurcationpoint}
\end{equation}

Finally, in Figure <ref>, $x_{k+1} = f(x_k;\,r)$ intersects $x_{k+1} = x_k$ at two points, so from Equation <ref>, the two fixed points are located at

\begin{equation}
\begin{split}
x_{s,\,1}\left(\frac{1}{2}\right) &= \frac{-1+\sqrt{1+4\left(\frac{1}{2}\right)}}{2} = \frac{-1+\sqrt{3}}{2}\\
x_{s,\,2}\left(\frac{1}{2}\right) &= \frac{-1-\sqrt{1+4\left(\frac{1}{2}\right)}}{2} = \frac{-1-\sqrt{3}}{2}
\end{split}
\end{equation}

where we define the functions

\begin{align}
x_{s,\,1}(r) &= \frac{-1+\sqrt{1+4r}}{2} \label{eq:quadratic-fixed-point-function1} \\
x_{s,\,2}(r) &= \frac{-1-\sqrt{1+4r}}{2} \label{eq:quadratic-fixed-point-function2}
\end{align}

as the locations of the two fixed points that exist for $r>-1/4$.

Bifurcations are also concerned with the stability of these fixed points, namely, whether they are attractors or repellers. Let us first consider $x_{s,\,2}(r)$. By the criteria for the attractiveness of fixed points in one dimension established in Section <ref> and Appendix <ref>, the attractiveness of $x_{s,\,2}(r)$ is dependent on $|f'(x_{s,\,2};\,r)|$. Since $f'(x;\,r) = -2x$,

\begin{equation}
\begin{split}
|f'(x_{s,\,2};\,r)| &= \left|-2\left(\frac{-1-\sqrt{1+4r}}{2}\right)\right| \\
&= 1+\sqrt{1+4r}
\end{split}
\label{eq:quadratic-fixedpoint2}
\end{equation}

This is obviously greater than 1 for all $r>-1/4$, so on this interval, $x_{s,\,2}$ is a repeller.

Now, let us consider $x_{s,\,1}$. By the criteria for the attractiveness of fixed points in one dimension, $x_{s,\,1}$ is an attractor when $|f'(x_{s,\,1};\,r)|<1$ and a repeller when $|f'(x_{s,\,1};\,r)|>1$. Let us first consider for what $r$ values $x_{s,\,1}$ is an attractor:

\begin{equation}
\begin{gathered}
|f'(x_{s,\,1};\,r)|<1 \\
-1 < f'(x_{s,\,1};\,r) < 1 \\
-1 < -2\left(\frac{-1+\sqrt{1+4r}}{2}\right) < 1 \\
-\frac{1}{4} < r < \frac{3}{4} \\
\end{gathered}
\label{eq:quadratic-fixedpoint1-attractor}
\end{equation}

Similarly, the values of $r$ where $x_{s,\,1}$ is a repeller are

\begin{equation}
\begin{gathered}
|f'(x_{s,\,1};\,r)|>1 \\
\begin{aligned}
f'(x_{s,\,1};\,r) &< -1 \hspace{20px}& f'(x_{s,\,1};\,r) &> 1 \\
1-\sqrt{1+4r} &< -1 & 1-\sqrt{1+4r} &> 1 \\
\sqrt{1+4r} &> 2 & \sqrt{1+4r} &< 0 \\
r &> \frac{3}{4}
\end{aligned}
\end{gathered}
\label{eq:quadratic-fixedpoint1-repeller}
\end{equation}

In summary, $x_{s,\,2}$ is a repeller for $r > -1/4$, while $x_{s,\,1}$ starts as an attractor for $-1/4 < r < 3/4$ then turns into a repeller for $r > 3/4$.

[->] (-0.75,-2.75)–(2.75,-2.75);
[->] (-0.75,-2.75)–(-0.75,1.75);
[thick, domain=-0.248:0.72]plot(,(-1+sqrt(1+4*))/2);
[thick, dashed, domain=0.78:2.25]plot(,(-1+sqrt(1+4*))/2);
[thick, dashed, domain=-0.248:2.25]plot(,(-1-sqrt(1+4*))/2);
[thick] (-1/4,-1/2) circle [radius=0.05];
[thick] (3/4,1/2) circle [radius=0.05];
[right] at (-1/4,-1/2) $\left(-1/4,\, -1/2\right)$;
[below right] at (3/4,1/2) $\left(3/4,\, 1/2\right)$;
[below] at (1,-3.25) $r$;
[left] at (-1.25,-0.5) $x$;
(0, -2.75)–(0, -2.85);
(1, -2.75)–(1, -2.85);
(2, -2.75)–(2, -2.85);
[below] at (0, -2.85) $0$;
[below] at (1, -2.85) $1$;
[below] at (2, -2.85) $2$;
(-0.75, -2)–(-0.85, -2);
(-0.75, -1)–(-0.85, -1);
(-0.75, 0)–(-0.85, 0);
(-0.75, 1)–(-0.85, 1);
[left] at (-0.85, -2) $-2$;
[left] at (-0.85, -1) $-1$;
[left] at (-0.85, 0) $0$;
[left] at (-0.85, 1) $1$;

Bifurcation diagram showing the fixed point attractors (solid curve) and repellers (dashed curves) for the quadratic map $f(x;\,r)=r-x^2$, with a saddle-node bifurcation at $(-1/4,\,-1/2)$

In Figure <ref>, we plot a bifurcation diagram of this system, which shows the locations of the attractors and repellers of the system as a function of $r$. Specifically, the solid curves represent fixed point attractors, while the dashed curves represent fixed point repellers. We make the bifurcation diagram by graphing Equation <ref> (the top half of the parabola) and Equation <ref> (the bottom half of the parabola). We know the saddle-node bifurcation is at $(-1/4,\,-1/2)$ by Equation <ref>, and we know where the fixed points are attracting and repelling by Equations <ref>, <ref>, and <ref>.

[<->] (0,2)–(0,-2);
[<->] (2,0)–(-2,0);
[thick, domain=-1.443:1.443, samples=100]plot(,-()^4+2*(1/2)*()^2-(1/2)^2+1/2);
[dashed, domain=-2:2]plot(,);
[above] at (0,2) $x_{k+2}$;
[right] at (2,0) $x_{k}$;

[<->] (0,2)–(0,-2);
[<->] (2,0)–(-2,0);
[thick, domain=-1.552:1.552, samples=100]plot(,-()^4+2*(3/4)*()^2-(3/4)^2+3/4);
[dashed, domain=-2:2]plot(,);
[above] at (0,2) $x_{k+2}$;
[right] at (2,0) $x_{k}$;

[<->] (0,2)–(0,-2);
[< ->] (2,0)–(-2,0);
[thick, domain=-1.653:1.653, samples=100]plot(,-()^4+2*(1)*()^2-(1)^2+1);
[dashed, domain=-2:2]plot(,);
[above] at (0,2) $x_{k+2}$;
[right] at (2,0) $x_{k}$;

Graphs of the second iterate of the quadratic map $x_{k+1} = r-x_k^2$, namely, the quartic function $x_{k+2} = -x_k^4+2rx_k^2-r^2+r$, for values of $r$ before, at, and after a period-doubling bifurcation, with fixed points and periodic points of period 2 at the intersection of the function with the line $x_{k+2} = x_k$

§.§.§ Period-Doubling Bifurcations

We can see from Figure <ref> that $x_{s,\,1}$ changes from an attractor to a repeller at $r=3/4$ since $f'(x_{s,\,1}(3/4);\,3/4) = -1$. This is clearly a qualitative change in the system's dynamics, so it must be another kind of bifurcation. To see what this bifurcation is, we need to explore the attractiveness of 2-cycles in this system. To do this, we consider the mapping function that iterates $f(x;\,r)$ twice:

\begin{equation}
\begin{split}
f^2(x;\,r) &= r-\left(r-x^2\right)^2 \\
&= -x^4+2rx^2-r^2+r
\end{split}
\label{eq:iterate-quadratic-twice}
\end{equation}

By Equation <ref>, periodic points of 2-cycles will follow

\begin{equation}
x_p = f^2(x_p;\,r) = -x_p^4+2rx_p^2-r^2+r
\end{equation}

Numerically, we can determine that the four solutions to this quartic equation are

\begin{align}
x_{p,\,{1}}(r) &= \frac{1+\sqrt{-3+4r}}{2} \label{eq:quadratic-periodicpoint-function1} \\
x_{p,\,{2}}(r) &= \frac{1-\sqrt{-3+4r}}{2} \label{eq:quadratic-periodicpoint-function2} \\
x_{p,\,{3}}(r) &= \frac{-1+\sqrt{1+4r}}{2} \label{eq:quadratic-periodicpoint-function3} \\
x_{p,\,{4}}(r) &= \frac{-1-\sqrt{1+4r}}{2} \label{eq:quadratic-periodicpoint-function4}
\end{align}

where $x_{p,\,1}(r)$ and $x_{p,\,2}(r)$ are two periodic points in a 2-cycle and $x_{p,\,3}(r)$ and $x_{p,\,4}(r)$ are stationary points of $f(x;\,r)$ since, by Equations <ref> and <ref>, $x_{p,\,{3}}(r) = x_{s,\,{1}}(r)$ and $x_{p,\,{4}}(r) = x_{s,\,{2}}(r)$. This makes sense because a stationary point of $f(x;\,r)$ is also a stationary point of $f^2(x;\,r)$.

These four periodic points can be seen graphically in a very similar manner to Figure <ref>: by graphing $x_{k+2} = f^2(x_k;\,r)$ and looking at where it intersects with $x_{k+2}=x_k$, the line of all periodic points of period 2. We show this in Figure <ref> by plotting the second iterate of the quadratic map for three values of $r$. In Figure <ref>, we plot before the second bifurcation shown in Figure <ref>. The two points of intersection are the fixed points $x_{p,\,{3}}(1/2)$ and $x_{p,\,{4}}(1/2)$, which match up graphically with $x_{s,\,{1}}(1/2)$ and $x_{s,\,{2}}(1/2)$ in Figure <ref>. In Figure <ref>, we plot at the bifurcation, and we can see that there is a qualitative change in $f^2(x_k;\,r)$ at the new point of tangency. Finally, after the bifurcation (Figure <ref>), we can see there are now four points of intersection, two of which are $x_{p,\,{1}}(1)$ and $x_{p,\,{2}}(1)$.

Using Equations <ref> and <ref>, we can determine that these points $x_{p,\,{1}}(1)$ and $x_{p,\,{2}}(1)$ are

\begin{equation}
\begin{split}
x_{p,\,1}(1) &= \frac{1+\sqrt{-3+4(1)}}{2} = 1 \\
x_{p,\,2}(1) &= \frac{1-\sqrt{-3+4(1)}}{2} = 0
\end{split}
\end{equation}

These points clearly make up a 2-cycle: $f\left(x_{p,\,1}(1);\,1\right) = 1-1^2 = 0 = x_{p,\,2}$ and $f\left(x_{p,\,2}(1);\,1\right) = 1-0^2 = 1 = x_{p,\,1}$. To determine whether or not this 2-cycle is an attractor or a repeller, we have to consider $\left|\left(f^2\right)'\left(x_{p,\,{1}}(1);\,1\right)\right|$ and $\left|\left(f^2\right)'\left(x_{p,\,{2}}(1);\,1\right)\right|$. Taking the derivative of Equation <ref>, we get that

\begin{equation}
\left(f^2\right)'(x;\,r) = -4x^3+4rx
\end{equation}

\begin{equation}
\begin{split}
\left(f^2\right)'\left(x_{p,\,1}(1);\,1\right) &= -4(1)^3+4(1)(1) = 0 \\
\left(f^2\right)'\left(x_{p,\,2}(1);\,1\right) &= -4(0)^3+4(1)(0) = 0 \\
\end{split}
\end{equation}

which is also clear graphically in Figure <ref>. Then, by the criteria for fixed point attractiveness established in Section <ref>,[We can use the criteria for fixed point attractiveness here because these periodic points are fixed points in $f^2(x;\,r)$.] $x_{p,\,{1}}(1)$ and $x_{p,\,{2}}(1)$ are both attractors in $f^2(x;\,r)$, meaning the 2-cycle $\{x_{p,\,1}(1),\,x_{p,\,2}(1)\}$ is an attractor.

Now, we would like to establish whether the 2-cycles for all values of $r$ are attractors or repellers. To make this easier, we realize that since $\left(f^2\right)'(x;\,r)$ is a polynomial, it changes continuously as $r$ changes. Therefore, instead of going through the inequality work we did for the map's fixed points, we can establish where the derivative's absolute value equals $1$, then check points on either side of these $r$ values for attractiveness. Considering $x_{p,\,1}(r)$ first, we can solve $\left|\left(f^2\right)'\left(x_{p,\,{1}}(r);\,r\right)\right| = 1$ by substituting in Equation <ref> and solving for $r$, which yields $r=3/4,\,5/4$. The same values of $r$ satisfy $\left|\left(f^2\right)'\left(x_{p,\,{2}}(r);\,r\right)\right| = 1$, which makes sense because a 2-cycle won't have one attracting fixed point and one repelling fixed point. Because the 2-cycle changes stability at $r=3/4$ and $r=5/4$, we want to consider their attractiveness on the intervals $r<3/4$, $3/4<r<5/4$, and $r>5/4$. First, we can see from Equations <ref> and <ref> that there are no 2-cycles when $r<3/4$. We have already done a test for a 2-cycle with an $r$ between $3/4$ and $5/4$, namely $r=1$, which we determined was an attractor. This indicates that 2-cycles for $3/4<r<5/4$ are attractors. Finally, since $\left|\left(f^2\right)'\left(x_{p,\,{1,\,2}}(r);\,r\right)\right|$ crosses $1$ at $r=5/4$, 2-cycles for $r>5/4$ must be repellers, which we can easily verify using any test value of $r>5/4$.

[->] (-0.75,-2.75)–(2.75,-2.75);
[->] (-0.75,-2.75)–(-0.75,1.75);
[thick, domain=-0.248:0.72, samples=100]plot(,(-1+sqrt(1+4*))/2);
[thick, dashed, domain=0.78:2, samples=100]plot(,(-1+sqrt(1+4*))/2);
[thick, dashed, domain=-0.248:2, samples=100]plot(,(-1-sqrt(1+4*))/2);
[thick, domain=0.753:1.2, samples=100]plot(,(1+sqrt(-3+4*))/2);
[thick, domain=0.753:1.2, samples=100]plot(,(1-sqrt(-3+4*))/2);
[thick, dashed, domain=1.3:2.25, samples=100]plot(,(1+sqrt(-3+4*))/2);
[thick, dashed, domain=1.3:2.25, samples=100]plot(,(1-sqrt(-3+4*))/2);
[thick] (-1/4,-1/2) circle [radius=0.05];
[thick] (3/4,1/2) circle [radius=0.05];
[thick] (5/4,1.207) circle [radius=0.05];
[thick] (5/4,-0.207) circle [radius=0.05];
[right] at (-1/4,-1/2) $\left(-1/4,\, -1/2\right)$;
[right] at (3/4,0.35) $\left(3/4,\, 1/2\right)$;
[right] at (5/4,1.15) $(5/4,\, (1+\sqrt{2})/2)$;
[right] at (5/4,-0.15) $(5/4,\, (1-\sqrt{2})/2)$;
[below] at (1,-3.25) $r$;
[left] at (-1.25,-0.5) $x$;
(0, -2.75)–(0, -2.85);
(1, -2.75)–(1, -2.85);
(2, -2.75)–(2, -2.85);
[below] at (0, -2.85) $0$;
[below] at (1, -2.85) $1$;
[below] at (2, -2.85) $2$;
(-0.75, -2)–(-0.85, -2);
(-0.75, -1)–(-0.85, -1);
(-0.75, 0)–(-0.85, 0);
(-0.75, 1)–(-0.85, 1);
[left] at (-0.85, -2) $-2$;
[left] at (-0.85, -1) $-1$;
[left] at (-0.85, 0) $0$;
[left] at (-0.85, 1) $1$;

Bifurcation diagram showing the fixed point and 2-cycle attractors (solid curves) and repellers (dashed curves) for the quadratic map $f(x;\,r)=r-x^2$, with a saddle-node bifurcation at $(-1/4,\,-1/2)$ and period-doubling bifurcations at $(3/4,\, 1/2)$, $(5/4,\, (1+\sqrt{2})/2)$, and $(5/4,\, (1-\sqrt{2})/2)$

Using this information, we can now extend our bifurcation diagram from Figure <ref> to Figure <ref>, which includes 2-cycles by also graphing the two functions in Equations <ref> and <ref> and indicating their stability. We can see now why the bifurcation point at $r=3/4$ with $f'(x_{s,\,1}(3/4);\,3/4) = -1$ is known as a period-doubling bifurcation: at $r=3/4$, the attractor changes from a fixed point attractor at $x=1/2$ to a periodic orbit with a point right above and right below $x=1/2$. Namely, the attractor's period doubles from 1 to 2. Similarly, at $r=5/4$, we have $\left(f^2\right)'\left(x_{p,\,{2}}(r);\,r\right)=-1$, which implies that the bifurcation of the 2-cycle from stable to unstable at $r=5/4$ is also a period-doubling bifurcation (from 2 to 4). We could prove this by continuing the algebra for $f^4(x;\,r)$, but instead, we will demonstrate a method to find these bifurcations numerically using a different, well-studied map.

§.§.§ The Logistic Map

The logistic map is another simple quadratic map defined by the following iteration function may:

\begin{equation}
f(x;\,r) = rx(1-x)
\end{equation}

It is probably the most famous example of complex dynamics emerging from a simple model. Despite its mathematical simplicity, rather than analytically analyzing the dynamics of the map as we did for the previous quadratic map, we will do it numerically.

To numerically plot a bifurcation diagram, we will use the method described by [p. 363]strogatz that plots the attractors of a system for many different values of $r$ and displays any saddle-node or period-doubling bifurcations.[These diagrams are more accurately called orbit diagrams [p. 368]strogatz because a true bifurcation diagram plots both the attractors and repellers of a system. However, since repellers are not important for this discussion, we will also refer to these diagrams that plot only attractors as bifurcation diagrams.] The method goes as follows: for a given value of $r$, start from an initial state $x_0$ and iterate some number of transients, say $100$ times, which allows the system to settle down to its attractor. Then, generate an orbit starting from $x_{100}$ for some number of iterations, say $100$ times again. This leaves us with the orbit $\{x_{100},\, x_{101},\,\hdots\,,\,x_{200}\}$. For our purposes, this orbit contains enough points in the attractor $A$ for a given value of $r$. Calculating these attractor orbits for many values of $r$, we get our bifurcation diagram by plotting these orbits on a graph of $x$ vs. $r$. The Python code in Appendix <ref> accomplishes this for $2\leq r\leq 4$, and the resulting graph is shown in Figure <ref>.

In this bifurcation diagram, we can see that the system starts with a fixed point attractor, but at around $r=3$, there is a period-doubling bifurcation, and the attractor turns into a 2-cycle. Near $r=3.45$, there is another period-doubling bifurcation to a 4-cycle, then another to an 8-cycle, and so on, the successive period-doubling bifurcations becoming faster and faster.[The factor by which the distance between successive bifurcations shrinks is a universal constant for systems approaching chaos by period-doubling known as the Feigenbaum constant [p. 523]layek. An interested reader is recommended to see the original paper by feigenbaum for details.]

Orbits of the logistic map for $0\leq t\leq 50$, graphed with the code in Appendix <ref>

Bifurcation diagram of the logistic map for $2\leq r\leq 4$, graphed with the code in Appendix <ref>

Lyapunov exponents of the logistic map for $2\leq r\leq 4$, graphed with the code in Appendix <ref>

To visualize these different periodic orbit attractors for different values of $r$, we can plot some points of the forward orbit $O^+(x_0)$ against time $t$ at different values of $r$. We use the code in Appendix <ref> to accomplish this for $x_0=0.1$ and the orbit $\{x_0,\,x_1,\,\hdots\,,\,x_{50}\}$, which is shown in Figure <ref>. In the figure, we draw lines between consecutive points $(k,\, x_k)$ and $(k+1,\, x_{k+1})$ in the orbit to make seeing the dynamics easier, but it is important to note that only the corners of these curves matter. In Figure <ref>, we plot at $r=2.75$, which is before the first period-doubling bifurcation. Therefore, the orbit settles down at a fixed point around $x=0.64$, corresponding to the attractor value shown in the bifurcation diagram (Figure <ref>) at $r=2.75$. After the first period-doubling bifurcation, we plot at $r=3.25$ (Figure <ref>), then after the second, we plot at $r=3.5$ (Figure <ref>). Comparing the 2-cycle and 4-cycle shown in Figures <ref> and <ref> to our bifurcation diagram in Figure <ref>, the periodic orbits correspond to their respective attractor sets.

As we can see in the bifurcation diagram in Figure <ref>, at around $r=3.57$, the successive bifurcations reach their limit, and the map approaches an attractor with essentially an infinite number of points. In other words, a forward orbit $O^+(x_0)$ for most values of $r>3.57$ will never settle down to a fixed point or periodic orbit; it will exhibit the aperiodic long-term behavior characteristic of chaotic systems. We can definitively show for which $r$ values the logistic map is chaotic by examining its sensitivity to initial conditions. Namely, we can calculate the Lyapunov exponent of the map for different values of $r$, which we accomplish using the code in Appendix <ref> and show a graph of in Figure <ref>. In this graph, we can see that the Lyapunov exponents stay less than or equal to 0 for all $r<3.57$. We should expect this because orbits approach a fixed point or periodic orbit attractor for these parameter values, meaning the system is not chaotic. However, when $r>3.57$, the Lyapunov exponents remain mainly positive, indicating that the system is chaotic. This matches with the chaotic-looking bifurcation diagram above it. Interestingly, there are some values of $r$ in the field of chaos above $r=3.57$ where the Lyapunov exponent is less than 0, most prominently where $r\approx3.83$. Matching this with the bifurcation diagram, we can see that when $r\approx3.83$, the logistic map briefly has a stable 3-cycle before dissolving into chaos again. These intervals of stability between regions of chaos displayed in Figure <ref> are known as periodic windows [p. 363]strogatz or islands of stability brandon.

Another interesting observation when comparing Figures <ref> and <ref> is that the period-doubling bifurcations always seem to occur when $\lambda = 0$. We can show this is the case by recalling from earlier in this section that a period-doubling bifurcation occurs when

\begin{equation}
\left(f^q\right)'(x_{p};\,r) = -1
\label{eq:period-doubling-bifurcation-requirement}
\end{equation}

where $x_{p}$ is a periodic point in an attracting $q$-cycle. Then, from Equation <ref>[We manipulate Equation <ref> by changing the chain rule step in Appendix <ref> (Equation <ref>) so that the derivatives are split by $q$ iterations rather than 1. In order for this to work, $t$ must be chosen so that $(t-1)/q\in\mathbb{N}$.] and choosing $x_0=x_p$,[If we didn't make this choice, the result would still be the same because $x_{qt}\to x_{p}$ for large $t$, and since $t\to\infty$, the factor $1/t$ will take care of all the non-zero terms. However, making this choice makes the math cleaner.]

\begin{equation}
\begin{split}
\lambda &= \lim_{t\to\infty}\frac{1}{t}\sum_{i=0}^{(t-1)/q}\ln\left|(f^q)'(x_{qi};\,r)\right| \\
&= \lim_{t\to\infty}\frac{1}{t}\sum_{i=0}^{(t-1)/q}\ln\left|(f^q)'(x_{p};\,r)\right| \\
&= \lim_{t\to\infty}\frac{1}{t}\sum_{i=0}^{(t-1)/q}\ln1 \\
&= 0
\end{split}
\end{equation}

We have now demonstrated a few connections between bifurcations and Lyapunov exponents using the logistic map example, which will be useful in our later analysis.

[<->] (0,2)–(0,-2);
[<->] (2,0)–(-2,0);
[above] at (0,2) $y$;
[right] at (2,0) $x$;
[thick, domain=-1.5:1.5, samples=100]plot(,0.25);
[thick, domain=-1.5:1.5, samples=100]plot(,-0.25);
[thick, domain=-0.25:0.25, variable=]  plot (-1.5, );
[thick, domain=-0.25:0.25, variable=]  plot (1.5, );

A rectangular region of states

[<->] (0,2)–(0,-2);
[<->] (2,0)–(-2,0);
[above] at (0,2) $y'$;
[right] at (2,0) $x'$;
[thick, domain=-1.5:1.5, samples=100]plot(,0.25+1-*);
[thick, domain=-1.5:1.5, samples=100]plot(,-0.25+1-*);
[thick, domain=-0.25:0.25, variable=]  plot (-1.5, +1-2.25);
[thick, domain=-0.25:0.25, variable=]  plot (1.5, +1-2.25);

The transformation of (a) by $T'$

[<->] (0,2)–(0,-2);
[<->] (2,0)–(-2,0);
[above] at (0,2) $y''$;
[right] at (2,0) $x''$;
[thick, domain=-1.5:1.5, samples=100]plot(0.5*,0.25+1-*);
[thick, domain=-1.5:1.5, samples=100]plot(0.5*,-0.25+1-*);
[thick, domain=-0.25:0.25, variable=]  plot (0.5*-1.5, +1-2.25);
[thick, domain=-0.25:0.25, variable=]  plot (0.5*1.5, +1-2.25);

The transformation of (b) by $T''$

[<->] (0,2)–(0,-2);
[<->] (2,0)–(-2,0);
[above] at (0,2) $y'''$;
[right] at (2,0) $x'''$;
[thick, domain=-1.5:1.5, samples=100]plot(0.25+1-*,0.5*);
[thick, domain=-1.5:1.5, samples=100]plot(-0.25+1-*,0.5*);
[thick, domain=-0.25:0.25, variable=]  plot (+1-2.25, 0.5*-1.5);
[thick, domain=-0.25:0.25, variable=]  plot (+1-2.25, 0.5*1.5);

The transformation of (c) by $T'''$

A chain of transformations that results in the Hénon map

§.§ Strange Attractors and Fractal Geometry

So far, our discussion of attractors has been limited to fixed points and periodic orbits, which are not chaotic. In this section, however, we will discuss attractors that do exhibit chaotic dynamics. According to ruelle, who coined the term, an attractor set $A$ is called strange if there are points in the attractor $\mathbf{x},\,\mathbf{a}\in A$ such that the distance between $\mathbf{f}^t(\mathbf{x})$ and $\mathbf{f}^t(\mathbf{a})$ grows exponentially with $t$ until it is bounded by the size of the attractor. In other words, strange attractors are chaotic, exhibiting sensitive dependence on initial conditions.[More generally, the term strange attractor can be used to describe any attractor with a fractal structure and not necessarily chaotic ones grebogi, but for the purposes of this paper, we will say that a strange attractor is one that is both chaotic and fractal.]

At first, the idea of an attractor being chaotic can seem counterintuitive. How can nearby trajectories get exponentially farther apart from each other while staying bounded in the attractor? Mathematically, this can be answered by looking at the Lyapunov exponents of the dynamics on a strange attractor. We show in Appendix <ref> that all the Lyapunov exponents of a fixed point attractor are negative because all nearby trajectories are attracted towards it. However, strange attractors have Lyapunov exponents that are both positive and negative; for example, a three-dimensional strange attractor has a Lyapunov spectrum of signs $\{+,\,0,\,-\}$ wolf. This is because a strange attractor is globally stable, but locally unstable, exhibiting both chaotic and attractive behavior. Infinitesimally close initial states grow further apart exponentially according to the $\lambda_1$ until the edge of the attractor prevents them from getting any further apart.

We can get an even better understanding of the concept of a strange attractor by looking at it geometrically. Strange attractors are generated through a “stretching and folding process” grassberger. As an intuitive example of how this works, we can imagine the stretching and folding of pastry dough. If we put a drop of food coloring in our pastry dough (representing a set of close initial states), then stretch and fold it a bunch of times, the food coloring will eventually spread throughout the entire thing. A map with a strange attractor will mathematically stretch and fold the attractor as it iterates, so a set of close initial conditions will eventually spread throughout the entire attractor. This gives us a good understanding of the seemingly counterintuitive nature of strange attractors: the stretching and folding property magnifies any initial perturbation in initial conditions, but by the nature of the attractor itself, trajectories stay bounded.

The Hénon attractor, graphed with the code in Appendix <ref>

§.§.§ The Hénon Map

To illustrate an example of a strange attractor, we will use a simple, well-studied map discovered by henon, who simulated the mathematical stretching and folding property of strange attractors by through a series of simple transformations in two-dimensional space. First, let us define the state vector $\mathbf{x}$ of the Hénon map as

\begin{equation}
\mathbf{x} = \begin{pmatrix}
x\e{1} \\
\end{pmatrix}
\begin{pmatrix}
x \\
\end{pmatrix}
\end{equation}

and two Hénon map parameters $a$ and $b$. Starting with a rectangular region of states (Figure <ref>), the first transformation $T'$ transforms $x$ into $x'$ and $y$ into $y'$ by stretching and folding the rectangle into a parabolic shape (Figure <ref>):

\begin{equation}
x' \\
\end{pmatrix}
\begin{pmatrix}
x \\
1+y-a x^2
\end{pmatrix}
\end{equation}

The next transformation $T''$ completes the folding process by contracting along the $x$-axis (Figure <ref>):

\begin{equation}
x'' \\
\end{pmatrix}
\begin{pmatrix}
b x' \\
\end{pmatrix}
\end{equation}

Finally, we complete the process with the transformation $T'''$ that reflects across the line $y=x$, which brings us back to our initial orientation with a stretched and folded set of states (Figure <ref>):

\begin{equation}
x''' \\
\end{pmatrix}
\begin{pmatrix}
y'' \\
\end{pmatrix}
\end{equation}

The composite transformation $T=T'''\,T''\,T'$ is the Hénon mapping:

\begin{equation}
x''' \\
\end{pmatrix}
\begin{pmatrix}
1 + y - a x^2 \\
b x
\end{pmatrix}
\end{equation}

Relabelling to our standard notation $\mathbf{x}_{k+1} = \mathbf{f}(\mathbf{x}_k)$,

\begin{equation}
\begin{pmatrix}
x_{k+1} \\
\end{pmatrix}
\begin{pmatrix}
1 + y_k - a x_k^2 \\
b x_k
\end{pmatrix}
\label{eq:henon-mapping}
\end{equation}

The Hénon attractor can be visualized computationally by using Equation <ref> to generate an orbit $O(\mathbf{x}_0)$. We accomplish this using the code in Appendix <ref>, which generates the orbit $O(\mathbf{x}_0) = \{ \mathbf{x}_{0},\, \mathbf{x}_{1},\, \hdots \,,\, \mathbf{x}_{250010} \}$ using $a = 1.4$, $b = 0.3$,[The choice of values of parameters $a$ and $b$ is the standard choice to ensure the proper amount of stretching and folding in order to generate the strange attractor. See the paper by henon for more insight into how these numbers were chosen.] and $\mathbf{x}_0 = \langle 0,\,0 \rangle$, then plots the orbit $\{ \mathbf{x}_{10},\, \mathbf{x}_{11},\, \hdots \,,\, \mathbf{x}_{250010}\}$. The result is shown in Figure <ref>.

To examine the strange chaotic nature of the Hénon attractor, specifically its global stability but local instability, let us examine its Lyapunov exponents. Recall from Section <ref> that the Lyapunov spectrum can be calculated using the eigenvalues of the matrix $J^{t\intercal}J^t$. From Equation <ref>, we know that the Hénon map is defined by the function

\begin{equation}
\mathbf{f}(\mathbf{x};\,a,\,b) =
\begin{pmatrix}
f^{[1]}(x,\,y;\,a,\,b) \\
\end{pmatrix}
\begin{pmatrix}
1+y-a x^2 \\
b x
\end{pmatrix}
\end{equation}

Therefore, by Equation <ref>, the Jacobian matrix $J(\mathbf{x})$ of the Hénon map is

\begin{equation}
J(\mathbf{x}) =
\begin{pmatrix}
\frac{\partial f^{[1]}}{\partial x} & \frac{\partial f^{[1]}}{\partial y} \\[4px]
\frac{\partial f^{[2]}}{\partial x} & \frac{\partial f^{[2]}}{\partial y}
\end{pmatrix}
\begin{pmatrix}
-2a x & 1 \\
b & 0
\end{pmatrix}
\end{equation}

which depends on $x$. Because we are finding the Lyapunov spectrum on the strange attractor, this $x$ value will keep changing and jumping around the attractor, so finding the matrix $J^{t\intercal}J^t$ and its eigenvalues analytically is impossible. Even numerically, approximating the spectrum using a large value of $t$ will require computationally expensive matrix multiplication, and the matrix product will likely overflow before reaching a sufficiently high value of $t$ for an accurate approximation. For this reason, we will use the QR factorization method for Lyapunov spectrum calculation that we derive in Appendix <ref>. Specifically, in the code in Appendix <ref>, we implement the algorithm described in Appendix <ref> and approximate the spectrum using a large value of $t$ in Equation <ref>:

\begin{equation}
\lambda_i = \lim_{t\to\infty}\frac{1}{t}\sum_{j=1}^t \ln \left|r^{(j)}_{ii}\right|
\end{equation}

This yields that

\begin{equation}
\lambda = \{\lambda_1,\,\lambda_2\} \approx \{0.419,\, -1.623\}
\label{eq:henon_lyap_spec}
\end{equation}

As we can see, the maximal Lyapunov exponent is positive, indicating the chaotic behavior on the strange attractor and its local instability. However, $\lambda_2$ is negative, indicating its global stability.[The values we calculate in Appendix <ref> and display in Equation <ref> agree with the values calculated in the thesis by asbroek, which we recommend an interested reader to see for a more in-depth analysis.]

$0.5\leq x\leq 0.75$, $0.16\leq y\leq 0.21$

$0.61\leq x\leq 0.635$, $0.187\leq y\leq 0.192$

Zooming in twice on the Hénon attractor, graphed by changing the $x$-limits and $y$-limits of the code in Appendix <ref>

§.§.§ Fractal Dimension and Measure

Now that we have quantified the chaotic dynamics on the Hénon attractor, we are ready to explore its geometry. To begin, let us zoom in on the Hénon attractor in Figure <ref> to examine its fine structure, or detail at a small scale. By altering the $x$ and $y$-limits of the code in Appendix <ref> to $0.5\leq x\leq 0.75$ and $0.16\leq y\leq 0.21$, we get the graph in Figure <ref>. As we can see, there appear to be six parallel curves: three bunched together at the top, two together below that, and one lone curve below that. However, if we zoom in again on the three curves at the top, which we show in Figure <ref>, we can see that it is actually composed of six curves grouped exactly as it was before: three, two, and one. If we were to zoom in more, we would see the same thing over and over no matter how far we zoom in henon. This property of the Hénon attractor's fine structure and detail not vanishing, but rather, remaining essentially unchanged as we zoom in to arbitrarily small scales is characteristic of geometric objects known as fractals mandelbrot1. In a way, fractals are a rebellion against calculus, which generally assumes that things eventually appear smooth if we look closely enough.

Fractal geometry is important for us to study because it is a defining feature of strange attractors and, more generally, multiple types of sensitivity to initial conditions. To understand why this is the case, recall that a strange attractor is generated through a process of stretching and folding. Once again using the Hénon map as an example, let us look at Figure <ref> again. We see that starting with a rectangle of initial conditions, iterating the map once gives us two layers. Doing it again would give us four, then eight, and so on. Taking this to the limit, we would have infinitely many infinitely thin layers, all separated by gaps of different sizes. This is essentially what the Hénon attractor is, and it is why strange attractors are fractals.

The geometry of fractal attractors can be complex and difficult to describe because of its departure from calculus and our standard mathematical tools. One way that we can quantitatively characterize fractal geometry is through the idea of a fractal dimension, which generalizes our standard idea of a dimension in a clever way.

In order to motivate the concept of a fractal dimension, we must first introduce the concept of measure. Measure is a way of generalizing the concepts of length in one dimension, area in two dimensions, and volume in three dimensions to any geometric object [p. 2]tao. Let us denote the measure of some set $S$ as $\mu(S)$. If $S$ is $n$-dimensional, meaning we can describe any point in it using $n$ numbers, we can denote its measure as $\mu^n(S)$, but we often omit this superscript when the dimension is clear from context. A valid measure $\mu(S)$ has the following properties [p. 250]alligood:

* For any set $S$,

\begin{equation}
\mu(S)\geq 0
\end{equation}

This makes sense intuitively because measure can't be negative.

* For a union set $S_1\cup S_2\cup\hdots = \bigcup_{i=1}^{\infty} S_i$ where no two sets $S_i$ intersect,

\begin{equation}
\mu\left(\bigcup_{i=1}^{\infty} S_i\right) = \sum_{i=1}^{\infty}\mu(S_i)
\end{equation}

In other words, the measure of a set is equal to the sum of the measures of all the distinct subsets that the set is made of. This property also works for a finite number of sets. Physically, this means that if we divide a three-dimensional object into a bunch of pieces, adding up the volumes of all the individual pieces will give the total volume of the object.

* If $\emptyset$ is the empty set, or the set with nothing in it, it is clear that

\begin{equation}
\mu(\emptyset) = 0
\end{equation}

* If $S$ is $n$-dimensional and $m>n$, then

\begin{equation}
\mu^{m}(S) = 0
\label{eq:measure-property-4}
\end{equation}

In other words, the length of a point, or the area of a line, or the volume of a sheet are all 0.

* For $S\subset\mathbb{R}^n$ and $\mathbf{x}\in\mathbb{R}^n$,

\begin{equation}
\mu(S+\mathbf{x}) = \mu(S)
\end{equation}

In other words, if we move something around in space, its measure doesn't change.

Now that we have defined measure, we can now begin to talk about dimension. First, let us consider the one-dimensional set $S=(0,\,a)$, where $(0,\,a)$ is the open interval from 0 to $a$. This line segment clearly has a length of $a$, so $\mu^1(S) = a$. Now, if we scale up this segment by a factor of 2, we get the set $2S = (0,\,2a)$. Splitting this up, we can see that $2S$ is the union of two segments with measure $a$, namely, the interval from 0 to $a$ and the interval from $a$ to $2a$ (by Property 5 of measures). Then, by Property 2 of measures, $\mu^1(2S) = 2\mu^1(S) = 2a$, which intuitively makes sense because the line segment is two times longer. In other words, scaling a one-dimensional object by 2 scales its measure by 2.

Let us now consider a two-dimensional set $S = (0,\,a)\times(0,\,a)$,[$S_1\times S_2$ is the Cartesian product of $S_1$ and $S_2$, which is the set of all vectors $\langle a,\,b \rangle$ such that $a\in S_1$ and $b\in S_2$.] which is a square with side length $a$. This square has an area of $a^2$, so $\mu^2(S) = a^2$. Scaling up this square by a factor of 2, we get $2S = (0,\,2a)\times(0,\,2a)$, which is a square of side length $2a$. This square has area $4a^2$, so $\mu^2(2S) = 4\mu^2(S)$, meaning scaling a two-dimensional object by 2 scales its measure by 4. For a three-dimensional set $S = (0,\,a)\times(0,\,a)\times(0,\,a)$, which is a cube with side length $a$ and volume $a^3$, scaling $S$ up by a factor of 2 will result in a cube with side length $2a$ and volume $8a^3$, so $\mu^3(2S) = 8\mu^3(S)$. In other words, scaling a three-dimensional object by 2 scales its measure by 8.

In summary, we have

\begin{equation}
\begin{split}
\mu^1(2S) &= 2\mu^1(S) \\
\mu^2(2S) &= 4\mu^2(S) \\
\mu^3(2S) &= 8\mu^3(S)
\end{split}
\end{equation}

The pattern here is clear: the dimension $n$ of a given object is the power to which a scaling factor is raised to scale its measure. In other words, if we scale an $n$-dimensional object $S$ by some scaling factor $\sigma$,

\begin{equation}
\mu(\sigma S) = \sigma^n\mu(S)
\label{eq:scaling-factor-measure-nonfractal}
\end{equation}

Then, if we have some set $S$ that defines a fractal object and it follows the equation

\begin{equation}
\mu(\sigma S) = \sigma^d\mu(S)
\label{eq:scaling-factor-measure}
\end{equation}

for any $\sigma$, we can interpret the exponent $d$ as the fractal dimension of $S$.

Covering the Hénon attractor in Figure <ref> with two-dimensional boxes of side length $\epsilon=1/4$

A clever way of practically determining this fractal dimension is through the method of box-counting. The process goes as follows: for a set $S\subset\mathbb{R}^n$, fill state space with $n$-dimensional boxes of side length $\epsilon$.[A one-dimensional box is an interval of length $\epsilon$, a two-dimensional box is a square of side length $\epsilon$, and so on.] Say that the number of boxes that contain at least one element of $S$ is $N(\epsilon)$, which is approximately proportional to the set's measure $\mu(S)$. In Figure <ref>, we cover the Hénon attractor with boxes of size $\epsilon=1/4$. Out of the 48 boxes shown, the attractor touches 28 of them, so $N(1/4)=28$.

Notice that scaling $S$ by some factor $\sigma$ is equivalent to scaling the box size $\epsilon$ by $1/\sigma$. For example, consider a rectangle in two-dimensional space. If we scale it up by 2, the rectangle will cover 4 times as much area, so $N(\epsilon)$ will scale by 4. However, if we cut $\epsilon$ in half, all the boxes will get 4 times smaller, so the same thing happens. Therefore, we can rewrite Equation <ref> as

\begin{equation}
N(\epsilon) \approx \left(\frac{1}{\epsilon^d}\right)c\mu(S)
\label{eq:box-counting-power-function}
\end{equation}

where $c$ is the constant of proportionality. Taking logs on both sides,[We use the natural logarithm, but of course, this will work for any convenient base.] we get

\begin{equation}
\ln N(\epsilon) \approx d\ln\left(\frac{1}{\epsilon}\right) + \ln(c\mu(S))
\label{eq:box-counting-line-regression}
\end{equation}

Solving for $d$,

\begin{equation}
d \approx \frac{\ln N(\epsilon) - \ln(c\mu(S))}{\ln(1/\epsilon)}
\label{eq:box-counting-approximation}
\end{equation}

Now, if we let $\epsilon$ approach 0, this box-counting method will capture all of the fine structure and detail of $S$, so in the limit, $N(\epsilon)$ is exactly proportional to the measure of the scaled set. Because $\ln N(\epsilon)$ and $\ln(1/\epsilon)$ grow without bound, the $\ln(c\mu(S))$ term becomes insignificant in the limit, so we can rewrite Equation <ref> as

\begin{equation}
d = \lim_{\epsilon\to 0}\frac{\ln N(\epsilon)}{\ln(1/\epsilon)}
\label{eq:box-counting-dimension-exact}
\end{equation}

This fractal dimension $d$ is known as the box-counting or Minkowski–Bouligand dimension [p. 70]ott.[There are many different other kinds of fractal dimensions that capture the geometry of fractals in different ways, like the information dimension, correlation dimension, Rényi dimension, and Hausdorff dimension, to name a few, but the box-counting dimension is a simple one that works well intuitively and computationally.]

To get a better understanding of box-counting, we will once again examine the Hénon attractor. Obviously, it is impractical to divide all of state space into an infinite number of infinitesimal boxes to determine the dimension of a geometric object, so instead, our strategy is to utilize Equation <ref>. Namely, if we numerically calculate a bunch of $N(\epsilon)$ values for different box sizes and plot the values on a graph of $N(\epsilon)$ vs. $1/\epsilon$, we should expect these points to approximately fall along a power function described by Equation <ref>. Then, by Equation <ref>, if we plot these points on a graph of $\ln N(\epsilon)$ vs. $\ln(1/\epsilon)$, we should see the points approximately fall along a line with a slope of $d$. Therefore, if we were to calculate the linear regression of these points $(\ln(1/\epsilon),\,\ln N(\epsilon))$, its slope would be approximately equal to the set's fractal dimension $d$.

$\epsilon$     $N(\epsilon)$

1/4     28

1/8     63

1/16     148

1/32     349

1/64     832

1/128     1921

1/256     4519

Some $N(\epsilon)$ values of the Hénon attractor, calculated using the code in Appendix <ref>

[->] (1.4,3.56)–(8.6,3.56);
[->] (1.4,3.56)–(1.4,13.39);
[thick, domain=1.5:8.5]plot(,1.227*+2.321);
[thick] (2,4.8074) circle [radius=0.075];
[thick] (3,5.9773) circle [radius=0.075];
[thick] (4,7.2095) circle [radius=0.075];
[thick] (5,8.4471) circle [radius=0.075];
[thick] (6,9.7004) circle [radius=0.075];
[thick] (7,10.908) circle [radius=0.075];
[thick] (8,12.142) circle [radius=0.075];
[below] at (5,3) $\log_2(1/\epsilon)$;
[left, rotate=90] at (0.55,9.5) $\log_2 N(\epsilon)$;
(2, 3.56)–(2, 3.46);
(3, 3.56)–(3, 3.46);
(4, 3.56)–(4, 3.46);
(5, 3.56)–(5, 3.46);
(6, 3.56)–(6, 3.46);
(7, 3.56)–(7, 3.46);
(8, 3.56)–(8, 3.46);
[below] at (2, 3.46) $2$;
[below] at (3, 3.46) $3$;
[below] at (4, 3.46) $4$;
[below] at (5, 3.46) $5$;
[below] at (6, 3.46) $6$;
[below] at (7, 3.46) $7$;
[below] at (8, 3.46) $8$;
(1.4, 4)–(1.3, 4);
(1.4, 5)–(1.3, 5);
(1.4, 6)–(1.3, 6);
(1.4, 7)–(1.3, 7);
(1.4, 8)–(1.3, 8);
(1.4, 9)–(1.3, 9);
(1.4, 10)–(1.3, 10);
(1.4, 11)–(1.3, 11);
(1.4, 12)–(1.3, 12);
(1.4, 13)–(1.3, 13);
[left] at (1.3, 4) $4$;
[left] at (1.3, 5) $5$;
[left] at (1.3, 6) $6$;
[left] at (1.3, 7) $7$;
[left] at (1.3, 8) $8$;
[left] at (1.3, 9) $9$;
[left] at (1.3, 10) $10$;
[left] at (1.3, 11) $11$;
[left] at (1.3, 12) $12$;
[left] at (1.3, 13) $13$;

A plot of the points in Table <ref> on a graph of $\log_2 N(\epsilon)$ vs. $\log_2(1/\epsilon)$ and their line of best-fit

We follow this method to calculate the fractal dimension of the Hénon map by first using the code in Appendix <ref> to generate values of $N(\epsilon)$ for various values of $\epsilon=2^{-k}$ for numerical simplicity. These values are shown in Table <ref>. Since we have $\epsilon$ values that are integer powers of 2, we take base 2 logs instead of natural logs and plot the points $(\log_2(1/\epsilon),\,\log_2 N(\epsilon))$ in Figure <ref>. Then, numerically determining the best-fit line of these points gives us

\begin{equation}
\log_2 N(\epsilon) = 1.23\log_2\left(\frac{1}{\epsilon}\right) + 2.32
\end{equation}

with an $R^2$ value of $0.9999$. We can therefore conclude that the Hénon attractor is approximately 1.23-dimensional.

Intuitively, this makes sense. Although the Hénon attractor lives in two-dimensional space, it is composed of curves that appear to be one-dimensional (see Figures <ref> and <ref>). Therefore, it should scale similarly to a one-dimensional object. However, it has a fine structure and self-similarity that touches more small-$\epsilon$ boxes than a simple curve would, so its dimension $d$ is slightly higher than 1. With this intuition, we are now prepared to give a formal definition of a fractal: a fractal is a geometric object with a non-integer dimension mandelbrot2.

This idea of a non-integer dimension captures our intuitive understanding of the detail and roughness of fractals, objects that don't have length, area, or volume, but rather, a positive measure of something in between. The “rougher” a geometric object is, the more boxes it needs to capture its fine detail at a small scale, so the higher its fractal dimension is. In fact, physical objects don't have integer dimensions because they aren't perfectly smooth geometric objects; they all have some roughness, and thus, they have non-integer dimensions and a fractal-like structure.[For more information about how fractals permeate the natural world, see the book by mandelbrot2 for more details.]

§.§.§ The Kaplan-Yorke Conjecture

In a paper by kaplan-yorke, a quantity now known as the Lyapunov dimension $d_l$ is introduced. To calculate this quantity, which we will soon discover has a strong connection to the fractal structure of attractors, recall from Section <ref> that a Lyapunov spectrum's exponents are ordered such that $\lambda_1\geq\lambda_2\geq\hdots\geq\lambda_n$. Now, let $\kappa$ be the largest index such that

\begin{equation}
\sum_{i=1}^{\kappa}\lambda_i\geq 0
\label{eq:kappa}
\end{equation}

Using this, we define the Lyapunov dimension $d_l$ as

\begin{equation}
d_l = \kappa + \frac{1}{|\lambda_{\kappa+1}|}\sum_{i=1}^{\kappa}\lambda_i
\label{eq:dl}
\end{equation}

Referring back to the Lyapunov spectrum of the Hénon map that we calculated earlier in this section, $\{\lambda_1,\,\lambda_2\}\approx\{0.419,\,-1.623\}$, we can see that the index $\kappa$ is equal to 1 since $\lambda_1+\lambda_2<0$, so the Lyapunov dimension $d_l$ of the Hénon attractor is

\begin{equation}
d_l = 1+\frac{\lambda_1}{|\lambda_2|} \approx 1+\frac{0.419}{|-1.623|} \approx 1.26
\end{equation}

Notice that this Lyapunov dimension $d_l$ is surprisingly close to the attractor's fractal dimension $d\approx 1.23$, which is off by less than 3%, well within the range of numerical error in calculating the fractal dimension and Lyapunov exponents. This is a result of the Kaplan-Yorke conjecture, which states that the Lyapunov dimension of an attractor is equal to its fractal dimension for some systems ott-attractor-dim.[More specifically, the Kaplan-Yorke conjecture states that the Lyapunov dimension is equal to an object's information dimension, which is calculated in a slightly different way to the box-counting dimension, but is, in almost all cases, close to or equal to it.] In other words,

\begin{equation}
d_l = d
\end{equation}

The Kaplan-Yorke conjecture is rather extraordinary because it relates the dynamics on an attractor to the attractor's geometry and structure, and it is a strong support of the fact that chaos and certain geometrical properties like fractal geometry are linked. Although some of the systems we examine in this paper do not actually obey the Kaplan-Yorke conjecture for various reasons, we will still use it as an indication of whether the attractor of a system is fractal.

§.§ Basins of Attraction

Recall from Section <ref> that Property 2 of attractors says that $A$ is an attractor if it attracts an open set of initial conditions $U$, where for all initial conditions $\mathbf{x}_0\in U$, this distance from $\mathbf{f}^t(\mathbf{x}_0)$ to $A$ goes to 0. The largest $U$ with this property is called the basin of attraction of $A$ [p. 332]strogatz, which we will denote as $\hat{A}$. Essentially, $\hat{A}$ is a region in state space composed of all the initial conditions that approach $A$ in the long term. The geometrical structure of different basins of attraction varies greatly among different dynamical systems, and depending on what type of basin an initial state belongs to, its long-term behavior can be fundamentally different ottboa. For this reason, our first priority concerning basins of attraction is establishing a method to classify different basins and quantify their size.

To do this, we will follow a similar approach to the one outlined by sprott. In this method, rather than dealing with the complex geometry of different attractors, we instead capture the general location and size of the attractor with two numbers. The first number is the attractor's mean or center of mass $\langle A\rangle$, which is the average of all the attractor's points. For a fractal attractor composed of an infinite number of points,[We concern our analysis with attractors like these because our equations can be trivially altered for attractors of a finite number of points.] let us say $A=\{\mathbf{a}_1,\,\mathbf{a}_2,\,\hdots\}$, where the states $\mathbf{a}_i$ are in no particular order. Then, the mean of $A$ is

\begin{equation}
\langle A\rangle = \lim_{j\to\infty}\frac{1}{j}\sum_{i=1}^j\mathbf{a}_i
\label{eq:mean}
\end{equation}

Our second number of interest is the attractor's standard deviation $\sigma_A$, or the variation of the attractors' points around its mean. From statistics [p. 56]stats, we know this is calculated by

\begin{equation}
\sigma_A = \sqrt{\lim_{j\to\infty}\frac{1}{j}\sum_{i=1}^j|\mathbf{a}_i-\langle A\rangle|^2}
\label{eq:standard-dev}
\end{equation}

Using this, we can define a quantity $\xi$ that we will call “normalized distance” from an attractor, which is the Euclidean distance between a state $\mathbf{x}$ and the attractor's mean, normalized by the attractor's standard deviation:

\begin{equation}
\xi = \frac{|\mathbf{x} - \langle A\rangle|}{\sigma_A}
\label{eq:normalized-distance}
\end{equation}

This notion of distance accounts for the fact that some attractors are bigger than others, and it will allow us to quantify basins relative to their attractor's size.

To motivate our classification method, we will introduce two more useful sets. First, let us say $S(\xi)$ is the set of all states that lie in an open $n$-dimensional ball[We define a one-dimensional ball to be a line segment, a two-dimensional ball to be a disk, a three-dimensional ball to be a filled-in sphere, and so on.] of radius $\xi$ centered at $\langle A\rangle$, where $n$ is the dimension of the state space $A$ lives in. In other words, $S(\xi)$ contains all the states less than $\xi$ away from $\langle A\rangle$. Mathematically, this definition is

\begin{equation}
S(\xi) = \{\mathbf{x}:\:|\mathbf{x}-\langle A\rangle|<\xi\}
\end{equation}

where $:$ means “such that.” Our second useful set is $\hat{A}(\xi)$, which is the subset of the basin of attraction $\hat{A}$ containing all the states in the $n$-dimensional ball of radius $\xi$ centered at $\langle A\rangle$ that get attracted to $A$. In other words, it is the intersection of $\hat{A}$ and $S(\xi)$:[The intersection of $S_1$ and $S_2$, denoted $S_1\cap S_2$, is the set of all the elements contained in both $S_1$ and $S_2$.]

\begin{equation}
\hat{A}(\xi) = \hat{A}\cap S(\xi)
\label{eq:a-hat-xi-def}
\end{equation}

Now, we are ready to introduce the function that our basin classification method is based around: $P(\xi)$, which is the probability that an initial state $\mathbf{x}_0\in S(\xi)$ is in the basin of attraction of $A$. In other words, it is the fraction of states in the $n$-dimensional ball with radius $\xi$ centered at $\langle A\rangle$ that get attracted to $A$. Using the measure function, we can write this as

\begin{equation}
P(\xi) = \frac{\mu(\hat{A}(\xi))}{\mu(S(\xi))}
\label{eq:p-xi-def}
\end{equation}

In the limit $\xi\to\infty$, $P(\xi)$ follows a power law sprott:

\begin{equation}
P(\xi) = \frac{P_0}{\xi^{\gamma}}
\label{eq:probability-function-power-law}
\end{equation}

Based on these parameters $P_0$ and $\gamma$, we can now divide basins of attraction into four distinct classes, ordered from largest to smallest:

* Class 1 basins have $P_0=1$ and $\gamma=0$. Then, the power law is $P(\xi)=1$, which means the probability of a state in any $n$-dimensional ball being in the basin of attraction is 100%. Therefore, these basins include all of state space (except perhaps a set of finite measure).

* Class 2 basins have $P_0<1$ and $\gamma=0$. Then, the power law is $P(\xi)=P_0$, which is independent of $\gamma$. Therefore, all $n$-dimensional balls centered at $\langle A\rangle$ contain the same fraction of states that go to $A$, meaning these basins must occupy a fixed fraction of state space.

* Class 3 basins have $0<\gamma<n$, where $n$ is the dimension of the state space the basin belongs to. Notice that, as we scale $\xi$, $P(\xi)$ drops off slower than $\mu(S(\xi))$ grows since as we recall from Section <ref> and our discussion of dimension, the measure of an $n$-dimensional ball scales according to $\xi^n$ while $P(\xi)$ scales according to $\xi^{-\gamma}$. Therefore, these basins extend to infinity in some directions but take up an increasingly small fraction of state space as we move out.

* Class 4 basins have $\gamma=n$. In this case, since $P(\xi)$ drops off at the same rate $\mu(S(\xi))$ grows, we can conclude that these basins occupy a finite region of state space and therefore have a well-defined size relative to their attractors. The linear, normalized size of a Class 4 basin is given by

\begin{equation}
\xi_0 = P_0^{1/n}
\end{equation}

It is worth noting that although using the normalized distance $\xi$ doesn't affect the basin class or the value of $\gamma$, it does affect $P_0$, and thus, the relative size of a Class 4 basin. See the paper by sprott for examples of all these kinds of basins in both two and three-dimensional state space.

Here, the system we will explore is once again the Hénon map, whose basin we will find and classify to demonstrate our numerical method of finding $P(\xi)$. Since points that do not end up in the Hénon attractor diverge to infinity henon, we can calculate $P(\xi)$ for the Hénon attractor easily by saying that if an initial state $\xi$ away from $\langle A\rangle$ eventually maps to a point greater than $10000\xi$ away from $\langle A\rangle$, it isn't in $\hat{A}$. Because $P(\xi)$ is a power function, we can find $\gamma$ using a similar method to the one we used to find the fractal dimension $d$ in Section <ref>. Namely, taking logs on both sides of Equation <ref>, we get that

\begin{equation}
\ln P(\xi) = -\gamma\ln\xi + \ln P_0
\label{eq:log-log-xi-relationship}
\end{equation}

where we purposefully take a $-\gamma$ out rather than leaving a $1/\xi$ inside the log because we are interested in the limit $\xi\to\infty$. Then, if we numerically find $P(\xi)$ for various $\xi = 2^k$, plotting $\log_2 P(\xi)$ vs. $\log_2 \xi$ will give us points along a line with slope $-\gamma$.

We can approximate $P(\xi)$ for a given $\xi$ by choosing a bunch of random states in $S(\xi)$ and iterating them to see if they end up in the basin. However, a problem arises in the case where there is a lot of space in $S(2^k)$ that doesn't contain the basin. Then, $P(2^k)$ may be very close to 0, so a numerical calculation testing only a finite number of random states may not be accurate. To counteract this issue, we detail a shell method for calculating values of $P(\xi)$ in Appendix <ref>. The concept behind this method is an iterative calculation of $P(2^{k+1})$ using $P(2^k)$ and $\Delta P(2^{k})$, which is the probability that a state lying in the shell centered at $\langle A\rangle$ with inner radius $\xi=2^k$ and outer radius $\xi=2^{k+1}$ is also in the basin of attraction $\hat{A}$. We show in Appendix <ref> that this can be accomplished using

\begin{equation}
P(2^{k+1}) = \frac{P(2^k)}{2^n} + \left(1-\frac{1}{2^n}\right)\Delta P(2^k)
\label{eq:iteration-of-shell-method}
\end{equation}

We can use this equation to generate a data set of $P(2^k)$ values using a Monte Carlo algorithm. First, we estimate $P(1)$ by first randomly picking initial states in the sphere of radius $\xi=1$ centered at $\langle A\rangle$ and checking whether or not they are in the basin by generating their orbits. Then, we get $P(1)$ by dividing the number of states that were in the basin of attraction by the total number of states we tested. Using a similar approach, we can estimate $\Delta P(1)$ by randomly picking states in the shell with inner radius $\xi=1$ and outer radius $\xi=2$ centered at $\langle A\rangle$. We can then calculate $P(2)$ using Equation <ref>:

\begin{equation}
P(2) = \frac{P(1)}{2^n} + \left(1-\frac{1}{2^n}\right)\Delta P(1)
\end{equation}

We can then numerically estimate $\Delta P(2)$ to calculate $P(4)$, and so on.

This seemingly roundabout way of finding $P(\xi)$ is useful because it allows us to distinguish $P(\xi)$ from 0 for large values of $\xi$. For example, even if $\Delta P(2^k)$ is 0 because the attractor lies entirely inside of the shell's inner radius $\xi=2^k$, calculating $P(2^{k+1})$ using Equation <ref> will give a non-zero result because we iterate up to it. However, if we were to calculate $P(2^{k+1})$ directly, our Monte Carlo algorithm may fail to give an accurate result due to the vanishingly small basin within $S(2^{k+1})$ that a computationally reasonable random selection of points wouldn't pick up on. Thus, our shell method allows us to find an accurate power function $P(\xi)$ in a computationally efficient way.

$\xi$     $P(\xi)$

$2^0=1$     0.8751

$2^1=2$     0.6976

$2^2=4$     0.2708

$2^3=8$     0.09465

$2^4=16$     0.03085

$2^5=32$     0.01000

$2^6=64$     0.003124

$\xi$     $P(\xi)$

$2^7=128$     9.220e-4

$2^8=256$     2.598e-4

$2^9=512$     7.469e-5

$2^{10}=1024$     2.167e-5

$2^{11}=2048$     6.918e-6

$2^{12}=4096$     1.730e-6

$2^{13}=8192$     4.324e-7

Some approximate $P(\xi)$ values of the Hénon basin, calculated using the code in Appendix <ref>

[->] (-0.5,3/2)–(13.5,3/2);
(-1,0)–(15, 0);
[->] (-0.5,3/2)–(-0.5,-24/2);
[thick, domain=-0.25:13.25]plot(,-1.7503/2*+1.9694/2);
[thick] (2,-1.885/2) circle [radius=0.15];
[thick] (3,-3.401/2) circle [radius=0.15];
[thick] (4,-5.019/2) circle [radius=0.15];
[thick] (5,-6.644/2) circle [radius=0.15];
[thick] (6,-8.322/2) circle [radius=0.15];
[thick] (7,-10.083/2) circle [radius=0.15];
[thick] (8,-11.91/2) circle [radius=0.15];
[thick] (9,-13.71/2) circle [radius=0.15];
[thick] (10,-15.49/2) circle [radius=0.15];
[thick] (11,-17.14/2) circle [radius=0.15];
[thick] (12,-19.14/2) circle [radius=0.15];
[thick] (13,-21.1411/2) circle [radius=0.15];
[thick] (0,-0.192/2) circle [radius=0.15];
[thick] (1,-0.52/2) circle [radius=0.15];
[above] at (6.25,2.5) $\log_2\xi$;
[left, rotate=90] at (-2.9,-3) $\log_2 P(\xi)$;
(0, 1.5)–(0, 1.7);
(1, 1.5)–(1, 1.7);
(2, 1.5)–(2, 1.7);
(3, 1.5)–(3, 1.7);
(4, 1.5)–(4, 1.7);
(5, 1.5)–(5, 1.7);
(6, 1.5)–(6, 1.7);
(7, 1.5)–(7, 1.7);
(8, 1.5)–(8, 1.7);
(9, 1.5)–(9, 1.7);
(10, 1.5)–(10, 1.7);
(11, 1.5)–(11, 1.7);
(12, 1.5)–(12, 1.7);
(13, 1.5)–(13, 1.7);
[above] at (0, 1.7) $0$;
[above] at (2, 1.7) $2$;
[above] at (4, 1.7) $4$;
[above] at (6, 1.7) $6$;
[above] at (8, 1.7) $8$;
[above] at (10, 1.7) $10$;
[above] at (12, 1.7) $12$;
(-0.5, 0)–(-0.7, 0);
(-0.5, 1/2)–(-0.7, 1/2);
(-0.5, 2/2)–(-0.7, 2/2);
(-0.5, -1/2)–(-0.7, -1/2);
(-0.5, -2/2)–(-0.7, -2/2);
(-0.5, -3/2)–(-0.7, -3/2);
(-0.5, -4/2)–(-0.7, -4/2);
(-0.5, -5/2)–(-0.7, -5/2);
(-0.5, -6/2)–(-0.7, -6/2);
(-0.5, -7/2)–(-0.7, -7/2);
(-0.5, -8/2)–(-0.7, -8/2);
(-0.5, -9/2)–(-0.7, -9/2);
(-0.5, -10/2)–(-0.7, -10/2);
(-0.5, -11/2)–(-0.7, -11/2);
(-0.5, -12/2)–(-0.7, -12/2);
(-0.5, -13/2)–(-0.7, -13/2);
(-0.5, -14/2)–(-0.7, -14/2);
(-0.5, -15/2)–(-0.7, -15/2);
(-0.5, -16/2)–(-0.7, -16/2);
(-0.5, -17/2)–(-0.7, -17/2);
(-0.5, -18/2)–(-0.7, -18/2);
(-0.5, -19/2)–(-0.7, -19/2);
(-0.5, -20/2)–(-0.7, -20/2);
(-0.5, -21/2)–(-0.7, -21/2);
(-0.5, -22/2)–(-0.7, -22/2);
(-0.5, -23/2)–(-0.7, -23/2);
[left] at (-0.7, 0) $0$;
[left] at (-0.7, -4/2) $-4$;
[left] at (-0.7, -8/2) $-8$;
[left] at (-0.7, -12/2) $-12$;
[left] at (-0.7, -16/2) $-16$;
[left] at (-0.7, -20/2) $-20$;

A plot of the points in Table <ref> on a graph of $\log_2 P(\xi)$ vs. $\log_2\xi$ and their line of best-fit

The code in Appendix <ref> implements this method to find $P(2^k)$ values for the basin of the Hénon attractor, which are displayed in Table <ref>. Using this code, we find the Hénon attractor's mean is approximately $\langle A\rangle\approx \langle 0.2571,\,0.0771 \rangle$ and its standard deviation is approximately $\sigma_A\approx 0.7526$, both of which are visually consistent with Figure <ref>. In Figure <ref>, we plot the values in Table <ref> on a $\log_2 P(\xi)$ vs. $\log_2\xi$ plot and take a linear regression,[We neglect the $P(1)$ and $P(2)$ values when taking the linear regression because we are interested in the limit $\xi\to\infty$.] which gives us,

\begin{equation}
\log_2P(\xi) = -1.750\log_2\xi + 1.969
\end{equation}

with an $R^2$ value of 0.999. By Equation <ref>, $\gamma\approx-1.750$ and $P_0\approx3.196$, which indicates that the Hénon attractor has a Class 3 basin.

The basin of the Hénon attractor for $-5\leq x\leq 5$ and $-5\leq y\leq 5$ using the parameters $a=1.4$ and $b=0.3$, graphed with the code in Appendix <ref>

In Figure <ref>, we visualize the Hénon basin using the code in Appendix <ref>. Here, points in the basin of the Hénon attractor are shown in white, and points that diverge to infinity are shown in black. In this regard, we can consider infinity the other attractor of the Hénon map, meaning the white region is the basin of the strange attractor in Figure <ref> and the black region is the basin of infinity. Visually, we can see that Figure <ref> aligns with the fact that the Hénon basin is Class 3: it has a quasi-parabolic shape that takes up less and less space as we extend outward.

The fact that the Hénon map has two attractors, the strange attractor and infinity, makes the Hénon map a bistable system, or a system with two distinct, coexisting attractors. Bistability is the simplest case of multistability, which is defined as the coexistence of several attractors for a given set of system parameters [pg. 6]pisarchikbook. Each attractor of a multistable system has an associated basin of attraction, and the set of points that separate two or more basins is called the boundary of the basins. As we can see, the basin boundary of the Hénon map with our standard parameters $a=1.4$ and $b=0.3$ is smooth, but by varying the parameters, we can get more interesting geometries. Many dynamical systems exhibit complex and counterintuitive geometrical properties resulting in sensitive dependence on initial conditions in their basins of attractions and basin boundaries. We will explore some of these geometrical properties in this section.

(0, 0)–(10, 0);
(10, 0)–(10, 4);
(10, 4)–(0, 4);
(0, 4)–(0, 0);
[very thick] (6, 4)–(4, 0);
[below right] at (0, 4) $\hat{A}$;
[below left] at (10, 4) $\hat{C}$;
(4.75, 1) circle [radius=1.5pt];
[thick] (4.75, 1) circle [radius=0.75];
[below] at (4.75, 1) $\mathbf{b}$;
[->, thick] (4.75, 1)–(4.75+0.866*0.75, 1+0.75/2);
at (4.75+0.866*0.375-0.15, 1+0.375/2+0.15) $\epsilon$;
(2, 2.5) circle [radius=1.5pt];s
[thick] (2, 2.5) circle [radius=0.75];
[below] at (2, 2.5) $\mathbf{a}$;
[->, thick] (2, 2.5)–(2+0.866*0.75, 2.5+0.75/2);
at (2+0.866*0.375-0.15, 2.5+0.375/2+0.15) $\epsilon$;
(8.5, 1.5) circle [radius=1.5pt];
[thick] (8.5, 1.5) circle [radius=0.75];
[below] at (8.5, 1.5) $\mathbf{c}$;
[->, thick] (8.5, 1.5)–(8.5+0.866*0.75, 1.5+0.75/2);
at (8.5+0.866*0.375-0.15, 1.5+0.375/2+0.15) $\epsilon$;
[-Stealth] (6.5, 3)–(5.5, 3);
[right] at (6.5, 3) $\Sigma$;

Two basins of attractions $\hat{A}$ and $\hat{C}$ divided by smooth basin boundary $\Sigma$ and example initial states $\mathbf{a}$, $\mathbf{b}$, and $\mathbf{c}$ with uncertainty $\epsilon$

§.§.§ Fractal Basin Boundaries

A given point $\mathbf{x}$ is on the boundary of a basin $\hat{A}$ if and only if every
neighborhood of $\mathbf{x}$ contains both points in $\hat{A}$ and points not in $\hat{A}$. Let us denote a basin boundary with $\Sigma$. Like before, it will be useful to quantify basin boundaries in a way that allows us to get a grasp on their complex geometry. To illustrate our method of quantifying basin boundaries, let us consider a region of the state space of a two-dimensional system with two attractors $A$ and $C$ (like the Hénon map). This region of state space contains a subset of the basins of both $A$ and $C$, which are divided by a basin boundary $\Sigma$. Therefore, if we were to choose a random initial state $\mathbf{x}_0$ in this region, it would lie in either basin $\hat{A}$ or $\hat{C}$ with a 100% probability ottboa. This is because the basin boundary $\Sigma$ has zero measure in two-dimensional state space as it is a one-dimensional object; from Property 4 of measures in Section <ref>, $\mu^2(\Sigma)=0$.

Now, let us say that this state $\mathbf{x}_0$ has some uncertainty $\epsilon$. Specifically, let us say that the actual state may be anywhere in a disk with radius $\epsilon$ centered at $\mathbf{x}_0$, or $|\mathbf{x}-\mathbf{x}_0|<\epsilon$.[For a general system, this would mean that the state may be anywhere in an $n$-dimensional ball with radius $\epsilon$ centered at $\mathbf{x}_0$.] In Figure <ref>, we display a schematic diagram of this region of state space with three initial states $\mathbf{a}$, $\mathbf{b}$, and $\mathbf{c}$. We can see that initial state the $\mathbf{a}$, despite having some uncertainty, will always be attracted to attractor $A$. Similarly, the initial state $\mathbf{c}$ will always be attracted to attractor $C$. However, we can see that the uncertainty $\epsilon$ causes the initial state $\mathbf{b}$ to have the possibility of being attracted to either $A$ or $C$. Picking a random state in our region of state space, we will now consider what the probability $\varrho(\epsilon)$ is that the uncertainty $\epsilon$ will cause us to make a mistake in predicting which attractor our random state will go to. As we can see in Figure <ref>, this probability $\varrho(\epsilon)$ is identical to the fraction of our region of state space that lies a distance $\epsilon$ from the basin boundary $\Sigma$. It is clear from Figure <ref> that this subset of states that lie a distance $\epsilon$ from $\Sigma$ is a strip of width $2\epsilon$ centered at $\Sigma$. Therefore, the fraction $\varrho(\epsilon)$ scales proportionally to $\epsilon$:

\begin{equation}
\varrho(\epsilon)\sim\epsilon
\label{eq:smooth-boundary-exponent}
\end{equation}

where $\sim$ means “is proportional to.”

However, Equation <ref> only applies to smooth boundaries, like the ones in Figures <ref> and <ref>. For fractal basin boundaries, namely, those with a non-integer dimension $d$, $\varrho(\epsilon)$ follows the power law

\begin{equation}
\varrho(\epsilon)\sim\epsilon^{\mathfrak{u}}
\label{eq:fractal-boundary-exponent}
\end{equation}

where $\mathfrak{u}$ is a number less than 1 known as the uncertainty exponent ottboa. mcdonald prove that the uncertainty exponent is related to a basin boundary's dimension $d$ by the following relation:

\begin{equation}
\mathfrak{u} = n-d
\label{eq:uncertainty-exp-frac-dim}
\end{equation}

where $n$ is the dimension of the system's state space. When $\mathfrak{u}=1$, Equation <ref> gives us $d=n-1$, but when $\mathfrak{u}<1$, the boundary is fractal because $d$ is then not an integer. When $\mathfrak{u}$ is substantially less than 1, a significant improvement in initial uncertainty $\epsilon$ causes only a marginal reduction in final state uncertainty, or which attractor a state will actually end up in grebogi-final-state. For example, if $\mathfrak{u}=0.1$, to reduce the final uncertainty $\varrho(\epsilon)$ by a factor of 10, we will have to reduce our initial uncertainty $\epsilon$ by a factor of $10^{10}$. Thus, an uncertainty exponent $\mathfrak{u}$ less than 1, indicating the presence of a fractal basin boundary, directly results in sensitive dependence on initial conditions [p. 151]ott, as it becomes increasingly difficult to predict where initial states will end up as we near the fractal boundary.

We are specifically interested in the behavior of $\varrho(\epsilon)$ when $\epsilon$ is small mcdonald, so to numerically calculate uncertainty exponents, we can find a bunch of $\varrho(\epsilon)$ values for $\epsilon = 2^{-k}$ and graph them on a plot of $\log_2\varrho(\epsilon)$ vs. $\log_2\epsilon$. Specifically, taking logs on both sides of Equation <ref>, we get

\begin{equation}
\ln\varrho(\epsilon) = \mathfrak{u}\ln\epsilon + c
\end{equation}

where $c$ is the natural log of the proportionality constant. We approximate a given $\varrho(2^{-k})$ value for the Hénon boundary using the code in Appendix <ref>, which implements a Monte Carlo algorithm. Namely, we choose a random initial state in $-5\leq x\leq 5$ and $-5\leq y\leq 5$ (shown in Figure <ref>), which we will denote as $\langle x_0,\, y_0 \rangle$, then test whether it belongs to the basin of the strange attractor or the basin of infinity by iterating it. Then, we test four perturbed states $\langle x_0+\epsilon,\, y_0 \rangle$, $\langle x_0-\epsilon,\, y_0 \rangle$, $\langle x_0,\, y_0+\epsilon \rangle$, and $\langle x_0,\, y_0-\epsilon \rangle$. If any of these states end up in a different attractor, then our initial state was uncertain. Using our standard parameters of the Hénon map ($a=1.4$ and $b=0.3$), the code in Appendix <ref> gives us $\mathfrak{u}=1$, which makes sense from Figure <ref> as the boundary appears smooth.

The basin of the Hénon attractor for $-5\leq x\leq 5$ and $-5\leq y\leq 5$ using the parameters $a=1.45$ and $b=0.3$, graphed with the code in Appendix <ref>

Step 1

Step 2

Two steps in the construction of the lakes of Wada, from pisarchikbook

However, if we consider the Hénon map with the parameters $a=1.45$ and $b=0.3$, we can see in Figure <ref> that its basin boundary is far more complex. Using the code in Appendix <ref> to calculate $\varrho(2^{-k})$ values and taking a linear regression (neglecting the first four points as we are interested in small values of $\epsilon$), we get that

\begin{equation}
\log_2\varrho(\epsilon) = 0.06954\log_2\epsilon - 3.050
\end{equation}

with an $R^2$ value of 0.984. Thus, $\mathfrak{u}\approx0.0695$ for our choice of parameters, indicating a large amount of final state uncertainty and “geometrical” chaos on the Hénon boundary. Visually, we can see in Figure <ref> why this is the case: any uncertainty near the boundary will significantly affect where an initial state will be attracted. From Equation <ref>, we can find that the dimension $d$ of the Hénon boundary for $a=1.45$ and $b=0.3$ is

\begin{equation}
d = n-\mathfrak{u} \approx 2-0.0695 = 1.9305
\end{equation}

which is significantly larger than 1. This indicates that the Hénon boundary for these parameters has a large amount of “roughness” and detail that significantly increases its fractal dimension from the $d=1$ dimension for the Hénon map's standard parameters, which is immediately obvious from Figures <ref> and <ref>.

§.§.§ Wada Basins

Let us consider an $n$-dimensional system with $w\geq3$ attractors $A_1,\,A_2,\,\hdots\,,\,A_w$ and associated basins $\hat{A}_1,\,\hat{A}_2,\,\hdots\,,\,\hat{A}_w$. In a system like this, it is possible for the basins to satisfy the counterintuitive Wada property, for which 3 or more basins $\hat{A}_1,\,\hat{A}_2,\,\hdots\,,\,\hat{A}_w$ share the same boundary $\Sigma$. Specifically, the basins are Wada if each boundary point of every basin is also a boundary point of every other basin kennedy. It is immediately evident that the presence of Wada basins is an even larger barrier to final state predictability than a standard fractal boundary because an initial uncertainty $d\mathbf{x}_0$ could potentially cause an initial state to end up in one of many distinct attractors. Thus, we want to be able to detect when basins exhibit the Wada property because it directly results in extreme sensitivity to initial conditions zhang.

To understand how it is possible for three basins to be Wada, we will use the lakes of Wada example described by kennedy. Imagine a white island surrounded by a red ocean. On the island, there are two lakes, one with blue water and one with green water. To construct the lakes of Wada, our first step, which is shown in Figure <ref>, is to dig a canal from the ocean into the island such that once the canal is dug, there are no points of land further than 1 unit away from red water. The second step, shown in Figure <ref>, is to dig a canal from the blue lake into the island such that each point in the remaining land is no further than $1/2$ units away from blue water. Similarly, the third step is to dig a canal from the green lake so that the remaining land is within $1/4$ units away from green water. After that, we return to the red ocean, digging a canal that makes the remaining land within $1/8$ units from red water, and so on, cycling through the three bodies of water. After an infinite number of steps, each point on the remaining land is arbitrarily close to all three bodies of water, so the remaining land is the boundary of all three. In other words, since every neighborhood of each land point contains points in all three bodies of water, the land is a Wada boundary of the three bodies.

This example clearly demonstrates how Wada boundaries are indicative of “geometrical” chaos due to extreme sensitivity to initial conditions. For this reason, we would like to have a method to detect the presence of the Wada property in any $n$-dimensional dynamical system with $w\geq3$ attractors. We will base our method of detecting Wada basins on the algorithm presented by daza, which uses a grid approach. This method goes as follows:

* Given an $n$-dimensional dynamical system with $w\geq3$ attractors, consider a “cubical” region of state space[We define a one-dimensional “cubical” region to be a line segment, a two-dimensional “cubical” region to be a square, a three-dimensional “cubical” region to be a standard cube, and so on.] that includes subsets of all the system's basins $\hat{A}_1,\,\hat{A}_2,\,\hdots\,,\,\hat{A}_w$.

* Cover the region with an $s\times s\times\hdots\times s$ (with $s$ being repeated $n$ times) grid $G$ composed of $s^n$ closed boxes $\Box_i$, where $i=1,\,2,\,\hdots\,,\,s^n$. It is beneficial to think of each box $\Box_i$ as an element of the set $G$. Mathematically, this means

\begin{equation}
\end{equation}

In a two-dimensional system, $G$ is an $s\times s$ grid, where the boxes $\Box_i$ are indexed from left to right, top to bottom.

* Define a function $C(\mathbf{x})$ that tells us which basin a given state belongs to. Specifically, for a state $\mathbf{x}\in\hat{A}_j$, define $C(\mathbf{x})=j$.

* For a given box $\Box_i$, define $C(\Box_i)=C(\mathbf{x})$, where $\mathbf{x}$ is the point at the center of $\Box_i$. Practically, it is useful to think of $C(\Box_i)$ as the color of $\Box_i$.

* Define $b(\Box_i)$ to be the set of boxes consisting of $\Box_i$ and all the boxes sharing at least 1 boundary point with $\Box_i$. Therefore, for a two-dimensional system, $b(\Box_i)$ is a $3\times3$ collection of boxes with $\Box_i$ as the central box:

\begin{equation}
\begin{split}
b(\Box_i) &= \bigcup_{j=-1}^1\bigcup_{k=-1}^1\Box_{i+js+k} \\
&= \Box_{i-s-1}\cup\Box_{i-s}\cup\Box_{i-s+1} \\
&\mathrel{\phantom{=}}\cup\,\Box_{i-1}\cup\Box_{i}\cup\Box_{i+1} \\
\end{split}
\end{equation}

* Define the function $K(\Box_i)$ to be the number of distinct colors $C(\Box)$ in $b(\Box_i)$ and calculate $K(\Box_i)$ for each box $\Box_i\in G$

* For every box with a $K(\Box_i)$ that is not 1 or $w$, draw a line segment from the center point of $\Box_i$ to the center point of a box in $b(\Box_i)$ that has a different color than $\Box_i$. Now, step $p=1$ is to calculate the color of the midpoint of this line segment. If this is a different color from all the other colors of the boxes in $b(\Box_i)$, there is now one more color in $b(\Box_i)$, so $K(\Box_i)$ increases by one. If $K(\Box_i)=w$ now, stop and move on to the next box with a $K(\Box_i)$ that is not 1 or $w$. Otherwise, move on to step $p=2$, which is to calculate the colors of two more points, namely, the points $1/4$ and $3/4$ up the line segment. Adjust $K(\Box_i)$ accordingly if either of these are new colors. Continue this process, step $p=3$ being dividing the line segment into eighths, step $p=3$ being dividing the line segment into sixteenths, and so on, until either $K(\Box_i)=w$ or the number of calculated points gets sufficiently large enough to assume that $K(\Box_i)<w$.

* For each step $p$ in 7. and for $k=1,\,2,\,\hdots\,,\,w$, define $G_k^p\subset G$ to be the set of all $\Box_i$ for which $K(\Box_i)=k$ immediately after step $p$. In other words, if $K(\Box_i)$ is calculated immediately after step $p$, then

\begin{equation}
G_k^p = \{\Box_i:\:K(\Box_i)=k\}
\end{equation}

After a given step $p$, this method sorts all the boxes in the grid $G$ into distinct sets:

\begin{equation}
G = \bigcup_{k=1}^w G_k^p
\end{equation}

By definition, the set $G_k^p$ that a given box $\Box_i$ belongs to tells us how many different colors $k$ are within $b(\Box_i)$ after $p$ steps. This means that, after an infinite number of steps, it indicates whether the box is inside a basin or on (or sufficiently near) a certain basin boundary. Specifically, $G_1^{\infty}$ is composed of boxes $\Box_i$ that are in the interior of a basin, as $k=1$ indicates that all the center points of all the boxes in a given $b(\Box_i)$ end up in the same attractor. $G_2^{\infty}$ is composed of boxes $\Box_i$ that are on or sufficiently near a boundary of two attractors since an infinite number of points in a neighborhood of the center point of $\Box_i$ led to two distinct attractors. Similarly, $G_3^{\infty}$ is composed of boxes that are on or sufficiently near a boundary of three attractors, and so on.

By definition, a system with $w$ attractors exhibits the Wada property if every basin boundary point is a boundary point of all $w$ basins of attraction. This means that for a system with Wada basins, every point in state space is either inside a basin or on a Wada boundary, meaning every box in the grid $G$ is either an element of $G_1^{\infty}$ or $G_w^{\infty}$. Mathematically, we can write this condition as

\begin{equation}
\lim_{p\to\infty}\sum_{k=2}^{w-1}\left|G_k^p\right|=0
\label{eq:wada-condition}
\end{equation}

where $\left|G_k^p\right|$ represents the number of boxes in, or the cardinality of, the set $G_k^p$.

We can label a system as Wada or not Wada by whether or not Equation <ref> is satisfied, but it will be useful to have a parameter to determine whether a system is partial Wada. A system with partially Wada basins contains some basin boundary points that are shared by all $w$ basins, but at least some boundary points that are not daza. To quantify this, we can define a parameter $W_k$ that tells us the fraction of boxes on or sufficiently near a boundary that belong to a boundary shared by $k$ basins:

\begin{equation}
\end{equation}

since all $G_j^{\infty}$ for $j\geq2$ contain only boxes on or sufficiently near a boundary. Clearly, if $W_w=1$, then the system is Wada, but if $0<W_w<1$, then the system is partial Wada by our established definition.

Attractors of the map in Equation <ref> with parameters $a=0.71$ and $b=0.9$, where the blue points are a 2-cycle attractor and the green points are a 6-cycle attractor, graphed with the code in Appendix <ref>

$-2.5\leq x\leq 2.5$, $-2.5\leq y\leq 2.5$

$0.6\leq x\leq 1.1$, $-0.75\leq y\leq -0.25$

Basins of attraction of the map in Equation <ref> with parameters $a=0.71$ and $b=0.9$, where the red area represents the basin of infinity, the blue area represents the basin of a 2-cycle attractor, and the green area represents the basin of a 6-cycle attractor, graphed with the code in Appendix <ref>

$p$     $G_1^p$     $G_2^p$     $G_3^p$

1     7186     1286     1528

2     7186     1117     1697

3     7186     892     1922

4     7186     732     2082

5     7186     581     2233

6     7186     512     2302

7     7186     475     2339

$p$     $G_1^p$     $G_2^p$     $G_3^p$
8     7186     446     2368

9     7186     410     2404

10     7186     385     2429

11     7186     347     2467

12     7186     324     2490

13     7186     280     2534

14     7186     250     2564

Some $G_k^p$ values for a $100\times 100$ grid over the region in Figure <ref>, calculated using the code in Appendix <ref>

[->] (0,0)–(0,15.5);
[->] (0,0)–(15,0);
(-2, 0)–(18, 0);
[thick] (0,15.04) circle [radius=0.15];
[thick] (1,12.86) circle [radius=0.15];
[thick] (2,11.17) circle [radius=0.15];
[thick] (3,8.92) circle [radius=0.15];
[thick] (4,7.32) circle [radius=0.15];
[thick] (5,5.81) circle [radius=0.15];
[thick] (6,5.12) circle [radius=0.15];
[thick] (7,4.75) circle [radius=0.15];
[thick] (8,4.46) circle [radius=0.15];
[thick] (9,4.10) circle [radius=0.15];
[thick] (10,3.85) circle [radius=0.15];
[thick] (11,3.47) circle [radius=0.15];
[thick] (12,3.24) circle [radius=0.15];
[thick] (13,2.80) circle [radius=0.15];
[thick] (14,2.50) circle [radius=0.15];

[below] at (8,-1.5) $p$;
[left, rotate=90] at (-2.25,7.85) $G_2^p$;
(0, 0)–(0, -0.2);
(1, 0)–(1, -0.2);
(2, 0)–(2, -0.2);
(3, 0)–(3, -0.2);
(4, 0)–(4, -0.2);
(5, 0)–(5, -0.2);
(6, 0)–(6, -0.2);
(7, 0)–(7, -0.2);
(8, 0)–(8, -0.2);
(9, 0)–(9, -0.2);
(10, 0)–(10, -0.2);
(11, 0)–(11, -0.2);
(12, 0)–(12, -0.2);
(13, 0)–(13, -0.2);
(14, 0)–(14, -0.2);

[below] at (0, -0.2) $0$;
[below] at (4, -0.2) $4$;
[below] at (8, -0.2) $8$;
[below] at (12, -0.2) $12$;

(0, 0)–(-0.2, 0);
(0, 1)–(-0.2, 1);
(0, 2)–(-0.2, 2);
(0, 3)–(-0.2, 3);
(0, 4)–(-0.2, 4);
(0, 5)–(-0.2, 5);
(0, 6)–(-0.2, 6);
(0, 7)–(-0.2, 7);
(0, 8)–(-0.2, 8);
(0, 9)–(-0.2, 9);
(0, 10)–(-0.2, 10);
(0, 11)–(-0.2, 11);
(0, 12)–(-0.2, 12);
(0, 13)–(-0.2, 13);
(0, 14)–(-0.2, 14);
(0, 15)–(-0.2, 15);
[left] at (-0.2, 0) $0$;
[left] at (-0.2, 2) $200$;
[left] at (-0.2, 4) $400$;
[left] at (-0.2, 6) $600$;
[left] at (-0.2, 8) $800$;
[left] at (-0.2, 10) $1000$;
[left] at (-0.2, 12) $1200$;
[left] at (-0.2, 14) $1400$;

A plot of some of the points in Table <ref> on a graph of $G_2^p$ vs. $p$, showing that $G_2^p\to0$ as $p\to\infty$

As an example of a dynamical system with three attractors, let us consider an alternative Hénon map defined by the function

\begin{equation}
\mathbf{f}(\mathbf{x};\,a,\,b) = \begin{pmatrix}
a-x^2+by \\
\end{pmatrix}
\label{eq:alternative-henon-map}
\end{equation}

with parameters $a=0.71$ and $b=0.9$. In this system, there are three attractors: a 2-cycle, a 6-cycle, and infinity nusse. The two periodic orbit attractors are shown in Figure <ref> using the code in Appendix <ref>, with the 2-cycle being represented by blue points and the 6-cycle being represented by green points. In Figure <ref>, we visualize the basins of this map using the code in Appendix <ref>, with the red area being the basin of infinity, the blue area being the basin of the 2-cycle, and the green area being the basin of the 6-cycle. Zooming in (Figure <ref>), we can see complex geometry on the boundary of the basins, not only being fractal, but also red, blue, and green points being intermingled with each other. For this reason, we suspect the basins might be Wada, so we implement our grid algorithm for detecting the Wada property in Appendix <ref> using a $100\times 100$ grid layered over the region in Figure <ref>. The results are shown in Table <ref>, and a plot $G_2^p$ against $p$ is shown in Figure <ref>. We can see that $G_2^p$ is tending to 0, so we can assume that Equation <ref> is satisfied. Therefore, all boundary points are Wada points, meaning this system exhibits the Wada property.

§.§.§ Milnor Attractors and Riddled Basins

In order to discuss the concept of riddled basins of attraction, we will first need to establish Milnor's definition of an attractor milnor, which is different from the definition we established in Section <ref>. We will consider a subset of state space $\mathcal{A}\subset\mathbb{R}^n$ to be a Milnor attractor if it satisfies the following two properties:

* The set of all points that $\mathcal{A}$ attracts, or its basin of attraction $\hat{\mathcal{A}}$, has a strictly positive measure. That is, $\mu(\hat{\mathcal{A}})>0$.

* There is no $\mathcal{A}'\subset\mathcal{A}$ such that $\hat{\mathcal{A}}'=\hat{\mathcal{A}}$.

Recall from Section <ref> that we required a standard attractor $A$ to attract an open set of initial conditions $U$. In Property 1 of Milnor attractors, we loosen this requirement to say that a Milnor attractor $\mathcal{A}$ attracts a region with positive measure, meaning the basin of $\mathcal{A}$ takes up a non-zero “volume” of state space but is not necessarily an open set. Another way of viewing this is that if $\mathcal{A}$ is a Milnor attractor, there exists a bounded region of state space such that if a state is randomly chosen from the region, there is a non-zero probability that the state will be attracted to $\mathcal{A}$ ott-riddled. Property 2 just ensures that every part of $\mathcal{A}$ plays a role in the attractor, meaning we can't take out a part of $\mathcal{A}$ and still attract the same set of initial conditions milnor-wiki.

When a certain kind of symmetry exists in a chaotic system, Milnor attractors of the system may be associated with basins that are riddled. Riddled basins of attraction occur in systems with, because of some kind of symmetry, a smooth invariant manifold, which is essentially a smooth surface in state space with the property that any initial condition in the surface will have a forward orbit that stays in the surface [p. 15]pisarchikbook. A riddled basin of attraction $\hat{\mathcal{A}}$ has the property that every point $\mathbf{a}\in\hat{\mathcal{A}}$ has points arbitrarily close to it belonging to another basin $\hat{\mathcal{C}}$ ott-riddled. Specifically, for any $\epsilon>0$ and any $\mathbf{a}\in\hat{\mathcal{A}}$, an $n$-dimensional ball centered at $\mathbf{a}$ with radius $\epsilon$, namely, $|\mathbf{x}-\mathbf{a}|<\epsilon$, will always contain a positive measure of points that get attracted to $\mathcal{C}$ instead of $\mathcal{A}$. In other words, even though $\hat{\mathcal{A}}$ has a positive measure, the set $\hat{\mathcal{A}}$ and its boundary set are the same because every neighborhood of every point in $\hat{\mathcal{A}}$ contains points both inside and outside of the set $\hat{\mathcal{A}}$ ottboa. It is now clear why riddled basins can only be associated with a Milnor attractor and an attractor by our previous definition. Since $\hat{\mathcal{A}}$ is not an open set, it doesn't qualify as the basin of a standard attractor $A$.

The practical implication of the riddled basin property is that we can never be sure whether a point $\mathbf{a}\in\hat{\mathcal{A}}$ will be attracted to $\mathcal{A}$ or some other attractor $\mathcal{C}$ because no matter how small we make the error $\delta\mathbf{a}$, there will always be a chance $\mathbf{a} + \delta\mathbf{a}$ lies in $\hat{\mathcal{C}}$ rather than $\hat{\mathcal{A}}$. For this reason, the uncertainty exponent of a riddled basin is $\mathfrak{u}=0$, so the presence of a riddled basin is the most extreme version of sensitive dependence on initial conditions resulting from a geometrical property of a chaotic system that we have discussed so far.

However, an even stronger geometrical property leading to sensitive dependence on initial conditions in a system is the presence of intermingled basins. Intermingled basins occur in $n$-dimensional systems with $i\geq 2$ Milnor attractors $\mathcal{A}_1,\, \mathcal{A}_2,\,\hdots\,,\, \mathcal{A}_i$. If all the system's basins of attractions $\hat{\mathcal{A}}_1,\, \hat{\mathcal{A}}_2,\,\hdots\,,\, \hat{\mathcal{A}}_i$ are riddled with each other, then we say that the basins are intermingled. More rigorously, this means that for any point $\mathbf{a}\in\mathbb{R}^n$, any $n$-dimensional ball centered around $\mathbf{a}$ will contain a positive measure of points in each of the basins $\hat{\mathcal{A}}_1,\, \hat{\mathcal{A}}_2,\,\hdots\,,\, \hat{\mathcal{A}}_i$, meaning, for each attractor $\mathcal{A}_1,\, \mathcal{A}_2,\,\hdots\,,\, \mathcal{A}_i$, there is a non-zero probability that a random initial state in the ball $|\mathbf{x}-\mathbf{a}|<\epsilon$ will go to that attractor alexander.

Proving that a given basin $\hat{\mathcal{A}}$ is riddled requires two steps alexander:

* Show that enough points are attracted to $\mathcal{A}$ to say that the probability of that happening is positive. In the paper by alexander, a method is described to do this using Lyapunov exponents.

* Show that sufficiently many points are repelled from $\mathcal{A}$. According to alexander, in systems with riddled basins, this is usually implied by symmetry.

Riddled basins can be found in models for a variety of different systems in many different fields. Some of these include electronic circuits, coupled chaotic oscillators, learning and mechanical systems, and interacting populations [p. 17]pisarchikbook. However, most ordinary systems do not contain smooth invariant manifolds, which are the requirement for riddled basins to exist, so riddled basins do not usually appear in practice ottboa. For this reason, we will not discuss in detail how to implement the two steps above, nor the mathematical theory behind smooth invariant manifolds. Our main purpose in discussing the basics behind riddled basins was simply to introduce the most extreme version of sensitivity to initial conditions emerging from geometrical aspects of dynamical systems. However, they will not be relevant in our analysis and discussion of the Rulkov maps.

§ SLOW-FAST SYSTEMS AND DYNAMICS

Slow-fast systems are dynamical systems that involve several variables evolving on different time scales zheng. Dynamical systems with separated slow and fast variables are invaluable for modeling many real-world applications, including lasers, chemical reactions, optoelectronic systems, and ecological systems omelchenko, zheng. For our purposes, however, we are interested in using slow-fast systems to model the behavior of biological neurons.

We will first discuss the continuous-time case of a two-variable slow-fast system. From Equation <ref>, a two-dimensional continuous-time dynamical system can be represented as

\begin{equation}
\begin{pmatrix}
\dot x \\
\dot y
\end{pmatrix}
\begin{pmatrix}
\frac{dx}{dt} \\[3pt]
\frac{dy}{dt}
\end{pmatrix}
\begin{pmatrix}
g^{[1]}(x,\,y) \\
g^{[2]}(x,\,y) \\
\end{pmatrix}
\end{equation}

This system can be transformed into the standard form of a slow-fast system as shown by [p. 70]ginoux by defining $g^{[1]}(x,\,y) = \chi(x,\,y)$ and $g^{[2]}(x,\,y) = \eta \omega(x,\,y)$, where $0<\eta\ll 1$. Then,

\begin{equation}
\begin{pmatrix}
\dot x \\
\dot y
\end{pmatrix}
\begin{pmatrix}
\chi(x,\,y) \\
\eta \omega(x,\,y) \\
\end{pmatrix}
\label{eq:slow-fast-def-fast-time}
\end{equation}

Assuming $\chi$ and $\omega$ are approximately on the same time scale, since $\eta$ is a small number, $\dot y$ is small at any given time, meaning $y$ changes slowly. For this reason, we call $x$ the fast variable and $y$ the slow variable. Often, we can approximate the dynamics of a slow-fast system by analyzing the fast and slow dynamics separately, treating $y$ as a slowly drifting parameter of $\frac{dx}{dt}$. However, to accurately describe the system's dynamics, we must analyze the slow and fast dynamics together as a two-dimensional map.

While slow-fast dynamics in continuous-time dynamical systems have been studied extensively, discrete-time slow-fast systems have received far less attention. In fact, an agreed-upon form of this type of system has not yet even been established. However, we will define a slow-fast discrete-time dynamical system as one governed by the iteration function

\begin{equation}
\begin{pmatrix}
x_{k+1} \\
\end{pmatrix}
\begin{pmatrix}
f(x_k,\,y_k) \\
y_k + \eta g(x_k,\,y_k)
\end{pmatrix}
\end{equation}

where $0<\eta\ll 1$. In this case, $x$ is the fast variable and $y$ is the slow variable, as each iteration changes $y$ by only a small amount. In this section, we will explore the dynamics and relevance of both continuous-time and discrete-time slow-fast systems, with a specific application to the dynamics and properties of neuronal models.

§.§ Neimark-Sacker Bifurcations