Border-collision bifurcations from stable fixed points to any number of coexisting chaotic attractors.

Piecewise-smooth maps have different functional forms in different parts of phase space. As a parameter of a piecewise-smooth map is varied, a bifurcation occurs when a fixed point collides with a switching manifold, where the functional form of the map changes. This type of bifurcation is termed a border-collision bifurcation (BCB).

BCBs have been identified in diverse applications. Classical examples include power converters, where BCBs can cause the internal dynamics of the converter to suddenly become quasi-periodic or chaotic [1, 2], and mechanical systems with friction, where BCBs can induce recurring transitions between sticking and slipping motion [3, 4]. More recently in [5] it is shown how BCBs can explain changes to the frequency of annual influenza outbreaks.

After the popularisation of BCBs by Nusse and Yorke in [6], it was quickly realised that BCBs often represent remarkably complicated transitions, equivalent to the amalgamation of several (even infinitely many) smooth bifurcations [7, 8, 9, 10]. It is natural to then ask, how complicated can the transition be? At a BCB the local attractor of a system can change from a stable fixed point to a higher period solution, a quasiperiodic solution, or a chaotic solution [11]. It can also split into several attractors. These attractors are created simultaneously and grow out of a single point. If the parameter that affects the bifurcation is varied dynamically, then in the presence of arbitrarily small noise it cannot be known a priori which attractor the system will transition to, although the size of the basin of attraction of an attractor can be expected to correlate positively with the likelihood that it will be selected [12, 13].

Examples of BCBs creating multiple attractors have been described by many authors [8, 14, 15, 16]. Examples of BCBs creating arbitrarily many or infinitely many attractors was described in [17, 18, 19]. The first numerical example of a BCB creating multiple chaotic attractors is possibly that given in Section 7 of Avrutin et. al. [20]. More recently Pumariño et. al. [21] showed that two-dimensional, piecewise-linear maps can exhibit $2^{m}$ coexisting chaotic attractors for any $m\geq 1$, and under a coordinate transformation their maps are equivalent to members of the border-collision normal form.

The purpose of this paper is to demonstrate these complexities further. Our approach extends that of Glendinning [22] who studied perturbations of the $n$-dimensional border-collision normal form for which the $n^{\rm th}$ iterate is a direct product of identical skew tent maps. But whereas Glendinning [22] used skew tent maps that have a chaotic attractor consisting of one interval, here we use skew tent maps that have a chaotic attractor consisting of $k$ disjoint intervals. This novelty generates multiple attractors. A similar strategy was employed by Wong and Yang [23] to get two chaotic attractors in the two-dimensional case.

Unlike previous works we take the extra step of proving that the bifurcation phenomenon is not an artifact of the piecewise-linear nature of the border-collision normal form. We show the phenomenon persists when nonlinear terms are added to the normal form. For a generic piecewise-smooth map (representing a mathematical model), near a BCB the map is conjugate to a member of the normal form plus nonlinear terms [11].

To prove chaos we show some iterate of the map is piecewise-smooth and expanding. Immediately this implies every Lyapunov exponent is positive, but stronger results have been obtained in the context of ergodic theory. Piecewise-$C^{2}$ expanding maps generically (i.e. on an open, dense subset within the space of all such maps) have at least one invariant measure that is absolutely continuous with respect to the Lebesgue measure [24]. The requirement that each piece is $C^{2}$ is not satisfied for BCBs that correspond to grazing-sliding bifurcations (where the quadratic tangency of the grazing trajectory of the underlying system of differential equations induces an order-$\frac{3}{2}$ error term). In this case one can look to [25] for more general results regarding invariant measures. In the two-dimensional case, if the map is piecewise-analytic the genericity condition is not needed [26, 27]. Analogous results for piecewise-linear maps are described in [28, 29].

The remainder of the paper is organised as follows. First in §2 we formally state the main results: Theorem 2.1 for two-dimensional maps and Theorem 2.2 for $n$-dimensional maps, for any $n\geq 2$. In §3 we introduce a simple form for $n$-dimensional maps about which perturbations will be performed, and show how nonlinear terms can be accommodated. In §4 we perform the straight-forward task of demonstrating that our class of perturbed maps have an asymptotically stable fixed point on one side of the BCB. Then in §5 we show that on the other side of the BCB the $n^{\rm th}$ iterate of the simple form is a direct product of identical skew tent maps.

