Solitons unveilings and modulation instability analysis for sixth-order coupled nonlinear Schrödinger equations in fiber bragg gratings

1.   Introduction

A Solitary wave is a particular kind of long-wave which is not dispersive and travels with constant speed, keeping its original shape. This wave appears among the solutions of nonlinear partial differential equations (NLPDEs) ^[1,2,3]. Unlike normal waves, solitary waves don't merge and when colliding with each other, and they keep their original shapes, sizes, and speeds. This particle-like behavior was the reason for calling these solitary waves "solitons" ^[4,5]. Gap solitons appear inside finite gaps within the domain of continuous systems. They differ from regular solitons in the group velocity dispersion of photonic band structure ^[6]. Recent studies have focused on the dynamics of high-order solitons in optical fibers with engineered dispersion profiles. These solitons exhibit complex interactions and splitting behaviors, which are crucial for applications in ultrafast optics and wavelength-division multiplexing (WDM) systems ^[7]. Advances in understanding the soliton self-frequency shift (SSFS), particularly in photonic crystal fibers (PCFs), have enabled the generation of tunable frequency combs and supercontinuum sources. These developments are significant for applications in metrology and spectroscopy ^[8]. Recent work has explored the use of solitons for nonlinear pulse compression in fibers, achieving few-cycle pulses with high peak powers. This has implications for ultrafast laser systems and attosecond science ^[9]. Recent research has investigated the behavior of gap solitons in chirped fiber Bragg gratings (FBGs), where the grating period varies along the fiber length ^[10,11]. These solitons exhibit unique trapping and slowing properties, making them promising for optical delay lines and buffering applications ^[12,13,14]. In optical fiber, information is sent as packets like solitons which provide high speed connectivity. This is because they travel in optical fiber at the speed of light ^[15]. Optical fiber is constructed by a core coated with a translucent cladding material having lower refractive index. Soliton's transmission through optical fiber depends on the equilibrium of chromatic dispersion with fiber nonlinearity. The decrease in chromatic dispersion restricts wave propagation through the fiber, which poses a significant problem ^[16]. Bragg gratings technology is used to solve this problem by compensating the low chromatic dispersion. In this technology, an intense ultraviolet source is applied to create repetitive refractive index differences in the core of the optical fibe, which is called FBGs. When the light travels through the FBGs, the refractive index differs. This causes a small light quantity to be reflected on boundaries. If the periods of the light wavelength and the grating match, a positive backward reflection is achieved as power couples from forward to backward directions. Light with other wavelengths can't travel through fiber. Gap solitons arise from the resonant reflection of light on Bragg gratings ^[6,17,18].

The applications of FBGs are numerous in optical communications and sensors. FBG sensors are used in geothermal engineering field observation, communication systems as multiplexers and demultiplexers with optical circulators or optical add-drop multiplexers, health monitoring systems for aircraft, underwater acoustic sensors, biomedical applications, robotics, and other fields. FBG sensors offer many advantages that make them widely used in real-life applications, such as being lightweight, stable, resistant to electromagnetic interference, fast-responding, corrosion-resistant, capable of multi-point sensing, affordable, and compact ^[17,19,20]. The use of FBGs in artificial intelligence and machine learning is expected to increase the abilities of optical fiber sensors in the future. This will introduce new levels of accuracy of data and real time monitoring ^[21,22].

These important applications encouraged intensive studies of solitons in FBGs with variant formulas of nonlinear refractive index like Kerr nonlinearity law ^[23], parabolic nonlinearity law ^[24], polynomial nonlinearity law ^[25], quadratic-cubic nonlinearity law ^[26], parabolic-nonlocal combo nonlinearity law^[27], and many others in the literature ^[28,29].

