Homoclinic puzzles and chaos in a nonlinear laser model

Of special interest in this paper is the role(s) of cod-2 inclination flip bifurcations (labeled as $IF$) in the parameter space and how they shape and organize the global bifurcation unfolding of the OPL-model, see Figs. 4,4,15. The IF bifurcation occurs when an orientable stable or unstable manifold transitions to a non-orientable one, when followed along the homoclinic orbit. Unlike the Shimizu-Morioka model [18], this bifurcation in the OPL model involves the primary homoclinic butterfly and not the double homoclinic loops. As mentioned previously, L.P. Shilnikov pointed out the conditions under which the IF bifurcation gives rise to the onset of the Lorenz attractor in $\mathbf{Z_{2}}$-symmetric systems. These conditions helped determine and verify the regions of existence of the Lorenz attractor, as well as enabled computer assisted proofs of chaotic dynamics in the canonical Lorenz model, without stable orbits and homoclinic tangencies [57, 58, 59, 60]. One such condition showed that the inclination-switch bifurcation can also lead to the onset of stable orbits nearby in the phase space, even in the case of an expanding saddle whose saddle index satisfies the condition $1/2<\nu<1$. Below, we will highlight the essence of the inclination-switch bifurcation. Its in-depth analysis is given in [4, 5].

3D images in Figure 6 illustrate the concept of an inclination flip bifurcation that transforms an orientable homoclinic loop of the saddle (Figure 6a) into a non-orientable one (Figure 6b). This is depicted by a single twist of the surface traced out by a small, local section of the strongly stable direction, due to the smaller eigenvalue $\lambda_{3}$ ($\dots<\lambda_{3}<\lambda_{2}<0<\lambda_{1}$) of the saddle. Let us explore the global return map $\mathbf{T}$ that takes a 2D cross-section, $\pi=\pi_{1}\cup\pi_{2}$, transverse to the stable manifold $W^{s}$ of the saddle with a saddle index $\nu<1$, onto itself along the homoclinic loop. The 1D outgoing separatrix of the saddle ( here $\Gamma_{1}$) comes back along the leading direction (vertical, due to $\lambda_{2}$). Before the inclination-flip (Figure 6a), the section $\pi_{1}$ (the right half, relative to $W^{s}$) is mapped back on to the original $\pi_{1}$ (darker wedge), along the orientable primary homoclinic loop $\Gamma_{1}$ (in black). A nearby trajectory (red) is also shown. Post the inclination flip (Figure 6b), in the non-orientable or twisted case, the section $\pi_{1}$ is mapped onto the opposite section $\pi_{2}$ (dark yellow). This makes the surface “spanned” by the unstable (due to $\lambda_{1}$) and the non-leading stable (due to $\lambda_{3}$) eigenvectors look like a Möbius band.

Figure 6c depicts the action near an IF bifurcation using a 2D local Poincaré map. The map takes a small interval $d_{1}\ll 1$ on $\pi$ into $d_{2}\sim d_{1}^{\nu}>d_{1}$. For $\nu<1$, the local map near the saddle is typically an expansion, stretching phase areas or volumes. Let us consider a small portion of the local manifold $\mathfrak{M}$ tangent to the span of the leading stable (due to $\lambda_{2}$) and the leading unstable eigenvectors of the saddle $O$. Near an IF bifurcation, as we take $\mathfrak{M}$ along the outgoing separatrix loop in a succession, it bends and hits the cross-section $\pi$ as a hook, narrowly squeezed due to the strongly stable eigenvalue $\lambda_{3}$. This hook $\mathbf{T}\pi_{1}$ of the pre-image $\pi_{1}$ looks like a Smale horseshoe on the cross-section (see Fig. 6c inset). This bending makes the image of $d_{2}$ shorter than the original $d_{1}$, making the global map $\mathbf{T}$ a contraction ($\mathbf{T}d_{1}<d_{1}$) as it overcomes the local expansion near the saddle ($d_{2}>d_{1}$). Note that the global map $\mathbf{T}$ is an expansion prior to the IF bifurcation. Figure 6d illustrates the onset of one of two “one-sided” chaotic attractors in the OPL-model near the inclination flip point on the primary homoclinic loop (at $a=3.7,b=0.48415447,\sigma=1.5$, next to the $IF_{1}$-point in Fig. 4) due to the presence of such a hook/Smale horseshoe (red), filled in with the points of intersection of transverse segments of the outgoing separatrix $\Gamma_{1}$ of the saddle $O$ on a nearby cross-section. The inset depicts a simulated analogue of the bending manifold $\mathfrak{M}$ sketched in Fig. 6c.

