Complexity and Chimera States in a Network of Fractional-Order Laser Systems

By applying the Adams-Bashforth-Moulton method (ABM), this paper explores the complexity and synchronization of a fractional-order laser dynamical model. The dynamics under the variance of derivative order q and parameters of the system have examined using the multiscale complexity algorithm and the bifurcation diagram. Numerical simulation outcomes demonstrate that the system generates chaos with the decreasing of q. Moreover, this paper designs the coupled fractional-order network of laser systems and subsequently obtains its numerical solution using ABM. These solutions have demonstrated chimera states of the proposed fractional-order laser network.


Introduction
The fractional-order calculus has solved various important problems, which eventually explained some of the natural phenomena [1][2][3]. In recent years [4][5][6], its application to chaotic systems-based communications and encryption has attracted lots of attention, and this is due to the high efficiency of fractional-order chaotic models [7,8]. Nevertheless, the literature has only a few studies about the complex dynamics of fractional-order laser chaotic systems. Therefore, this study discovers the complex dynamics of a fractional-order laser system for broadening the applications of such models.
Complex dynamics analysis is crucial for both integer-order and fractional-order chaotic systems, and this is due to its importance for determining proper parameters. The most used methods are the phase portraits, bifurcation diagram, 0-1 test [9], and maximum Lyapunov exponents. The complexity of the system's time series is considered the most efficient technique for extracting the chaotic system's properties [10][11][12][13]. For the time being, the permutation entropy (PE) [14], the statistical complexity measure (fuzzy entropy (FuzzyEn) [15], and SCM [16]), the spectral entropy (SE) [17], and the C 0 algorithm [18] are the most widely used for estimating the complexity of the time series. Meanwhile, the SE [17] and C 0 [18] algorithms are more precise methods to evaluate the complexity of chaotic time series. Recently, Coast et al. [19] proposed multiscale SampEn (MSampEn) to estimate the complexity of time series more reliably and accurately. By employing the idea of multiscale, the C 0 and the multiscale SE algorithms are employed to examine the complexity performance of the fractional-order chaotic laser model [20].
Therefore, investigating the complexity of an isolated laser system is important since it provides a basis for the study of the dynamical behaviors of a laser network. As a matter of fact, a complex network coupled by nonlinear systems such as neurons has become a hot research topic. Structure, dynamics, synchronization, and the stability of different networks have been widely investigated [21][22][23]. As a result, coupling is widely used in the oscillators and nonlinear systems for the synchronization control due to its simplicity and efficiency. Recently, chimera states in the network of coupled oscillators show emergent patterns, which has attracted much attention of scholars [24][25][26][27]. Moreover, different kinds of fractional-order calculuses have been introduced to the nonlinear systems [28][29][30], and it is found that those systems have more complex dynamics and the derivative order is also a bifurcation parameter. As a result, introducing the fractional derivative to the complex networks becomes necessary.
Specifically, researchers also keep abreast of research hotspots of dynamics and synchronization control in the laser networks [31][32][33]. For instance, Gonzélez et al. [31] investigated synchronization emerges in a small network of delay-coupled lasers anḑ Safak et al. [32] investigated the performance of a synchronous multi-color laser network with daily sub-femtosecond timing drift. It should be noted that the coupled laser networks are widely investigated due to the simple controls. Since the fractional-order calculuses provide a more effective way for modeling of nonlinear systems [28][29][30]34], it is crucial to study how the fractional derivative order, the coupling strength and the structure affect the complex dynamics of the network, especially to analyze different chimera states in such networks. In this paper, we construct a ring network of the fractional-order laser system and analyze its dynamical behaviors, where the system parameters are chosen according to the dynamical analysis results. Moreover, the appearance of a chimera state is explored by varying the network coupling strength, the connection number of each node and the fractional-order derivative order.
The remainder of this paper is arranged as follows. Section 2 introduces the system of the fractional-order chaotic laser model and obtains its numerical solution by applying Adams-Bashforth-Moulton (ABM). In Section 3, complexity in the fractional-order laser system is analyzed. In Section 4, the coupled ring network of the fractional-order laser systems is examined, and the simulation and analysis of the network are carried out. The results are summarized in Section 5.

Dynamics of the Fractional-Order Laser System
Here, we briefly describe the fractional-order of a laser model. Subsequently, we obtain its solution by the ABM.