In §6 we review attractors of skew tent maps focussing on chaotic attractors that are comprised of $k\geq 2$ disjoint intervals. In §7 we fatten these intervals to obtain trapping regions for the skew tent maps. Then in §8 we take Cartesian products of intervals to obtain boxes, where unions of the boxes form trapping regions for the $n$-dimensional maps. In §9 we use Burnside’s lemma and combinatorical arguments to derive an explicit formula for the number of trapping regions that result from this construction as a function of $n$ and $k$. We believe that each trapping region contains a unique attractor for sufficiently small perturbations, but a proof of this remains for future work. In the case $n=2$ we can obtain any number of trapping regions.

Then in §10 we prove that some iterate of the map is piecewise-smooth and expanding, and in §11 collate the results to prove Theorems 2.1 and 2.2. Lastly §12 provides some final remarks.

In this section we state our main results on the genericity of BCBs where a stable fixed point bifurcates into several chaotic attractors. Throughout the paper ${\rm int}(\cdot)$ denotes the interior of a set.

For any trapping region $\Omega$, the set $\bigcap_{i\geq 0}f^{i}(\Omega)$ is an attracting set, by definition [30]. A topological attractor is then a subset of an attracting set that is dynamically indivisible, in some sense. Different authors use different definitions for this indivisibility constraint (e.g. there exists a dense orbit) [31]. For this paper we just require the fact that any trapping region contains at least one attractor.

To motivate the statement of the next definition, note that in Definition 2.1 the restriction of $f$ to $\Omega$ is a map $f:\Omega\to\Omega$.

The dynamics near any generic BCB in two dimensions can be described by a map of the form

$$x\mapsto\begin{cases}\begin{bmatrix}\tau_{L}x_{1}+x_{2}+\mu\\ -\delta_{L}x_{1}\end{bmatrix}+E_{L}(x;\mu),&x_{1}\leq 0,\\ \begin{bmatrix}\tau_{R}x_{1}+x_{2}+\mu\\ -\delta_{R}x_{1}\end{bmatrix}+E_{R}(x;\mu),&x_{1}\geq 0,\end{cases}$$

(2.1)

where $x=(x_{1},x_{2})$ is the state variable, $\mu\in\mathbb{R}$ is the bifurcation parameter, and $E_{L}$ and $E_{R}$ are $C^{1}$ and $\mathpzc{o}\mathopen{}\mathclose{{}\left(\|x\|+|\mu|}\right)$ [11]. The last condition means $\frac{\|E_{L}(x;\mu)\|}{\|x\|+|\mu|}\to 0$ as $(x;\mu)\to(0;0)$ (regardless of how the limit is taken), and similarly for $E_{R}$. Since BCBs are local, only one switching condition is relevant and coordinates have been chosen so that it is $x_{1}=0$. The BCB of (2.1) occurs at the origin $x=0$ when $\mu=0$. The nonlinear terms $E_{L}$ and $E_{R}$ often have no qualitative effect on the bifurcation (as shown below for our setting) and by dropping these terms we obtain a piecewise-linear family of maps known as the two-dimensional border-collision normal form. Notice this form has four parameters $\tau_{L},\delta_{L},\tau_{R},\delta_{R}\in\mathbb{R}$, in addition to $\mu$.

Theorem 2.1 is proved in §11. Fig. 1 shows a typical example with $N=2$. This figure is for (2.1) with

$\displaystyle\tau_{L}$

$\displaystyle=-0.02,$

$\displaystyle\delta_{L}$

$\displaystyle=-0.62,$

$\displaystyle\tau_{R}$

$\displaystyle=-0.02,$

$\displaystyle\delta_{R}$

$\displaystyle=3,$

(2.2)

and

$\displaystyle E_{L}(x;\mu)$

$\displaystyle=\begin{bmatrix}-x_{1}^{2}\\ 0\end{bmatrix},$

$\displaystyle E_{R}(x;\mu)$

$\displaystyle=\begin{bmatrix}0\\ 0\end{bmatrix}.$

