1. Introduction

Polarization dynamics and its manipulation of innovative vertical-cavity surface-emitting lasers (VCSELs) have been one hot topic in recent years. Of particular interest is to learn how to extract new theoretical insight from the knowledge of their dynamical characteristics. Unlike traditional edge-emitting lasers, free-running VCSELs can produce a plethora of fascinating nonlinear dynamical behaviors, such as polarization switching [1,2], polarization bistability [3], and even polarization chaos [4] or hyperchaos [5]. All of these benefit from their specific geometric structure, which results in the competition between the two preferred orthogonally and linearly polarized modes. These compelling dynamics may offer the potential for engineering applications ranging from microwave signal generation [6,7] to secure communication [8–10]. For example, the nonlinear period-one (P-1) dynamics can be employed for the photonic generation of microwave and millimeter-wave signals [6], while the polarized chaos can act as the broadband carrier for embedding and transmitting information in secure communication [8,9]. Therefore, knowledge and control of the polarization dynamics offered by VCSELs become of paramount interest in polarization-sensitive applications.

On the one hand, much research effort has been paid to understand how VCSELs work in practice. For the first step, it is necessary to combine experimental studies of VCSELs with numerical simulations of their mathematical models. From an optical-electrical-mathematical point of view, modeling the polarization dynamics of a VCSEL is one of the main problems. Considerable efforts have been devoted to developing the model. In 1995, San and his co-workers first presented the well-known four-level Spin-flip model (SFM) to account for the polarization dynamics of the two heavy-hole-related circularly polarized transitions [11,12]. Since then, this model and its extended counterparts, such as the models of VCSELs with optical injection [13,14], optical feedback [15,16], and a saturable absorber [17–20], are adopted successfully to reproduce the aforementioned polarized dynamic behaviors observed in experiments. So far, the majority of previous investigations on the polarization dynamics of VCSELs have been performed in the framework of the SFM without anisotropies misalignment. However, in many practical situations, the anisotropies in the cavity may not be always perfectly aligned [21,22]. Virte et al. demonstrated that asymmetric polarization dynamic behaviors observed in the experiment are consistent with numerical simulations by incorporating a small misalignment between the anisotropies of the cavity into the framework SFM model [3,23]. Moreover, this case also is considered in the spin-VCSEL by Adams et al. [24].

On the other hand, up to now, most previous investigations on polarization dynamics based on SFM models were restricted to a few specific parameter values or intervals, i.e., these polarization characteristics were frequently presented through several isolated time-response plots and one-dimensional bifurcation diagrams, in which means that only partial information on the dynamics was revealed. Capturing the global behavior of the polarization characteristics of a VCSEL is an interesting but challenging task, especially for SFM models involved with multiple parameters. The indisputable fact is that, in practice, there exists a degree of mutual trade-offs and constraints between the parameters of the system, which together contribute to the evolution of the system dynamics. Therefore, then, questions naturally arise: how do the SFM model used in the Refs. [3,23,25,26] would tend to evolve in terms of laser dynamics under the influence of multiple system parameters or in the case of a wider range of system parameters? Are there new dynamics emerging?

Motivated by the facts mentioned above, this contribution aims to offer a systematic characterization of the dynamics of a free-running VCSEL in two-dimensional diagrams of distinct parameter combinations of the modified model. Specifically, we report two kinds of high-resolution stability diagrams, i.e., isospike diagrams and modified Lyapunov stability diagrams, which help us show a rather detailed classification of periodic and chaotic stable oscillations observed in two-parameter space of several parameter combinations. Interestingly, our stability diagrams reveal a plethora of unanticipated regularities (self-similar periodic structures) for a free-running VCSEL such as Arnold tongues and shrimp-shaped structures. Furthermore, some bifurcation curves such as Hopf bifurcations (HB) and torus bifurcations (TB) are obtained by the continuation technology, which helps us better understand the underlying dynamics of the foregoing structures.

The paper is organized as follows. In Section 2, the model and its dimensionality reduction counterpart of the SFM are briefly described. Then Section 3 explains how to obtain stability phase diagrams and the corresponding numerical results with different parameter combinations are shown. Next, the multistability is investigated briefly through the basin of attraction. At last, our conclusions and discussions are presented in Section 4.

2. System modeling

In this section, the theoretical model of the SFM including a small misalignment between the amplitude and phase anisotropies is illustrated. The modified rate equations are written as [3,23,25]:

