Higher-Order Interactions in the Network of Erbium-Doped Fiber Lasers

1. Introduction

Coupled systems are traditionally modeled as networks of elements, or nodes, connected in pairs by links or edges [1]. This pairwise paradigm has been remarkably successful in describing a wide range of dynamical phenomena. However, in many real-world systems, nodes interact in groups rather than exclusively in pairs. Considering only dyadic interactions fails to capture collective interactions among multiple elements simultaneously. To address this limitation, the framework of higher-order interactions (HOI) has been introduced [2,3,4,5,6].

The principal mathematical tools for representing HOI in complex networks are hypergraphs and simplicial complexes. In hypergraphs, a hyperedge generalizes the concept of an edge by connecting an arbitrary number of nodes. Simplicial complexes offer a more structured approach, where the topology is built from simplices of different dimensions (nodes, edges, triangles, tetrahedra, etc.), thereby providing a natural classification of first-order (pairwise), second-order (three-body), and HOI. Although the HOI framework has gained substantial traction in recent years across diverse fields, ranging from social systems and neuroscience to epidemic spreading and ecological dynamics [7,8,9,10,11], only a limited number of studies have incorporated HOI into laser models [6]. This scarcity is striking, as laser systems are known to exhibit complex dynamics and synchronization states that may be profoundly influenced by collective, multi-element couplings. Consequently, the current lack of research on HOI in optical systems severely restricts our understanding of the emergent collective dynamics and the range of possible behaviors that HOI may unlock.

The aim of this work is to investigate the dynamics and synchronization of a network of three erbium-doped fiber lasers (EDFLs) using a simplicial complex framework that explicitly accounts for both first-order (pairwise) and second-order (three-body) interactions. EDFLs hold significant technological importance, particularly in optical fiber communications, as they generate stable emission at the low-loss wavelength of $1.55\mu \mathrm{m}$ with high energy efficiency [12,13]. Since the EDFL is a highly complex multistable system, it is particularly convenient for investigating HOI. Nevertheless, the collective dynamics of EDFLs in higher-order coupling configurations remains entirely unexplored. While the individual behavior of EDFLs has been thoroughly characterized [14,15], and a recent study examined pairwise-coupled configurations [16], no work to date has addressed how second-order (three-body) interactions shape the dynamics and synchronization of EDFL networks.

Here, we consider a simplest nontrivial network of three EDFLs, modeled via linear diffusive couplings that enable both pairwise and collective (second-order) interactions. The remainder of this paper is organized as follows. Section 2 introduces the laser model and its normalized version used in numerical simulations. Section 3 presents the mathematical formalism of simplicial complexes and the generalized Laplacian matrix associated with the network topology. Using these tools, we formulate the network structure and define the metrics employed to characterize its dynamics and synchronization. Section 4 analyzes the obtained results and discusses the effects of introducing HOI on the laser dynamics and network synchronization. Finally, Section 5 summarizes the main findings, offers concluding remarks, and outlines directions for future research.

2. Numerical Model

2.1. Laser Model

The EDFL is modeled by two differential equations for the intracavity power P and the average population inversion I [14]. The model also accounts for relevant effects such as excited-state absorption (ESA) at the laser wavelength and pump depletion during propagation along the fiber:

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dP}{dt}}& =& {\displaystyle \frac{2}{L{T}_r}}P\left\{r_w{\alpha }_0\left[I(\xi -\eta )-1\right]-{\alpha }_{th}\right\}+P_{sp}\hfill \end{array}$$

(1)

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dI}{dt}}& =& {\displaystyle \frac{{\sigma }_{12}{r}_{w}P}{\pi r_0^2}}(I\xi -1)-{\displaystyle \frac{I}{\tau }}+P_{pump},\hfill \end{array}$$

(2)

where L is the length of the erbium-doped fiber and $T_r=2n_0(L+l_0)/c$ is the photon lifetime in the cavity ($l_0$ being the length of the Bragg grating couplers), $r_w=1+exp\left[2{(r_0/w_0)}^2\right]$ represents the overlap between the fundamental laser mode and the erbium-doped core of the active fiber, ${\alpha }_0=N_0{\sigma }_{12}$ is the small-signal absorption at the laser wavelength ($N_0=N_1+N_2$ being the total concentration of erbium ions), and ${\sigma }_{12}$ is the absorption cross section from the ground state 1 to the excited level 2.

We assume that the cross sections of the stimulated emission and return transition is the same (${\sigma }_{12}={\sigma }_{21}$), i.e., $\xi =({\sigma }_{12}+{\sigma }_{21})/{\sigma }_{12}=2$. The coefficient $\eta ={\sigma }_{23}/{\sigma }_{12}$ expresses the ratio between the ESA cross section and the fundamental absorption. The parameter ${\alpha }_{th}={\gamma }_0+(1/2L)ln(1/R)$ corresponds to the threshold loss inside the cavity, where ${\gamma }_0$ is the non-resonant fiber loss and R is the total reflection coefficient of the couplers.

The spontaneous emission in the fundamental laser mode is expressed as

$$P_{sp}={10}^{-3}{\displaystyle \frac{I}{\tau T_r}}{\displaystyle \frac{{\lambda }_g{w}_0^2}{{\left(r_0^2{\alpha }_{0}L\right)}^4}}{\displaystyle \frac{1}{4{\pi }^2{\sigma }_{12}}},$$

where ${\lambda }_g$ is the laser wavelength, $r_0$ is the radius of the erbium fiber core, and $w_0$ is the radius of the fundamental mode.

The average population inversion I is defined as

$$I={\displaystyle \frac{1}{n_{0}L}}{\int }_0^L{N}_2\left(z\right)dz,$$

where $N_2\left(z\right)$ represents the population of level 2 as a function of the position along the fiber, and $n_0$ is the “cold refractive index” of the doped core.

The pump power absorbed by the fiber is given as

$$P_{pump}={\displaystyle \frac{P_p\left[1-exp\left(-\beta {\alpha }_{0}L(1-I)\right)\right]}{N_0\pi r_0^{2}L}},$$

where $P_p$ is the input pump power to the erbium fiber and $\beta ={\alpha }_p/{\alpha }_0$ is a dimensionless coefficient representing the ratio between the absorption coefficients at the pump wavelength ${\lambda }_p$ (${\alpha }_p$) and the laser wavelength ${\lambda }_g$ (${\alpha }_0$).

If external modulation is included, the modulated pump power is modeled as

$$P_p^m=P_p\left[1+m_{0}sin\left(2\pi F_{m}t\right)\right],$$

(3)

where $m_0$ is the modulation depth and $F_m$ is the modulation frequency.

2.2. Normalized Model

To facilitate the numerical analysis of the erbium-doped laser system, the original Eqs. (1) and () are nondimensionalized, which yields the following equations:

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dx}{d\theta }}& =& axy-bx+c(y+e),\hfill \end{array}$$

(4)

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dy}{d\theta }}& =& -dxy-(y+e)+P_p^m\left\{1-exp\left[-f\left(1-{\displaystyle \frac{y+e}{g}}\right)\right]\right\}.\hfill \end{array}$$

(5)

The parameters used in the simulations are listed in Table 1 [16].

Subsequently, in order to match the numerical integration of the model with the sampling frequency of the real-time platform, an additional temporal rescaling is introduced.

Let $\theta$ be the original time variable and $\tau$ the scaled time variable related by

$$\theta =H\tau ,\mathrm{so}\mathrm{that}d\theta =Hd\tau ,$$

where H is a temporal scaling parameter that adjusts the time base of the simulated system. Applying the scaled time variable $t=H\tau$ to the original model yields the rescaled dimensionless system [17]:

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dx}{d\tau }}& =& H\left(axy-bx+c(y+e)\right),\hfill \end{array}$$

(6)

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dy}{d\tau }}& =& H\left(-dxy-(y+e)+P_p^m\left\{1-exp\left[-f\left(1-{\displaystyle \frac{y+e}{g}}\right)\right]\right\}\right),\hfill \end{array}$$

(7)

$$\begin{array}{ccc}\hfill P_p^m& =& P_p\left[1+m_{0}sin\left(H\left(2\pi F_{m}t\right)\right)\right].\hfill \end{array}$$

(8)

