Minimal universal laser network model: Synchronization, extreme events, and multistability

1 Introduction

Among the various collective dynamics that can emerge in coupled dynamical systems featured by complex networks – defined as a set of dynamical nodes connected through links in a specific structure and coupling scheme [1] – synchronization holds a prominent position. As a fundamental concept, synchronization denotes a state wherein all oscillators in a network exhibit identical temporal evolutions [2]. Consequently, the solution of the connected system is not only consistent in state space but also in the time domain. This specific state is referred to as full or complete synchronization, while synchronization encompasses a broader class of collective dynamics, indicating the emergence of coherence among the dynamical nodes [3]. Different synchronization patterns can be defined based on the level of coherence observed. For instance, cluster synchronization describes a state where the solutions of connected oscillators converge to different attractors, forming clusters with full coherence [4]. Another notable subgroup of synchronization is the chimera state, wherein the coexistence of one or a few coherent clusters with an incoherent one is observed [5,6]. Additionally, various other types of synchronization include phase synchronization [7], lag synchronization [8], solitary state [9], generalized synchronization [10], and relay synchronization [11]. Each of these synchronization patterns provides insight into the intricate dynamics of coupled dynamical systems on complex networks.

Synchronization is a concept extensively studied across various fields, with particular significance in neuroscience [12–14] and physics [15–17]. In neuroscience, modeling brain networks involves intricate neuronal models and synaptic pathways. These models are crucial for understanding how synchronization manifests within the brain’s complex network architecture [18–20]. In the realm of physics, certain physical systems play a vital role in the study of synchronization. Among these, laser systems are particularly noteworthy, representing a cornerstone of nonlinear dynamics and chaos theory in real-world applications. Numerous versions of laser systems have been proposed in the literature, each offering unique insights into the dynamics of synchronized systems. For instance, Haken introduced the single-mode laser model [21], while Agrawal developed a mathematical model for semiconductor lasers [22]. Additionally, Ciofini et al. introduced the four-level CO₂ laser model [23], while Meucci et al. proposed the minimal universal laser (MUL) model [24] and a modulated laser system [25]. As highlighted by Uchida et al., the synchronization of chaotic laser models is crucial as it facilitates information communication in both analog and digital forms [26]. Recognizing this significance, numerous studies have delved into examining the synchronization of networked laser systems by proposing physically meaningful coupling schemes. For instance, Sugawara et al. [27] and Mariño et al. [28] explored master-slave laser network configurations. DeShazer et al. [29] and Hillbrand et al. [30] focused on phase synchronization in coupled laser systems. Lag synchronization in an undirected ring network of laser systems was investigated by Mihana et al. [31]. Zhang et al. [32] documented chimera state and cluster synchronization exhibited by delay-coupled laser systems. Similarly, Chembo Kouomou and Woafo [33] found cluster synchronization in a two-dimensional laser network, while Röhm et al. [34] reported a chimera state when the connections of the laser systems involve delays. The phenomenon of cluster synchronization has also been observed in a nonlocally coupled semiconductor laser system, as noted by Kazakov et al. [35]. Additionally, Roy et al. [36] investigated the synchronization of hyperchaotic non-autonomous optically modulated CO₂ laser systems, focusing on the effects of ring and star coupling configurations on synchronization cost. Mehrabbeik et al. [15] also investigated the multistability of a modulated laser network and reported the emergence of cluster synchronization, chimera states, and solitary states within the network. These studies collectively contribute to our understanding of the synchronization phenomena in laser systems and their implications for information processing and communication.

Extreme events represent a category of dynamical properties characterized by the sudden emergence of high-amplitude oscillations or fluctuations within a system. These fluctuations are often undesired and unpredictable, posing significant challenges in their management [37,38]. Extreme events can manifest either in the dynamics of an isolated system [39,40] or in the behavior of coupled systems [41–43]. Multistability is another crucial dynamical property observed in both isolated [25,44] and coupled models [15,45]. It denotes a condition wherein the system lacks a unique solution, resulting in the network’s solution being contingent upon the selection of initial states [46]. Extreme multistability [47,48] and megastability [49,50] represent special cases of multistability wherein the system demonstrates an infinite number of solutions. In extreme multistability, these solutions are uncountable, while in megastability, they are countable.