A stability phase diagram is a graph that usually displays regions of varied dynamical behaviors, namely equilibrium points, periodicity, quasiperiodicity, chaos, and even hyperchaos. There are some tools able to characterize the dynamical behavior of a point in a given parameter-space diagram, such as diagrams based on Lyapunov exponents, and isospike diagrams [28]. And the Lyapunov exponents can be obtained through the QR decomposition method, see Ref. [29] for more details. In this work, we prefer the latter approach due to its higher computation efficiency and better dynamical discrimination, however, we will also utilize the former for comparison. Firstly, one covered the bi-parameter plane of interest with a mesh of 2000 × 2000 = 4 × 10⁶ equidistant points. Then, for each point of the mesh, the temporal evolution was obtained by integrating numerically Eqs. (4)–(8) using the standard fourth-order Runge-Kutta algorithm with a fixed time-step of h = 0.001. To avoid misjudgment of the dynamics by the transient behavior, the first 10⁶ steps were discarded and the results of the subsequent 10⁶ steps were taken as asymptotically stable solutions and their spikes (the local maxima) were recorded to count the number of peaks per period. Note that we referred to the number, $n({1,\textrm{}2,\textrm{}3, \cdots } )$, of repetitions of the peak per period in the RCP intensity oscillations, as period-n, or P-n simply. As done in Ref. [28], the periodic patterns with more than 19 were plotted “recycling colors mod 19”, i.e., by taking as their color index the remainder of the integer division of the number of peaks, namely $n({20,\textrm{}21,22, \cdots } )$, by 19. In addition, the fact is that so far there is no general method to assess the dependence of the stability phase on the initial conditions. As the difference between the stability phases obtained by scanning the parameters in different ways is small, herein, we constructed the stability diagrams by performing all integrations starting from an arbitrarily selected initial condition $({\mathrm{\Omega }(0 ),{Q_R}(0 ),{Q_I}(0 ),\textrm{}N(0 ),m(0 )} )= ({0.0975,\textrm{}0.2785,\textrm{}0.569,\textrm{}0.9575,\textrm{}0.9649} )$ by scanning parameters horizontally from left to right. Besides, the computation of high-resolution diagrams is a rather time-consuming task. To solve this problem, we used MATLAB and C language hybrid programming combined with parallel computation to obtain these stability diagrams, and all numerical simulations were carried out on a computer-owned two Intel Xeon CPU E5-2667.

3. Results and discussion

We now describe the broad picture of how the dynamics of the extended SFM vary as some key parameters are changed. The effects of some key parameters on the laser dynamics were preliminarily investigated by using two-dimensional maps of the 0–1 test and permutation entropy [14]. However, here, we concentrate on more subtle and intricate dynamical evolution in these parameter spaces by the combination of the informative isospike diagrams based on RCP temporal intensity and the Lyapunov stability diagrams for a free-running VCSEL where the anisotropy misalignment is taken into count.

3.1 Periodic Arnold tongue

In this section, to obtain a global view of the dynamics and the comprehensive effects of key aforementioned parameters, we first consider the evolutions of the polarization dynamical behaviors depending on the spin-flip relaxation rate ${\gamma _\textrm{s}}$ and the phase anisotropy ${\gamma _p}$. The results for the several different pump currents are displayed in Fig. 1.

 

Fig. 1. The isopike diagrams in the ${\gamma _s} \times {\gamma _p}$ plane for: (a) $\mu = 1.5$, (b) $\mu = 5$, (c) $\mu = 7$, (d) $\mu = 10$.

Download Full Size | PDF

As we can see, with a low value of $\mu $, such as $\mu = 1.5,\textrm{}$ three diverse dynamical regions are identified in the ${\gamma _s} \times {\gamma _p}$ plane, as shown in Fig. 1(a), which includes continuous-wave (CW) oscillations marked in blue, period-one oscillations in green, and a smaller zone of complex dynamics in black. Moreover, the periodic Arnold tongues exist in this scenario, which is investigated later in some detail. As the injection current is increased, e.g., $\mu = 5,\textrm{}$ the region with complex dynamics is expanded dramatically in size, that is, the chaotic or quasiperiodic regions grow in size. At the same time, more periodic islands with different peaks per period embedded in the chaotic region are encountered. Further increasing the injection current, as shown in Fig. 1(c) where $\mu = 7,$ the domain of the chaotic state continues to grow further, and the P-1 region located at the bottom is also extended. For $\mu = 10$, the parameter plane is almost covered by the chaotic region except for a few scattered areas of the various peaks per period, as shown in Fig. 1(d), which means that in this scenario there is a larger parameter region of pure chaos available suited for the chaos-related application. Meanwhile, another conspicuous feature observed in Fig. 1(d) as $\mu $ increases is the P-1 domain at the bottom indicates a well-defined horizontal V-type shape. From Fig. 1(b) to Fig. 1(d), it is not difficult to observe that an increase $\mu $ in also leads to a shift in the minima of the Hopf bifurcation towards greater values of the spin relaxation rate. In summary, the value of $\mu $ affects the shape, size, and type of the regions of nonlinear dynamics as reported in our previous work [14]. Additionally, at the bottom of the diagram, namely for regimes of the smaller linear birefringence, there are many fascinating self-similar structures. To further identify these intriguing and subtle structures, the regions inside white boxes labeled A, B, C, and D are investigated in detail in the following content.

