Basin Entropy and Shearless Barrier Breakup in Open Non-Twist Hamiltonian Systems

1. Introduction

The study of non-integrable Hamiltonian system is one of the main disciplines in the field of nonlinear dynamics [1]. A large number of Hamiltonian systems of physical interest can be analytically and numerically investigated through the use of area-preserving mappings in a Poincaré surface of section of the phase space. In particular, there have been most extensively studied those systems satisfying the so-called twist property for which, loosely speaking, there are no orbits with the same frequency in terms of the corresponding angle variables in phase space [2]. One outstanding example of area-preserving mappings satisfying the twist property is the Chirikov-Taylor map [3].

One major advantage of working with area-preserving twist maps is that many powerful results of Hamiltonian theory like KAM and Poincaré-Birkhoff theorems and Aubry-Mather theory are valid provided the system obeys the twist property [2]. However, in various Hamiltonian systems of physical interest, chiefly fluids and plasmas, the twist property fails to be satisfied, what has motivated the study of area-preserving non-twist maps [4]. One paradigmatic example of the latter category is the standard non-twist map introduced by del Castillo-Negrete and Morrison, for which the twist condition is locally violated [5].

The non-twist character of such maps has a profound influence on their dynamical properties. For example, due to the non-monotonic character of the frequency profile, there appear twin island chains failing to exhibit the well-known island overlapping. Instead, the islands suffer a kind of reconnection process and produce robust shearless transport barriers, that modify the transport properties displayed by non-twist maps [6].

Among the wide variety of physically relevant systems described by area-preserving non-twist maps we mention the magnetic field line structure in tokamaks and stellarators [7,8,9], planetary orbits [10,11], stellar pulsations [12], atomic physics [13,14], condensed matter [15] and sheared geostrophic flows [5,16]. In all these non-twist systems the existence of shearless transport barriers represent a local obstruction to chaotic diffusion of phase space trajectories. One of the relevant problems involving non-twist systems is how to characterize numerically the destruction of those shearless barriers [17]. This problem has been investigated in considerable detail for the standard non-twist map, thanks to a special property that is the existence of the so-called indicator points [18].

For general area-preserving non-twist maps, however, the absence of such indicator points makes the numerical task of determining the breakup of the shearless barrier a difficult one. Some methods have been proposed for this task. In the present paper, we propose an alternative method to determine the shearless curve breakup by using a definition of entropy applied to basins of behavior [19,20]. The latter is an extension, for general dynamical systems, of the usual basin of attraction concept. Since Hamiltonian systems do not possess attractors, we can define an analogous behavior by opening its domain and considering the escape of trajectories. In this sense the basin of escape is the set of initial conditions (in the Poincaré surface of section) generating trajectories that escape through that exit.

Due to the underlying structure of the dynamics in a chaotic orbit of a non-integrable Hamiltonian system, the structure of escape basins is highly fractal. The basin entropy quantifies the uncertainty related to the fractality of the escape basins and of their common boundary, and has been used in many Hamiltonian systems with this purpose. In the present work, however, we used the basin entropy specifically to determine for which value of the perturbation parameter (measuring, so as to speak, the strength of the system non-integrability) the shearless barrier suffers a breakup.

This paper is structured as follows: in Section 2 we introduce the specific area-preserving non-twist map used. In Section 3 we exhibit the escape basin structure and its characterization using an uncertainty exponent. Section 4 introduces the concepts of basin entropy and basin boundary entropy for the escape basins of an open Hamiltonian system. In Section 5 we show how to use these concepts to characterize the shearless barrier breakup. Our Conclusions are left to the final Section.

2. Area-preserving non-twist maps

Let us consider a Hamiltonian system with N degrees of freedom, characterized by canonical pairs $(p_i,q_i)$, $i=1,2,\dots N$. This is an integrable system if one can obtain a symplectic transformation to action-angle coordinates

$$(I,\theta )\to \left(q\right(I,\theta ),p(I,\theta \left)\right),$$

(1)