This article explores the synchronization dynamics of the MUL laser model, given its capability to exhibit chaos and adhere to laser dynamical properties within a low-dimensional framework. Following the introduction of the MUL laser model, Section 2 delves into the chaotic dynamics of the MUL model through bifurcation and Lyapunov analysis. In this section, we also introduce a suitable coupling scheme and global configuration alongside the formulation of the MUL network model. Moving to Section 3, we conduct a stability analysis of the synchronization state of the constructed network, employing both analytical methods and numerical approaches for validation. Furthermore, a comprehensive analysis is undertaken to elucidate the prominent synchronization patterns observed within the network. Finally, Section 4 provides a conclusion, summarizing the key findings and implications of the study.

2 Model

The CO₂ laser model’s unique characteristics were captured in a three-dimensional chaotic model introduced in the study by Meucci et al. [24], aligning closely with a laser system’s physics and inherent properties. This proposed laser model is mathematically described as follows:

where the variables x and y represent the fast (output intensity) and slow (population inversion) components, respectively, while z serves as the feedback variable. Notably, the feedback variable z functions as a low-pass filter, with its input comprising the fast variable x combined with a bias term denoted as B 0 . Consequently, it exerts a nonlinear influence on the fast variable, ultimately resulting in linear regulation. Furthermore, letting k = k 0 ( 1 + k 1 z 2 − y ) defines the decay rate of the fast variable, where k 0 represents the non-modulated cavity loss and k 1 denotes the modulation depth. The parameters γ and β describe the decay rates of the slow and feedback variables, respectively. p 0 signifies the normalized pump strength, α is a scaling factor for normalizing the fast variable, and R acts as a control parameter regulating the input of the feedback term.

Under the transformation y → p 0 y , t → t γ , ε 1 → k 0 γ , ε 2 → β γ , B 1 → R α = R γ 2 k 0 , the dimensionless form of System (1), commonly referred to as the MUL model, can be expressed as follows:

With B 0 serving as the control parameter, Figure 1 depicts the bifurcation and Lyapunov analysis of System (1), employing fixed values of ε 1 = 200 , ε 2 = 6 , k 1 = 12 , p 0 = 1.208 , and B 1 = 0.555 . The parameter values are set in accordance with the study by Meucci et al. [24]. The initial condition is set as ( x 0 , y 0 , z 0 ) = ( 0.5 , 1 , 0 ) , and the simulation spans 5,000 time units. Both bifurcation and Lyapunov exponents (LEs) diagrams are computed using the forward method, where the initial condition of the first simulation is selected manually, and the last solution value is employed as the initial condition for the subsequent simulation while systematically exploring the parameter value domain in ascending order.

Figure 1

Dynamical analysis of the MUL system as a function of the parameter B 0 : (a) the bifurcation diagram and (b) the LE spectra. Both diagrams are plotted using the forward method. Other parameters are set at ε 1 = 200 , ε 2 = 6 , k 1 = 12 , p 0 = 1.208 , and B 1 = 0.555 . To initialize the simulations, ( x 0 , y 0 , z 0 ) = ( 0.5 , 1 , 0 ) is considered, and each step of the simulation spans 5,000 time units. In the chaotic region, two distinct chaotic regimes are identified, separated by an interior crisis. In the first chaotic regime, the MUL system displays small-amplitude chaos, while in the second regime, it manifests homoclinic chaos.

Figure 1 illustrates that the MUL model transitions into the chaotic region (at B 0 ≈ 0.1239 ) via a period-doubling bifurcation and exits the chaotic region through a boundary crisis (at B 0 ≈ 0.1247 ). Upon closer examination of Figure 1, two distinct chaotic regimes emerge within the chaotic region. Initially, in the first chaotic regime ( 0.1239 < B 0 < 0.124224 ), the MUL model exhibits small amplitude chaos, as illustrated in Figure 2(a) and (b) for B 0 = 0.124 . Following an interior crisis at B 0 = 0.124224 , the second chaotic regime ensues and persists until B 0 = 0.1247 . Within this regime, the MUL model displays homoclinic chaos, as demonstrated in Figure 2(c) and (d) for B 0 = 0.1243 .

Figure 2

Trajectory (left panels) and temporal pattern (right panels) of the chaotic attractors of the MUL model for (a, b) B 0 = 0.124 , in which the system exhibits small-amplitude chaos, and (c, d) B 0 = 0.1243 , in which the system exhibits homoclinic chaos. Other parameters are set at ε 1 = 200 , ε 2 = 6 , k 1 = 12 , p 0 = 1.208 , and B 1 = 0.555 . To initialize the simulations, ( x 0 , y 0 , z 0 ) = ( 0.5 , 1 , 0 ) is considered, and each step of the simulation spans 10,000 time units. In the chaotic region, the MUL system demonstrates homoclinic chaos for higher values of the parameter B 0 .