System Model
Maxwell-Bloch equations are partial differential equations that formulate the chaotic scenario of a laser model. These equations describe the dynamics of a two-state quantum interacting with the electromagnetic mode of an optical resonator. Using some transformation on Maxwell-Bloch equations, we introduce a set of ODEs as follows [35,36] where Q 1 , Q 2 , a, b, c, and d are control parameters. Chaos and multi-periodic orbits of the system (1) appear for different sets of parameters [35].
Suppose that q ∈ (0, 1] is a fractional derivational order. Thus, the corresponding fractional-order laser system is defined by where, D q t 0 is known as the Caputo fractional-order derivative [37], which is given by

Solution Based on the Adams-Bashforth-Moulton Algorithm
For obtaining the solution of the system (2), we use the Adams-Bashforth-Moulton (ABM) algorithm [38], which is given by 0 represent the initial condition. Now, let us introduce the Volterra integral equation, which is given by where, is the fractional-order integral. It is crucial to state that the Volterra integral equation is equivalent to the system (4). However, consider h = T/N, t j = jh (j = 0, 1, · · · , N ∈ Z + ), then the discrete solution of system (4) is denoted as in which Now, we employ the function "FDE12.m" (https://www.mathworks.com/matlabcentral/ fileexchange/32918-predictor-corrector-pece-method-for-fractional-differential-equations (accessed on 20 January 2021)) under Matlab software to implement the ABM algorithm. However, Figure 1 shows the phase diagrams of system (2) for different values of the derivative order q. As shown in this figure, a periodic circle and chaotic attractor are observed with derivative order q = 0.98 and 0.998, respectively.

Bifurcation Analysis
To investigate the dynamics of system (2), three different cases are considered as follows. It is crucial to state here that the initial values of the system in these three cases are set as x 0 = [0.1, 0.2, 0.3, 0.4, 0.5], and the sequence length is N = 10 4 with removing the first 10 4 points of data.
(1) Fix a = 2.0, b = 0.8, c = 70, d = 0.01, Q 1 = 0.005, Q 2 = 0.01 and vary q from 0.98 to 1 with a step size of 8 × 10 −5 . Figure 2a,b depict the dynamics of system (2) using two significant analyses, which are the bifurcation diagram and Lyapunov exponents, respectively. Multiperiodic or periodic appears when q takes smaller values, but it becomes to be chaotic when q > 0.993. As a result, with the derivative order, system (2) owns bright dynamics.
(2) Fix b = 0.8, c = 70, d = 0.01, Q 1 = 0.005, Q 2 = 0.01, q = 0.999 and vary a from 1 to 5 with a step size of 0.016. Figure 2b,e shows that the system is totally chaotic with different values of parameter a, and the largest Lyapunov exponents increase with the growth of parameter a. It means that the complexity of the system increases with the growing values of a.  Figure 2c,f, the system becomes too chaotic when c > 148. However, starting from q = 0.98, the maximum Lyapunov exponent is close to zero. According to the bifurcation diagram results in the above three cases, the system is chaotic with a relatively large derivative order q, the parameters a and c. Vary q from 0.95 to 1 with a step size of 5 × 10 −4 , and vary a from 1 to 5 with a step size of 0.04. The largest Lyapunov exponents contour plot of the system is illustrated in Figure 3a. Obviously, the system has relatively large Lyapunov exponents in the top right corner of the parameter plane. Meanwhile, the minimum order for chaos is about q = 0.96 according to Figure 3a. Moreover, we vary the derivative order 0.97 to 1 with a step size of 3 × 10 −4 , and increase c from 70 to 200 with a step size of 1.3. The largest Lyapunov exponents contour plot is depicted in Figure 3b. It is clear that the higher complexity values of the system appear with larger values of q and the smaller values of c.

Multiscale Complexity Analysis
Suppose that the time series {x(n), n = 0, 1, 2, 3, · · · , N − 1}, and its corresponding discrete Fourier transformation is given by where k = 0, 1, · · · , N − 1 and j = √ −1 is the imaginary unit. If the power of a discrete power spectrum with the kth frequency is |X(k)| 2 , then the probability of this frequency is defined as Using the discrete Fourier transformation, the summation runs from k = 0 to k = N/2 − 1. The entropy of this probability distribution with the same summation limits in Equation (8) is [17] Assume the square value of {X(k), k = 0, 1, 2, · · · , N − 1} is In addition, supposeX where ξ is a control parameter, then the C 0 complexity algorithm can be defined as [18] C 0 (x, r) = Now, the consecutive coarse-graining of a one-dimensional discrete time series is given by [19] y τ (j) = 1 where 1 ≤ j ≤ [N/τ], τ represents the non-overlapping window length. Distinctly, for τ = 1, the sequence y τ represents the prime sequence {x(n), n = 0, 1, 2, · · · , N − 1}. However, the multiscale complexity in this paper is defined as [20] and MSE and MC 0 are employed to estimate the complexity of system (2), where the maximum scale factor τ max = 25. Complexity results with q varying are plotted in Figure 4a Figure 5. Firstly, it is obvious that the complexity results coincide with the bifurcation diagram and Lyapunov exponents. Secondly, it means that the system exhibits high complexity for real applications. Here, we vary a from 1 to 5 with a step size of 0.04, and increase c from 70 to 200 with a step size of 1.3. The complexity contour plots with q = 0.98, 0.99 and 0.999 are shown in Figure 6. It can be observed that the wider region of high complexity values appear when q takes larger values. Meanwhile, the system exhibits larger complexity values when parameter c is larger, and the system has wider range of high complexity when the parameter a is relatively larger. In conclusion, q, a, and c determines the complexity performance of the system.

Building of the Laser Network
In this section, the coupled ring network of different identical fractional-order laser systems is proposed. The equations of the network are defined as where κ is the coupled strength, 2K is the number of the nearby nodes, Nodes contains the index of the nearby nodes, and there are N nodes (i = 1, 2, · · · , N) in the network. In this paper, each node contains a fractional-order laser system. As a matter of fact, the node i is connected with its nearby 2K nodes. Figure 7 shows the networks of 50 nodes with different values of K. When K = 1, each node is connected with its two neighbors. The dynamics of this kind of network are investigated by He et al. [39] but with different systems. When K > 1, there are more connections, and the network becomes more and more complicated with the increase of K. To find all the connected nodes (Nodes), Algorithm 1 is designed. The input of variables are i, K and N. Here, we have 1 ≤ i ≤ N, 2K < N. For instance, when K = 1, i = 1 and N = 50, the nodes found are j = 50 and j = 2, while if i = 10, the nodes find are j = 9 and j = 11. The output of the function is a matrix, thus nodes found for i = 1 is 2]. Then the equations can be programmed, realized using the found Nodes, namely, j takes values from the matrix Nodes. (i, K, N).