The geometric model of the Lorenz attractor (presented in [57]) can be elaborated using a 2D return map near the primary homoclinic butterfly of the two separatrix loops of the saddle. For a system to possess a genuine chaotic attractor devoid of homoclinic tangencies or stable orbits, this map should satisfy a few analytical conditions. The boundary of its existence region is marked by a violation of these criteria. The 2D map near the IF bifurcations can be further reduced to a simplified 1D map (Fig. 7) in the following form [5]:

where $1/2<\nu={|\lambda_{2}|}/{\lambda_{1}}<1$ is the saddle index, $A$ is the separatrix value, and $\mu$ is the splitting parameter of the homoclinic loop [12]. Near the IF bifurcations when $0\sim|A|\ll 1$, the term $o(|x_{n}|^{2\nu})$ is no longer negligible. Fig. 7 demonstrates the bent shape of such a return map, derived from the simulated 2D map, shown in Fig. 6d, for the “one-sided” chaotic attractor. As $\mu$ is decreased further, both “one-sided” (left/right) non-oriented attractors merge into a symmetric one. The map illustrates possible dynamics near the primary IF bifurcation due to period-doubling bifurcations in the 1D Poincaré map, leading to the emergence of multiple secondary homoclinic curves from this cod-2 point, as well as multiple stable periodic orbits (Figs. 8,9,15,16) when the return map starts bending again and again to resemble the return map of a Shilnikov saddle-focus but with a finite number of Smale horse-shoes [18].

A typical signature for Lorenz-like systems is the complex universal and self-similar characteristic spirals in the parameter plane, organized around a central point called a Bykov terminal point (T-point) as seen in Figs. 4,4 and detailed magnifications are presented in Fig. 8 [30]. Each characteristic spiral around a T-point in the parameter plane corresponds to a homoclinic loop of the saddle $O$ in the phase space, such that with each turn of the spiral approaching the T-point, the outgoing separatrix of saddle makes increasing number of loops around a saddle-focus $P^{\pm}$, before finally forming a closed heteroclinic connection between the saddle and the saddle focus at the T-point. Fig. 2d shows a one-way heteroclinic connection between the saddle and the saddle focus for parameter values near the T-point $T_{1}$ in Fig. 8b, with $(a,b)\sim(3.68179,0.3517)$. Here, the unstable separatrix $\Gamma_{1}$ (red) of the saddle $O$ makes one loop around $P^{+}$, followed by another loop around $P^{-}$, then merges with the incoming separatrix of $P^{+}$, thus effectively making infinite loops around $P^{+}$ before emerging out. The symmetric heteroclinic connection along $\Gamma_{2}$ from $O$ to $P^{-}$, is shown in blue.

Near a T-point, there exist countably many secondary T-points with increasing complexity of heteroclinic connections in the phase space and with similar bifurcation structures as the central T-point in the parameter plane, although on a smaller scale (see Figs. 8, 9). Multiple inclination flip (IF) homoclinic bifurcations of the saddle occurring along the characteristic spirals of the T-points (detected with MatCont and shown as white dots in Figs. 8a,9) give rise to saddle-node and period-doubling bifurcations of periodic orbits [48, 5]. Fine organization of the structure of the chaotic regions and stability windows near the T-point and surrounding IF points is revealed in greater detail in 9b. In addition, the unfolding of a T-point also includes other curves corresponding to the homoclinic connections of the saddle-focus satisfying the Shilnikov condition [24, 54, 56] that give rise to a denumerable set of saddle periodic orbits nearby [25], as well as those corresponding to heteroclinic connections between both saddle-foci [19, 46, 61].