It is worth mentioning that within the chaotic region, various periodic windows can be discerned. The most expansive ones, located within the first and second chaotic regimes, are extended between 0.124104 < B 0 < 0.124134 and 0.12444 < B 0 < 0.124572 , respectively.

To delve into the synchronization dynamics of the MUL model, we must examine the coupled models within a predetermined network framework. Coupling the MUL models involves integrating the linear diffusive function of the fast variable from neighboring nodes into the feedback mechanism. Consequently, the structure of the MUL network model can be outlined as follows:

For every node indexed by i , X i = [ x i , y i , z i ] T represents a vector encapsulating the dynamical states of the MUL model. Meanwhile, [ f ( X i ) , g ( X i ) , h ( X i ) ] T governs the nodal dynamics. The parameter σ i j denotes the coupling strength, dictating the connection’s intensity between nodes i and j . The network structure of N coupled MUL models, outlining the neighboring nodes for each node i , is represented as G i j in Network (3). By defining L i j as the Laplacian matrix constructed from G i j , where L i j = G i j = − 1 , L i i = ∑ j = 1 N G i j = N , and σ i j = σ for all pairs of i and j , the network can describe the globally coupled MUL model with a uniform coupling strength of σ . Consequently, it can be simplified into the following system:

In the forthcoming section, we explore the synchronization state of Network (4) under homoclinic chaos conditions (with parameters set as identical to those in Figure 2(c) and (d)), leveraging insights from both analytical approaches and numerical simulations.

3 Synchronization analysis

In the synchronization state, all nodes exhibit an identical simultaneous temporal pattern, meaning X 1 = X 2 = ⋯ = X N = X s . Consequently, the temporal behavior of MUL models becomes correlated and coordinated across the network. According to Pecora and Carroll [51], the concept of the master stability function (MSF) suggests that the synchronization state remains stable if all nodes exhibit resilience to locally added perturbations to their synchronized dynamics, expressed as X i = X s + ζ i , where ζ i represents additive perturbations to the dynamics of node i . Here, since the diffusive term becomes zero in the synchronization state, in a linear diffusive coupling function, as in our case, we have X s = X , indicating that the dynamics of each coupled MUL model become akin to the dynamics of an uncoupled MUL model. In simpler terms, if the solution of the perturbation system converges to ζ ˙ = [ ζ ˙ x , ζ ˙ y , ζ ˙ z ] T = [ 0 , 0 , 0 ] T as a stable fixed point, the synchronization manifold can be stable at a specific value of the connection strength σ . In Lyapunov analysis, a negative largest Lyapunov exponent (LLE) signifies fixed-point dynamics. Therefore, the values of σ for which the LLE of the perturbation system turns negative offer the necessary conditions for synchronizing the MUL network. According to the MSF technique, the perturbation system of the MUL network (Network (4)) can be obtained as follows:

where K = σ λ i is the normalized coupling parameter, and λ = [ λ 1 , λ 2 , … , λ N ] T is the vector of eigenvalues of the Laplacian matrix L . This normalization allows the results to be generalized across networks of varying sizes and structures. In our case (global coupling structure), λ 1 , λ 2 , … , λ N = N , and thus, K = σ N . It is worth noting that F v ( . ) represents the derivative of the function F with respect to the variable v . Figure 3(a) depicts the LLE of the perturbation system ( Ψ ) plotted against the parameter K for N = 100 globally coupled MUL systems. As no region of synchronization ( Ψ < 0 ) is discernible in Figure 3(a), we can infer that Network (4) is unable to attain a complete synchronization state. The asynchronizability can also be assessed through numerical analysis by calculating the time-averaged synchronization error, which can be obtained using the following formula:

Figure 3

Analytical (left panel) and numerical (right panel) analysis of the MUL network synchronization with global coupling scheme and N = 100 in terms of (a) the LLE of the perturbation system (MSF analysis) and (b) the time-averaged synchronization error of MUL network (error analysis). Parameters are consistent with those in Figure 2(c) and (d), with initial conditions set randomly around ( x 0 , y 0 , z 0 ) = ( 0.5 , 1 , 0 ) . The MUL network fails to achieve complete synchronization within the specified coupling scheme and configuration.