Synchronization and Chimera States in the Network
According to the complexity results and phase diagrams, which have been shown in Figure 1, we choose the system parameters a = 2.0, b = 0.8, c = 70, d = 0.01, Q 1 = 0.005 and Q 2 = 0.05. When q = 0.98, the system is periodic, while q = 0.998, the system is chaotic. Meanwhile, the number of the nodes in the network is 500, namely, N = 500. The initialization of network is set using random numbers. As a result, in the simulations, we choose it as x 0 = random(5 * N, 1). Meanwhile, let us define the network error with respect to time as It shows the difference between different nodes. If Error(t → ∞) → 0, the network is synchronized. Otherwise, there is no synchronization. Meanwhile, if there are coexisting zero and nonzero values, chimera states can be observed.
The dynamics of the network with different values of q and number of connections K are investigated. Figure 8 demonstrates the spatiotemporal patterns of the ring network of fractional-order laser systems with 500 nodes and couple strength κ = 5. When q = 0.98, the spatiotemporal patterns of the network are shown in Figure 8a-c, and the errors are shown in Figure 8d-f. Meanwhile, q = 0.998, the simulation results are shown in Figure 8g-l. When K = 10, the networks with q = 0.98 and 0.998 are not synchronized. However, when q = 0.98 and K takes larger values, the network is synchronized. Moreover, when q = 0.998 and K equals 50 and 100, it can be seen in Figure 8g-l that the spatiotemporal behaviors demonstrate the coherent and incoherent states, thus there are chimera states. As a result, when q = 0.998 and κ = 5, there is no synchronization with different K. We fix q = 0.998 and K = 100, the spatiotemporal patterns of the network are obtained after three simulations and the results are shown in Figure 9. Since the initial conditions of each experiment is different due to the randomness, it shows that the network is also sensitive to the initial values.
Let the coupling strength κ = 5 and derivative order q =0.98 and 0.998, and vary the connection number K =15, 55, 95, 125 and 155. Each experiment is carried out two times with different initial conditions, and the network errors with respect to time are shown in Figure 10a-d. It shows that the network can be synchronized when q = 0.98 and K > 15. However, when q = 0.998, it shows in Figure 10c,d that the network cannot be synchronized and non-stationary chimera states are observed, where the positions of coherent and incoherent oscillators vary in time. Fix the connection number K = 50 and derivative order q =0.98 and 0.998, and the coupling strength κ is set as 5, 15, 25, 35 and 45. The analysis results are shown in Figure 10e-h. For each q, the experiments are carried two times and they are obtained with different initial conditions. Obviously, the network is synchronized when q = 0.98 and the chimera state is observed when q = 0.998. Generally, larger coupling strength κ and more connections K cannot make the network synchronized. Figure 8. Spatiotemporal patterns and errors of the ring network of fractional-order laser systems with 500 nodes and the couple strength κ = 5: (a-c) Spatiotemporal patterns with q = 0.98 and K equal to 10, 50, 100, respectively; (d-f) errors with q = 0.98 and K equal to 10, 50, 100, respectively; (g-i) spatiotemporal patterns with q = 0.998 and K equal to 10, 50, 100, respectively; (j-l) errors with q = 0.998 and K equal to 10, 50, 100, respectively. Therefore, it can be concluded that the network is hard to be synchronized when the fractional-order laser system is chaotic (q = 0.998). Furthermore, the network is synchronized with proper coupling strength and connection numbers when the laser system is periodic (q = 0.98). Since each node contains a fractional-order laser system, the network is a fractional-order network, and its dynamics are determined by the derivative order q. Meanwhile, states of the network are determined by the initial conditions, and the chimera states are observed in the network. We hold the opinion that the laser network can have complex dynamical behaviors.