Topologically identical regions in the parametric plane can be identified by defining a formal convergent power series P [29, 62] for a sequence $\{k_{i}{\}}_{i=p}^{q}$ defined as:

For example, the value of this power series for the sequence $\{10100101...\}$ with $p=0$ and $q=7$ can be computed as : $P=\dfrac{1}{2^{8}}+\dfrac{0}{2^{7}}+\dfrac{1}{2^{6}}+\dfrac{0}{2^{5}}+\dfrac{0}{2^{4}}+\dfrac{1}{2^{3}}+\dfrac{0}{2^{2}}+\dfrac{1}{2}=0.64453125.$ The value of the sum $P$ lies between 0 (for the sequence $\{\overline{0}\}$) and 1 (for the sequence $\{\overline{1}\}$ in the limit as $(q-p)\rightarrow\infty$). Two different sets of parameter values that show topologically identical behavior in their trajectories, result in identical sequences $\{k_{i}\}$, and thus, have the same $P$ values. Therefore, $P$ serves as a dynamic invariant, which we can use in the construction of biparametric scans such as Fig. 4, as follows. The parameter $\sigma$ is kept constant while we vary the parameters $a$ and $b$ (at a resolution of 2000x2000). For each set of values of these parameters, we follow the right unstable separatrix $\Gamma_{1}$ of $O$ to obtain the sequence $\{k_{i}\}$ and the power series $P$. We then use a colormap to project the $P$ values onto the biparametric plane, where topological similarity between regions is identified by their equivalent colors, and the boundaries between adjacent regions represent homoclinic bifurcation curves. The colormap is constructed by discretizing the range of $P$, i.e., [0,1], into $2^{8}$ bins, and assigning them RGB-color values from 0 through 1, for each of Red, Green and Blue colors in decreasing, random, and increasing fashion, respectively. Greater weight is assigned to the symbols towards the end of the sequences $\{k_{i}\}$ in the computation of $P$, to maximize the color contrast between neighboring regions in a scan, that differ in the last symbol, separated by a homoclinic curve at their boundary. Numerical integration is performed using the classic Runge-Kutta method (RK4) with fixed step size $dt=0.01$. The computation of these trajectories is massively parallelized by running on separate GPU threads using CUDA. This allows for the construction of bi-parametric scans such as Fig. 12, in as little as 8 seconds on an Nvidia Tesla K40 GPU. Visualizations are done in Python. In order to construct complex bifurcation structures in the 3D parametric space, such as Figs. 17a,19,21, we obtain a large number of biparametric scans in the $(a,b)$ parametric plane, as we continuously vary the third parameter $\sigma$. Such arrays of scans are then rendered together using the open source volume exploration tool Drishti [69], which performed the best with our huge datasets, compared to other available tools for 3D rendering.