In Eq. (6), ⟨ ⋅ ⟩ denotes the averaging function over a specific time period t . In a complete synchronization state, where all nodes exhibit identical temporal patterns ( X i = X j = X s ), E = 0 is expected, while E > 0 indicates an asynchronous state. Figure 3(b) illustrates the time-averaged synchronization error plotted against the parameter K . The results confirm that Network (4) fails to achieve synchrony, as no region with E ≈ 0 is detected.

Although the MUL network model cannot achieve complete synchrony, our comprehensive analysis of the collective dynamics of Network (4) reveals intriguing monostable or multistable synchronization patterns. For instance, at lower values of K , the two-cluster synchronization state emerges as the most discernible synchronization pattern. Figures 4 and 5 illustrate the two-cluster synchronization state for K = 0.1 and K = 0.12 , respectively. In both figures, these synchronous clusters also exhibit anti-phase synchronization. However, the trajectories indicate that for K = 0.1 , the network is monostable as the solutions are nearly identical, while for K = 0.12 , the network demonstrates bistability as the solutions differ in their periodic bursting modes.

Figure 4

Illustration of the monostable two-cluster synchronization state in the MUL network model with global coupling structure for K = 0.1 , in terms of (a) spatiotemporal pattern, (b) nodal trajectories, (c) temporal evolutions, and (d) the snapshots of the last time sample. Parameters are consistent with those in Figure 2(c) and (d), with initial conditions set randomly around ( x 0 , y 0 , z 0 ) = ( 0.5 , 1 , 0 ) . The MUL models in one cluster also exhibit anti-phase synchronization with the ones in the other cluster.

Figure 5

Illustration of the multistable two-cluster synchronization state in the MUL network model with global coupling structure for K = 0.1 , in terms of (a) spatiotemporal pattern, (b) nodal trajectories, (c) temporal evolutions, and (d) the snapshots of the last time sample. Parameters are consistent with those in Figure 2(c) and (d), with initial conditions set randomly around ( x 0 , y 0 , z 0 ) = ( 0.5 , 1 , 0 ) . The MUL models in one cluster also exhibit anti-phase synchronization with the ones in the other cluster. The dynamics of the MUL models in each cluster differ slightly from those in the other cluster.

Figures 6 and 7 depict how the rise in the value of K induces multistable chimera states. Moreover, a noteworthy observation is the occurrence of extreme events stemming from the presence of high-amplitude oscillations amidst lower-amplitude ones. To explore the potential occurrence of extreme event scenarios, a qualifying threshold is established as follows:

where μ x j and σ x j denote the mean and standard deviation of the local maxima of the time evolution of the fast variable x of the node j , respectively. Here, the coefficient of the standard deviation σ is assigned a high value to ensure that only events with a significant deviation from the mean value μ are considered. The threshold (Thresh), demonstrated by a black dashed line, is computed over a significant period of 5,000 time units, equivalent to 1 0 6 iterations or samples. Any spike that surpasses the threshold can be identified as an extreme event. Additionally, Figure 7 illustrates that a greater increase in the value of K results in a few high-amplitude oscillations.

Figure 6

Illustration of the multistable chimera state with extreme events in the MUL network model with global coupling structure for K = 0.24 , in terms of (a) spatiotemporal pattern, (b) nodal trajectories, (c) temporal evolutions, and (d) the snapshots of the last time sample. Parameters are consistent with those in Figure 2(c) and (d), with initial conditions set randomly around ( x 0 , y 0 , z 0 ) = ( 0.5 , 1 , 0 ) . Extreme events can be observed in the behavior of both coherent and incoherent MUL models.

Figure 7

Illustration of the multistable chimera state with extreme events in the MUL network model with global coupling structure for K = 0.25 , in terms of (a) spatiotemporal pattern, (b) nodal trajectories, (c) temporal evolutions, and (d) the snapshots of the last time sample. Parameters are consistent with those in Figure 2(c) and (d), with initial conditions set randomly around ( x 0 , y 0 , z 0 ) = ( 0.5 , 1 , 0 ) . Extreme events can be observed in the behavior of both coherent and incoherent MUL models. As the parameter K gradually increases, the number of extreme events decreases.