Next, we will focus on the periodic Arnold tongue structures. Box A in Fig. 1(a) is magnified in Fig. 2(a), from which, one can observe three distinct domains consisting of the P-1 oscillations marked in green, non-zero CW laser intensities painted in blue as well as the mixed structure with multiple periodic oscillations and complex dynamics. In the mixed structure, a set of periodic structures numbered as 7, 8, 9, 10, 11, 12, 13, and so on, with the number being related to the period of the representative structure. These structures are similar to the Arnold tongues encountered in both discrete systems [30–32] and continuous systems [33–35]. Thus, this indicates the existence of a period-adding sequence for larger tongues with periodicity increased by the unity one, which is the typical characteristic of periodic accumulation structure. Besides, another distinctive feature of this structure is that each tail of the tongue is hierarchically split into multiple periodic regions with a distinct number of peaks through a cascade of period-doubling bifurcations and finally evolves into chaotic domains with the decrease of the phase anisotropy ${\gamma _p}$. Meanwhile, as the spin-flip relaxation rate ${\gamma _s}$ increases, the distance between the separate tongue becomes smaller and smaller, making it difficult to distinguish the isolated periodic structures with different numbers of peaks unless a more refined local magnification view is used. Since the isospike diagram cannot distinguish quasiperiodic, chaotic as well as hyperchaotic solutions, and the effective quantifier for discriminating the foregoing dynamics is the Lyapunov exponents, an improved Lyapunov stability diagram is obtained by comparing the signs of the first three largest Lyapunov exponents, as shown in Fig. 2(b). According to the instructions of the color bar, it is not difficult to identify the quasiperiodic dynamics and other behaviors that cannot be distinguished in the isospike diagram. Concretely, the panel ${\mathrm{\gamma }_s} \times {\mathrm{\gamma }_p}$ is divided roughly into four regions, which are the region of polarized laser oscillations with steady-state (ST) marked in grey color, chaotic states in red, quasiperiodic states in blue, as well as the green periodic oscillation region in the upper left corner and in the middle of the diagram where the periodic Arnold’s tongues (green) embedded in the quasiperiodic (blue) area. Even though Figs. 2(a) and 2(b) are essentially two very different diagrams, the boundaries of the chaotic and periodic zones delineated by the Lyapunov exponent are in excellent agreement with the one obtained by computing the periodicity. Additionally, the Lyapunov exponent stability diagram gives good insight into the distribution skeleton of various dynamical behaviors in parameter space. We also provide some bifurcation curves, as illustrated in Fig. 2(a), that have been computed by employing AUTO [36] and MATCONT [37]. The panel is overlaid with white curves that correspond to saddle-node bifurcations (SNB). Another curve coded in red is the TB, which is the demarcation line between P-1 oscillations and quasi-periodic oscillations. Specifically, with the increase (decrease) of ${\mathrm{\gamma }_s}$ (${\mathrm{\gamma }_p}$), the periodic polarized laser oscillations with one spike per period disappear, giving rise to an invariant torus featuring quasiperiodic motion. Then the tips of the larger Arnold tongue embedded in the quasiperiodic region are connected with the TB curve. The bifurcation curves present quite good agreement with the isospike diagram and reveal the type of bifurcations that occur in the VCSEL system given by Eqs. (4)–(8) for specific parameter-planes.

 

Fig. 2. The characterization of the Arnold tongue structure: (a) isospike diagram, (b) Lyapunov stability diagram, (c) the first three Lyapunov exponents spectrum and one-dimensional bifurcation diagram, the maximum of the RCP intensity time series, varying ${\gamma _s} \in [{22,32} ]$ with ${\gamma _p} = 10$, (d) Poincaré section diagram of seven representative trajectories showing period-7 ∼ 13 orbits, respectively.

Download Full Size | PDF