The significance of our DCP tool is that not only can it reveal the short term transient dynamics and the underlying homoclinic, heteroclinic, saddle, and T-point spiral structures, but it can also be employed to reveal the long term behavior of the system (Fig. 14). When looking at the long term behavior of a trajectory after omitting a long transient, one should be aware of the shift-symmetry of periodic orbits. Depending on the length of the initial transient omitted, the same symmetric 8-shaped periodic orbit (Fig. 2b - blue) can be represented by either $\{\overline{01}\}$) or $\{\overline{10}\}$), and hence, needs to be normalized. In addition, if there is an asymmetric periodic orbit such as $\{\overline{011}\}$), the symmetry of the system implies the existence of another periodic orbit which is its symmetric counterpart, $\{\overline{100}\}$). This implies that all six of the periodic sequences $\{\overline{011}\}$), $\{\overline{110}\}$), $\{\overline{101}\}$), $\{\overline{100}\}$), $\{\overline{001}\}$) and $\{\overline{010}\}$) need to be normalized. We accomplish this by means of an algorithm we call Periodicity Correction (PC), where we identify the periodic structure within a sequence, along with its symmetric counterpart, and then normalize the sequence to the smallest valued circular permutation. Thus, both $\{\overline{01}\}$) and $\{\overline{10}\}$) are normalized to $\{\overline{01}\}$) and all six of the asymmetric periodic orbits described previously are normalized to the lowest valued sequence $\{\overline{001}\}$). PC helps in identifying regions of simple, Morse-Smale dynamics of stable fixed points and periodic orbits by eliminating artifacts arising from symmetries of the system or the symbolic sequence. Structurally unstable, chaotic regions are marked by the lack of periodicity in their symbolic representations. Alternatively, we can also employ a compression algorithm, like Lempel-Ziv-76 (LZ) [37], to find the complexity of a binary sequence. LZ scans a sequence from left to right, and adds a new word to the vocabulary every time a previously unencountered sub-string is detected. Since the circular permutations of a periodic orbit will have identical complexity, this approach can also detect windows of stability amidst structurally unstable chaotic regions.

While PC can detect stable periodic orbits efficiently even at short sequence lengths, LZ requires very long sequences and is not as effective. Within chaotic regions, on the other hand, PC cannot distinguish between different sequences while LZ separates very complex strings from not-so-complex strings, akin to Lyapunov exponents. To harness the best of both worlds, we therefore combine both PC and LZ. For a long term biparametric sweep, we first run the symbolic sequences through PC to detect the existence of, and to normalize, periodic orbits. The normalized sequences are then used to compute the power series sum $P$ and colored with the colormap as previously described. We then run the aperiodic sequences through LZ to detect their complexity. We color the LZ-complexity values in shades of gray, with darker gray representing greater complexity. Together, this results in long term bi-parametric sweeps such as Figs. 9,14,15b,16b,, where darker shades of gray represent greater structural instability and chaotic dynamics, while all other solid colors represent simple Morse-Smale dynamics of stable fixed points or periodic orbits. On the whole, we coin this collection of symbolic tools the ‘Deterministic Chaos Prospector (DCP)’. The open source software developed is available at https://bitbucket.org/pusuluri_krishna/deterministicchaosprospector/.

We start this section with $\sigma=1.5$. Fig. 4 shows the bifurcation diagram in $(a,b)$-parametric plane at $\sigma=1.5$ using symbolic sequences $\{k_{i}{\}}_{i=5}^{12}$. Two sets of parameters that result in topologically identical dynamics for the symbolic sequence under consideration, are colored identically in the bi-parametric sweep. Boundaries between such regions represent a shift in the dynamics due to homoclinic bifurcations. The corresponding bifurcation diagram obtained using continuation is shown in Fig. 4. Blue curves represent homoclinic bifurcation curves while heteroclinic bifurcation curves are shown in red.