Discussion
For the neural process in nature, there are many potential applications regarding the chimera states. For instance, the bump states in the neural system localized regions of coherent oscillation surrounded by incoherence deduced the chimera states. Furthermore, various types of neuronal diseases, including schizophrenia, Parkinson's disease, Alzheimer's disease, and epileptic seizures connect to the sorting of co-existence of synchronization and desynchronization [40]. However, as far as we know, although the optical chimera state is not a new topic and it has already been reported [41,42], there are few reports on the chimera state of the fractional-order laser network. As a matter of fact, an optical chimera state was observed in the experiments and simulations, and it shows that the amplitude-phase coupling is necessary for the formation of the chimera states.
Here, the network investigated is a ring coupled network with limited numbered neighbor links. In real situations, the network can be with other structures. However, the analysis results can be concluded that the fractional-order laser network can have complex dynamics and even show chimera states. The main reason comes from the complex dynamics of the laser system the amplitude-phase coupling. When the original laser system is chaotic, chimera states are observed. Otherwise, when the original laser system is periodic, the network can be synchronized.
However, we still hope there will be more experiments, which can point out the existence of fractional-order laser systems or can better explain the physical significance. Meanwhile, the existence of chimera states in the laser is observed. However, applications of chimera states in the laser system should be further investigated.

Conclusions
This paper has explored the complexity and synchronization of a fractional-order laser system. Multiscale SE and C 0 algorithms have been applied to analyze the complexity performance of the time series of the system. Meanwhile, the ring network of the fractional-order laser models has been built with the presence of periodic and chaotic states. The following conclusions are stressed.
(1) Symmetric chaotic attractors and periodic circuits are observed, and chaos is found in the system.
(2) It can be observed that the fractional-order laser system has a wide range of high complexity regions in the parameter planes.
(3) It exhibits rich dynamics with the variance of q, a, and c, which was verified by a bifurcation diagram, phase diagram, MSE complexity and MC 0 complexity.
(4) Simulation results show that the chimera states can be observed only when the laser systems exhibit chaotic behaviors. Meanwhile, the network is synchronized with the periodic behaviors of the systems.
In fact, our analysis has demonstrated the complex behaviors of the proposed fractionalorder laser network, especially the existence of different kinds of chimera states. Our next work will focus on the applications of chimera states in the fractional-order laser network.  Data Availability Statement: All datasets generated for this study are included in the article.

Acknowledgments:
The authors would like to thank the four anonymous reviewers for their constructive comments and insightful suggestions.

Conflicts of Interest:
The authors declare no conflict of interest.