(2.3)

The bifurcation diagram, panel (a), shows that as the value of $\mu$ passes through $0$, a stable fixed point turns into two coexisting attractors (coloured orange and black). Panel (b) shows these attractors in phase space for one value of $\mu$. Each attractor appears to be two-dimensional, like in [23], but the orange attractor has three connected components while the black attractor has six connected components. Numerically we observe the orange attractor is destroyed in a crisis [32] at $\mu\approx 0.01$ after which all forward orbits converge to the black attractor. This example is explained further in §12.

We now consider BCBs in maps with more than two dimensions. The form (2.1) generalises from $\mathbb{R}^{2}$ to $\mathbb{R}^{n}$ ($n\geq 2$) as

$$x\mapsto f(x;\mu)=\begin{cases}C_{L}x+e_{1}\mu+E_{L}(x;\mu),&x_{1}\leq 0,\\ C_{R}x+e_{1}\mu+E_{R}(x;\mu),&x_{1}\geq 0.\end{cases}$$

(2.4)

where again $E_{L}$ and $E_{R}$ are $C^{1}$ and $\mathpzc{o}\mathopen{}\mathclose{{}\left(\|x\|+|\mu|}\right)$, and

$\displaystyle C_{L}$

$\displaystyle=\begin{bmatrix}c^{L}_{1}&1\\ c^{L}_{2}&&\ddots\\ \vdots&&&1\\ c^{L}_{n}\end{bmatrix},$

$\displaystyle C_{R}$

$\displaystyle=\begin{bmatrix}c^{R}_{1}&1\\ c^{R}_{2}&&\ddots\\ \vdots&&&1\\ c^{R}_{n}\end{bmatrix}$

(2.5)

are companion matrices whose first columns are vectors $c^{L},c^{R}\in\mathbb{R}^{n}$. In (2.4) and throughout the paper we write $e_{j}$ for the $j^{\rm th}$ standard basis vector of $\mathbb{R}^{n}$. We now generalise Theorem 2.1 from $\mathbb{R}^{2}$ to $\mathbb{R}^{n}$.

Theorem 2.2 is proved in §11. In (2.6) the sum is over all divisors of $a$, and $\varphi(d)$ is (by definition) the number of positive integers that are less than or equal to $d$ and coprime to $d$. For example with $k=10$ and $n=6$, we have $a=3$ whose divisors are $1$ and $3$. Since $\varphi(1)=1$ and $\varphi(3)=2$, we have

$$N[10,6]=\frac{1}{60}\mathopen{}\mathclose{{}\left(1\cdot 10^{6}+2\cdot 10^{2}}\right)=16670.$$

Table 1 lists the values of $N[k,n]$ for small values of $k$ and $n$. If $k$ is a multiple of $n$, then $a=1$ and $N=\frac{k^{n-1}}{n}$. If $n$ is prime and $k$ is not a multiple of $n$, then $a=n$ has only two divisors: $1$ and $n$. Here $\varphi(n)=n-1$ so $N=\frac{k^{n-1}+n-1}{n}=\mathopen{}\mathclose{{}\left\lceil\frac{k^{n-1}}{n}}\right\rceil$. Interestingly this implies $k^{n-1}-1$ is a multiple of $n$ which is Fermat’s little theorem; compare [33, 34]. Also $N[k,2]=\mathopen{}\mathclose{{}\left\lceil\frac{k}{2}}\right\rceil$ and $N[2,n]$ is Sloane’s integer sequence $A000016$ [35] which arises in coding problems [36]; see [37] for its occurrence in another dynamical systems setting.

In the absence of the nonlinear terms $E_{L}$ and $E_{R}$, (2.4) is the $n$-dimensional border-collision normal form, first considered in [38]. In this section we introduce a two-parameter subfamily of the normal form about which perturbations will be taken.

Given parameters $a_{L},a_{R}\in\mathbb{R}$ and $\sigma=\pm 1$, we form the $n$-dimensional map

$$g(y;a_{L},a_{R},\sigma)=\begin{cases}A_{L}y+\sigma e_{1}\,,&y_{1}\leq 0,\\ A_{R}y+\sigma e_{1}\,,&y_{1}\geq 0,\end{cases}$$