The introduction of the scaling factor H enables the temporal normalization of the model dynamics, ensuring that the numerical integration in MATLAB is performed with a step size compatible with the sampling frequency of the physical system, while preserving the nonlinear dynamics of the laser.

3. Higher-Order Interactions (HOI)

3.1. Definitions

A wide range of natural and engineered systems have been successfully modeled as networks of interacting units. Traditionally, such models have relied exclusively on pairwise (dyadic) interactions. However, this approach often fails to capture the full complexity of real-world systems, where collective interactions involving more than two nodes frequently occur and can fundamentally alter the emergent dynamics [18].

In recent years, HOI have gained substantial relevance in the study of complex systems, enabling the analysis of multi-element interactions in biological, neurological, physical, and social systems, among others [7,11]. Mathematical frameworks such as simplicial complexes and hypergraphs have revealed that system dynamics can change considerably when HOI are taken into account, compared to purely pairwise descriptions [6,7,8,9,10].

Among these frameworks, simplicial complexes provide a particularly structured and mathematically rigorous approach to modeling HOI. In essence, a network with higher-order couplings can be represented by a simplicial complex of order d, whose formal definition is given below [3,6].

The general dynamics of a node $X_i$ in a network with interactions up to order d is modeled as follows:

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{X}_i}{dt}}& =F\left(X_i\right)+{\sigma }_1{\displaystyle \sum _{j_1=1}^N}{a}_{i{j}_1}^{\left(1\right)}{H}^{\left(1\right)}(X_i,X_{j_1})+{\sigma }_2{\displaystyle \sum _{j_1=1}^N}{\displaystyle \sum _{j_2=1}^N}{a}_{i{j}_1{j}_2}^{\left(2\right)}{H}^{\left(2\right)}(X_i,X_{j_1},X_{j_2})\hfill \\ \hfill & +{\sigma }_d{\displaystyle \sum _{j_1=1}^N}\dots {\displaystyle \sum _{j_d=1}^N}{a}_{i{j}_1\dots j_d}^{\left(d\right)}{H}^{\left(d\right)}(X_i,X_{j_1},\dots ,X_{j_d}),\hfill \end{array}$$

(9)

where $F\left(X_i\right)$ represents the intrinsic dynamics of node i, ${\sigma }_d$ is the coupling strength of order d, $a_{i{j}_1\dots j_d}^{\left(d\right)}$ corresponds to the d-th order adjacency tensor, $H^{\left(d\right)}$ denotes the d-th order coupling function, and N is the total number of nodes in the network.

3.2. Simplicial Complexes

Simplicial complexes are topological structures composed of simplices of different dimensions (e.g., nodes, edges, triangles, tetrahedra), which capture many-body interactions among the elements of a system [3]. A simplex represents an interaction between two or more nodes, and its topological representation is shown in Figure 1, where an n-simplex consists of $n+1$ nodes.

A simplicial complex K is a finite collection of simplices $S\left(K\right)$ that satisfies the following properties:

1.

If $\alpha \in S\left(K\right)$, then every face of $\alpha$ also belongs to $S\left(K\right)$.

2.

If $\alpha ,\beta \in S\left(K\right)$ and $\beta \subseteq \alpha$, then $\beta$ is a face of $\alpha$.

Each n-simplex $\alpha \in S\left(K\right)$ is denoted as

$$\alpha =[v_0,v_1,\dots ,v_n],$$

where $v_0,v_1,\dots ,v_n$ are distinct vertices of K, and the ordering of the vertices determines an orientation.

The faces of a simplex $\alpha$ are all simplices formed by nonempty proper subsets of its vertices. In particular, the $(n-1)$-dimensional faces of an n-simplex are obtained by removing exactly one vertex:

$$\partial \alpha =\left\{[v_0,\dots ,\widehat{v_i},\dots ,v_n]\right|i=0,\dots ,n\},$$

where the symbol $\widehat{v_i}$ indicates that the vertex $v_i$ is omitted.

Figure 2. Decomposition of a simplicial complex of order 2 with its lower-order faces. If $\alpha =[v_0,v_1,v_2]$ is a 2-simplex (oriented triangle), its 0-dimensional faces (vertices) are $\left[v_0\right],\left[v_1\right],\left[v_2\right]$, its 1-dimensional faces (edges) are $[v_0,v_1],[v_1,v_2],[v_2,v_1]$, and its 2-dimensional face (triangle) is $[v_0,v_1,v_2]$.

Figure 2. Decomposition of a simplicial complex of order 2 with its lower-order faces. If $\alpha =[v_0,v_1,v_2]$ is a 2-simplex (oriented triangle), its 0-dimensional faces (vertices) are $\left[v_0\right],\left[v_1\right],\left[v_2\right]$, its 1-dimensional faces (edges) are $[v_0,v_1],[v_1,v_2],[v_2,v_1]$, and its 2-dimensional face (triangle) is $[v_0,v_1,v_2]$.

To algebraically capture the relationship between simplices of dimension n and $n-1$, the chain space $C_n\left(K\right)$ is introduced, defined as the vector space generated by all n-simplices of the simplicial complex K. The boundary operator is then defined as a linear transformation

$${\partial }_n:C_n\left(K\right)⟶C_{n-1}\left(K\right),$$

which acts on an oriented n-simplex $[v_0,\dots ,v_n]$ as

$${\partial }_n[v_0,\dots ,v_n]={\displaystyle \sum _{i=0}^n}{(-1)}^i[v_0,\dots ,\widehat{v_i},\dots ,v_n],$$

where the symbol $\widehat{v_i}$ indicates that the vertex $v_i$ is omitted.

By fixing an ordering for all simplices of dimension n and $n-1$, the boundary operator can be represented in matrix form through boundary matrix $B_n$. The columns of $B_n$ are indexed by the n-simplices of the complex, while the rows are indexed by the $(n-1)$-simplices. Thus, $B_n$ is a matrix of size $|{\mathcal{S}}_{n-1}|\times |{\mathcal{S}}_n|$, where ${\mathcal{S}}_k$ denotes the set of k-simplices of K.

The entries of the matrix $B_n$ are given by