where $I=\{I_1,I_2,\dots I_n\}\in B\subset \mathbb{N}$ (B is an open set) and $\theta =\{({\varphi }_1,{\varphi }_2,\dots {\varphi }_N)mod2\pi \}$, such that ${\varphi }_i$ parameterize the motion on a N-dimensional torus. In terms of these action-angle coordinates the Hamiltonian becomes

$$H\left(q\right(I,\theta ),p(I,\theta \left)\right)=K\left(I\right),$$

(2)

and the corresponding Hamilton´s equations are

$$\begin{array}{cc}\hfill \frac{dI}{dt}& =-\frac{\partial K\left(I\right)}{\partial \theta }=0,\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill \frac{d\theta }{dt}& =\frac{\partial K\left(I\right)}{\partial I}=\omega \left(I\right),\hfill \end{array}$$

(4)

where $\omega \in \left\{({\omega }_1,{\omega }_2,\dots {\omega }_N)\right\}$ are the frequencies corresponding to each irreducible circuit on the N-torus.

Since a $N=1$ degree of freedom system with time-independent Hamiltonian is always integrable, the simplest non-integrable systems have $N=2$. In addition, let us consider that the non-integrability comes from a weak perturbation of an integrable system, in the standard form

$$H(I_1,I_2;{\theta }_1,{\theta }_2)=H_0(I_1,I_2)+\epsilon H_1(I_1,I_2;{\theta }_1,{\theta }_2),$$

(5)

where $\epsilon \ll 1$ for a quasi-integrable system. The integrable system is characterized by two frequencies ${\omega }_1(I_1,I_2)$ and ${\omega }_2(I_1,I_2)$. The Hamiltonian $H_0$ is said to satisfy the so-called twist condition if

$$\frac{\partial {\omega }_i}{\partial I_j}=\frac{\partial H_0}{\partial I_j\partial I_i}\ne 0,$$

(6)

i.e. the integrable system does not have two phase-space trajectories with the same frequency. Conversely, if this condition fails to be satisfied at any value of the action, the system is said to be non-twist.

Many Hamiltonian system of physical interest satisfy the twist property. Moreover, for this class of systems many powerful results of Hamiltonian theory are valid, such as Kolmogorov-Arnold-Moser theorem, Poincaré-Birkhoff theorem, Aubry-Mather theory, and so on. On the other hand, there are non-twist systems of interest mainly in hydrodynamics and plasma physics, as commented in the Introduction. For non-twist systems there are different dynamical properties that have been investigated over the last years.

The energy $H=E$ is a constant of the motion inasmuch the Hamiltonian does not depend explicitly on time. Hence the motion, which occurs in a 4-dimensional phase space, actually is limited to a 3-dimensional energy surface $H=H(I_1,E;{\theta }_1,{\theta }_2)$. Moreover, considering a Poincaré surface of section ${\theta }_2=const.$$mod2\pi$ we can reduce the continuous-time flow generated by solving (3)-(4) to a discrete-time mapping in the plane $I_1\times {\theta }_1$, with the general form

$$\begin{array}{cc}\hfill I_{n+1}& =I_n+\epsilon h({\theta }_n,I_{n+1}),\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill {\theta }_{n+1}& ={\theta }_n+f\left(I_{n+1}\right)+\epsilon g({\theta }_n,I_{n+1})mod2\pi ,\hfill \end{array}$$

(8)

where $f\left(I\right)$ is the so-called winding number, and h and g represent the effects of the perturbation term in the Poincaré map. The twist condition (6) reads

$$\frac{d{\theta }_{n+1}}{d{I}_n}\ne 0.$$

(9)

The Poincaré map preserves the symplectic area in the surface of section if

$$\frac{\partial g({\theta }_n,I_{n+1})}{\partial {\theta }_n}+\frac{\partial h({\theta }_n,I_{n+1})}{\partial I_{n+1}}=0.$$

(10)