(3.7)

where

$\displaystyle A_{L}$

$\displaystyle=\begin{bmatrix}0&1\\ \vdots&&\ddots\\ 0&&&1\\ a_{L}\end{bmatrix},$

$\displaystyle A_{R}$

$\displaystyle=\begin{bmatrix}0&1\\ \vdots&&\ddots\\ 0&&&1\\ a_{R}\end{bmatrix}.$

(3.8)

The matrices $A_{L}$ and $A_{R}$ are companion matrices (2.5) for which $c^{L}$ and $c^{R}$ are scalar multiples of $e_{n}$. Specifically $c^{L}=d^{L}$ and $c^{R}=d^{R}$ where

$\displaystyle d^{L}$

$\displaystyle=a_{L}e_{n}\,,$

$\displaystyle d^{R}$

$\displaystyle=a_{R}e_{n}\,.$

Let us now explain the utility of (3.7). To a map of the form (2.4), we perform the spatial scaling $y=\frac{x}{|\mu|}$, assuming $\mu\neq 0$, to produce the map

$$y\mapsto\frac{f(|\mu|y;\mu)}{|\mu|}=\tilde{f}(y;\mu).$$

(3.9)

This map can be written as

$$\tilde{f}(y;\mu)=g(y;a_{L},a_{R},{\rm sgn}(\mu))+\begin{cases}(c^{L}-d^{L})y_{1}+\frac{E_{L}(|\mu|y;\mu)}{|\mu|},&y_{1}\leq 0,\\ (c^{R}-d^{R})y_{1}+\frac{E_{R}(|\mu|y;\mu)}{|\mu|},&y_{1}\geq 0.\end{cases}$$

(3.10)

By choosing $\mu$ small and $c^{L}$ and $c^{R}$ close to $d^{L}$ and $d^{R}$, we can make the difference between (2.4), in scaled coordinates, as close to the simple form (3.7) as we like. Formally we have the following result that follows immediately from the assumption that $E_{L}$ and $E_{R}$ are $C^{1}$ and $\mathpzc{o}\mathopen{}\mathclose{{}\left(\|x\|+|\mu|}\right)$.

In this section we generalise Section 2 of Glendinning [22] to accommodate higher order terms.

For the remainder of the paper we consider $g$ with $\sigma=1$. Here we show that the $n^{\rm th}$ iterate of $g$ is conjugate to a direct product of skew tent maps. To this end, we let

$$h(z;a_{L},a_{R})=\begin{cases}a_{L}z+1,&z\leq 0,\\ a_{R}z+1,&z\geq 0,\end{cases}$$

(5.12)

be a family of skew tent maps.

A detailed analysis of the skew tent map family (5.12) was done in Ito et. al. [39], see also Maistrenko et. al. [40]. For any $a_{L},a_{R}\in\mathbb{R}$, (5.12) has non-negative Schwarzian derivative almost everywhere so has at most one attractor [41]. Fig. 2 catagorises this attractor throughout the $(a_{L},a_{R})$-parameter plane. It contains regions $P_{k}$, for $k\geq 1$, where there exists a stable period-$k$ solution. It also contains regions $Q_{k}$, for $k=1$ and $k=2^{\ell}$ for $\ell\geq 3$, regions $R_{k}$, for even $k\geq 4$, and regions $S_{k}$, for $k\geq 2$, where there exists a chaotic attractor consisting of $k$ disjoint closed intervals.

As evident in Fig. 2, each $P_{k}$ is situated below $P_{k-1}$. For all $k\geq 3$, $R_{2k}$ is narrow and located immediately to the right of $P_{k}$, while $S_{k}$ is similarly narrow and located immediately to the right of $R_{2k}$.

In this paper we perturb about instances of (3.7) for which the pair $(a_{L},a_{R})$ belongs to some region $S_{k}$. In comparison Wong and Yang [23] use $(a_{L},a_{R})=(0.5,-2)$, which lies on the boundary of $P_{2}$ and $R_{4}$, while Glendinning allows any $\frac{6}{11}<a_{L}<\frac{10}{11}$ and $a_{R}=-2$, which includes points in $R_{4}$, $S_{2}$, and $Q_{1}$. All constructions use $|a_{L}|<1$ to ensure a stable fixed point for small $\mu<0$ by Lemma 4.1.