A homoclinic bifurcation curve of the trivial equilibrium $O$ emerges from $BT$, which will be of primary interest as we vary $\sigma$. As we follow this primary curve we encounter a codim 2 point $DRS$ with two real stable eigenvalues and the leading stable direction changes. Here, a double, twisted homoclinic loop emerges. As we move along the primary homoclinic bifurcation curve $H_{0}$ from $DRS$ towards $IF_{1}$ in the parametric plane (Fig. 4), in the phase space (Fig. 5a) the primary homoclinic loop is oriented and approaches $O$ from the right (light red curve). Moving on along $H_{0}$, we find a first inclination flip bifurcation $IF_{1}$ of complex type, where the primary homoclinic orbit approaches $O$ along its stable manifold (Fig. 5a - dark red). As we move further along $H_{0}$, a twisted single loop approaches $O$ from the left (green), and $\nu_{1}$ and $\nu_{2}$ increase again when we find a second inclination flip point $IF_{2}$. Beyond $IF_{2}$, the homoclinic orbit changes its shape from a single loop around one non-trivial equilibrium to a double loop configuration around both non-trivial equilibria, that gradually grows bigger (light blue $\rightarrow$ dark blue). Finally the primary curve gets close to the DRS point where the continuation terminates. This corroborates theoretical expectations although verification that the exact end point of this primary homoclinic curve is at the DRS point, is still open.

From the eigenvalues (see Table 2), we find that the inclination flip at $IF_{1}$ and other cases are all of complex type so that infinitely many homoclinic bifurcation curves emerge from this point $IF_{1}$. The homoclinic bifurcation curves all spiral to the T-points such as $T_{0}$ (marked by the heteroclinic connection $\{1\overline{0}\}$). From the T-point $T_{0}$, two heteroclinic bifurcation curves also emerge connecting the saddle-foci $P^{\pm}$. These curves are rather straight and it may be appreciated that the homoclinic spirals to secondary T-points end up close to these heteroclinic curves. The inclination flip points, T-points, saddles($S$) and closed loops are further explored in Sections 4.3,4.4.

The $(a,b)$-parametric sweep for $\sigma=2.0$ using symbolic sequences $\{k_{i}{\}}_{i=3}^{10}$ and continuation is shown in Fig. 10. We find results are globally quite similar for the pitchfork and AH curves. The primary homoclinic bifurcation curve $H_{0}$ has some marked differences though. First, the codimension-two neutral saddle point has disappeared simultaneously altering the nature of the BT point. Second, the primary homoclinic bifurcation curve no longer ends at the DRS point but spirals towards the T-point. In this case, as we sample parameters beyond $IF_{2}$, in the phase space the twisted loop starts making more and more revolutions around the stable manifold of $P^{-}$, until it forms a heteroclinic connection to $P^{-}$ at the T-point. As shown in Fig. 5b, this transfiguration of the primary homoclinic orbit into a heteroclinic connection to $P^{-}$ happens through the orbits light red $\rightarrow$ dark red $\rightarrow$ green $\rightarrow$ light blue as we move along $H_{0}$. Detailed transition of the behavior of $H_{0}$ between $\sigma=1.5$ and $\sigma=2.0$ is explored further in Section.4.4. We also note that the second leading stable eigenvalue has changed from complex at $\sigma=1.5$ to real for $\sigma=2.0$. This transition occurs at $\sigma\approx 1.90$.

Moving on to $\sigma=10.0$, the bifurcation diagram using symbolic sequences $\{k_{i}{\}}_{i=2}^{8}$ is shown in Fig. 12. The corresponding diagram using parametric continuation is shown in Fig. 12. The inclination flip $IF_{2}$ on the primary homoclinic is of a different type (type B), see Table 2 for the eigenvalues. This involves small regions with stable limit cycles. Indeed, we find a limit point of cycles bifurcation (LPC) indicating a small region with stable periodic orbits. These LPCs are only connected to the primary homoclinic bifurcation. This may not be apparent from Figure 12 as the primary curve does not seem to make a loop towards the T-point. Detailed calculations though suggest the primary curve makes a turn near $\nu_{1}=0.3$, but accurate continuation did not succeed in that parameter region. All other inclination flip points are of type C. This implies that only the primary homoclinic bifurcation results in the emergence of additional stable periodic orbits nearby. The label $NS$ stands for the neutral homoclinic saddle with a zero saddle value. Here, $T_{0}$ ($\{1\overline{0}\}$), $T_{1}$ ($\{10\overline{1}\}$), $T_{2}^{1}$ ($\{11\overline{0}\}$) and $T_{2}^{2}$ ($\{11\overline{0}\}$) represent various T-points. Note that the T-points $T_{2}^{1}$ and $T_{2}^{2}$, above and below the saddle $S$, have identical construction and the same heteroclinic connection $\{11\overline{0}\}$. $T_{g}$ with the heteroclinic connection ($\{\overline{1}\}$) is a special case. Although $T_{g}$ resembles a T-point, it falls in the region of the periodic orbit ($\{\overline{01}\}$) which devours the spiral structure of $T_{g}$, leaving behind a ghost of the T-point. $T_{g}$ and the semi-annular isolas near $C$ are detailed in Fig. 23. Further details of the special organizing structures are provided in the following sections.

