Mapping Chaos: Bifurcation Patterns and Shrimp Structures in the Ikeda Map

Dynamical systems are crucial for understanding the complex behaviors that arise in various natural and engineered processes. Among the most extensively studied dynamical systems are the Lorenz system Lorenz (1963), the Hénon map Gallas (1993), the logistic map May (1976); Strogatz (2018), and the Duffing oscillator Guckenheimer and Holmes (2013); Moon and Holmes (1979). These nonlinear systems, along with many others, display a rich spectrum of behaviors ranging from regular and predictable to chaotic and unpredictable dynamics. The phase space of such systems can generally be divided into three regions: regular Kamphorst and de Carvalho (1999); Oliveira and Leonel (2010); Oliveira and Robnik (2012), chaotic Pustyl’nikov (1995); Sinai (1970); Bunimovich (1979); Robnik (1983); Oliveira, Vollmer, and Leonel (2011); Oliveira and Leonel (2012a), and mixed regions Leonel and McClintock (2005); Oliveira, Bizao, and Leonel (2009); Oliveira, Leonel, and Robnik (2011); Oliveira, Silva, and Leonel (2015); Page et al. (2020).

Regular regions are characterized by periodic or quasi-periodic trajectories that repeat over time, leading to stable and predictable behavior. Chaotic regions, in contrast, exhibit sensitive dependence on initial conditions, where even minute differences in starting points result in vastly divergent trajectories, leading to an unpredictable and complex phase space. Mixed regions are particularly intriguing as they encompass both regular and chaotic dynamics, with stable islands of regularity embedded within chaotic seas, making the overall behavior of the system highly intricate.

In this paper, we explore the properties of a dynamics of the Ikeda map Ikeda (1979); Ikeda, Daido, and Akimoto (1980); Watanabe and Strogatz (1994). Introduced by Kensuke Ikeda in the late 1970s, the Ikeda map models the behavior of light in a nonlinear optical cavity. It has since become a quintessential example in the study of chaos and complex systems. The Ikeda map is defined by a set of iterative equations that describe the evolution of a point in the complex plane. Despite its relatively simple form, the Ikeda map exhibits a rich variety of dynamical behaviors, including fixed points, periodic orbits, and chaotic attractors. The map is particularly notable for its sensitivity to parameter changes, which can lead to sudden transitions between regular and chaotic dynamics as we will show below. Additionally, as we will demonstrate, increasing the dissipation parameter leads to a period-doubling bifurcation cascade Linsay (1981), where Feigenbaum’s constant, $\delta$ Feigenbaum (1978, 1979); Hanias, Avgerinos, and Tombras (2009); Chen et al. (2012), which quantifies the rate of these bifurcations, can be calculated numerically. We then examine the two-dimensional parameter space, specifically focusing on the dissipation parameters associated with the real and imaginary components of the map, and the model reveals the presence of self-similar structures known as “shrimps".

Initial investigations by Gaspard et al. Gaspard, Kapral, and Nicolis (1984) in 1984, Rössler et al. Rössler et al. (1989) in 1989, and Komuro et al. Komuro et al. (1991) in 1991 laid the groundwork for understanding self-similar periodic structures in two-dimensional mappings. Gaspard and colleagues focused on Chua’s system, revealing complex bifurcation patterns, while Rössler’s study on the logistic map and Komuro’s exploration of the Double Scroll circuit further illustrated the ubiquity of these structures in nonlinear dynamical systems.