A further simplification consists in choosing $g({\theta }_n,I_{n+1})=0$ and $h({\theta }_n,I_{n+1})=sin{\theta }_n$, a choice that fulfills the symplectic condition (10). In this case, the twist condition (9) reduces to $df/dI\ne 0$, i.e. the winding number profile should be monotonic over the range of the action variable.

An example of area-preserving non-twist map, where $f\left(I\right)=k(I^2-1)$, was introduced by Weiss in the context of advection of passive scalars (see Appendix A)

$$\begin{array}{cc}\hfill I_{n+1}& =I_n-ksin\left({\theta }_n\right)\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill {\theta }_{n+1}& ={\theta }_n+k\left(I_{n+1}^2-1\right)mod2\pi ,\hfill \end{array}$$

(12)

where k is a parameter representing the non-integrable perturbation. Since $df/dI=2kI$ the twist condition is not satisfied at $I=0$. Indeed, non-twist maps usually have non-monotonic winding number profiles. A shorthand notation for this map is $\mathbf{x}\mapsto \mathbf{M}\left(\mathbf{x}\right)$, where $\mathbf{x}=(I,\theta )$ and $\mathbf{M}$ is given by (11)-(12). Since this is a Hamiltonian system, there exists an inverse map ${\mathbf{M}}^{-1}$.