We conclude this section by presenting a summary of the bifurcation diagrams for varying $\sigma$-values. This also serves as a precursor for the construction of 3D parametric sweeps of Sec.4.4. Fig. 13(a-h) reveal variations in the bifurcation structure using $\{k_{i}{\}}_{i=4}^{11}$ for $\sigma$-values $1.2$, $1.32$, $1.5$, $1.72$, $1.76$, $2.0$, $6.0$, $7.375$ and $10.0$. See also supplementary Movie M1 that reveals these variations in further detail. Between $\sigma=1.72$ (Fig. 13d) and $\sigma=1.76$ (Fig. 13e), the primary homoclinic bifurcation curve branches and starts spiraling on to the primary T-point $T_{0}$. As we increase $\sigma$ further, the inclination flip $IF_{1}$ and the homoclinic bifurcation curves from $IF_{1}$ to $T_{0}$ change their orientation. Between $\sigma=6$ (Fig. 13g) and $\sigma=7.375$ (Fig. 13h), $H_{0}$ branches again to separate from $T_{0}$.

For all three sigma values of $1.5$ (Fig. 4), $2.0$ (Fig. 10) and $10.0$ (Fig. 12), the primary T-point $T_{0}$ is marked by the symbolic sequence $\{1\overline{0}\}$. The phase trajectory at $T_{0}$ is as shown by the light blue curve of Fig. 5b, where $\Gamma_{1}$ makes one loop around $P^{+}$, followed by a heteroclinic connection to $P^{-}$. The trajectory effectively makes an infinite number of loops around $P^{-}$, before leaving the neighborhood of $P^{-}$. Fig. 8a shows a magnification near $T_{0}$ for $\sigma=1.5$ using $\{k_{i}{\}}_{i=5}^{12}$. It also shows multiple secondary T-points surrounding the primary T-point. Fig. 8b with $\{k_{i}{\}}_{i=3}^{10}$ shows a zoom of a small region with an analogous pair of T-points $T_{1}$ and $T_{2}$, above and below $H_{0}$, respectively, for $\sigma=1.5$ (see Fig. 4 white box). The heteroclinic connection at $T_{1}$ is shown in Fig. 2d. Here $\Gamma_{1}$ makes one loop around $P^{+}$, followed by one loop around $P^{-}$, and then forms a heteroclinic connection to $P^{+}$. Hence, the symbolic representation of $\{10\overline{1}\}$. Similarly, the heteroclinic connection at $T_{2}$ is given by $\{11\overline{0}\}$. Such T-points always come in pairs on either side of a homoclinic bifurcation curve. Such a T-point pair can also be seen for $\sigma=10.0$ (Fig. 12), with $T_{1}$ ($\{10\overline{1}\}$) above $H_{0}$, and $T_{2}^{1}$ ($\{11\overline{0}\}$ )below it. The T-point $T_{2}^{2}$ below the saddle $S$, also has the same structure in the parametric plane and the same heteroclinic connection ($\{11\overline{0}\}$) as $T_{2}^{1}$. This is explored further in Sec.4.4.