Each $S_{k}$ is bounded by three smooth curves. The upper boundary curve is where $h^{k}(0)=0$ and is given by

$$a_{R}=-\frac{1-a_{L}^{k-1}}{(1-a_{L})a_{L}^{k-2}}.$$

(6.14)

The left and right boundary curves of $S_{k}$ are given by

	$\displaystyle a_{L}^{2k-2}a_{R}^{3}+a_{L}-a_{R}$	$\displaystyle=0,$		(6.15)
	$\displaystyle a_{L}^{k-1}a_{R}^{2}+a_{R}-a_{L}$	$\displaystyle=0,$		(6.16)

respectively. Some explanation for these is provided below.

The ordering (6.17) is illustrated for $k=3$ in Fig. 3. Lemma 6.17 is a consequence of calculations done in [39, 40] and, although a little tedious, is not difficult to derive directly from (6.14)–(6.16). Note the ordering (6.17) also holds throughout $R_{2k}$.

We now use the forward orbit of $0$ to define intervals

$$\begin{split}I_{0}&=\mathopen{}\mathclose{{}\left[h^{k+1}(0),h(0)}\right],\\ I_{i}&=\mathopen{}\mathclose{{}\left[h^{i+1}(0),h^{i+k+1}(0)}\right],~{}\text{for all $i=1,\ldots,k-1$},\end{split}$$

(6.18)

see again Fig. 3. The next result follows immediately from (6.17).

By Lemma 6.2 orbits cycle through the intervals $I_{i}$, one of which ($I_{k-1}$) contains the critical point $z=0$. Thus the restriction of $h^{k}$ to any of these intervals is a continuous piecewise-linear map with two pieces, see Fig. 4. That is, $h^{k}$ is conjugate to an instance of the skew tent map family (5.12). Specifically, $h^{k}$ is conjugate to $h(z;\tilde{a}_{L},\tilde{a}_{R})$, where $\tilde{a}_{L}=a_{L}^{k-2}a_{R}^{2}$ and $\tilde{a}_{R}=a_{L}^{k-1}a_{R}$, because in each cycle orbits undergo either one or two iterations under the right piece of $h$, and the remaining iterations under the left piece of $h$.

If $(\tilde{a}_{L},\tilde{a}_{R})\in Q_{1}$ then $h(z;\tilde{a}_{L},\tilde{a}_{R})$ is transitive on $\tilde{I}=[h^{2}(0),h(0)]$. Refer to [42] (Lemmas 2.1 and 2.2) for a simple proof of this that assumes $\tilde{a}_{L}>1$, as is the case here. This implies $h(z;a_{L},a_{R})$ is transitive on $I_{0}\cup\cdots\cup I_{k-1}$. Since the value of $\tilde{a}_{L}$ is relatively large, to be in $Q_{1}$ we just need $(\tilde{a}_{L},\tilde{a}_{R})$ to lie between the two boundaries of $Q_{1}$ that are labelled in Fig. 2. This is because below $a_{L}a_{R}+a_{L}-a_{R}=0$ the interval $\tilde{I}$ is not forward invariant, while above $a_{L}a_{R}^{2}+a_{R}-a_{L}=0$ transitivity fails because typical orbits cannot reach a neighbourhood of the fixed point in $\tilde{I}$.

By replacing $a_{L}$ and $a_{R}$ with $\tilde{a}_{L}$ and $\tilde{a}_{R}$ in the formulas for these boundaries we produce (6.15) and (6.16). This explains the particular formulas (6.15) and (6.16) and gives the following result.

We now fatten the intervals $I_{i}$ (6.18) to form a trapping region for the attractor $I_{0}\cup\cdots\cup I_{k-1}$ of the one-dimensional map $h$.

To motivate the next construction, recall from Lemma 5.13 that $g^{n}$ is conjugate (via a translation) to a direct product of $n$ copies of $h$. So by Lemma 7.20 to obtain trapping regions for $g^{n}$ we can take unions of Cartesian products of the $J_{i}$, some shifted by $-1$ to account for the translation in (5.13). But to obtain trapping regions for $g$ we have to work a bit harder.