It is well known that these periodic structures immersed in quasi-periodic regions feature self-similar properties with pattern repeatability at different length scales. Next, to illustrate the self-similar structure more visually, by increasing ${\gamma _s}$ from 22 to 32 along the cyan line in Fig. 2(a) where ${\gamma _p} = 10$, Fig. 2(c) shows the spectrum of the first three Lyapunov exponents (top) and the bifurcation diagram (bottom) by recording the maximum value of the RCP intensity series. Obviously, these two kinds of diagrams are consistent with each other. As can be seen, it is easy to distinguish the periodic, quasiperiodic as well as chaotic states by observing the signs of the first and third Lyapunov exponents, while the bifurcation diagram illustrates period-adding cascades with the characteristic alternation of periodicity and quasiperiodicity. In the bifurcation figure, there are seven vertical lines labeled from 7 to 13, which refer to the number of peaks presented in one period of the RCP intensity. To further manifest their periodicity, the corresponding Poincaré sections in the case of the parameters defined by these vertical lines are shown in Fig. 2(d). The numbers of the points on the Poincaré sections, marked with different colors, represent periods, which also verify the correctness of the one-parameter bifurcation (Fig. 2(c)) and the isospike diagram. At the same time, by increasing ${\gamma _s}$ from 22 along the cyan line in Fig. 2(a), the RCP laser intensity stays in nonchaotic regions until about ${\gamma _s} = 29$, being characterized by the fact that the largest Lyapunov is approximately equal to zero. Then further increasing ${\gamma _s}$, some inconsecutive chaotic windows appear in the parameter interval of the periodic or quasi-periodic states. Up to now, the Arnold tongue topology of the free-running VCSEL governed by the Eqs. (4)–(8) has been visualized by using the isospike diagram, Lyapunov stability diagram, bifurcation diagram, and Poincaré section diagram.

It should be noted that the dynamics described above also exist in other cases, which are verified by numerous numerical simulations. To avoid duplication, next, we only continue to discuss a more detailed and informative investigation of the self-organization of the Arnold tongues that exist in another situation. Herein, we take the parameter region inside the bottom white rectangular box C in Fig. 1(c) as an example, and the zoomed-in diagram is plotted in Fig. 3(a). As for the rest of the cases, the fine dynamics, for instance, enlargements of the parameter regions located in the white boxes B and D of Fig. 1(b) as well as Fig. 1(d), respectively, are given in Fig. 4.

 

Fig. 3. Two complementary kinds of stability diagrams illustrating the infinite Arnold tongues in the ${\gamma _s} \times {\gamma _p}$ plane: (a) the magnification of the box in Fig. 1(c), (b) the corresponding Lyapunov stability diagram, (c) the magnification of the box in Fig. 3(a), (d) the magnification of the box in Fig. 3(c).

Download Full Size | PDF

 

Fig. 4. A magnification of the rectangle highlighted by white color in (a) Fig. 1(c) and (b) Fig. 1(d), respectively; (c) – (d) the corresponding Lyapunov stability diagram, (e) magnification of the inner region of the box in Fig. 4(b); (f) time series of the RCP intensity for labeled in dots in Fig. 4(e), from the top to the bottom row represents the time series corresponding to the clockwise marked points A1, B1, C1 in Fig. 4(e) respectively. Peak refers to the number of spikes, while T is the oscillation period.

Download Full Size | PDF

By looking back to the colored tongues in Figs. 3(a) and Fig. 3(b), we again verify the existence of the Arnold tongue structures, and these periodic structures immersed in large quasiperiodic regions feature the following characteristics: (i) these larger tongues are spaced at P-1, and their periodicity increases sequentially as ${\gamma _s},{\gamma _p}$ decrease simultaneously to form a period-adding sequence $4 \to 5 \to 6 \to 7 \to 8 \to \cdots $, which is different from the formation mechanism of tongues in Fig. 2 (in the case of the appropriate parameters, their periodicity can be tuned by sole modulation of the spin-flip relaxation rate ${\gamma _s}$, namely, the choice of the material used in the active region of the VCSEL). Similarly, their tips are located on the curve of the TB. (ii) These tongues organize themselves in a Fare treelike sequence which can be generated by adding two rationals $p/q$ and $p^{\prime}/q^{\prime}$ corresponding to any two larger Arnold tongues to give another rational $(p + {p^{\prime}})/(q + {q^{\prime}})$ locating between them, which results in a period $(q + {q^{\prime}})$ structure in between the period q and $q^{\prime}$ structures. Figure 3(c) offers an enlargement of the region inside the white box in Fig. 3(a), which shows the organized structures between the tongues of period-6 and 9. One can easily identify the period-13 tongue with rotation number 2/13 that emerged in between the period-6 and period-7 islands with rotation numbers 1/6 and 1/7, respectively. In this similar manner, there exists a period-15 tongue with a rotation number of 2/15 between the periodic structures of rotation numbers 1/7 and 1/8, and so on. Thus, it is not difficult to obtain an alternative organization of the periodic structure with the following sequence of numbers:

The sequence in Eq. (9) belongs to the incomplete part of the Fare tree containing all rational numbers between 0 and 1. A similar finding was reported in another system by Rao et al. [33]. (iii) To further reveal the more fine-scale dynamic structure between the two tongues, the white rectangle box of Fig. 3(c) is magnified in Fig. 3(d) to have a good view, which shows that the largest structure between period 7 and period-13 structures, is a period-20 structure. Note that 20 is the sum of 7 and 13. Similarly, the largest structure between period-7 and period-20 structures is a period-27 structure, and 27 is the result of the sum of 7 and 20; the largest structure between period-7 and period-27 structures is a period-34 structure, and 34 is the result of the sum of 7 and 27. This behavior repeats itself, which results in the birth of the sequence of numbers ${b_1} = 7,\textrm{}{b_2} = 13,\textrm{}{b_3} = 20,{b_4} = 27,{b_5} = 34,{b_6} = 41,\textrm{}{b_7} = 48,{b_8} = 55, \cdots $. The corresponding closed form can be described as

It should be noted that this sequence is essentially a subset of the foregoing sequence $(q + {q^{\prime}})$. (iv) Finally, a well-known Fibonacci-like sequence is found. Specifically, the periodicities of the tongues are increasing in the following pattern: ${c_1} = 6 \to {c_2} = 7 \to {c_3} = 13 \to {c_4} = 20 \to {c_5} = 33 \to {c_6} = 53 \to \cdots $, where the sequence ${\{{{c_n}} \}_{n \in \mathrm{{\rm N}}}}$ obeys the recurrence relation ${c_{n + 2}} = {c_n} + {c_{n + 1}}$. A crucial characteristic of the sequence is the ratio of consecutive terms, namely, ${c_{n + 1}}/{c_n}$ converges asymptotically to the golden ratio $\mathrm{\Phi } = \left( {1 + \sqrt 5 } \right)/2 \approx 1.61803398$ . Thus, each periodic cycle transits to the next one by increasing its periodicity according to the Fibonacci sequence, with the respective ratio reaching a dense state, i.e., the maximum or final periodicity of the quasi-periodic cycle. Meanwhile, with standard continuation technology, some bifurcation curves are overlapped on the isospike diagram. In addition to the three bifurcation curves mentioned above, the period-double bifurcation (PDB) marked with a black curve, which results in P-1 oscillation move to the period-2 one, is given in Fig. 3(a). Among them, two SNBs are connected with a codimension-2 bifurcation point, namely, the cusp bifurcation labeled by the black star, while the TBs are associated with the birth of the two families of the Arnold tongues. At the same time, with the increase of the line birefringence rate, CW oscillations are unstabilized after undergoing the HB, then evolved into the P-1 one.

In Fig. 4, these self-similar structures can also be unambiguously observed. Particularly, a cascade of Arnold tongues with distinct periodicity is observed [see Figs. 4(a) and 4(b)], among which some tongues cascade covered by more than one different number of peaks. To the best of our knowledge, this was not reported in previous works. Moreover, the corresponding Lyapunov stability diagrams delineate chaotic, periodic, quasi-periodic, and steady-state regions, as shown in Figs. 4(c) and 4(d). To visualize this elaboration, Fig. 4(e) provides a zoomed-in view of the region highlighted by white in Fig. 4(b). In Fig. 4(e), the main body of the largest tongue is composed of five separate domains with a distinct number of peaks, which are the regions characterized by the regular modes with 19, 20, 21, 39, and 40 spikes per period, respectively. Interestingly, these tongues are split into multiple periodic parts, but they still follow similar rules (i) – (iv) mentioned above. It is worth mentioning that, unlike in Fig. 3, the tail region of the tongues is no longer following the route of period-doubling bifurcation to chaos as the parameters change. Furthermore, Fig. 4(f) shows the temporal evolution of dots labeled A1, B1, and C1 marked in Fig. 4(e), which helps us to better understand the distribution of the periodic structures. The red dots indicate the peaks of the time series, which are used to visualize the period length. In Fig. 4(f), peak refers to the number of spikes, and T to the oscillation period. This diagram, in addition to revealing the number of repetitions of the peak within a period, has the very distinctive feature that the duration of each period is exactly identical, i.e., the period length is always equal to 1.182, whereas the number of repeated peaks per period is different. At this point, the self-similar Arnold tongues leading to different period-adding sequences are elaborated in detail.