In the $k\to 0$ limit we have an identity map $(I\mapsto I,\theta \mapsto \theta$. For relatively small values of k the system becomes non-integrable and one can observe quasiperiodic orbits spanning the entire interval $[-\pi ,\pi ]$ and invariant curves inside islands centered at a periodic orbit of the map (11)-(12). This is a consequence of the non-monotonicity of the winding number profile, i.e. there will be two orbits with the same winding number (a phenomenon also known as degeneracy, and which appears only for non-twist systems). For a larger value of the perturbation parameter k we observe such a collision of periodic orbits, involving a reconnection of the islands’ separatrices. This is actually a global bifurcation changing the topology of the orbits as some parameter is varied through a critical value.

In the Figure 1(a)-(d) we show phase portraits generated using the Weiss map (11)-(12) for different values of the parameter k. There can be observed chaotic orbits near the former islands’ separatrices, that no longer exist due to the homoclinic tangle formed after reconnection. A distinguished feature of non-twist maps is the presence of a robust shearless barrier between the two islands [depicted in blue in Figure 1(a)]. Such a shearless barrier corresponds to a local extremum of the winding number profile for the map (11)-(12). This barrier prevents global transport related to the chaotic orbits, the shearless barrier clearly separates the chaotic regions near the separatrices.

An increase of k, however, causes the shearless barrier breakup and the mixing of the chaotic regions associated with each island [Figure 1(b)]. The latter, on its turn, occupies a wider fraction of the phase space as k is further increased [Figure 1(c)-(d)]. The shearless barrier breakup occurs for a critical value of k between $0.50$ and $0.55$ but a precise determination is difficult to make only from these phase portraits. In the following we will consider a systematic way to accomplish this task.

In order to distinguish between chaotic and non-chaotic orbits we have computed the Lyapunov exponents for this map

$${\lambda }_{1,2}=\underset{n\to \infty }{lim}\frac{1}{n}log\left|\right|{\mathbf{DM}}^n\left(\mathbf{x}\right)·{\mathbf{u}}_{1,2}\left|\right|,$$

(13)

where $\mathbf{DM}$ is the tangent map corresponding to Eqs. (11)-(12) and ${\mathbf{u}}_{1,2}$ are mutually orthogonal eigendirections. Due to the area-preserving nature of the Weiss’s map, it follows that ${\lambda }_1+{\lambda }_2=0$, such that it suffices to present results for the largest Lyapunov exponent ${\lambda }_1$. A color map for the latter is showed in Figure 2 for different values of the parameter k, and corresponding to the same values used in the phase portraits of Figure 1.

The islands’ interior, comprising quasi-periodic closed orbits, is related to vanishing Lyapunov exponents, whereas the chaotic region near the islands’ separatrices have positive values of ${\lambda }_1$ [Figure 2(a)]. Moreover, the existence of a shearless barrier is clearly separates the local islands’ separatrices. The same pattern is observed for higher k, but the chaotic orbits have a considerably higher value of ${\lambda }_1$, an almost tenfold increase [Figure 2(b)]. By the same way, the shearless barrier breakup can be observed by the Lyapunov exponent colormap [Figure 2(c)]. The value of ${\lambda }_1$ also increase for higher k [Figure 2(d)].

3. Escape basins

The Hamiltonian system given by the map equations (11)-(12) is opened by considering that the particles can escape by one or more exits in the ($x,y$) phase space [24,25]. The sets of initial conditions, that reach each one of the exists, after a given number of map iterations, form their corresponding basins of escape. If the exits are placed at regions with non-chaotic orbits their basins are relatively simple. On the other hand, if the exits are placed in chaotic orbits their corresponding escape basins are fractal, with a fractal basin boundary. This results from the properties of the chaotic saddle, an invariant non-attracting chaotic set formed by the intersections of the unstable and stable manifolds of unstable periodic orbits embedded in the chaotic region [31].

We will consider two possible square exits of length $0.2$, centred in ($0.0,-1.1$) and ($\pi -0.1,1.0$), denoted by ${\mathcal{B}}_1$ and ${\mathcal{B}}_2$, respectively. Given a initial condition $(x_0,y_0)$. For each iteration of the map (11)-(12), we make the following test: if $x_n,y_n$ inside one of the squares we stop iterations and save the values of the initial condition.

In the Figure 3 we showed the escape basins for different values of the parameter k, the purple pixels correspond to the initial conditions that escape through the square located in ${\mathcal{B}}_1$ and the orange one escapes thorough the exit ${\mathcal{B}}_2$, and the white ones correspond to initial conditions that do not escape in our computation time ${10}^5$, these points are trapped inside islands. For the case $k=0.5$ Figure 3(a) the basins are separated by invariant curves, in Figure 3(b) the basins are mixed, however the initial conditions that are close to the central island tends to belong to ${\mathcal{B}}_1$, while the initial conditions close to the the points $x=-\pi$ and $x=\pi$ tends to the the ${\mathcal{B}}_2$ exit. The intermixing of the basins increases with the increment of k, as showed in the Figure 3(c) and (d). The Figure 4 shows the fractal structure form by the basins in a fine scale.

The escape basins are mixed at arbitrarily fine scales, as also the escape time, i.e. the number of map iterations that an orbit takes to reach one of the openings, has a complicated distribution in the phase space. The time that a initial condition takes to leave the system is showed in Figure 5 (in a color bar), as a function of the initial condition. Reddish colors correspond to higher escape times this occurs around the islands and the invariant curves Figure 5(a), bluish colors correspond to a small escape time, and white pixels are orbits that do not escape. It is clear the formation of escape channels, paths to each the initial conditions escapes for very small times.

4. Basin entropy

In order to, quantify the final state uncertainty produced by the fractality, we apply the concept of basin entropy, developed by Daza and coworkers [19,20]. It was originally developed for basins of attraction and their boundaries, but it can be extended for a more general setting, which is basins of behavior. For open systems, for which the desired behavior is the escape of orbits, we can work with the corresponding escape basins and their boundaries. We have applied this methods in a variety of problems involving magnetic field lines in Tokamaks [28,29], drift motion of charged plasma particles [27], and light scattering through black holes [30].

Let us consider a bounded phase space region $\mathcal{A}$ which includes a part of the escape basin boundary, and cover this region into boxes by using a mesh of $M\times M$ points. We assign to each mesh point a random variable, whose values characterize each different escape. The basin entropy is obtained by applying the information entropy definition to this set. For open systems, we consider a number $N_A$ of exits through which orbits can escape. We assign to each mesh point (which stands for an initial condition), an integer (pseudo-)random variable (called color) labeled from 1 to $N_A$.

The region $\mathcal{A}$ is covered with a regular grid of N boxes with sidelength $\epsilon =n/M$, where $n\in \mathbb{N}$. Let $p_{i,j}$ denote the probability that the jth color is assigned to the ith box, where $i=1,2,\dots$. The fraction of points $p_{ij}$ belonging to a basin inside a box i is computed for each box, considering that the colors inside a box are equiprobable, i.e. there is statistical independence. The information (Gibbs) entropy of the ith box is

$$S_i=-\sum _{j=1}^{m_i}{p}_{ij}log{p}_{ij}.$$

(14)

where $m_i\in [1,N_A]$ is the number of colors for the ith box. The total entropy for the mesh covering the region $\mathcal{A}$ results from summing over the N boxes, or $S={\sum }_{i=1}^N{S}_i$. The basin entropy is defined as the total entropy divided by the number of boxes

$$S_b=\frac{S}{N}=-\frac{1}{N}\sum _{i=1}^N\sum _{j=1}^{N_A}{p}_{ij}log{p}_{ij}.$$

(15)

The system considered in the present work has two exits, named 1 and 2, with the corresponding escape basins, as described in the previous Section. The corresponding probabilities $p_{i,j}$ satisfy $p_{i,1}+p_{i,2}=1$ for each i, such that the basin entropy reads

$$\begin{array}{cc}\hfill S_b& =-\frac{1}{N}\sum _{i=1}^N\left(p_{i,1}log\left(p_{i,1}\right)+p_{i,2}log\left(p_{i,2}\right)\right),\hfill \\ \hfill & =-\frac{1}{N}\sum _{i=1}^N\left\{p_{i,1}\left[log\left(p_{i,1}\right)-log(1-p_{i,1})\right]+log(1-p_{i,1})\right\}=-\frac{1}{N}\sum _{i=1}^N{S}_i.\hfill \end{array}$$

(16)

From the computational point of view, the escape basin are discretized into pixels with equal size, such that each square box contains $N_p^2$ pixels, where $N_p$ is the number of pixels contained by the box with sidelength $ϵ$. For any given box i the corresponding probability $p_{i,1}$ takes on a discrete value out of the following set

$$p_{i,1}\in \left[0,\frac{1}{N_p^2},\frac{2}{N_p^2},\cdots ,1-\frac{1}{N_p^2},1\right].$$

Notice that those boxes for which $p_{i,1}=0$ or 1 do not contribute to the computation of the basin entropy $S_b$ because $S_i=0$ for such cases. What remains is the contribution of the boxes at the escape basin boundary, namely those containing pixels of both escape basins. Hence the possible values of the probabilities $p_{i,1}$ for the remaining $N_b$ boxes are given by

$$p_m=\frac{m}{N_p^2},(m=1,2,\dots N_p^2-1).$$

(17)

Considering that a given fraction $q_m$ of the $N_b$ boxes has a basin probability given by (17), the basin entropy (16) becomes

$$\begin{array}{cc}\hfill S_b& =-\frac{1}{N}\sum _{m=1}^{N_p^2-1}{q}_m{N}_b{S}_m\hfill \\ \hfill & =-\frac{N_b}{N}\sum _{m=1}^{N_p^2-1}{q}_m{S}_m=-C\frac{N_b}{N},\hfill \end{array}$$

(18)

where C is a constant that depends on the distribution of the quantity $q_m$. For fractal basins boundaries, which is just the case of the escape basins here investigated, the values of $q_m$ are concentrated around a mean value with a Gaussian-like distribution.

Let us denote by d and D the box-counting dimensions of the escape basin and its corresponding basin boundary, respectively. Considering that it takes a number N of boxes with sidelength $ϵ$ in the phase space region $\mathcal{A}$, it follows that $N\sim \tilde{N}{ϵ}^{-d}$ for small enough $ϵ$, where $\tilde{N}$ is a constant. By the same token, since it takes a number $N_b$ of those boxes to cover the corresponding basin boundary, then $N_b\sim {\tilde{N}}_b{ϵ}^{-D}$, where ${\tilde{N}}_b$ is another constant and $ϵ$ is also small enough. Substituting both expressions into Eq. (18) we obtain a relation between the basin entropy and the box-counting dimensions of the basin and its corresponding boundary:

$$S_b=-C\frac{{\tilde{N}}_b}{\tilde{N}}{ϵ}^{d-D}.$$

(19)

This equation can be further transformed by using the concept of uncertainty exponent. Since each initial condition is determined, in the two-dimensional phase space, up to a given uncertainty $ϵ$, it can be represented by a disk of radius $ϵ$ centered at that initial condition. If this disk does not intercept the basin boundary, all its interior points converge to the same escape and the initial condition is $ϵ$-certain. Otherwise, if the disk intercepts the basin boundary, it is called $ϵ$-uncertain. The fraction of $ϵ$-uncertain disks is known to scale with the uncertainty $\delta$ as a power-law: $f\left(ϵ\right)\sim {ϵ}^{\alpha }$, where $\alpha =d-D$ is called uncertainty exponent. Substituting into (19) there results

$$ln{S}_b\left(ϵ\right)=\alpha lnϵ+ln\left(\frac{{\tilde{N}}_b}{\tilde{N}}C\right).$$

(20)

We used the linear relation Eq. 20 to estimate the uncertainty exponent $\alpha$ for the escape basins of the Weiss map considered in Section II. For each value of the box sidelength $ϵ$ we computed the basin entropy using (16), and we repeated the procedure for a number of values of $ϵ$ with $M=1000$ and n in the interval $[15,35]$, the results being depicted in Figure 6. We have used a least squares fit to obtain a value of $\alpha =0.0066\pm 0.0003$ in (a) and $\alpha =0.0054\pm 0.0002$ in (b), where the numerical error arises from the fitting.

The fractality is quantify with the aid of the basin and boundary basin entropy’s, we consider a grid of boxes, so that, each box contains 25 initial conditions, from Daza et al. [19] this values produces the optimum results of the basin entropy. In Figure 7(a), we showed the entropy as a function of the parameter k. The basin entropy is zero until k reach a critical values $k_c$ where the shearless curve is broken leading to the mix of the two basins. We estimate the values of $k_c$ as the values of k that produced $S_b\ne 0$, meaning that the basins are mixed. This values is $k_c=0.535$, the entropy sets close to the maximum value $ln2$ meaning that there is a great uncertainty of the final state, caused by the fractal structures of the system. In Figure 7(b) is presented a magnification of the basins for $k_c$, in the region were that was the invariant curves.

5. Conclusions

In this work, we investigated the escape of chaotic orbits in a non twist map called Weiss map. Considering opening in the phase space, we can calculate the escape basins, they present a fractal structures given by the underlying dynamical structure of the chaotic orbits. The escape time also showed fractal structures with presence of paths where the escape is very fast.

In order to, quantify the fractality, we used the concept of basin entropy, a quantity of the uncertainty of the final state caused by the fractal structures. Moreover, we showed a way to compute the uncertainty exponent with the basin entropy, for $k=0.6$, $alpha=0.0066$ what indicates that the basin boundary is extremely involved. The system exhibits a collection of invariant curves that act as boundaries between chaotic regions surrounding two main islands. However, as k increases, these curves become broken, until only one curve remains the shearless curve, this is broken when $k=k_c$, so that the two chaotic regions are connected. The value of $k_c$ were estimated using the entropy basin concept, to which we find the value of $k_c=0.53$.

The basin entropy concept we used in this work is based on the Boltzmann-Gibbs-Shannon-von Neumann version of the entropy, which is an extensive quantity by definition. However, in dynamical systems exhibiting complex behavior, including coexistence among periodic, quasiperiodic, and chaotic orbits like the Weiss map we considered in the present paper, there has been argued that a non-extensive entropy would be better suited, like Tsallis entropy [32,33]. This is particularly interesting when the chaotic transport is characterized by anomalous diffusion, which is heavily influenced by stickiness and other dynamical features of the chaotic orbit. A further extension of the present approach would be therefore to adapt the basin entropy concept vis-a-vis of the Tsallis non-extensive entropy.