Next, we show how the inclination flip points serve not only as the confluence of infinitely many homoclinic bifurcation curves, but also as the originators of multiple windows of stability with respect to long term dynamics in the system. Fig. 8a ($\sigma=1.5$) shows several inclination flips (white dots), starting from $IF_{2}$ and leading all the way beyond $T_{0}$, from which the homoclinic bifurcation curves arise leading up to the secondary T-points surrounding $T_{0}$. The specific values of IF bifurcation points are obtained with continuation. Fig. 15 shows the behavior near the inclination flip point $IF_{1}$ at $\sigma=1.5$ for both transient (compare Fig. 4) and long term dynamics (compare Fig. 14b). Fig. 15a ($\{k_{i}{\}}_{i=5}^{12}$) reveals multiple homoclinic bifurcation curves that converge to $IF_{1}$, while Fig. 15b ($\{k_{i}{\}}_{i=999}^{1999}$ with DCP) shows the convergence of multiple stability windows (solid colors) towards $IF_{1}$, interspersed within regions of structural instability and chaos (darker gray implies greater LZ-complexity and less compressible binary sequences thus generated from the model equations). Similar transient and long term behavior near $IF_{2}$ for $\sigma=10.0$ (compare Fig. 12,14d) is summarized in Fig. 16a ($\{k_{i}{\}}_{i=16}^{23}$) and Fig. 16b ($\{k_{i}{\}}_{i=999}^{1999}$ with DCP), respectively.

Finally, we magnify the behavior at $\sigma=1.2$ which nicely summarizes the behavior of the inclination flips and secondary T-points for short term (Fig. 9a, $\{k_{i}{\}}_{i=2}^{9}$) and long term dynamics (Fig. 9b, $\{k_{i}{\}}_{i=1000}^{1999}$ with DCP). Several inclination flips (white dots) close to the primary T-point $T_{0}$ are clearly seen as the congregation centers of homoclinic curves leading up to secondary T-points vis-à-vis transient dynamics, while those same inclination flip points also give rise to stability windows in the long term, amidst regions of structurally unstable and chaotic dynamics. See Fig. 13a,14a for the corresponding global dynamics.

In this section, we describe a saddle in the parameter space that causes branch switching between homoclinic bifurcation curves. As seen in Fig. 4 and Fig. 10, between $\sigma=1.5$ and $\sigma=2.0$, the primary homoclinic bifurcation curve shifts from a complete loop towards the codimension-2 point $DRS$, to a spiral organization around the primary T-point $T_{0}$. This implies the existence of a saddle in between the two $\sigma$ values that branches the primary homoclinic bifurcation curve $H_{0}$, as seen in Fig. 13d,e. Fig. 17(left) shows a detailed 3D $(a,b,\sigma)$ bifurcation diagram close to this saddle. It is constructed using 100 bi-parametric sweeps in the $(a,b)$-plane with $1.7418$ (top surface) $\leq\sigma\leq 1.7439$ (bottom surface), using $\{k_{i}{\}}_{i=0}^{7}$ and 3D volume is rendered with Drishti software [69]. This branching of $H_{0}$ occurs at $\sigma\approx 1.7428$ and is marked by $P$. Detailed zooms close to $P$, shown in Fig. 19, further unravel this 3D organization of the homoclinic bifurcation surface at the homoclinic branching saddle $P$. A similar 3D bifurcation surface with a saddle in the Shimizu-Morioka system is presented in [12] and shown in Fig. 19. Fig. 17(right) shows the homoclinic bifurcation curves close to the branching saddle $P$ obtained with continuation at $\sigma$ values of $1.73$ (top), $1.74$ (middle) and $1.75$ (bottom).