Gallas’s groundbreaking work in 1993 Gallas (1993) represented a pivotal moment in the study of dynamical systems. Through a detailed exploration of the parameter space of the Hénon map, Gallas identified the presence of complex shrimp-shaped domains—regions characterized by periodic behavior amidst chaotic dynamics. This discovery not only highlighted the significance of these structures but also spurred a wave of subsequent research focused on uncovering similar patterns in other models. As noted by Vitolo, Glendinning, and Gallas (2011), "Shrimps are formed by a regular set of adjacent windows centered around the main pair of intersecting superstable parabolic arcs. A shrimp is a doubly infinite mosaic of stability domains, comprising an innermost main domain and all adjacent stability domains arising from two period-doubling cascades and their corresponding chaos regions. It is important to distinguish these shrimp from their innermost main domain of periodicity." Since Gallas’s seminal work, shrimp-shaped domains have been recognized in a wide array of theoretical models, including but not limited to those examined in studies by Gallas himself Gallas (1995, 1994), as well as in investigations by Hunt et al. Hunt et al. (1999), Bonatto et al. Bonatto, Garreau, and Gallas (2005), and Oliveira et al. Oliveira and Leonel (2013, 2012b, 2014). Furthermore, in addition to theoretical explorations, E. N. Lorenz Lorenz (2008), one of the foremost pioneers of chaos theory, extended the investigation of shrimp-shaped domains into new territories, focusing on their appearance in complex systems. Lorenz’s work highlighted the richness and complexity of these structures, solidifying their importance in the study of chaotic dynamics. These theoretical studies have contributed to a deeper understanding of the conditions under which these domains emerge and their implications for the dynamics of the systems in which they appear. More recently, the existence of shrimp-shaped domains has been experimentally verified in physical systems, most notably in electronic circuit Maranhao et al. (2008); Stoop, Benner, and Uwate (2010); Stoop et al. (2012); Viana et al. (2010). These experimental observations have bridged the gap between theory and practice, demonstrating that these intricate structures are not merely mathematical curiosities but are also observable in real-world systems.

To classify regions in the parameter space as either regular or chaotic, we primarily use the Lyapunov exponent. Although other quantities, such as particle velocity, can also be employed—like in the well-known Fermi-Ulam model Diego and Leonel (2011)—they may not effectively distinguish between regular and chaotic behavior. In this study, we focus on the Lyapunov exponent. Our approach involves starting with a fixed initial condition, allowing for an extended transient period, and then computing the Lyapunov exponent. For each parameter combination of the dissipation paramenters, we assign a color based on the computed Lyapunov exponent. We then increment the parameters, using the final values of the dynamical variables before the increment as the new initial condition, ensuring that we remain within the basin of the same attractor. This methodology reveals well-organized, self-similar shrimp-shaped structures embedded within a broader region of chaotic attractors.

The Ikeda map, originally introduced to model the dynamics of a laser system, is described by the following 1D equation involving a complex variable:

$$z_{n+1}=A+Bz_{n}e^{i\left(\theta-\frac{\phi}{|z_{n}|^{2}+1}\right)},$$

(1)

where the parameters $\theta$ and $\phi$ typically represent phase angles that control the behavior of the system, $z_{n}=x_{n}+iy_{n}$ is a complex number with real part $x_{n}$ and imaginary part $y_{n}$, and $|z_{n}|^{2}=x_{n}^{2}+y_{n}^{2}$ is the squared magnitude of $z_{n}$. In order to obtain the two-dimensional map the describes the dynamics of the system, first, one compute the exponential term in Eq. 1.Thus, by using Euler’s formula, we have:

$$e^{i\left(\theta-\frac{\phi}{|z_{n}|^{2}+1}\right)}=\cos\left(\theta-\frac{\phi}{|z_{n}|^{2}+1}\right)+i\sin\left(\theta-\frac{\phi}{|z_{n}|^{2}+1}\right).$$

(2)

Substitute $z_{n}=x_{n}+iy_{n}$ and the exponential term into the map:

$$z_{n+1}=A+B(x_{n}+iy_{n})\left[\cos\left(\theta-\frac{\phi}{x_{n}^{2}+y_{n}^{2}+1}\right)+i\sin\left(\theta-\frac{\phi}{x_{n}^{2}+y_{n}^{2}+1}\right)\right]e^{C}.$$

(3)

Expanding and separating the real and imaginary parts, we obtain:

$\displaystyle z_{n+1}=A$

$\displaystyle+$

$\displaystyle Be^{C}\left\{\left[x_{n}\cos\left(\theta-\frac{\phi}{x_{n}^{2}+y_{n}^{2}+1}\right)-y_{n}\sin\left(\theta-\frac{\phi}{x_{n}^{2}+y_{n}^{2}+1}\right)\right]+i\left[x_{n}\sin\left(\theta-\frac{\phi}{x_{n}^{2}+y_{n}^{2}+1}\right)+y_{n}\cos\left(\theta-\frac{\phi}{x_{n}^{2}+y_{n}^{2}+1}\right)\right]\right\}.$

(4)

Assuming that $A=1$ and $Be^{C}=u_{i}$, with $i=x$ for the real component or $i=y$ for the imaginary componedt of the mapping and defining $t_{n}=\theta-\frac{\phi}{x_{n}^{2}+y_{n}^{2}+1}$. Then the map can be simplified as:

$$z_{n+1}=1+u_{i}\left[(x_{n}\cos t_{n}-y_{n}\sin t_{n})+i(x_{n}\sin t_{n}+y_{n}\cos t_{n})\right].$$

(5)

Finally, by separating the real and imaginary parts of $z_{n+1}$, we obtain the discrete two-dimensional non-linear mapping that describes the system’s dynamics:

$\displaystyle S:\left\{\begin{array}[]{ll}x_{n+1}&=\text{Re}(z_{n+1})=1+u_{x}(x_{n}\cos t_{n}-y_{n}\sin t_{n}),\\ y_{n+1}&=\text{Im}(z_{n+1})=u_{y}(x_{n}\sin t_{n}+y_{n}\cos t_{n}),\end{array}~{},\right.$

(8)

where $t_{n}=\theta-\frac{\phi}{x_{n}^{2}+y_{n}^{2}+1}$ and $u_{x}\in[0,1]$ and $u_{y}\in[0,1]$ are the dissipation factors for the real and imaginary components of the map. From now on, we will fix $\theta=0.4$ and $\phi=0.6$. These factors modulate the transformation of the real and imaginary parts of the complex variable through each iteration. This transformation effectively maps the original complex map into a system of coupled real-valued equations, which can be analyzed to study the dynamical behavior of the system. It is important to mention that if $u_{x}=u_{y}=1$ all the results for the Hamiltonian area-preserving map are recovered. From the map $S$ [see Eqs. (8)], one can easily obtain the Jacobian matrix, $J$, which is defined as

$$J=\begin{pmatrix}\frac{\partial x_{n+1}}{\partial x_{n}}&\frac{\partial x_{n+1}}{\partial y_{n}}\\ \frac{\partial y_{n+1}}{\partial x_{n}}&\frac{\partial y_{n+1}}{\partial y_{n}}\end{pmatrix}$$

(9)

with coefficients given by the following expressions

	$\displaystyle\frac{\partial x_{n+1}}{\partial x_{n}}$	$\displaystyle=$	$\displaystyle u_{x}\left[\cos t_{n}-(x_{n}\sin t_{n}+y_{n}\cos t_{n})\frac{\partial t_{n}}{\partial x_{n}}\right]$		(10)
	$\displaystyle\frac{\partial x_{n+1}}{\partial y_{n}}$	$\displaystyle=$	$\displaystyle u_{x}\left[-\sin t_{n}-(x_{n}\sin t_{n}+y_{n}\cos t_{n})\frac{\partial t_{n}}{\partial y_{n}}\right]$		(11)
	$\displaystyle\frac{\partial y_{n+1}}{\partial x_{n}}$	$\displaystyle=$	$\displaystyle u_{y}\left[\sin t_{n}+(x_{n}\cos t_{n}-y_{n}\sin t_{n})\frac{\partial t_{n}}{\partial x_{n}}\right]$		(12)
	$\displaystyle\frac{\partial y_{n+1}}{\partial y_{n}}$	$\displaystyle=$	$\displaystyle u_{y}\left[\cos t_{n}+(x_{n}\cos t_{n}-y_{n}\sin t_{n})\frac{\partial t_{n}}{\partial y_{n}}\right]$		(13)

where $\frac{\partial t_{n}}{\partial x_{n}}$ and $\frac{\partial t_{n}}{\partial y_{n}}$ are given by

$$\frac{\partial t_{n}}{\partial x_{n}}=\frac{2\phi x_{n}}{(1+x_{n}^{2}+y_{n}^{2})^{2}},\quad\frac{\partial t_{n}}{\partial y_{n}}=\frac{2\phi y_{n}}{(1+x_{n}^{2}+y_{n}^{2})^{2}}$$

After some calculation one can show that the map is area preserving only when $u_{x}=u_{y}=1$ since the determinant of the Jacobian matrix is $\det(J)=u_{x}u_{y}$. As an ilustration, Figure 1 illustrates the structure of the phase space for the map 8 with varying values of $u_{x}$ and $u_{y}$. Considering $u_{x}=u_{y}=1$ [ Fig. 1(a)], we observe a chaotic sea interspersed with a set of Kolmogorov-Arnold-Moser (KAM) islands. In this scenario, the phase space exhibits regions of chaotic behavior, where trajectories appear to move unpredictably and ergodically. Amidst this chaotic sea, KAM islands emerge as stable, quasi-periodic regions where the motion is regular and confined. These islands represent the remnants of the invariant tori that survive the perturbation introduced by the system’s non-linearity. The coexistence of chaotic regions and KAM islands highlights the complex and rich structure of the phase space, illustrating the intricate interplay between order and chaos in the dynamical system. As dissipation is introduced, the phase space structure undergoes significant changes. The system exhibits a complex interplay between chaotic and periodic dynamics, as demonstrated in Fig. 1(b), where both chaotic and periodic orbits coexist. Specifically, Fig. 1(b) illustrates an attractive fixed point (indicated by a $\times$) alongside a chaotic attractor, each with its own basin of attraction, as further detailed in Fig. 2. As the system parameters vary, the dynamics undergo significant changes. For instance, in Fig. 1(c), a chaotic attractor emerges, signaling a transition from periodic to chaotic behavior. However, with further decreases in the dissipation parameters, the chaotic attractor is eventually replaced by a set of attracting fixed points, as shown in Fig. 1(d). This progression underscores the system’s intricate and diverse behavior as it navigates different regions of parameter space. Figure 2 presents the basins of attraction for both the attracting fixed point (orange) and the chaotic attractor (blue) depicted in Fig. 1(b). To generate this figure, the ranges $u_{x}\in[-11,30]$ and $u_{y}\in[-25,25]$ were partitioned into grids of 2000 intervals each, resulting in a total of $4\times 10^{6}$ distinct initial conditions. Each initial condition was iterated up to $n=5\times 10^{5}$. While other attractors may exist, their basins of attraction, if present, are either too small to detect or lie outside the scope of the initial conditions considered in this study.

Furthermore, as we have shown, by slightly changing the control parameter, the behavior of the initial conditions can shift from chaotic to periodic. In Fig. 3 (a-b), we see that for lower values of $u_{x}$, the system tends to exhibit periodic behavior, characterized by regular oscillations of $x$ and $y$. As $u_{x}$ increases, these periodic regions become interspersed with chaotic intervals, where the trajectories of $x$ and $y$ become irregular and sensitive to initial conditions. The transition between these behaviors is marked by bifurcations, indicating changes in the system’s stability and the emergence of new dynamical regimes. We can observe the transition from regular to chaotic behavior by examining the bifurcation diagram. In this analysis, we consider the case where $u_{y}=0.6$. To explore typical behaviors, specifically the bifurcation diagrams as the control parameter $u_{x}$ varies, we use the initial conditions $x_{0}=0.2$ and $y_{0}=0.3$. Figure 3(a) shows the behavior of $x$ plotted against the control parameter $u_{x}$, where a sequence of period-doubling bifurcations is evident. A similar sequence is observed for the asymptotic variable $y$, as shown in Figure 3(b). It is important to note that the bifurcations of the same period in both (a) and (b) occur for the same values of the control parameter $u_{x}$. Feigenbaum Feigenbaum (1978, 1979) was the first to observe a “universal” feature in the behavior of bifurcations in dynamical systems. Specifically, he noticed that as a system transitions to chaos through a series of period-doubling bifurcations, the ratios of the differences in control parameter values at which these bifurcations occur converge geometrically at a constant rate, denoted as $\delta$. This discovery indicates that there is a universal behavior in a wide range of systems approaching chaos.

The procedure to determine the Feigenbaum constant $\delta$ is methodical: (a) Identify the bifurcation points: Let $u_{x}(1)$ be the control parameter value at which a period-1 orbit (a single stable cycle) bifurcates into a period-2 orbit (a stable cycle that repeats every two periods); (b) Continue the process: Let $u_{x}(2)$ be the value where the period-2 orbit bifurcates into a period-4 orbit, and $u_{x}(3)$ where the period-4 orbit bifurcates into a period-8 orbit, and so on; (c) Generalize the parameter values: In general, the parameter $u_{x}(n)$ corresponds to the control parameter value at which a period-$2^{n}$ orbit is born. The Feigenbaum constant $\delta$ is then defined as the limit of the ratio of successive differences between these control parameter values as $n$ approaches infinity. Mathematically, it is expressed as:

$$\delta=\lim_{n\rightarrow\infty}\frac{u_{x}(n)-u_{x}(n-1)}{u_{x}(n+1)-u_{x}(n)}.$$

(14)

This constant $\delta$ captures the rate at which the bifurcations occur and is found to be the same for a wide variety of dynamical systems, highlighting the universality of this behavior. The theoretical value of the Feigenbaum constant $\delta$ is approximately $4.669201609\ldots$. This value has been confirmed through both numerical calculations and experimental observations, and it plays a crucial role in the understanding of the transition to chaos in nonlinear dynamical systems. How can we determine the parameter where bifurcation occurs? One effective tool is the Lyapunov exponent. As discussed by Eckmann and Ruelle Eckmann and Ruelle (1985), the Lyapunov exponents are defined as:

$$\lambda_{j}=\lim_{n\rightarrow\infty}{1\over{n}}\ln|\Lambda_{j}|~{}~{},~{}~{}j=1,2~{}~{},$$

(15)

where $\Lambda_{j}$ are the eigenvalues of $M=\prod_{i=1}^{n}J_{i}(x_{i},y_{i})$ and $J_{i}$ is the Jacobian matrix evaluated over the orbit $(x_{i},y_{i})$. However, a direct implementation of a computational algorithm to evaluate Eq. (15) has a severe limitation to obtain $M$. For the limit of short $n$, the components of $M$ can assume different orders of magnitude for chaotic orbits and periodic attractors, making the implementation of the algorithm impracticable. To avoid such a problem, $J$ can be written as $J=\Theta T$ where $\Theta$ is an orthogonal matrix and $T$ is a right up triangular matrix. $M$ is rewritten as $M=J_{n}J_{n-1}\ldots J_{2}\Theta_{1}\Theta_{1}^{-1}J_{1}$, where $T_{1}=\Theta_{1}^{-1}J_{1}$. A product of $J_{2}\Theta_{1}$ defines a new $J_{2}^{\prime}$. In a next step, one can show that $M=J_{n}J_{n-1}\ldots J_{3}\Theta_{2}\Theta_{2}^{-1}J_{2}^{\prime}T_{1}$. The same procedure can be used to obtain $T_{2}=\Theta_{2}^{-1}J_{2}^{\prime}$ and so on. Using this procedure, the problem is reduced to evaluate the diagonal elements of $T_{i}:T_{11}^{i},T_{22}^{i}$. Finally, the Lyapunov exponents are given by

$$\lambda_{j}=\lim_{n\rightarrow\infty}{1\over{n}}\sum_{i=1}^{n}\ln|T_{jj}^{i}|~{}~{},~{}~{}j=1,2~{}~{}.$$

(16)

If at least one of the $\lambda_{j}$ is positive, the orbit is said to be chaotic. Figure (3)(c) shows the behavior of the Lyapunov exponents corresponding to Fig. 3(a-b). It is evident that when bifurcations occur, the exponent $\lambda_{j}$ vanishes. Based on the numerical data obtained through the calculation of Lyapunov exponents, Feigenbaum’s $\delta$ for the Ikeda map is found to be $\delta=4.669248396257327\ldots$. Here, we have considered bifurcations up to the tenth order (see Table I), and our result is in good agreement with Feigenbaum’s $\delta$ up to $10^{-4}$.

Furthermore, the detailed examination of these bifurcation diagrams reveals the underlying mechanisms driving the system towards chaos. Overall, these bifurcation diagrams serve as a crucial tool for visualizing and understanding the complex dynamics of the system as it responds to changes in the control parameter $u_{x}$. The clear sequence of period-doubling bifurcations leading to chaos underscores the rich and intricate nature of the system’s behavior under the given conditions.

To obtain a more comprehensive understanding of the model’s dynamics, we investigate the parameter space where both components of the dissipation parameters, $u_{x}$ and $u_{y}$, vary. By systematically varying $u_{x}$ and $u_{y}$, we can map out the regions of the parameter space that correspond to different dynamical behaviors. This includes identifying zones of periodic motion, chaotic behavior, and the transitions between them. Such a parametric study provides insights into the robustness of the observed phenomena and helps to uncover the underlying structure of the phase space.

The detailed exploration of the parameter space also reveals how the interplay between the dissipation parameters $u_{x}$ and $u_{y}$ influences the overall dynamics. For example, increasing $u_{x}$ might enhance the stability of periodic orbits, while varying $u_{y}$ could affect the onset of chaos. By understanding these relationships, we can better predict the system’s response to changes in the control parameters and potentially develop strategies to control or exploit the dynamics for specific applications.

To thoroughly investigate the parameter space of the system described in equation 8, we systematically varied the dissipation parameters $u_{x}$ and $u_{y}$. For each combination of these parameters, after discarding a significant transient period, we calculated the Lyapunov exponent, which is a crucial measure of the system’s sensitivity to initial conditions. The value of the Lyapunov exponent determines whether the system exhibits chaotic or regular behavior.

Once the Lyapunov exponent was computed for each pair $(u_{x},u_{y})$, we assigned a corresponding color to visually represent the stability or chaos of the system within the parameter space. Specifically, Figure 4 provides a detailed view of the parameter space for the Ikeda map, where a shrimp-shaped structure becomes evident. The color scale used in this figure is carefully designed to distinguish between different types of behavior: regions exhibiting regular, stable dynamics are colored in shades of red to yellow, while regions displaying chaotic behavior are colored in shades of green to blue. This shrimp-shaped structure is a well-known feature in the study of dynamical systems and indicates the presence of complex bifurcation patterns.

Our findings in this parameter space are in agreement with previous studies, particularly with the results presented in Refs. Gallas (1993); Baptista and Caldas (1997); Rössler et al. (1989); Mackay and Tresser (1986); Celestino et al. (2014). To construct the figure, we divided the range of both $u_{x}$ and $u_{y}$ into 2000 equal intervals, resulting in a grid of $4\times 10^{6}$ distinct parameter combinations. For each combination, we began with initial conditions $x_{0}=0.2$ and $y_{0}=0.3$, and then iteratively followed the attractor as the parameters $u_{x}$ and $u_{y}$ were incrementally adjusted. After each increment, the final state of the previous combination served as the new initial condition, ensuring continuity in the exploration of the attractor’s evolution.

However, it is important to note that this approach can inadvertently omit information about other potential attractors. This omission occurs because the chosen initial conditions might lie within the basin of attraction of a specific attractor, thereby excluding the exploration of other possible attractors that might exist for different initial conditions.

Figure 4(b) offers a magnified view of the main structure in Figure 4(a), revealing the intricate bifurcation patterns characteristic of each shrimp-shaped structure. Each shrimp features a central body, generally representing stable dynamics, followed by an infinite sequence of bifurcations that adhere to the pattern $k\times 2^{n}$, where $k$ denotes the period of the central body. These bifurcations give rise to increasingly complex and chaotic dynamics, showcasing the rich and multifaceted behavior of the Ikeda map as the parameters are varied. The detailed bifurcation sequences provide valuable insights into the system’s transition from regular to chaotic dynamics, highlighting the underlying complexity and sensitivity of the parameter space.

This study explores the rich dynamical behavior of the Ikeda map, a nonlinear system originally designed to model light in a nonlinear optical cavity. Despite its seemingly simple form, the Ikeda map exhibits a complex array of dynamical phenomena, including periodic orbits, chaotic attractors, and intricate bifurcation structures.

Through our analysis, we have demonstrated that the Ikeda map features a diverse range of dynamical regimes depending on the values of the dissipation parameters $u_{x}$ and $u_{y}$. The phase space analysis reveals a striking transition from chaotic to periodic behavior as these parameters vary. Specifically, the map exhibits a period-doubling bifurcation cascade, characteristic of systems approaching chaos. We have successfully quantified the Feigenbaum constant $\delta$ associated with these bifurcations, finding it to be consistent with the theoretically established value of approximately $4.669201609...$. This confirms the universal nature of the bifurcation cascade observed in the Ikeda map and aligns with the behaviors seen in other dynamical systems.

Furthermore, the identification of shrimp-shaped structures within the parameter space provides a compelling visualization of the map’s dynamical complexity. These structures represent intricate domains of stability interspersed with chaotic regions, underscoring the interplay between order and chaos in the map’s dynamics. Our use of Lyapunov exponents as a tool for distinguishing chaotic from regular regions has proven effective, offering a clear method for analyzing the map’s behavior.

Overall, this work highlights the Ikeda map’s utility as a model for studying chaotic systems and bifurcation phenomena. The observed transitions between periodic and chaotic behaviors, along with the rich parameter space structures, illustrate the map’s capability to capture the essence of dynamical complexity. Future research could extend these findings by exploring other parameter ranges or by examining the influence of additional nonlinearities, further enhancing our understanding of chaos and stability in nonlinear dynamical systems.

I would like to dedicate this work to the memory of Prof. Dr. Jason Alfredo Carlson Gallas, with whom I had the privilege of meeting during my time as a Postdoctoral Researcher at the Institute of Multiscale Simulation, Friedrich-Alexander-Universität Erlangen-Nürnberg, under the supervision of Prof. Dr. Thorsten Pöschel in 2012.