Write

	$\displaystyle J_{i}$	$\displaystyle=[p_{i},q_{i}],$
	$\displaystyle h\mathopen{}\mathclose{{}\left(J_{i-1\,\text{mod}\,k}}\right)$	$\displaystyle=[r_{i},s_{i}],$

for all $i=0,1,\ldots,k-1$. Observe $p_{i}<r_{i}<s_{i}<q_{i}$ for all $i$ by (7.20). Now let

$\displaystyle t_{i,j}$

$\displaystyle=p_{i}+\frac{(r_{i}-p_{i})(j-1)}{n},$

$\displaystyle u_{i,j}$

$\displaystyle=q_{i}+\frac{(s_{i}-q_{i})(j-1)}{n},$

(7.21)

and

$$K_{i,j}=\mathopen{}\mathclose{{}\left[t_{i,j}-1,u_{i,j}-1}\right],$$

(7.22)

for all $i=0,1,\ldots,k-1$ and $j=2,3,\ldots,n$. The next result uses the following notation: given $Z\subset\mathbb{R}$ and $a\in\mathbb{R}$, we write $Z+a$ to abbreviate $\{z+a\,|\,z\in Z\}$. Equation (7.23) is an immediate consequence of the ordering illustrated in Fig. 5.

We write $\mathbb{Z}_{k}$ for the set $\{0,1,\ldots,k-1\}$ with addition taken modulo $k$. Given a vector $v\in\mathbb{Z}_{k}^{n}$ (so $v_{j}\in\mathbb{Z}_{k}$ for each $j=1,2,\ldots,n$), we let

$$\Phi_{v}=J_{v_{1}}\times K_{v_{2},2}\times K_{v_{3},3}\times\cdots\times K_{v_{n},n}\,,$$

(8.24)

be a box in $\mathbb{R}^{n}$. By Lemma 7.2 and the definition of $g$, each such box maps under $g$ to the interior of another such box. Specifically $\Phi_{v}$ maps to the interior of $\Phi_{\psi(v)}$, where the map $\psi:\mathbb{Z}_{k}^{n}\to\mathbb{Z}_{k}^{n}$ is defined by

$$\psi(v)=\mathopen{}\mathclose{{}\left(v_{2},v_{3},\ldots,v_{n},v_{1}+1}\right).$$

(8.25)

Formally we have the following result.

We now use Lemma 3.1 to extend this result to maps of the form (2.4). Here we use the following notation: given $\Omega\subset\mathbb{R}^{n}$, we write $\mu\Omega$ to abbreviate $\{\mu y\,|\,y\in\Omega\}$.

The orbit of $v\in\mathbb{Z}_{k}^{n}$ under $\psi$ is the set

$${\rm orb}(v)=\mathopen{}\mathclose{{}\left\{\psi^{i}(v)\,\middle|\,i\geq 0}\right\}.$$

(9.28)

Given $v\in\mathbb{Z}_{k}^{n}$ let

$$T_{v}=\bigcup_{w\in{\rm orb}(v)}\Phi_{w}\,.$$

(9.29)

Then the scaled set $\mu T_{v}$ is a trapping region for any map $f$ that satisfies the conditions of Lemma 8.2. The number of mutually disjoint trapping regions given by this construction is equal to the number of orbits of $\psi$. The purpose of this section is to prove the following result.

Fig. 6 shows an example with $k=4$ and $n=3$. This figure is for the simple form (3.7) in three dimensions with $(a_{L},a_{R})=(0.47,-10)\in S_{4}$. By Proposition 9.1 the number of trapping regions (9.29) is $N[4,3]=6$. Numerically we observe each trapping region has a single three-dimensional attractor. Five of the trapping regions are comprised of $12$ boxes. The sixth (black in Fig. 6) is comprised of four boxes and corresponds to the orbit

$$\begin{split}v&=(2,1,0),\\ \psi(v)&=(1,0,3),\\ \psi^{2}(v)&=(0,3,2),\\ \psi^{3}(v)&=(3,2,1),\\ \psi^{4}(v)&=v.\end{split}$$

To prove Proposition 9.1 we first establish three lemmas.