We now focus on another kind of saddle ubiquitous to parametric sweeps, serving as a bridge between two identical T-points, marked $S$ in Figs.4,12. The interesting feature of such a saddle is that the T-point spirals above and below have identical construction and the same heteroclinic connections, as described previously in Fig. 12 for $\sigma=10$, with $T_{2}^{1}$ above and $T_{2}^{2}$ below $S$. In order to see how such T-points on either side of a saddle are related to each other, we run multiple sweeps with closely varying $\sigma$ values and observe the structural changes in the spiral organization. Fig. 20 shows such chaotic mixing close to the saddle $S$ as $\sigma$ values are varied from (a) $\sigma=1.372$, (b) $\sigma=1.352$, (c) $\sigma=1.288$, and through (d) $\sigma=1.264$. As $\sigma$ changes, the T-points above and below the saddle appear to be merging together, as depicted in the transitions: Fig. 20a$\rightarrow$b and Fig. 20c$\rightarrow$d.

The 3D scrolling structure of the bifurcation surface around the hyperbolic saddle $S$ for $1.344\leq\sigma\leq 1.374985$ constructed using 2000 bi-parametric sweeps using $\{k_{i}{\}}_{i=4}^{11}$ and rendered with open-source scientific visualization software Drishti ²²2Drishti is available at https://github.com/nci/drishti is shown in Fig. 21. Fig. 21a,b,c reveal with gradually increasing depth, the continuous structural connections between the T-points on either side of this bridging saddle $S$. As we move along the T-point curve by varying $\sigma$, the curve undergoes a change in orientation and re-enters the bifurcation planes of previous T-points, giving birth to the saddles in between them. Note that each sweep has the 2000x2000 resolution. All 2000 such produced slices imply a total computation of $2000^{3}=8\times 10^{9}$ trajectories. Despite being computationally heavy, this was achieved in just about 8 hours on a single Nvidia Tesla-K40 GPU.

As seen in Fig. 20, as the T-points on either side of a bridging saddle merge together with varying $\sigma$, it results in the formation of isolas of homoclinic bifurcation curves resembling concentric annular structures. This is due to the fact that the T-point curve changes its orientation with respect to the bifurcation planes and the T-point becomes non-transverse to this plane, which was described as a codimension-two-plus-one event [70, 71]. This also becomes clear in Fig. 22 which shows the top and bottom surfaces (rotated through $90^{\circ}$) of the 3D T-point curve in Fig. 21. As we slice several $(a,b)$ bi-parametric sweeps for changing $\sigma$-values, moving towards the bottom surface in Fig. 22b, the T-points appear to be merging as they become non-transverse to the bifurcation plane and we can only observe isolas made of homoclinic bifurcation curves.

We end this section with a description of isolas of a different nature with semi-annular structure as seen near the point $C$ of Fig. 12. The T-point $T_{g}$ (with the heteroclinic connection ) next to $C$ is of particular interest with respect to these isolas. To understand the relation between these isolas and $T_{g}$, we vary $\sigma$ in Fig. 23 through (a) 11.1253 (b) 11.2459 (c) 11.4501 and (d) 11.5522. As $\sigma$ increases, $T_{g}$ starts crossing the boundary between the stable periodic orbit (green region) and structurally unstable chaotic regions in the parametric plane, the semi-annular isolas start merging around $T_{g}$ to gradually reveal the full spiral organization of homoclinic bifurcation curves around the T-point in Fig. 23(d). This explains how in Fig. 12, when the T-point falls inside the region of the periodic orbit, the spiral structure of $T_{g}$ is devoured by the periodic orbit, leaving behind the remnants of the spirals surrounding the ghost of the T-point $T_{g}$, in the form of semi-annular isolas near the boundary. Note that when we say periodic orbit $\{\overline{01}\}$, for the sake of simplicity, we only mean an orbit that produces a periodic sequence in its symbolic representation. Its trajectory may still be chaotic but follow a periodic symbolic sequence as also seen in Fig. 6d for a one sided chaotic trajectory, whose symbolic representation would still be periodic $\{\overline{1}\}$. In fact, the region surrounding $T_{g}$ seems to be the two sided chaotic trajectory following the periodic sequence $\{\overline{01}\}$.