According to Figures 8 and 9, as K gradually increases, the MUL network model once again demonstrates cluster synchronization states. However, in this scenario, some MUL models enter an oscillation death state, while others synchronize into a few periodic orbits. For example, in Figure 8, the coexistence of a fixed point cluster with three period-2 clusters can be detected for K = 0.26 . On the other hand, for K = 0.3 , Figure 9 shows that the number of periodic clusters reduced to one period-1 cluster. The greater increase in the K values makes the periodic cluster become smaller in amplitude, and finally, at K = 1 , two fixed point clusters coexist, as shown in Figure 10. In this condition, all MUL models achieve two different oscillation death states.

Figure 8

Illustration of the multistable four-cluster synchronization state in the MUL network model with global coupling structure for K = 0.26 , in terms of (a) spatiotemporal pattern, (b) nodal trajectories, (c) temporal evolutions, and (d) the snapshots of the last time sample. Parameters are consistent with those in Figure 2(c) and (d), with initial conditions set randomly around ( x 0 , y 0 , z 0 ) = ( 0.5 , 1 , 0 ) . Here, three period-2 clusters coexist with a fixed point cluster.

Figure 9

Illustration of the multistable two-cluster synchronization state in the MUL network model with global coupling structure for K = 0.3 , in terms of (a) spatiotemporal pattern, (b) nodal trajectories, (c) temporal evolutions, and (d) the snapshots of the last time sample. Parameters are consistent with those in Figure 2(c) and (d), with initial conditions set randomly around ( x 0 , y 0 , z 0 ) = ( 0.5 , 1 , 0 ) . Here, a period-1 cluster coexists with a fixed point cluster.

Figure 10

Illustration of the multistable two-cluster synchronization state in the MUL network model with global coupling structure for K = 1 , in terms of (a) spatiotemporal pattern, (b) nodal trajectories, (c) temporal evolutions, and (d) the snapshots of the last time sample. Parameters are consistent with those in Figure 2(c) and (d), with initial conditions set randomly around ( x 0 , y 0 , z 0 ) = ( 0.5 , 1 , 0 ) . Here, two fixed-point clusters coexist.

We also observe that higher values of K cause the fixed points to approach each other gradually, yet the network never reaches the same solution. This observation aligns with the outputs from the MSF analysis and elucidates the transitions observed throughout the error and MSF diagrams. We also discover that the identified synchronization patterns are nearly robust to the selection of the initial condition around ( x 0 , y 0 , z 0 ) = ( 0.5 , 1 , 0 ) .

4 Conclusion

In this study, we introduced a network model of MUL systems by incorporating the linear diffusive function of neighboring fast variables into the feedback term. Subsequently, we investigated the synchronization dynamics of this constructed networked model, employing a global coupling scheme. Leveraging the MSF techniques, we provided analytical insights to assess the stability of the synchronization state of the MUL network model. Our results revealed that the coupled MUL models fail to achieve a complete synchronization state, a finding corroborated through numerical analysis by computing the time-averaged synchronization error values. Moving forward, we conducted comprehensive numerical analyses to explore the synchronization patterns within the network. Our findings unveiled that MUL models tend to synchronize into different clusters under varying coupling parameter values. Specifically, in both low and high coupling parameter regimes, the emergence of synchronization clusters was observed. However, under weak connections, two clusters exhibiting homoclinic periodicity coexist, evolving into an anti-phase synchronization state. As connections strengthen, the coexistence of small-amplitude periodic orbits alongside a fixed-point cluster is detected, indicating that some MUL models fail to oscillate and converge to an oscillation death state. With the further strengthening of connections, all MUL models cease oscillating, achieving two distinct synchronous oscillation death states. Despite the increasing proximity of fixed-point solutions, synchronous collapse is never achieved, and the network solution remains almost unchanged until the MUL models become unstable. Conversely, at intermediate coupling strengths, chimera states emerge, accompanied by extreme events observed in MUL model synchronization. Across most observed synchronization patterns, multistability is evident, with solutions demonstrating robustness to variations in the selection of initial conditions around a working point. The results obtained from the coupled MUL network model have practical implications for optimizing synchronization in laser arrays, optical communication systems, and other networked systems. These findings contribute to designing more efficient and stable systems by providing insights into how network size, coupling strength, and structure affect synchronization behavior.

Introducing fractional-order dynamics could allow for a more comprehensive understanding of systems where historical states influence present dynamics [52,53]. Although the current model assumes classical local time derivatives to capture the short-term dynamics of the CO₂ laser system, future work could explore the generalization of this model using fractional calculus. Fractional-order derivatives can account for memory effects and nonlocal temporal dependencies, which may enhance the accuracy of the model when long-term or delayed feedback influences are present.