$${\left(B_n\right)}_{\alpha ,\beta }=\left\{\begin{array}{cc}+1,\hfill & \mathrm{if}\mathrm{the}(n-1)-\mathrm{simplex}\alpha \mathrm{is}\mathrm{a}\mathrm{face}\mathrm{of}\beta \mathrm{and}\mathrm{the}\mathrm{orientation}\mathrm{coincides},\hfill \\ -1,\hfill & \mathrm{if}\alpha \mathrm{is}\mathrm{a}\mathrm{face}\mathrm{of}\beta \mathrm{but}\mathrm{with}\mathrm{opposite}\mathrm{orientation},\hfill \\ 0,\hfill & \mathrm{if}\alpha \mathrm{is}\mathrm{not}\mathrm{a}\mathrm{face}\mathrm{of}\mathrm{the}n-\mathrm{simplex}\beta .\hfill \end{array}\right.$$

In this way, each column of $B_n$ encodes how an n-simplex decomposes into its $(n-1)$-dimensional faces, while each row indicates in which higher-order simplices a given $(n-1)$-simplex participates. The matrix $B_n$ enables the computation of topological properties of the complex, such as the number of connected components or cavities [3].

3.3. Higher-Order Laplacian Matrix

From the boundary matrices it is possible to construct the simplicial Laplacian of order n, which generalizes the classical concept of the graph Laplacian matrix [3,5,18]. Let $B_n$ be the boundary matrix that relates simplices of dimension n with those of dimension $n-1$, and $B_{n+1}$ the matrix that relates simplices of dimension $n+1$ with those of dimension n. The Hodge Laplacian of order n is defined as:

$$L_n=B_n^{\top }{B}_n+B_{n+1}{B}_{n+1}^{\top }$$

(10)

where $B_n^{\top }{B}_n$ is known as the down Laplacian and captures the interactions induced by $(n-1)$-dimensional faces and $B_{n+1}{B}_{n+1}^{\top }$ is known as the up Laplacian and captures the interactions induced by $(n+1)$-dimensional simplices that contain each n-simplex.

From an algebraic viewpoint, the Hodge Laplacian $L_n$ acts on n-chains and combines the topological information coming from lower-dimensional faces and higher-dimensional simplices that contain each n-simplex.

Figure 3. Relationship between simplices of different dimensions through the boundary matrices. The matrix $B_1$ connects vertices with edges, while $B_2$ connects edges with triangles. Their transposes $B_1^{\top }$ and $B_2^{\top }$ allow reconstruction of the combinations of edges forming each triangle, or vertices forming each edge.

Figure 3. Relationship between simplices of different dimensions through the boundary matrices. The matrix $B_1$ connects vertices with edges, while $B_2$ connects edges with triangles. Their transposes $B_1^{\top }$ and $B_2^{\top }$ allow reconstruction of the combinations of edges forming each triangle, or vertices forming each edge.

In the case $n=0$, the term $B_0$ is not defined (since there are no simplices of dimension $-1$), therefore

$$L_0=B_1{B}_1^{\top },$$

(11)

that coincides with the Laplacian matrix of a traditional graph.

Thus, the simplicial Laplacian provides a unified tool for studying dynamics on nodes ($n=0$), edges ($n=1$), and higher-order simplices ($n\ge 2$), enabling the modeling of HOI such as those present in simplicial networks.

3.4. Network Construction

A three-node network with clockwise orientation is constructed (Figure 4), where each node represents an EDFL.

Linear diffusive coupling is used, implemented in the pump modulation of the system, yielding the following dynamical equations for each node:

$$\begin{array}{ccc}\hfill {\displaystyle \frac{d{x}_i}{d\tau }}& =& H\left(a{x}_i{y}_i-b{x}_i+c(y_i+e)\right)\hfill \end{array}$$

(12)

$$\begin{array}{ccc}\hfill {\displaystyle \frac{d{y}_i}{d\tau }}& =& H\left(-d{x}_i{y}_i-(y_i+e)+P_p^m\left\{1-exp\left[-f\left(1-{\displaystyle \frac{y_i+e}{g}}\right)\right]\right\}\right),\hfill \end{array}$$

(13)

where the modulated pump term is given as follows:

$$P_p^m=P_p\left[1+{\sigma }_1{\displaystyle \sum _{j=1}^N}{L}_0\left(x_j-x_i\right)+{\sigma }_2{\displaystyle \sum _{j=1}^N}{\displaystyle \sum _{k=1}^N}{L}_1\left(x_j+x_k-2x_i\right)\right].$$

(14)

Here, ${\sigma }_1$ and ${\sigma }_2$ are the first- and second-order coupling strength parameters, respectively; $L_0$ and $L_1$ correspond to the higher-order Laplacian matrices [19]; N is the number of nodes in the network; and the summation terms represent the corresponding first- and second-order diffusive coupling contributions.

From this network, the boundary matrices capturing the connections between nodes (1-simplices) in Table 2 and HOI (2-simplices) in Table 3 are constructed.

Using these matrices, the higher-order Laplacians are obtained:

$$L_0=\left(\begin{array}{ccc}2& -1& -1\\ -1& 2& -1\\ -1& -1& 2\end{array}\right)$$

$$L_1=\left(\begin{array}{ccc}3& 0& 0\\ 0& 3& 0\\ 0& 0& 3\end{array}\right)$$

Here, the matrix $L_0$ corresponds to the first-order Laplacian, and $L_1$ represents the second-order Laplacian matrix. With this formulation, the model is completely defined in terms of the dynamics of each laser and the first- and second-order interactions that modulate the pump.

3.5. Dynamical and Synchronization Measures

To study the three-node EDFL network with linear diffusive coupling, the system of differential equations is solved numerically using a fourth-order Runge–Kutta (RK4) scheme with a fixed step size of $h=0.1$ and a total of $2^{20}$ time steps. Analysis of the resulting time series yields the following dynamical and synchronization measures:

• Bifurcation diagrams. For each combination of the coupling parameters $({\sigma }_1,{\sigma }_2)$, the local maxima of the laser intensity are extracted from the time series. The last 35 maxima are retained to ensure that transient behavior has been discarded. These maxima are then plotted as a function of the coupling parameter to characterize transitions between periodic and chaotic dynamics.
• Lyapunov exponents. The maximum Lyapunov exponent is computed for each pair $({\sigma }_1,{\sigma }_2)$ using the Wolf algorithm [20]. The calculation is performed on the numerically generated time series after discarding initial transients, considering only the asymptotic regime. The resulting exponent is plotted as a function of the coupling parameter to identify regions of stability ($\lambda \le 0$) and chaos ($\lambda >0$), thereby providing a quantitative confirmation of the dynamical transitions observed in the bifurcation diagrams.
• Fourier spectrum. To characterize the frequency content of the numerically obtained time series, the Fourier spectrum of the laser intensity is computed in the asymptotic regime. The discrete Fourier transform (DFT) is employed, defined as:

  $$X\left(k\right)={\displaystyle \sum _{n=0}^{N-1}}x\left(n\right)e^{-i2\pi kn/N},$$

  (15)

  where $x\left(n\right)$ represents the laser intensity time series, N is the total number of sampled points, and k is the frequency index. The power spectrum is obtained from the modulus of the transform, enabling identification of the dominant frequency components of the system. A finite number of well-defined peaks indicates periodic behavior, whereas a broad or continuous spectrum with multiple components is typically associated with quasiperiodic or chaotic dynamics. This analysis complements the information obtained from bifurcation diagrams and time series [21].
• Average synchronization error. The instantaneous deviation between two lasers i and j in the $(x,y)$ phase space is defined as [22]:

  $$|e_{ij}\left(t\right)|=\sqrt{{(x_j\left(t\right)-x_i\left(t\right))}^2+{(y_j\left(t\right)-y_i\left(t\right))}^2},$$

  (16)

  where x and y denote appropriate state variables (e.g., intracavity power and population inversion). From this instantaneous measure, the average synchronization error over all laser pairs is computed by averaging over time and over all unique pairs:

  $$\langle E\rangle ={\displaystyle \frac{1}{T}}{\int }_0^T{\displaystyle \frac{1}{M}}\sum _{i<j}\left|e_{ij}\left(t\right)\right|dt,$$

  (17)

  where $M=\left({\displaystyle \genfrac{}{}{0pt}{}{N_{\mathrm{lasers}}}{2}}\right)$ is the total number of distinct laser pairs (for three lasers, $M=3$). This metric quantifies the global synchronization level of the network over time, with $\langle E\rangle =0$ indicating perfect synchronization.
• Average phase difference. The instantaneous phase of each laser is obtained as [22]:

  $${\varphi }_i\left(t\right)=arctan\left({\displaystyle \frac{y_i\left(t\right)-y_{i0}}{x_i\left(t\right)-x_{i0}}}\right),$$

  (18)

  where $(x_{i0},y_{i0})$ denotes a reference center of the trajectory (e.g., the centroid of the attractor or an equilibrium point). The phase difference between lasers i and j is then computed as:

  $$\Delta {\varphi }_{ij}\left(t\right)={\varphi }_j\left(t\right)-{\varphi }_i\left(t\right).$$

  (19)
  This measure allows the identification of in-phase ($\Delta {\varphi }_{ij}\approx 0$), anti-phase ($\Delta {\varphi }_{ij}\approx \pi$), or drifting phase relationships.
• Normalized cross-correlation. The normalized cross-correlation between lasers i and j is defined as [22]:

  $$C_{ji}\left(\tau \right)={\displaystyle \frac{〈\left(x_j\left(t\right)-{\overline{x}}_j\right)\left(x_i(t+\tau )-{\overline{x}}_i\right)〉}{\sqrt{〈{\left(x_j\left(t\right)-{\overline{x}}_j\right)}^2〉〈{\left(x_i(t+\tau )-{\overline{x}}_i\right)}^2〉}}},$$

  (20)

  where ${\overline{x}}_i$ and ${\overline{x}}_j$ denote the temporal averages of the signals, and $\langle ·\rangle$ represents time averaging. The maximum value of $C_{ji}\left(\tau \right)$ indicates the strongest correlation between the two signals, while the lag ${\tau }_{max}$ at which this maximum occurs quantifies the time delay that maximizes their similarity.
• Similarity measure. This quantity evaluates the structural correspondence between two time series, even in the presence of delays, and is defined as [22]:

  $$S_{ji}\left(\tau \right)={\displaystyle \frac{〈{\left[x_j(t+\tau )-x_i\left(t\right)\right]}^2〉}{\sqrt{〈x_i{\left(t\right)}^2〉+〈x_j{\left(t\right)}^2〉}}},$$

  (21)

  where the numerator quantifies the mean squared difference between the delayed signal $x_j$ and the reference signal $x_i$, and the denominator normalizes by the root-sum-square of the signal powers. Lower values of $S_{ji}\left(\tau \right)$ indicate greater similarity between the two time series.
• Global synchronization error. This measure quantifies the average pairwise deviation among all laser states across the entire system [23]:

  $$e\left(t\right)=\sqrt{{\displaystyle \frac{1}{N(N-1)}}\sum _{i=1}^N\sum _{j\ne i}^N{\parallel {\mathbf{x}}_i\left(t\right)-{\mathbf{x}}_j\left(t\right)\parallel }^2},$$

  (22)

  where ${\mathbf{x}}_i\left(t\right)=(x_i\left(t\right),y_i\left(t\right))$ represents the state vector of laser i at time t, and $\parallel ·\parallel$ denotes the Euclidean norm. This metric approaches zero when the lasers are highly synchronized, providing a scalar measure of global coherence.
• Kuramoto order parameter. This quantity describes the degree of phase synchronization among the lasers [23]:

  $$r\left(t\right)={\displaystyle \frac{1}{N}}\left|{\displaystyle \sum _{i=1}^N}{e}^{j{\varphi }_i\left(t\right)}\right|,$$

  (23)

  where ${\varphi }_i\left(t\right)$ is the instantaneous phase of laser i, $j=\sqrt{-1}$ (the imaginary unit), and $N=3$ is the total number of lasers in the network. The order parameter satisfies $0\le r\left(t\right)\le 1$, with $r\left(t\right)\approx 1$ indicating high phase synchronization (all phases aligned) and $r\left(t\right)\approx 0$ reflecting complete phase desynchronization (phases uniformly distributed on the circle).

4. Results

This section presents the dynamical behavior of the three-node EDFL network under different coupling configurations. We systematically investigate the influence of HOI by comparing four distinct cases, which differ in the presence and relative strength of first-order (pairwise) and second-order (three-body) couplings. Unless otherwise stated, all simulations are performed using the same integration protocol and parameter ranges to ensure fair comparison across cases. The specific coupling scenarios are described below.

4.1. Parameter Cases Studied

• Case 1: Simplicial complex of order 1 (pairwise coupling only). The first-order coupling strength ${\sigma }_1$ is varied in the range $0\le {\sigma }_1\le 4$, while the second-order coupling is set to ${\sigma }_2=0$. This case serves as a reference baseline, representing dynamics arising from purely pairwise interactions.
• Case 2: Second-order coupling only (${\sigma }_1=0$, varying ${\sigma }_2$). Here, ${\sigma }_2$ is varied while ${\sigma }_1=0$. Although this configuration does not formally constitute a simplicial complex (as lower-order faces are absent), it is examined to assess the dynamical effects of second-order interactions in isolation.
• Case 3: Simplicial complex of order 2 with fixed ${\sigma }_2$. The first-order coupling strength ${\sigma }_1$ is varied in the range $0\le {\sigma }_1\le 4$, while the second-order coupling is held constant at ${\sigma }_2=1.5$. This configuration incorporates both pairwise and three-body interactions. The value ${\sigma }_2=1.5$ was selected after preliminary simulations across a range of ${\sigma }_2$ values, as it yields the richest dynamical behavior, exhibiting multiple transitions.
• Case 4: Simplicial complex of order 2 with fixed ${\sigma }_1$. The second-order coupling strength ${\sigma }_2$ is varied in the range $0\le {\sigma }_2\le 4$, while the first-order coupling is held constant at ${\sigma }_1=1.5$. This value is chosen for consistency with Case 3, facilitating direct comparison between the two scenarios.

For Cases 1, 2, and 3, the last 300,000 data points of each time series are analyzed to ensure that transient dynamics were excluded. For Case 4, the last 600,000 data points are used to achieve clearer visualization of the observed behaviors. These datasets serve as the basis for computing global and pairwise synchronization metrics, phase differences, and the Kuramoto order parameter.

Bifurcation diagrams are constructed from the last 35 local maxima of the intensity time series for each laser. For the similarity and cross-correlation analyses, the interval from 600,000 to 800,000 time steps is selected. This choice avoids artifacts arising from time shifts introduced during the evaluation, which preclude the direct use of the final portion of the time series.

Each case is evaluated by considering the different faces of the simplicial complex, with analysis tools tailored to each interaction level:

• Face 0 (0-simplices: individual nodes). For the dynamics of isolated nodes (Face 0), we analyze bifurcation diagrams, maximum Lyapunov exponents, time series, phase space reconstructions, and Fourier spectra. These tools characterize the local dynamical behavior of each laser under the influence of network coupling.
• Face 1 (1-simplices: pairwise edges). For each edge connecting a pair of nodes, we evaluate the average synchronization error, phase difference, similarity measure, and normalized cross-correlation. These metrics quantify the degree and quality of pairwise synchronization between lasers.
• Face 2 (2-simplices: triangles, corresponding to second-order interactions). For the three-node triangle (Face 2), we study the global synchronization of the entire system using the global synchronization error and the Kuramoto order parameter. These measures capture collective behavior emerging from higher-order interactions that cannot be reduced to pairwise descriptions.

This hierarchical approach associates each level of interaction with the most appropriate analytical tools, enabling a comprehensive characterization of the complex dynamics arising from higher-order couplings in the EDFL network.

4.1.1. Face 0: Individual Node Dynamics

To characterize the local dynamics of individual nodes (Face 0 of the simplicial complex), bifurcation diagrams of the laser intensity for each laser are constructed. Figure 5 shows these results for laser $x_1$, where the last 35 local maxima of the time series are plotted after discarding the transient regime.

For very small coupling values, specifically in the range $0\le {\sigma }_1\le 0.007$, the lasers generate a stable emission. As the coupling strength is increased, the laser intensity becomes oscillate periodically in the supercritical Hopf bifurcation with the interval $0.007\le {\sigma }_1\le 0.4$. Then, in the torus bifurcation at ${\sigma }_1=0.4$ the oscillations become quasiperiodic, and at $0.6\le {\sigma }_1\le 0.718$ a well-defined periodic window appears, showing multiple intensity peaks within each cycle with a phase shift between the three lasers. Finally, for ${\sigma }_1>0.718$, the oscillations are chaotic. Notably, in the subinterval $1.712\le {\sigma }_1\le 1.812$, a significant bistable window is observed, where the system alternates between chaotic and periodic behaviors.

To corroborate the bifurcation diagram results, the maximum Lyapunov exponent was computed from the time series using the Wolf algorithm [20]. For $0\le {\sigma }_1\le 0.718$, the exponent is close to zero, confirming the existence of a stable fixed point, limit cycle, or quasiperiodicity. As the coupling strength is further increased, the exponent becomes positive, indicating a chaotic regime. Within this region, the bistable zone observed in the bifurcation diagram ($1.712\le {\sigma }_1\le 1.812$) is reflected in the Lyapunov exponent, which oscillates between values near zero (periodic) and approximately $0.75$ (chaotic), confirming the coexistence of attractors.

To further characterize the dynamical regimes observed in the bifurcation diagram, Figure 6 illustrates the most representative behaviors, showing time series (left column), phase space portraits (middle column), and power spectra (right column) for selected values of ${\sigma }_1$.

Periodic regime (${\sigma }_1=0.3$).Figure 6(a–c) correspond to ${\sigma }_1=0.3$. The time series in Figure 6(a) shows periodic oscillations with constant amplitude. The corresponding phase-space trajectory, presented in Figure 6(b), is a stable limit cycle. The Fourier spectrum in Figure 6(c) reveals a pulse repetition frequency $f_1$ and its higher harmonics at $n{f}_1$ ($n=1,2,\dots$), a behavior typical of weakly coupled laser systems.

Quasiperiodic regime (${\sigma }_1=0.45$). Increasing the coupling strength induces a torus bifurcation, after which the motion becomes quasiperiodic (Figure 6(d–f)). As seen in the time series (Figure 6(d)), the maximum laser intensity reaches values similar to those observed under weaker coupling; however, the local maxima now exhibit varying amplitudes. A constant phase shift among the oscillations of the three lasers gives rise to a rotating phase wave along the laser ring [16]. The phase-space trajectory assumes a toroidal shape (Figure 6(e)), characteristic of quasiperiodic dynamics driven by two incommensurate frequencies ($f_1=2.3$ kHz and $f_2=10.0$ kHz), as shown in the power spectrum (Figure 6(f)). The interaction between these frequencies generates idler frequencies and linear combinations of the form $m{f}_1\pm n{f}_2$ (with integers $m,n$), a hallmark of two-frequency quasiperiodicity.

Complex periodic regime (${\sigma }_1=0.66$). For ${\sigma }_1=0.66$ (Figure 6(g–i)), the system exhibits more complex periodic behavior. The time series (Figure 6(g)) reveals multiple intensity peaks per cycle (period 3), suggesting amplitude modulation or the emergence of a secondary bifurcation of the limit cycle. Despite this increased complexity, overall periodicity is preserved. The phase-space trajectory in Figure 6(h) remains a closed loop but becomes wider, reflecting an increase in oscillation amplitude and, consequently, in the intracavity energy of the laser system. The power spectrum in Figure 6(i) exhibits well-defined discrete peaks at subharmonics and multiple harmonics of the fundamental frequency. Nevertheless, a noticeable broadening of the harmonics is also observed, indicating a partial loss of phase coherence at higher frequencies — a phenomenon commonly associated with routes to chaos in coupled laser systems.

Bistable regime (${\sigma }_1=1.75$). Finally, at ${\sigma }_1=1.75$ (Figure 6(j–l)), the system displays bistability, wherein the asymptotic dynamics depends sensitively on the initial conditions. The superimposed time series in Figure 6(j) and the phase-space trajectories in Figure 6(k) reveal the coexistence of periodic (green) and chaotic (black) regimes. The power spectrum reflects both behaviors: a period-3 state yields well-defined discrete peaks (harmonics and subharmonics), whereas the chaotic state produces peak broadening — a typical signature of chaotic dynamics in laser systems [24,25]. The observed subharmonics coincide with idler frequencies, namely $2f_1/3=f_1-f_1/3$ and $4f_1/3=f_1+f_1/3$.

4.1.2. Face 1: Pairwise Synchronization Analysis

Let us now examine the Face 1 of the simplicial complex, corresponding to pairwise interactions along each link of the network. Synchronization is characterized using several complementary measures.

Average synchronization error.Figure 7(a) shows the average synchronization error computed via Equation (17) (see also [22]), which quantifies the mean squared difference between the intensities of paired lasers.

For very small values of ${\sigma }_1$, the error remains near zero. As the coupling strength is increased, the error grows, reaching a maximum at ${\sigma }_1\approx 2.25$, a value that lies within the regime of periodic single-peak oscillations. Beyond this point, small fluctuations appear, reflecting complex dynamical transitions. As ${\sigma }_1$ continues to increase, the error gradually decreases, approaching ${\sigma }_1\approx 1.5$, though it never reaches zero, indicating the absence of complete synchronization.

Maximum synchronization error.Figure 7(b) presents the maximum synchronization error for each value of ${\sigma }_1$. In contrast to the average error, the maximum error exhibits an approximately exponential increase before saturating at a nearly constant value. To better understand this behavior, we examine the time series of the synchronization error between two specific lasers (e.g., $E_{12}$) at two representative coupling strengths. At ${\sigma }_1=0.5$ (Figure 7(c)), the error shows intermittent bursts. At ${\sigma }_1=3.5$ (Figure 7(d)), the intervals during which the system remains synchronized ($E=0$) become longer and more frequent. However, during brief episodes of desynchronization, the error spikes sharply, accounting for the large maximum error values observed in Figure 7(b).

Similarity and correlation measures.Figure 8(a) and Figure 8(b) present the similarity measure $S_{ij}$ and the normalized cross-correlation $C_{ij}$, respectively. To ensure coherent comparison across different dynamical regimes, both metrics are computed over time intervals corresponding to one complete oscillation cycle of the system for each value of ${\sigma }_1$, i.e., using the number of points spanning a full waveform.

In Figure 8(a), the similarity measures for the three node pairs ($S_{12}$, $S_{23}$, and $S_{31}$) are shown. In parameter intervals corresponding to periodic dynamics, $S\approx 0$, indicating that the signals are identical (up to a possible time shift). In the first bifurcation region and within the quadiperiodicity window, the similarity gradually increases, reflecting the growing phase shift between lasers. In the chaotic region, the similarity values form a scattered cloud without apparent structure. Only a few isolated points approach zero, corresponding to narrow periodic windows that briefly interrupt the chaotic behavior.

Figure 8(b) presents the cross-correlation for the three node pairs ($C_{12}$, $C_{23}$, and $C_{31}$), which exhibits qualitatively similar behavior. Within periodic windows, $C\approx 1$, indicating strong correlation. In the first bifurcation region and the quasiperiodicity window, C gradually decreases. In the chaotic region, the correlation values are scattered approximately between $0.2$ and $0.7$, though several points remain close to unity.

The inverse relationship between similarity and correlation confirms that even when the system does not achieve complete synchronization, episodes of delayed or advanced synchronization occur. Specifically, during periodic intervals, the signals maintain structural coherence with a temporal displacement. Similar behavior has been reported in previous studies, where partial synchronization emerges in relation to the structure of a simplicial complex [3,5,8].

Phase difference analysis.Figure 9(a) shows the average phase difference between node pairs. For very small coupling values, the phase difference is large because the interaction between lasers is weak. As the coupling strength is slightly increased, the phase difference decays, exhibiting a small plateau around ${\sigma }_1\approx 0.2$–$0.3$. In this region, the interaction is sufficiently strong to reduce the phase difference but not strong enough to achieve complete synchronization. With further increases in coupling, the average phase difference decreases with small fluctuations, settling near approximately $0.01$, but never reaching zero, indicating partial phase synchronization.

To gain deeper insight, three representative coupling values (low, intermediate, and high) are selected to visualize the time series of the phase difference between laser 1 and laser 2. At ${\sigma }_1=0.3$ (Figure 9(b)), the phase difference exhibits rapid and frequent fluctuations, with brief intervals where it momentarily reaches zero. Increasing the coupling to ${\sigma }_1=1.7$ (Figure 9(c)), the intervals during which ${\varphi }_{12}=0$ become longer, indicating more stable synchronization. At the highest coupling shown, ${\sigma }_1=3.8$ (Figure 9(d)), these synchronized intervals extend further. Across all three cases, the phase difference remains bounded by a maximum value of approximately $1.6$ radians. As the coupling strength is increased, the fraction of time during which the lasers remain synchronized also increases.

4.1.3. Face 2: Global Synchronization Analysis

Concluding the analysis of case 1, we now examine Face 2 of the simplicial complex, which captures the collective dynamics of the complete three-node network. Figure 10 presents two complementary metrics for global synchronization.

Global synchronization error.Figure 10(a) shows the global average synchronization error E as a function of ${\sigma }_1$. For very small coupling strengths, $E\approx 0$, indicating near-perfect synchronization. As ${\sigma }_1$ increases to approximately $0.1$, the error rises sharply to a value of about $2.3$, after which it remains relatively stable until ${\sigma }_1\approx 0.5$, where perturbations begin to appear. As ${\sigma }_1$ increases further beyond $2.0$, the error gradually decreases, reaching approximately $1.6$ at ${\sigma }_1=4.0$. Notably, the lasers do not achieve complete synchronization ($E=0$) for any ${\sigma }_1>0.1$. This behavior is qualitatively similar to that observed for pairwise synchronization (Face 1), although the global error values are slightly higher (compare with Figure 7(b)). Nevertheless, the overall consistency between Face 1 and Face 2 metrics suggests that the network responds uniformly to changes in the coupling strength.

Kuramoto order parameter.Figure 10(b) presents the Kuramoto order parameter R, which quantifies the degree of phase synchronization among the three lasers. For weak coupling (${\sigma }_1$ near zero), $R\approx 0.84$, indicating moderate phase coherence. As ${\sigma }_1$ increases, R rises sharply to approximately $1.0$, signifying strong phase synchronization. Subsequently, R decreases slightly to about $0.94$, then increases again to approximately $0.98$, before stabilizing around $R\approx 0.94$ for larger ${\sigma }_1$. At the highest coupling values examined, R approaches $1.0$, indicating robust phase synchronization across the network.

The sustained high values of R (typically $R>0.9$) demonstrate that the network remains phase-synchronized over most of the parameter range, even though complete synchronization (i.e., identical trajectories) is not achieved. This indicates that the phase differences and individual trajectories evolve coherently in phase space. This is a signature of phase synchronization in coupled oscillatory systems [22].

4.2. Case 2: Second-Order Coupling Only (${\sigma }_1=0$, varying ${\sigma }_2$)

We now examine the dynamical effects of second-order coupling in isolation by setting the first-order coupling strength to zero (${\sigma }_1=0$) and varying ${\sigma }_2$ in the range $0\le {\sigma }_2\le 100$. Strictly speaking, this configuration does not constitute a complete simplicial complex, as the lower-order faces (pairwise interactions) are absent. Nevertheless, studying this scenario allows us to isolate the influence of second-order interactions on the laser dynamics.

For coupling strengths in the range $0<{\sigma }_2<40$, the system exhibits no oscillations; the dynamics remains at a stable fixed point. As ${\sigma }_2$ is increased beyond 40, a stable periodic orbit emerges, as shown in the bifurcation diagram (Figure 11(a)). This limit cycle coexists with the fixed point. Representative time series and phase portrait at ${\sigma }_2=47$ are presented in Figure 11(b) and (c), respectively. Notably, for $47<{\sigma }_2<100$, the system exhibits bistability, characterized by the coexistence of two distinct attractors: the original stable fixed point and a stable limit cycle. In this regime, the asymptotic dynamics depends sensitively on the initial conditions.

The observed dynamics is not determined solely by the coupling strength ${\sigma }_2$; rather, it also depends on the topological structure defined by the connections among simplices of different dimensions. The bistability observed here can be interpreted as arising from the fact that increasing ${\sigma }_2$ acts as a strong nonlinear perturbation rather than as a canonical representation of HOI as typically formulated in the theoretical framework [3,5]. In other words, when first-order couplings are absent, the second-order term does not simply refine or extend pairwise dynamics but instead introduces fundamentally new dynamical features, including multistability.

4.3. Case 3: Second-Order Simplicial Complex with Fixed ${\sigma }_2$

We now investigate the influence of HOI by fixing the second-order coupling strength at ${\sigma }_2=1.5$ while varying the first-order coupling in the range $0\le {\sigma }_1\le 4$. As in Case 1, we analyze the dynamics hierarchically: Face 0 (individual node dynamics), Face 1 (pairwise synchronization), and Face 2 (global synchronization).

4.3.1. Face 0: Individual Node Dynamics

We begin with Face 0 of Case 3, which characterizes the local dynamics of individual lasers under the combined influence of first- and second-order couplings.

Figure 12(a) presents the bifurcation diagram of the peak intensity of laser 1. For $0\le {\sigma }_1\le 0.02$, the system converges to a stable fixed point. At ${\sigma }_1\approx 0.02$, a supercritical Hopf bifurcation occurs, giving rise to periodic oscillations with amplitude that increases with ${\sigma }_1$ up to ${\sigma }_1\approx 0.71$, where a second Hopf bifurcation leads to quasiperiodic dynamics.

In the interval $0.711\le {\sigma }_1\le 1.21$, the system exhibits quasiperiodic behavior. For stronger coupling, $1.211\le {\sigma }_1\le 1.26$, a narrow periodic window appears, characterized by the coexistence of several periodic attractors (multistability). In the range $1.261\le {\sigma }_1\le 1.559$, quasiperiodicity reemerges, though now with irregular oscillations, as discussed below.

As ${\sigma }_1$ is further increased, a wide periodic window appears in the interval $1.6\le {\sigma }_1\le 1.76$, again exhibiting multistability. This is followed by a chaotic region for $1.761\le {\sigma }_1\le 2.903$. Subsequently, a third periodic window emerges in the range $2.904\le {\sigma }_1\le 3.554$, also displaying multistability, before the system finally transitions into a chaotic regime for larger ${\sigma }_1$.

To complement the bifurcation analysis, the maximum Lyapunov exponent is computed using the Wolf algorithm [20], as shown in Figure 12(b). For coupling strengths in the ranges $0\le {\sigma }_1\le 1.761$ and $2.904\le {\sigma }_1\le 3.554$, the Lyapunov exponent remains near zero, exhibiting only small fluctuations, consistent with periodic or quasiperiodic dynamics. In contrast, within the parameter intervals identified as chaotic in the bifurcation diagram, the Lyapunov exponent takes positive values, confirming strong sensitivity to initial conditions.

Figure 13 shows four representative dynamical regimes observed in Case 3, showing time series (left column), phase space portraits (middle column), and Fourier spectra (right column).

Figure 13(a–c) correspond to ${\sigma }_1=1.1$. The time series shows quasiperiodic oscillations with a periodic envelope, characteristic of a higher-dimensional torus in phase space. The phase-space trajectory forms a well-defined toroidal structure, while the Fourier spectrum reveals three incommensurate fundamental frequencies: $f_1=10.87$ kHz, $f_2=2.18$ kHz, and $f_3=34.79$ kHz. Additional peaks appear at linear combinations of the form $n{f}_1\pm m{f}_2\pm p{f}_3$, a hallmark of quasiperiodicity.

At ${\sigma }_1=1.4$ (Figure 13(d–f)), the system remains quasiperiodic, with the envelope exhibiting irregular oscillations. This behavior suggests a loss of temporal coherence arising from nonlinear interactions among the internal modes of the network. The power spectrum shows dominant peaks near multiples of a fundamental frequency ($f_1\approx 10.376$ kHz) and its higher harmonic. The additional small peak appears at $f_2\approx 2.0$ kHz, resulting a low-frequency envelope in the time series. Since these frequencies do not yet satisfy an exact harmonic relationship, the system operates near a phase-locking regime. In this regime, frequencies tend to synchronize but retain small differences associated with residual quasiperiodic dynamics [5].

For ${\sigma }_1=1.66$ (Figure 13(g–i)), the system exhibits multistability, where the asymptotic dynamics depends on the initial conditions. The phase space reveals three coexisting limit cycles, each with distinct amplitudes. The power spectra for these coexisting attractors shares the same dominant fundamental frequency but differs in their spectral power distributions, reflecting the different amplitudes of the oscillations.

Finally, at ${\sigma }_1=2.5$ (Figure 13(j–l)), the system operates in a chaotic regime. The time series exhibits significant variations in peak intensity, while the phase-space trajectory displays a complex, non-repeating structure. The power spectrum exhibits a spectral broadening, characteristic signatures of chaotic dynamics in laser systems.

4.3.2. Face 1: Pairwise Synchronization Analysis

We now examine Face 1 of the simplicial complex for Case 3, characterizing pairwise synchronization along each network link using multiple complementary metrics.

Average synchronization error.Figure 14(a) shows the average synchronization error E as a function of ${\sigma }_1$.

Similar to Case 1, the error increases from $E=0$ to $E\approx 2.25$ as the coupling strength is increased from ${\sigma }_1=0$ to ${\sigma }_1\approx 0.1$. For $0.1\lesssim {\sigma }_1\lesssim 1.0$, the error remains approximately constant, and the lasers oscillate in a periodic regime. As the coupling strength is further increased, the error undergoes small fluctuations and gradually decreases to $E\approx 1.6$ because the dynamics becomes more complex. As in Case 1, the error never reaches zero, meaning that complete synchronization is never achieved.

Maximum error.Figure 14(b) shows that the maximum synchronization error increases monotonically with ${\sigma }_1$. The temporal evolution of the error is presented in Figure 14(c) and Figure 14(d) for ${\sigma }_1=0.5$ and ${\sigma }_1=3.5$, respectively. The error oscillates periodically at a frequency corresponding to the laser power oscillations. As the coupling strength is increased, the frequency of these error oscillations decreases. Consequently, although the peak error amplitude increases, the average error decreases because the intervals during which $E=0$ become longer. In other words, increasing ${\sigma }_1$ improves the average synchronization, even though brief episodes of desynchronization persist. A similar coexistence of synchronization and desynchronization in oscillatory systems with higher-order interactions has been reported by Majhi et al. [7].

Similarity and correlation measures.Figure 15(a) and Figure 15(b) present the similarity measure S and the normalized cross-correlation C for the three laser pairs $(1,2)$, $(2,3)$, and $(3,1)$ as functions of ${\sigma }_1$.

For small coupling values (${\sigma }_1<0.75$), the similarity is close to $S\approx 0$ and the correlation is close to $C\approx 1$. This indicates that the laser dynamics are nearly identical. Under weak coupling, the lasers oscillate periodically and become synchronized with a slight time delay [22,23].

As ${\sigma }_1$ is increased, S increases while C decreases, reflecting a progressive loss of synchronization. The abrupt rise and fall of S (and the corresponding inverse behavior in C) are associated with the emergence of quasiperiodic dynamics. In the interval $1.5<{\sigma }_1<1.8$, three distinct branches coexist, indicating multistability in which three limit cycles coexist.

For stronger coupling, in the chaotic region, we observe a total dispersion of the measures: S becomes high and C becomes low, indicating a significant loss of similarity between the laser dynamics. Within the chaotic region, a periodic window appears in the range $2.7<{\sigma }_1<3.2$, where five coexisting attractors are observed. In this region, the structural relationships between the signals exhibit partial coherence [7,9].

Phase difference analysis.Figure 16(a) shows the average phase difference as a function of ${\sigma }_1$ for the three laser pairs $(1,2)$, $(2,3)$, and $(3,1)$. The three curves are nearly identical, indicating symmetric coupling among the lasers. As ${\sigma }_1$ increases, the average phase difference decreases, approaching $\varphi \approx 0.02$, which indicates that the lasers become almost in-phase synchronized. Qualitatively similar behavior has been observed in other coupled systems with HOI [6,9,10].

The time series of the phase difference ${\varphi }_{12}$ are displayed in Figure 16(b–d) for three representative coupling strengths: ${\sigma }_1=0.3$ (weak coupling), ${\sigma }_1=1.7$ (intermediate coupling), and ${\sigma }_1=3.8$ (strong coupling). In all cases, the phase difference represents periodic spikes, i.e., the intervals of a synchronous dynamics ${\varphi }_{12}=0$ interrupted by asynchronous dynamics (${\varphi }_{12}>0$, with a frequency corresponding to the laser oscillations.

The time series of the phase difference at weak coupling (${\sigma }_1=0.3$) is shown in Figure 16(b). It represents periodic spikes of the same amplitude ${\varphi }_{12}\approx 1.57$ rad ($\pi /2$), indicating a consistent phase shift between the lasers during desynchronized episodes.

At intermediate coupling (${\sigma }_1=1.7$), the time series in Figure 16(c) shows that the phase difference represents period-doubling oscillations in the phase difference envelope.

At strong coupling (${\sigma }_1=3.8$), the time series in Figure 16(d) represents the irregularity in the peak phase difference, indicating a chaotic modulation of the phase difference envelope.

Notably, as ${\sigma }_1$ is increased, the intervals between successive phase difference switches become longer, meaning that the lasers remain synchronized for extended periods. Consequently, the average phase difference decreases with increasing coupling strength, consistent with the trend observed in Figure 16(a).

4.3.3. Face 2: Global Synchronization Analysis

Let us now present the results of the global analysis corresponding to Face 2 of the simplicial complex, namely, the global average synchronization error and the Kuramoto order parameter.

Global synchronization error.Figure 17(a) shows the global synchronization error averaged over all coupling values. As in the previous cases, for very weak coupling (${\sigma }_1<0.2$), the error increases with ${\sigma }_1$, indicating that the lasers become desynchronized. As the coupling strength is further increased up to ${\sigma }_1\approx 1.0$, the error remains high and relatively stable. For ${\sigma }_1>1.0$, the error begins to oscillate and gradually decreases; however, it never reaches zero. Thus, complete synchronization is never achieved in this configuration.

Kuramoto order parameter. In contrast, as shown in Figure 17(b), the Kuramoto order parameter R first decreases from $R=1.0$ to approximately $R\approx 0.84$ as ${\sigma }_1$ increases from 0 to ${\sigma }_1\approx 0.1$, indicating partial desynchronization. As the coupling strength is increased further to ${\sigma }_1\approx 0.2$, the order parameter rises to $R\approx 0.92$, a value it maintains up to ${\sigma }_1\approx 1.0$. Beyond this point, R increases further and eventually saturates at $R=1.0$, indicating that the three lasers exhibit in-phase synchronization [23].

The inverse relationship between the pairwise phase difference (Figure 16(a)) and the order parameter confirms that as the phase discrepancy between individual laser pairs decreases, the global phase coherence of the entire network increases.

4.4. Case 4: Second-Order Simplicial Complex with Fixed ${\sigma }_1$

Let us now investigate the influence of HOI by fixing the first-order coupling strength at ${\sigma }_1=1.5$ while varying the second-order coupling in the range $0\le {\sigma }_2\le 4$. This represents the complementary scenario to Case 3. Similar to the previous cases, we analyze the dynamics hierarchically: Face 0 (individual node dynamics), Face 1 (pairwise synchronization), and Face 2 (global synchronization).

4.4.1. Face 0: Individual Node Dynamics

We begin with the analysis of Face 0, focusing on the local dynamics of individual lasers. Figure 18 presents the bifurcation diagram of the laser peak intensity and the corresponding maximum Lyapunov exponent as functions of the coupling strength ${\sigma }_2$.

In contrast to the previous cases, the system dynamics in Case 4 is considerably more complex. For $0<{\sigma }_2<1.514$, the lasers exhibit chaotic dynamics interspersed with two narrow periodic windows: one near ${\sigma }_2\approx 0.9$ and another near ${\sigma }_2\approx 1.25$, where the maximum Lyapunov exponent ${\lambda }_{\max }\approx 0$. In the interval $1.515\le {\sigma }_2\le 1.615$, we observe the coexistence of three stable periodic orbits. At ${\sigma }_2\approx 1.615$, the system undergoes a torus bifurcation and becomes quasiperiodic.

Overall, the dynamics in this case encompasses chaos, periodicity, quasiperiodicity, and multistability, making it the richest among the four cases studied.

Figure 19 illustrates representative time series (left column), phase portraits (middle column), and power spectra (right column) for selected values of ${\sigma }_2$.

Chaotic regime:Figure 19(a–c) illustrate the chaotic regime at ${\sigma }_2=0.3$. The time series exhibits irregular fluctuations, the phase portrait reveals a complex, non-repeating structure, and the power spectrum displays a broad, continuous distribution—all characteristic signatures of chaos.

Multistable regime:Figure 19(d–f) show the multistable regime at ${\sigma }_2=1.56$, where three distinct periodic attractors coexist for different initial conditions. These attractors have the same dominant frequencies but differ in their phase-space portraits, reflecting distinct oscillation amplitudes.

Quasiperiodic regimes:Figure 19(g–i) and Figure 19(j–l) illustrate quasiperiodic regimes at ${\sigma }_2=2.0$ and ${\sigma }_2=3.5$, respectively. In both cases, the dynamics are characterized by two incommensurate frequencies, giving rise to distinct envelope modulations.

For ${\sigma }_2=2.0$ (Figure 19(g–i)), the two fundamental frequencies are $f_1=10.64$ kHz and $f_2=16.10$ kHz. Their irrational ratio confirms the quasiperiodic nature of the dynamics, and the power spectrum exhibits additional peaks at linear combinations of these frequencies.

For ${\sigma }_2=3.5$ (Figure 19(j–l)), the characteristic frequencies are $f_1=10.53$ kHz and $f_2=38.00$ kHz. The different frequency ratios in the two quasiperiodic regimes give rise to qualitatively different envelope structures, as evident in the time series and phase portraits.

4.4.2. Face 1: Pairwise Synchronization Analysis

We now examine Face 1 of the simplicial complex for Case 4, characterizing pairwise synchronization along each network link.

Figure 20(a) shows the average pairwise synchronization error E as a function of ${\sigma }_2$. Unlike in previous cases, the average error remains relatively constant across the entire range of ${\sigma }_2$, varying between approximately $2.0$ and $2.25$. As in the previous cases, complete synchronization ($E=0$) is never achieved.

In contrast, as shown in Figure 20(b), the maximum synchronization error decreases monotonically as ${\sigma }_2$ increases. This indicates that while the average level of desynchronization remains stable, the worst-case deviations between lasers become less severe with stronger second-order coupling.

Figure 20(c) and Figure 20(d) present representative time series of the synchronization error for link (1,2) at ${\sigma }_2=0.5$ (weak second-order coupling) and ${\sigma }_2=3.5$ (strong second-order coupling), respectively.

At ${\sigma }_2=0.5$ (Figure 20(c)), the error exhibits large-amplitude peaks interspersed with intervals of near-zero error, indicating episodes of brief synchronization followed by sudden desynchronization.

At ${\sigma }_2=3.5$ (Figure 20(d)), the peak amplitude of the error is considerably smaller, and the intervals between successive error peaks are shorter. This suggests that stronger second-order coupling suppresses large desynchronization events, leading to more regular and smaller-amplitude fluctuations in the synchronization error.

Figure 21 presents the similarity S and correlation C as functions of the coupling strength ${\sigma }_2$.

As ${\sigma }_2$ is increased, S progressively decreases while C increases. The regions where $S\approx 0$ and $C\approx 1$ correspond to periodic windows in which the laser dynamics are highly similar. The intervals where multiple distinct branches appear (approximately three) indicate regions of multistability, where several attractors coexist.

For stronger couplings (${\sigma }_2>1.7$), where quasiperiodic motion occurs, both S and C exhibit alternating increases and decreases. This behavior indicates that the system in this regime alternates between synchronous and asynchronous dynamics, reflecting intermittent loss and recovery of coherence.

Figure 22(a) shows the average phase difference $\varphi$ between laser oscillations as a function of ${\sigma }_2$. In contrast to Case 3, where $\varphi$ decreased with increasing coupling, here $\varphi$ progressively increases, indicating a loss of phase coherence as second-order coupling strengthens.

Figure 22(b–d) present the corresponding time series of ${\varphi }_{12}$ for three representative coupling strengths. At ${\sigma }_2=0.3$ (Figure 22(b)), intervals where $\varphi \approx 0$ are observed, corresponding to nearly complete synchronization. At ${\sigma }_2=1.7$ (Figure 22(c)), the intervals of phase coincidence become shorter, while the frequency of $\varphi$ oscillations increases. At ${\sigma }_2=3.8$ (Figure 22(d)), both the amplitude and frequency of the $\varphi$ fluctuations are higher, indicating a significant degradation of phase synchronization. These observations confirm that increasing ${\sigma }_2$ systematically undermines phase coherence in the network.

4.4.3. Face 2: Global Synchronization Analysis

Finally, to study global synchronization corresponding to Face 2 of the simplicial complex, we analyze the global average synchronization error and the Kuramoto order parameter.

Figure 23(a) shows the global average synchronization error E. As ${\sigma }_2$ increases, E increases slightly, indicating a gradual loss of global coherence in the system. The fluctuations observed in the average error suggest the presence of intermittent dynamics, likely arising from the coexistence of partially synchronized and desynchronized regimes.

Figure 23(b) presents the global order parameter R, defined within the Kuramoto framework as a quantitative measure of the degree of phase synchronization in the ensemble [26]. As ${\sigma }_2$ increases, R gradually decreases from values close to unity, indicating strong, though not complete, synchronization, to lower values, signifying a reduction in network synchronization coherence [5,27].

The inverse relationship between the global synchronization error (increasing) and the order parameter (decreasing) confirms that stronger second-order coupling in this configuration (fixed ${\sigma }_1=1.5$) systematically degrades both amplitude and phase coherence across the network.

5. Discussion and Conclusion

We investigate synchronization dynamics of a network of three EDFLs with first- and second-order interactions, modeled within a simplicial complex framework. Four coupling configurations are analyzed, corresponding to different combinations of pairwise (${\sigma }_1$) and higher-order (${\sigma }_2$) coupling strengths.

5.1. Summary of Dynamical Regimes

The analysis of individual node dynamics (Face 0) using bifurcation diagrams, maximum Lyapunov exponents, time series, phase-space representations, and Fourier spectra demonstrates that the incorporation of HOI significantly enriches collective dynamics of the network.

In Case 1 (${\sigma }_2=0$, pairwise coupling only), the system exhibits predominantly chaotic behavior interspersed with narrow periodic and quasiperiodic windows. This case serves as a baseline for purely first-order interactions.

In Case 2 (${\sigma }_1=0$, second-order coupling only), the dynamics are considerably simpler, as this configuration does not constitute a complete simplicial complex. Consistent with the topological framework and the use of the higher-order Laplacian matrix, a limited diversity of dynamical behaviors is observed.

In Case 3 (${\sigma }_2=1.5$, fixed second-order coupling), the system exhibits more complex dynamics, including multistability and quasiperiodicity, while chaotic intervals are reduced compared to Case 1.

In Case 4 (${\sigma }_1=1.5$, fixed first-order coupling), the number of periodic and quasiperiodic windows in which multistability is observed increases, making this the richest dynamical scenario among the four cases.

These findings are consistent with previous studies reporting that HOI increase system complexity [5,6].

5.2. Pairwise Synchronization (Face 1)

The analysis of local synchronization between pairs of nodes (Face 1) reveals that, in all cases, regions of periodic behavior exhibit delayed or anticipated synchronization. This is demonstrated through similarity and cross-correlation measures. Such synchronization characteristics are common in coupled oscillatory systems where trajectories maintain a functional relationship with a constant time shift.

In quasiperiodic regimes, delayed synchronization persists but not uniformly across all parameter values, as reflected by oscillations in the similarity and correlation measures, which alternate between high and low values. In contrast, in multistable periodic and chaotic windows, delayed or anticipated synchronization is not clearly identified.

The pairwise synchronization error analysis complements these findings, showing that the error never approaches zero. This indicates that complete synchronization is not achieved in any of the cases studied.

5.3. Phase Difference Analysis

The study of phase differences between nodes reveals that, although the average phase differences are small, they never reach zero, as confirmed by time series analysis. Partial phase synchronization exists, but complete synchronization is not achieved, consistent with the standard classification of synchronization regimes in coupled nonlinear systems [5].

Notably, for Cases 1 and 3, the average phase difference decreases as the coupling strength is increased. In contrast, for Case 4, the average phase difference increases with coupling. This opposite behavior highlights that the local coherence of the system depends critically on the type of interaction applied (first-order versus second-order coupling).

5.4. Global Synchronization (Face 2)

The global synchronization analysis employs two complementary metrics: the global synchronization error and the Kuramoto order parameter. The global synchronization error exhibits a qualitatively similar trend to the local synchronization error, albeit with larger absolute values. This confirms that complete synchronization is absent across the entire network.

The Kuramoto order parameter R shows an inverse relationship to the local phase difference analysis. While $R=1$ would indicate a state of global phase synchronization [26], this unit value is never attained in any of the four cases, indicating that the system does not achieve complete phase synchronization. Instead, partial phase coherence is maintained across the network.

5.5. Comparative Analysis and General Implications

A comparison of the four cases demonstrates that the overall changes in dynamics are not always obvious when each case is analyzed separately. However, the results collectively show that second-order interactions enrich the dynamical repertoire of the system and modify its synchronization properties. Conversely, the absence of either coupling mechanism (i.e., the lack of a simplex) restricts the dynamics and reduces complexity.

In the analysis of a nonlinear system such as the EDFL with first- and second-order interactions, a wide variety of dynamical regimes is observed, including periodic, quasiperiodic, multistable, and chaotic behavior, depending on both the coupling strength and its configuration.

Regarding local synchronization, coherence exists between pairs of laser nodes, with delayed or anticipated synchronization manifesting in periodic regimes and partial synchronization in quasiperiodic intervals. Synchronization is not observed in chaotic or multistable zones. Globally, the network does not reach complete synchronization, as evidenced by the global synchronization error, but the order parameter indicates that the network remains partially phase-synchronized.

5.6. Final Remarks

In summary, this work deepens our understanding of HOI in coupled nonlinear systems. Our analysis is especially relevant to higher-order networks, where collective coherence and synchronization are key to emergent behavior.

Comparing four case studies reveals that HOI fundamentally shape collective dynamics by increasing system complexity. Additionally, synchronization is shown to depend on both coupling strength and interaction type.

These findings highlight the necessity of including higher-order effects in future models of complex networked systems.