3.2 Self-similar shrimp structure

We now investigate the newly built organized structure, i.e., the well-known shrimp structure. Since then, this eye-catching dynamical feature has been reported in a variety of nonlinear dynamical systems in diverse fields ranging from biology [38–43] to electronics [44–46] as well as others [47,48]. A magnification of the marked parameter region by the white box in Fig. 2(a) is depicted in Fig. 5, showing the emergence of shrimp-shaped structures of different periods immersed in the chaotic domain. The distribution of the shrimps obeys the law of period-adding sequence of $25 \to 27 \to 29 \to 31 \to 33 \to 35 \to 37 \to 39 \to \cdots $, with the period being increased by 2 (the direction from left to right), and their structures share self-similar features. Meanwhile, the trajectory of the distribution of the period-adding structure exhibits a parabolic shape. Besides, Fig. 5 displays how shrimps organize themselves regularly along very specific directions in the ${\gamma _s} \times {\gamma _p}$ parameter space. All these shrimps are immersed in the domain of the chaotic regions, which in turn are surrounded by the periodic regions. Combing with Fig. 2, the above numerical simulation results also corroborate the fact that the Arnold tongues are intersectional away from the TB curve. More recently, Varga et al. and other groups have proved that the appearance of the patterns is the result of the coalescence and interaction of two pairs of period-doubling and saddle-node codimension-two bifurcation curves [49–51]. Therefore, we speculate that the shrimp-shaped structures observed in this paper share the same formation mechanism as the aforementioned literature.

 

Fig. 5. A zoom of the region highlighted by white in Fig. 2(a), shows the appearance and distribution of the self-similar shrimp-shaped structures.

Download Full Size | PDF

3.3 Dynamics in the plane of other parameter combinations

The evolution of the dynamics in the parameter space has been analyzed in detail in the previous sections. It is also interesting to check dynamics in the other parameter plane, To this end, Figs. 6 and 7 show the isospike diagrams obtained by varying ${\gamma _s} \times \mu $ and ${\gamma _p} \times \mu $ parameter plane, respectively.

 

Fig. 6. The isospike diagrams in the ${\gamma _s} \times \mu $ plane for (a) ${\gamma _p} = 5$, (b) ${\gamma _p} = 25$, (c) ${\gamma _p} = 50$, (d) ${\gamma _p} = 100$.

Download Full Size | PDF

 

Fig. 7. The isopike diagrams in the ${\gamma _p} \times \mu $ plane for (a) ${\gamma _s} = 5$, (b) ${\gamma _s} = 25$, (c) ${\gamma _s} = 50$, (d) ${\gamma _s} = 100$.

Download Full Size | PDF

From these diagrams in Fig. 6, we see that the phase anisotropy ${\gamma _p}$ has a significant influence on the size of the non-periodic region located in the parameter plane. Specifically, when ${\gamma _p} = 5,$ as shown in Fig. 6(a), it can be clearly observed that multiple periodic cascades of self-similar structures are embedded in a complex dynamical region characterized by aperiodic properties in the middle of the diagram. This parametric region surrounded by a white rectangle is enlarged in Fig. 8(a) to allow a closer sight. Increasing ${\gamma _p}$ to 25, the region in the upper right of Fig. 6(a) consisting of CW and P-1 oscillations evolves into non-periodic states except for a small region allowing for the periodic state, as shown in Fig. 6(b). By further increasing ${\gamma _p}$, Fig. 6(c) demonstrates that the range of non-periodic regions begins to shrink, with the part of the CW oscillation region located in the bottom-right entering the P-1 oscillation. Finally, in the case of ${\gamma _p} = 100\; $[Fig. 6(d)], the quasiperiodic and chaotic regions are suppressed and only smaller areas of chaos are preserved. Note that there is no periodic island immersed in the chaotic region marked in black in Fig. 6(d), which may be ideal candidates for chaos-based applications, such as chaotic communication and random number generation.

 

Fig. 8. Magnification of box (a) in Fig. 6(a) and (b) Fig. 7(c), respectively.

Download Full Size | PDF