The standard model of light propagation in nonlinear FBGs is based on a system of coupled-mode nonlinear Schrödinger equations for the right (forward) and left (backward) traveling waves ^[6,30]. Many different forms of coupled-mode systems were suggested and investigated earlier in the literature with higher dispersion orders and different nonlinearity laws. For example, in ^[31], the auxiliary equation method was used to investigate a coupled system of second order dispersion with Kudryashov's law of self-phase modulation. The Lie symmetry analysis is applied to a different version of the coupled system of fourth order dispersion with Kerr law of nonlinearity in ^[23]. The Jacobi elliptic function approach is used in ^[32] to investigate another coupled system of second order dispersion and parabolic law of nonlinearity. In ^[33], three integration schemes, the unified Riccati equation method, new extended auxiliary equation method, and unified auxiliary equation method, are applied on a different version of the coupled system of second order dispersion with parabolic nonlocal combo nonlinearity. In ^[34], a different version of the coupled system of second order dispersion and anti-cubic nonlinearity is investigated by the extended auxiliary equation approach. The extended trial function approach is applied to a coupled system of second order with parabolic nonlocal combo nonlinearity in ^[27] and many others ^[35,36].

Following the same motivation, this paper investigates a coupled system of sixth-order, incorporating both even and odd dispersion orders, describing the propagation of light waves in FBGs with Kerr and parabolic nonlocal combo laws of nonlinearity by the modified Sardar sub-equation (MSSE) technique. Rezazadeh et al. ^[37] introduced an efficient method, the Sardar sub-equation method (SSE), to solve the Benjamin-Bona-Mahony equation. SSE is a general method from which, under specific conditions, many techniques such as the functional variable method and first integral techniques can be derived^[37]. Akinyemi et al. ^[38] proposed a modification to SSE, called MSSE, from which many techniques, such as the $(G' /G)$ - expansion method, generalized auxiliary equation method, and extended tanh-function method, can be derived. Subsequently, various forms of MSSE were introduced, involving different substitutions or wave transformations. These forms either simplify the sub-equation or introduce a novel procedure to solve it. The advantages of MSSE over other techniques for solving NLPDEs include computational efficiency, fast equation solving, and high accuracy, particularly in cases of high nonlinearity ^[39,40]. One observed limitation of the MSSE technique in this research is its inability to generate dark solitons.

In this work, the MSSE technique is applied to solve the following coupled system of highly dispersive perturbed nonlinear Schrödinger equations in FBGs with parabolic nonlocal combo nonlinear refractive index ^[41]:

and

which models the light wave propagation in nonlinear FBGs with $Q\equiv Q(x, \, t)$ and $R \equiv R(x, \, t)$ as forward and backward light waves. The coordinate along the fiber is $x$ and $t$ is time. $a_{1}, \, b_{1}, \, \alpha_{1}, \, \alpha_{2}$ and $a_{2}, \, b_{2}$ are the inter-modal dispersion and the chromatic dispersion coefficients, respectively. $a_{n}, \, b_{n}\, (n = 3, 4, 5, 6)$ are coefficients of $n^{th}$ order dispersion. $c_{n}, \, e_{n}\, (n = 1, \, 2)$ are self-phase modulation coefficients. $d_{n}, \, g_{n}\, (n = 1, \, 2)$ are coefficients of cross phase modulation. $f_{n}\, (n = 1, \, 2)$ are nonlinear terms coefficients. $l_{n}, \, m_{n}\, (n = 1, \, 2)$ are nonlocal nonlinearity terms coefficients. $\beta_{n}$, $\gamma_{n}\,$ and $\sigma_{n}\, (n = 1, 2)$ are detuning parameters, self-steeping terms, and four wave mixing parameters coefficients, respectively. $\theta_{n}, \, \mu_{n}\, (n = 1, 2)$ are nonlinear dispersion terms coefficients^[41]. This coupled system was suggested and investigated by the extended auxiliary equation approach earlier in ^[41], producing only bright gap solitons and singular gap solitons. Different versions of this coupled system were studied earlier in the literature. For example, a modified version of Eqs (1.1) and (1.2) including a conformable fractional derivative is studied in ^[42] using the modified extended direct algebraic method producing bright and singular solitons, as well as hyperbolic and trigonometric solutions. Meanwhile in ^[43], a different version of Eqs (1.1) and (1.2) with only nonlocal law of self-phase modulation and differential group delay was studied using SSE approach, yielding bright and singular solitons. Also, in ^[44], the dynamics of highly dispersive gap solitons within the framework of the Kundu-Eckhaus equation with a new addition of multiplicative white noise were investigated by two methods, the extended simplest equation approach and the generalized Riccati equation mapping scheme, yielding bright, singular, and dark-singular straddled solitons.

The motivation and novelty of the present study is to derive novel solutions for the coupled system of Eqs (1.1) and (1.2). In addition, to apply the MSSE technique for the first time to the coupled system of Eqs (1.1) and (1.2) is to get benefit of the great advantages of the MSSE technique in dealing with this highly dispersive system with high nonlinearity. The MSSE technique succeeded to derive various types of solutions such as bright gap solitons, singular gap solitons, which were previously obtained in ^[41], as well as hyperbolic, exponential, singular periodic, and rational solutions, which are novel solutions that were not previously obtained in the literature. The obtained solutions are verified by direct substitution in the system under study. Furthermore, the modulation instability analysis for the governing system of Eqs (1.1) and (1.2) is investigated using linear stability analysis, and the analytical expression for the modulation instability gain spectrum is derived for the first time in the literature, to the best of our knowledge. The modulation instability can cause the appearance of optical rogue wave in the process of supercontinuum generation ^[45]. Furthermore, visual illustrations of the modulation instability gain spectrum for the system under study and are presented in this work for the first time. The modulation instability phenomena suggests that any small deviations from the steady-state solutions will amplify and lead to diverging behavior. That is why the modulation instability analysis is important for higher-order nonlinear models. These models exhibit an instability that requires the investigation of steady-state modulation due to the interaction of nonlinear and dispersive effects ^[46,47,48].

The structure of the rest of the paper is: Section 2 explains the MSSE technique algorithm. Then, the application of the MSSE technique algorithm step by step to the coupled system of Eqs (1.1) and (1.2) and the obtained analytic solutions, using the Mathematica software, is in Section 3. In Section 4, the modulation instability analysis for the governing system using linear stability analysis is presented. The derivation of analytical expression for the modulation instability gain spectrum is presented also in Section 4. After that, the graphical illustrations of a sample of derived solutions, with 2D and 3D graphs, as well as the visual illustrations of the modulation instability gain spectrum, with 2D and 3D graphs, are presented in Section 5. The conclusion is in Section 6.

2.   The algorithm of the suggested technique

The suggested MSSE technique is:

For the NLPDE with the unknown wave function $f(x, \, t)$:

① Let

with $Y(\eta)$, $\rho$, $\, r\,$, $\lambda$, and $v$ as the wave's amplitude, number, frequency, length, and velocity, respectively.

Hence, Eq (2.1) reduces to:

$H$ is a function in $Y, \, Y' = \frac{dY}{d\, \eta}, \, Y'' = \frac {d^{2}Y}{d\, \eta^{2}}, ...etc$.

② The suggested solution of Eq (2.3) is:

with $\varphi(\eta)$ satisfying:

where M is an integer, $K_{j}\, (j = 0, \, 1, \, ..., M)$, and $w_{0}, \, w_{1}$, and $\, w_{2}$ are constants.

③ M is known from Eq (2.3) by balancing the highest nonlinearity and dispersion orders.

④ Inserting Eqs (2.4) and (2.5) in Eq (2.3) produces a polynomial in $\varphi(\eta)$. The polynomial's coefficients are then set to vanish, developing a system of equations in $K_{j}\, (j = 0, \, 1, \, ..., M)$, $w_{0}, \, w_{1}$, $\, w_{2}$, $\, \rho, \, r, \, \lambda$, and $v$. Software like Mathematica produces this system's solution.

⑤ The general analytical solutions of Eq (2.5), with k as any arbitrary constant, are:

Case Ⅰ: For $w_{0} = 0, \, w_{1} > 0$, and $w_{2} \ne 0$,

Case Ⅱ: For $w_{0} = 0, \, w_{1} > 0, \,$ and $\, w_{2} = \pm4 B_{1}B_{2}$, with $B_{1}$ and $B_{2}$ as two arbitrary parameters,

Case Ⅲ: For $w_{0} = \frac{w_{1}^{2}}{4\, w_{2}}, \, w_{1} < 0, \,$ and $\, w_{2} > 0$, with $B_{1}$ and $B_{2}$ as two arbitrary parameters,

Case Ⅳ: For $w_{0} = 0, \, w_{1} < 0$, and $w_{2} \ne 0$,

Case Ⅴ: For $w_{0} = \frac{w_{1}^{2}}{4\, w_{2}}, \, w_{1} > 0, \,$ and $\, w_{2} > 0$, with $B_{1}^{2}- B_{2}^{2} > 0$ for any two arbitrary parameters $B_{1}$ and $B_{2}$,

Case Ⅵ: For $w_{0} = 0$ and $w_{1} > 0$,

Case Ⅶ: For $w_{0} = w_{1} = 0$ and $w_{2} > 0$,

These listed steps of the MSSE algorithm will be applied to the system under study in the next section in detail.

3.   MSSE technique implementation

3.1.   Methodology

The first step of implementation of the MSSE technique, explained in Section 2, is to assume the analytic solutions to Eqs (1.1) and (1.2) as:

where $P_{1}(\eta)$ and $P_{2}(\eta)$ are the forward and backward waves' amplitudes, respectively. The insertion of Eq (3.1) in Eqs (1.1) and (1.2) generates the following real and imaginary parts:

and

Now, under the assumption:

where the constant is $\psi\in \mathbb{R}$$-\left\{{0, 1}\right\}$, Eq (3.2) and Eq (3.5) are reduced to:

Equating the coefficients of Eqs (3.9) and (3.10) to zero gives:

Both of the Eqs (3.7) and (3.8) are the same and can be written as:

where

From Eqs (3.11) and (3.13), the following relations are deduced:

The next step of the MSSE technique in Section 2 is the balancing of the dispersion term $P_{1}^{(6)}(\eta)$ and the nonlinear term $P_{1}^{5} (\eta)$ in Eq (3.12), $i.e.,$ $5\, M = (M+6)$, and, hence, $M = 3/2$. In order to get an integer value of $M$, the following transformation is introduced:

which, when inserted in Eq (3.12), gives:

Now, balancing the nonlinear term $U^{12}(\eta)$ and the dispersion term $U^{5}(\eta)\, U^{(6)}(\eta)$ of Eq (3.16), $i.e.,$ $12 M = 5\, M+(M+6)$, and, hence, $M = 1$, yields:

The following system is generated from the insertion of Eqs (2.5) and (3.17) in Eq (3.16), followed by forcing the coefficients of $\varphi(\eta)$ to vanish:

The following result is found for the system in Eq (3.18) using the Mathematica software (14):

Result:

Solution verification:

The obtained result is verified by direct substitution as follows:

First, the insertion of parameters from Eq (3.19) into Eqs (2.5), (3.15), (3.17), and (3.1) yields:

and, hence,

After that, the solutions in Eq (3.22) are substituted into Eqs (1.1) and (1.2), followed by the separation of real and imaginary parts. Under the conditions in Eq (3.11), the imaginary part vanishes and the real part in Eq (3.12) reduces to:

Finally, the substitution of $P_1$ from Eq (3.21) and its derivatives, with the aid of Eq (3.20), into Eq (3.23), yields an identity using the Mathematica software where all coefficients of $\varphi$ reduce to zero.

3.2.   Derived analytic solutions

The generated result's parameters are inserted in Eq (3.17), then in Eq (3.1). Then, using Eq (2.6) through Eq (2.23), the analytic solutions for forward and backward traveling waves in the coupled system of Eqs (1.1) and (1.2) are:

Form 1: With $w_{0} = 0, \, w_{1} > 0$, and $w_{2} \ne 0$, the bright ($w_{2} < 0$) and singular ($w_{2} > 0$) solitons are:

and

Form 2: With $w_{0} = 0, \, w_{1} > 0\,$, and $\, w_{2} = \pm4 B_{1}B_{2}$, $B_{1}$ and $B_{2}$ are arbitrary parameters, and the hyperbolic solutions are:

Form 3: With $w_{0} = 0, \, w_{1} < 0$, and $w_{2} > 0$, the singular periodic solutions are:

and

Form 4: With $w_{0} = 0$, $\, w_{1} > 0$, and $w_{2} < 0$, the exponential solutions are:

and

Form 5: For $w_{0} = w_{1} = 0$ and $w_{2} > 0$, the rational solutions are:

The visual illustrations of some of these obtained solutions will be presented in Section 5 after investigating the modulation instability of the system under study in the next section.

4.   Modulation instability analysis

The analysis of modulation instability (MI) is a crucial aspect for understanding the behavior of various physical systems. This type of analysis focuses on the conditions under which modulation instability occurs and its implications for the dynamics of wave patterns. The standard linear stability technique is employed to explore the MI of the system under study in this section. The analysis starts by the formation of a steady-state solution followed by introducing perturbations to the steady-state solution, then ends by finding an analytical expression for the modulation instability gain spectrum ^[49,51].

Assume that Eqs (1.1) and (1.2) have the following steady-state solutions:

where $\sqrt{\mathcal{A}}$ is a real constant amplitude (normalized power), $B$ is the phase shift, and $\psi\in \mathbb{R}$$-\left\{{0, 1}\right\}$ as in Eq (3.6). The substitution of the steady-state solutions in Eqs (4.1) and (4.2) into Eqs (1.1) and (1.2) gives:

To perform a stability analysis, small perturbations are introduced in Eqs (4.1) and (4.2), then the perturbed steady-state solutions take the following form:

$F(x, t)$ is the perturbation complex amplitude and $\mid F(x, t)\mid < < \sqrt{\mathcal{A}}$. By incorporating Eqs (4.4) and (4.5) into Eqs (1.1) and (1.2) then linearizing the equations with respect to $F(x, t)$ and its conjugate $F^*(x, t)$, both equations reduce to the following equation:

where the coefficients obey Eq (3.14). Now, assume the solution of Eq (4.6) is in the form:

where $r_1, \, r_2$ are constants, $\omega$ is the perturbation frequency, and $L$ is the perturbation wave number. The substitution of Eq (4.7) into Eq (4.6), followed by the separation of the coefficients of $\ e^{i (L\ x-\omega\ t)}$ and $\ e^{-i (L\ x-\omega\ t)}$, yields:

which can be rearranged in matrix form as:

where

$T_1 = \omega -\alpha _1 L+\left(-a_6 L^5 -a_5 L^4 +a_4 L^3 +a_3 L^2 -a_2 L-a_1 \right) L \psi + 2 A^2 \left(e_1+ f_1 \psi ^2+ g_1 \psi ^4 \right)$

$+A \left(c_1+ \left(\gamma _1 + \mu _1 \right) L+ \left(d_1 + \sigma _1 \right) \psi ^2 \right)$,

$T_2 = A \left(\psi ^2 \left(2 A \left(f_1+g_1 \psi ^2\right)+d_1+\sigma _1\right)+2 A e_1+c_1\right)$,

$T_3 = A\left(\psi ^2 \left(2 A \left(f_1+g_1 \psi ^2\right)+d_1+\sigma _1\right)+2 A e_1+c_1\right)$,

and

$T_4 = \alpha _1 L-\omega +\left(-a_6 L^5 +a_5 L^4 +a_4 L^3 -a_3 L^2 -a_2 L +a_1 \right) L \psi +2 A^2 \left(e_1+ f_1 \psi ^2+ g_1 \psi ^4 \right)$

$+ A \left(c_1- \left(\gamma _1 + \mu _1 \right) L+ \left(d_1 +\sigma _1 \right)\psi ^2 \right)$.

The system of equations in Eq (4.9) has a nontrivial solution if: $\begin{vmatrix} T_1 & T_2 \\ T_3 & T_4 \\ \end{vmatrix} = 0,$ and, hence, the dispersion relation is:

where $S = L^2 \psi \left(a_6 L^4-a_4 L^2+a_2\right)$. The stability of the perturbed steady-state solutions in Eqs (4.4) and (4.5) depends on $\omega$ in Eq (4.10). If $\omega$ is complex, the steady-state solution becomes unstable as the perturbations will grow exponentially in Eqs (4.4) and (4.5) this is known as MI. The MI suggests that any small deviations from the steady-state will amplify and lead to diverging behavior. Meanwhile, if $\omega$ is real, the steady-state solution is stable and the MI does not occur ^[49,50]. Therefore, from Eq (4.10), the MI occurs if

$S\left(S -2A\left(2 A \left(e_1+ f_1 \psi ^2+ g_1 \psi ^4 \right)+ \left(c_1+ \left(d_1 + \sigma _1 \right) \psi ^2\right)\right)\right) < 0$.

Otherwise, the steady-state solution is stable. Consequently, the MI gain spectrum $G(A, L)$, which quantifies the degree of instability, can be written as:

From Eq (4.11), the MI gain spectrum is a function of the normalized power, the perturbation wave number, the dispersion coefficients, the phase modulation coefficients, and the nonlinearity coefficients. The visual illustrations of the MI gain spectrum in Eq (4.11) are presented in the next section.

5.   Results

5.1.   3D and 2D graphics of derived solutions

The MSSE technique is implemented for the coupled system of Eqs (1.1) and (1.2), modeling the forward and backward light waves propagation through FBGs with Kerr and parabolic nonlocal combo laws of nonlinearity. The proposed MSSE technique succeeded to derive various types of solutions such as bright gap solitons, singular gap solitons, which were previously obtained in ^[41], as well as hyperbolic, exponential, singular periodic, and rational solutions, which are novel solutions that were not previously obtained in the literature. Graphs (3D and 2D) of a sample of achieved solutions are shown in Figures 1–4.

The graphs of bright solitons $|Q_{1, 1}|^2$ and $|R_{1, 1} |^2$ from Eqs (3.24) and (3.25) are presented in Figure 1. Bright solitons are characterized by a centered peak of maximum intensity ^[52]. These solitons are crucial in optical communication systems because they can propagate over long distances without distortion, supporting high-speed data transmission efficiency ^[42]. The 2D and 3D plots illustrate the stable propagation of bright solitons, emphasizing their role in maintaining signal integrity and enabling high-performance optical communication networks.

The graphs of singular solitons \(|Q_{1, 2}|^2\) and \(|R_{1, 2}|^2\) from Eqs (3.26) and (3.27) are shown in Figure 2. Singular solitons feature a centered infinite amplitude and are valuable for studying extreme wave phenomena, aiding in the development of realistic models ^[52]. The plots highlight the unique properties of singular solitons, which can be used to model extreme wave behavior and enhance the understanding of nonlinear dynamics in FBGs.

Figure 3 presents the 2D and 3D graphs of singular periodic solutions \(|Q_{3, 1}|^2\) and \(|R_{3, 1}|^2\) from Eqs (3.30) and (3.31). Singular periodic solutions have significant applications in advanced signal processing techniques ^[52]. The graphical representations demonstrate the oscillatory nature of these solutions, which can be exploited for pulse shaping, wavelength conversion, and other signal processing applications in optical communication systems.

The dynamics of rational solutions \(|Q_5|^2\) and \(|R_5|^2\) from Eqs (3.38) and (3.39) are depicted in Figure 4. The 2D and 3D plots illustrate the behavior of rational solutions, providing insights into their potential applications in modeling complex wave interactions and optimizing signal transmission in FBGs.

5.2.   3D and 2D graphics of MI gain spectrum

In this section, the graphical illustrations of the MI gain spectrum $G(A, L)$ derived in Eq (4.11) versus different parameters are discussed and plotted in Figures 5–7.

Figure 5 (a) presents the 3D graph of the MI gain spectrum $G(A, L)$ derived in Eq (4.11). Figure 5 (b) presents the 2D graph of $G(A, L)$ versus the perturbation wave number $L$ for three different values of normalized power $A = \{0.5, 1, 1.5\}$. Figure 5 (c) presents the 2D graph of $G(A, L)$ versus the normalized power $A$ for three different values of perturbation wave number $L = \{1, 1.5, 2\}$. The graph of $G(A, L)$ versus $L$ in Figure 5 (b) shows that the gain does not depend on the sign of $L$. Also, Figure 5 (b) shows that the gain has two local maxima. If the normalized power $A$ increases, these two local maxima occur at higher values of $L$, and the instability region width increases. The graph of $G(A, L)$ versus $A$ in Figure 5 (c) shows that as the wave number $L$ increases, the onset of the instability gain occurs at higher values of the normalized power $A$.

Figure 6 (a) and (b) present the MI gain spectrum $G(A, 1.5)$ derived in Eq (4.11) versus the normalized power $A$ and $G(0.5, L)$ versus the perturbation wave number $L$, respectively, for $f_{1} = \{0.5, 1, 1.5\}$. Figure 6 (c) and (d) present the MI gain spectrum $G(A, 1.5)$ derived in Eq (4.11) versus the normalized power $A$ and $G(0.5, L)$ versus the perturbation wave number $L$, respectively, for $e_{1} = \{-0.5, 0.5, 1.5\}$. Figure 6 (e) and (f) present the MI gain spectrum $G(A, 1.5)$ derived in Eq (4.11) versus the normalized power $A$ and $G(0.5, L)$ versus the perturbation wave number $L$, respectively, for $g_{1} = \{0.5, 1, 1.5\}$.

Figure 7 (a) and (b) present the MI gain spectrum $G(A, 1.5)$ derived in Eq (4.11) versus the normalized power $A$ and $G(0.5, L)$ versus the perturbation wave number $L$, respectively, for $l_{1} = \{0.6, 0.8, 1\}$. Figure 7 (c) and (d) present the MI gain spectrum $G(A, 1.5)$ derived in Eq (4.11) versus the normalized power $A$ and $G(0.5, L)$ versus the perturbation wave number $L$, respectively, for $a_{2} = \{0.6, 1, 1.5\}$. Figure 7 (e) and (f) present the MI gain spectrum $G(A, 1.5)$ derived in Eq (4.11) versus the normalized power $A$ and $G(0.5, L)$ versus the perturbation wave number $L$, respectively, for $a_{6} = \{0.6, 1, 1.5\}$.

For constant value of wave number $L$, Figure 6 (a), (c), (e) and Figure 7 (a), (c), (e), show that at a certain value of the normalized power $A$, any increase in the values of the coefficients $f_1, \, e_1, \, g_1, \, l_1, \, a_2, \,$ and $a_6$ causes an increase in the gain. Figure 6 (b), (d), (f) and Figure 7 (b) show that when the values of the coefficients $f_1, \, e_1, \, g_1,$ and $\, l_1$ increase, the maximum of the gain spectrum has higher value and occurs at higher value of the perturbation wave number $L$. Meanwhile, for a constant value of the normalized power $A$ in Figure 7 (d) and (f), if the values of the coefficients $a_2$ and $\, a_6$ increase, the maximum of the gain spectrum has the same value and occurs at slightly lower perturbation wave number $L$.

6.   Conclusions

The novelty of this work includes several achievements: First, the application of the MSSE technique for the first time to solve the coupled system of Eqs (1.1) and Eq (1.2). This coupled system of highly dispersive (sixth-order with all even and odd orders present) perturbed nonlinear Schrödinger equations represent the propagation of forward and backward light waves in FBGs with Kerr and parabolic nonlocal combo nonlinearity laws. The study of FBGs is vital due to their extensive utilization in optical communication systems and sensors. The MSSE technique has many advantages such as low computational cost, high consistency, and simple calculations. One observed limitation of the MSSE technique in this research is its inability to generate dark solitons. This encourages the search of other techniques capable of generating dark solitons. Second, the derivation of various types of solutions, using Mathematica software, such as bright gap solitons and singular gap solitons, which were previously obtained in ^[41], as well as hyperbolic, exponential, singular periodic, and rational solutions, which are novel solutions that were not previously obtained in the literature. The dynamics of a sample of the obtained analytic solutions are illustrated with 2D and 3D graphs for parameter values that satisfy the limiting conditions. The obtained solutions are verified by direct substitution in the system under study. During this research, the $\phi^6$ method was investigated for the coupled system under study. The $\phi^6$ did not yield any solutions, consistent with ^[41], which states that this coupled system has no Jacobi elliptic solutions. Third, the derivation of the MI analysis for the governing system by the linear stability technique and obtaining the analytical expression for the MI gain spectrum for the first time in the literature. The FBG system studied in this work combines nonlocal nonlinearity with a periodic grating structure, leading to a unique MI gain spectrum. The derived spectrum shows bandgap-like features due to the grating, as well as modifications due to the nonlocal nonlinearity. The inclusion of Kerr nonlinearity and parabolic nonlocal terms in the FBG system lead to a richer MI gain spectrum compared to simpler systems, with multiple peaks or suppression regions depending on the parameter values. The inclusion of sixth-order dispersion in this work represents a significant extension of previous studies. Higher-order dispersion terms can introduce additional instability bands or modify the existing ones, leading to a more complex MI gain spectrum. The derived MI gain spectrum in this work has additional features due to the grating-induced dispersion and the nonlocal nonlinearity. This asymmetry would highlight the unique role of the grating in modifying the MI dynamics. The inclusion of sixth-order dispersion and parabolic nonlocal nonlinearity in the FBG system further distinguishes it from conventional fibers, where lower-order dispersion and local nonlinearity dominate. The MI gain spectrum derived in this work provides critical insights into the behavior of nonlinear waves in FBGs and has significant physical and practical implications for optical communication systems. MI can lead to the exponential growth of perturbations in the optical signal, causing signal degradation or amplification depending on the system parameters. The derived MI gain spectrum helps identify the conditions under which MI occurs, such as specific power levels, dispersion coefficients, and nonlinearity parameters. This knowledge is crucial for designing FBG-based systems to either mitigate MI-induced signal degradation or exploit it for signal amplification in controlled scenarios. MI is often associated with the breakup of continuous waves into solitons or pulse trains. The MI gain spectrum provides a quantitative measure of the instability regions, which can be used to predict the formation of solitons and their stability. This is particularly important in FBGs, where solitons are used to counteract chromatic dispersion and maintain signal integrity over long distances. Fourth, visual illustrations of the modulation instability gain spectrum for the system under study is presented in this work for the first time. The modulation instability analysis enhances the understanding of optical pulse behavior in FBG systems with nonlinear refractive indices within turbulent nonlinear Schrödinger equations (NLSEs). Dispersion management, in conjunction with modulation techniques and the NLSEs, shapes, evolves, and stabilizes these pulses, thereby enhancing their practical applications in optical communications and signal processing ^[53,54,55].

Author contributions

Noha M. Kamel: Formal analysis, Software, Methodology, Writing-original draft, Visualization; Hamdy M. Ahmed: Validation, Supervision, Conceptualization, Writing-review & editing; Wafaa B. Rabie: Investigation, Writing-review & editing, Visualization.

Use of AI tools declaration

The authors declare that they have not used Artificial Intelligence tools in the creation of this article.

Conflict of interest

All authors declare no conflicts of interest in this paper.