In Fig. 7, we explore the evolution of the dynamics in the ${\gamma _p} \times \mu $ parameter plane for various spin-flip rates ${\gamma _s}$. It is worth noting that these two parameters can be tunable in experiments, with the injection current being easier to control, followed by the birefringence, which, for example, can be manipulated by applying in-plane anisotropic strain [52]. First, for a small value of ${\gamma _s}$, larger regions of chaos can be obtained, and VCSELs can produce chaotic polarized light when the injected current exceeds the threshold by a small fraction [see Fig. 7(a)]. As ${\gamma _s}$ increases, the size of the chaotic region is reduced, and multiple narrow periodic bands occur within the region. Moreover, the regions located on the left of Fig. 7(a), originally exhibiting non-periodic states, evolve into P-1 oscillation, as shown in Fig. 7(b). With the ${\gamma _s}$ is further increased, not only is the size of the chaotic region shrinking, but more complex periodic structures appear in this region [see Figs. 7(c) and 7(d)]. Intriguingly, the left side of Fig. 7(c) again exhibits multiple V-shaped structures with different periodic states, whose local enlargement is shown in Fig. 8(b). When ${\gamma _s} = 100$, self-similar periodic structures can also be observed. In contrast to Fig. 7(c), the tongues become larger and have been reduced in number (periodicity) in Fig. 7(d). Besides, from Figs. 7(c) and 7(d), we can speculate that the self-similar structures with periodicity can be readily found in the lower birefringence regime.

3.4 Multistability and basin of attraction

So far, in this work, the Arnold tongue, shrimp-shaped structures have been visualized by using high-resolution stability diagrams, bifurcation diagrams, temporal diagrams, and Poincaré section diagrams. Now, of interest is the investigation of multistability, which refers to the coexistence of multiple topologically inequivalent attractors exhibited by a system with fixed parameters and changing initial conditions. Although previous works, both experimental and numerical simulations, have confirmed the existence of bistability or multistability in laser systems [53–55]. A large number of facts confirm the existence of coexisting attractors by scanning the bifurcation parameter both forward and backward in one-dimensional bifurcations, i.e., the attractor following technique. Following this methodology, it is relatively piecemeal to obtain the types of attractor coexistence. How can a more complete classification of attractor coexistence be obtained? This is one of the aims of the present contribution, namely to provide a systematic analysis of the influence of initial conditions on the dynamics of the system, in other words, to perform a detailed classification of the coexistence of multiple attractors.

Figure 9 displays the results of a numerical simulation designed to illustrate the appearance of multistability by changing the way of scanning parameters, for example, following the attractor from left to right [Fig. 9(a)], and from right to left [Fig. 9(b)], from bottom to top [Fig. 9(c)] and from top to bottom [Fig. 9(d)]. In all cases, the initial condition of the first point in the integration is $({\mathrm{\Omega }(0 ),{Q_R}(0 ),{Q_I}(0 ),\textrm{}N(0 ),m(0 )} )= ({0.0975,\textrm{}0.2785,\textrm{}0.569,\textrm{}0.9575,\textrm{}0.9649} ).$

 

Fig. 9. An amplified view of the small red rectangle in Fig. 3(a). It is characterized by the crossover of the periodic structure, leading to bi-stability between multiple attractor pairs. Four pairs of time-evolving diagrams corresponding to four annotated points selected from different overlapping regions show the coexistence of topologically non-equivalent attractor pairs. The colors of the time-domain waveforms are based on their associated periodicity as mentioned in the color bar at the top of the isospike diagram.

Download Full Size | PDF

From the four subfigures in Fig. 9, one can easily distinguish some of the prominent multistability regions. Here, we can clearly observe the overlapping of multiple pairs of coexisting attractors of distinct periodicities. To be accurate, we select four different $({\gamma _s},\textrm{}{\gamma _p})$ pairs, belonging to four overlapping regions, and mark them with blue dots. For each of these four points, we draw two temporal diagrams corresponding to two slightly different initial values, as shown in Fig. 10. Concretely, ${Q_R}(0 )$ is set to 0.2785 and 0.279, respectively, and the remaining initial conditions are fixed as $({\mathrm{\Omega }(0 ),{Q_I}(0 ),\textrm{}N(0 ),m(0 )} )= ({0.0975,\textrm{}0.569,\textrm{}0.9575,\textrm{}0.9649} )$. For ease of reference, the color of each time-domain waveform is matched to the one associated with its periodicity, as used in the isospike diagram. When $({\gamma _s},\textrm{}{\gamma _p}) = ({11.84,\textrm{}4.606} )$, the RCP intensity exhibits period-12 oscillation in the case of ${Q_R}(0 )= 0.2785$, while a different initial state, namely, ${Q_R}(0 )= 0.279$, leads to period-14 oscillation, as shown in Fig. 10(a). As for other points, the RCP intensity tends to oscillate with distinct periods when starting the integration from different initial conditions [see Figs. 10(b) to 10(d)], in other words, a slight change in the initial state leads to completely different attractors in the phase space, which is also an inherent property of non-linear systems. To obtain a global view of the dynamics concerning the initial condition, the investigation on the basins of attraction for the SFM model is a good choice.

 

Fig. 10. Time traces of the RCP intensity corresponding to the points labeled (a) A2, (b) B2, (c) C2, (d) D2 in Fig. 9.

Download Full Size | PDF

By fixing the four sets of parameters in Fig. 9 and the initial conditions $({{Q_I}(0 ),N(0 ),m(0 )} )= ({0.569,\textrm{}0.9575,\textrm{}0.9649} )$, Fig. 11 shows the basins of attraction located on the $\mathrm{\Omega }(0 )\times {Q_I}(0 )$ initial condition plane, which corroborates the coexistence of multiple attractors. For $({\gamma _s},\textrm{}{\gamma _p}) = ({11.9,\textrm{}4.606} )$, Fig. 11(a) displays the coexistence of period-12, period-14, and period-18 as well as chaotic attractors. It is evident that the bi-parameter regions with overlapping colors in Fig. 9 are not only coexisting with attractors of the corresponding state, but also with other types of attractors. Next, we will report two self-similar structures, namely, the “eye” of the anti-chaos [56] and periodicity. The former is an attractive nested ring structure, more recently presented by Rao et al. in a memristor-based Shinriki’ circuit, which is different from the “eye” of the chaos formed by the accumulation of a sequence of forward and reverse period-doubling bifurcation [57], and its central region is periodic rather than chaotic motion. There exists a similar eye of the anti-chaos in Figs. 11(b) and 11(c) for $({\gamma _s},\textrm{}{\gamma _p}) = ({12.42,\textrm{}4.616} )$ and $\textrm{}({\gamma _s},\textrm{}{\gamma _p}) = ({12.45,\textrm{}4.56} )$, respectively. It should be noted that the structure of the eye of anti-chaos encountered in the present contribution is more like that of an eye, for example, the outermost ring of chaos marked in black forms the outer outline of the eye, while the inner ring of chaos and the periodic structure contained within it construct the pupil of the eye. Furthermore, another self-organization, as shown in Fig. 11(d), i.e., the eye of the periodicity, to the best of our knowledge, is observed for the first time in a free-running VCSEL. Specifically, there are three types of periodic attractors coexisting in the basin of attraction, namely the period-6, period-7, and period-16 attractors, where the regions of period-16, marked in bright green, is distributed on the top and bottom. Meanwhile, the attraction domains of period-6 and period-7 shaped ring are alternately nested in the initial condition plane, and these structures are characterized by fractal properties. Finally, the pure area of period-6 is located at the center of the nested ring structure, thus resulting in the formation of the eye of the periodicity.

 

Fig. 11. Four basins of attraction for the SFM model related to the four annotated points in the overlap zone of the ${\gamma _s} \times {\gamma _p}$ isospike diagram presented in Fig. 9. Periodicity is also indicated by the respective color bars. For ease of reference, the colors associated with each periodicity are kept the same as those used in the associated isospike diagram.

Download Full Size | PDF

4. Conclusion

In summary, we have explored numerically the global stabilities in a free-running VCSEL by considering the appearance of the misalignment between the linear phase and amplitude anisotropies. We have reported two self-similar structures in multiple parameter spaces: one is the self-similar Arnold tongues with period-adding sequences immersed in the quasiperiodic region, and the other is the shrimp-shaped structures embedded in chaotic regimes. We have offered the possible underlying mechanisms behind the Arnold tongues with the standard continuation techniques, that is, the birth of the Arnold tongues is associated with TB. By further amplifying these self-similar structures, we have found that the Arnold tongues obey the law of period-adding, which is confirmed by the single parameter bifurcation diagram and Poincaré section diagram. Moreover, the sequences of periodicities are similar to the Fibonacci sequence, the part of the Fare tree; there exist some Arnold tongue cascades covered by distinct numbers of peaks, which are verified by the time series diagram. Finally, we have systematically explored the coexisting of the multiple attractors through four isospike diagrams and basin attraction. Of interest are the two self-similar structures we observed in the basin of attraction, namely, the eye of the anti-chaos and periodicity, respectively. It should be emphasized that the occurrence of the latter, to the best of our knowledge, is the first observation in laser systems. These two kinds of fractal basins confirm the existence of multistability in the free-running VCSEL. We hope that these findings will contribute to the theoretical search for a more accurate description of the hierarchical structure of complex non-linear phenomena and provide guidance for the preparation of materials of VCSELs as well as for engineering applications.