Optical solitons, qualitative analysis and multistability response to study the dynamical behaviour of light wave promulgation

Abstract

Our aim is to explore novel optical solitons within the context of an extended higher-order nonlinear Schrödinger equation that governs the behaviour of propagating light waves. Primarily, this research finds abundant types of solitons, such as singular, kink and trigonometric functions subject to certain constraint conditions. We have utilized $G^{\mathrm{\prime }}/(b{G}^{\mathrm{\prime }}+G+a)$-expansion method, which has proved effective in retrieving these solitons and presenting them for further analysis. The model's dynamic behaviour is then investigated through bifurcation, quasi-periodic oscillations, chaotic behaviour, and sensitivity. These include methods like phase portrait rendering, time series scrutiny, Lyapunov exponents calculation, and the assessment of multi-stability. Finally, sensitivity analysis is conducted at three distinct initial conditions, revealing that the model displays a high level of sensitivity, with substantial alterations occurring in response to minor changes in the initial conditions. The results of this study are revolutionary, intriguing, and possess crucial theoretical importance in evolution disruptions.

Keywords:

• Generalized higher-order NLSE
• $G^{\mathrm{\prime }}/(b{G}^{\mathrm{\prime }}+G+a)$-expansion method
• chaotic behaviour
• bifurcation analysis
• sensitivity analysis

1. Introduction

Solitons represent an exceptional category of waves capable of traversing extensive distances without distortion and retaining their original form even upon encountering other waveforms. Solitons have been observed in various fields, including plasma physics, elastic media, biology, chemical physics, fluid dynamics, high-energy physics and especially in optical fibres [1–5]. Solitons in the optical domain are particularly interesting because of their potential applications in communication and information handling. Therefore, exploration in this area has been growing rapidly. In recent times, long-haul optical fibres have become the standard for transmitting data for several purposes ranging from cross-oceanic to cross-continental.

The electronic, scientific and industrial domains have undergone a revolutionary transformation in the current age. Significant scientific advancements have surfaced due to mankind's curious nature to unveil various physical phenomena and to construct mathematical tools to comprehend their complexities. However, these advancements in scientific disciplines are not enough and have opened new horizons of innovation and cultivation. To keep up with the pace of the current era, an instrument is required that allows humanity to construct intricate structures with exacting features and to utilize highly precise apparatus for analysis and innovation. These capabilities eventually led to the creation of femtosecond lasers. The laser processing strategies have produced a range of sub-micrometer resolution devices, including optical waveguides, three-dimensional photonic crystals and specialized optical instruments for ultra-short pulses. Furthermore, laser processing has made a significant contribution to the field of tissue engineering as well. In the last thirty years, researchers globally have dedicated their endeavours to examining extremely brief pulses lasting picoseconds or even shorter durations across diverse domains, in practical applications.

Several methodologies have been used to attain soliton solutions of different non-linear models in previous literature. A few to mention are; trial equation approach [6], Hirota bilinear transformation technique [7], tanh-function method [8], modified simple equation method [9], new extended direct algebraic approach [10], modified exp-function method [11], Jacobi elliptic function method [12] and unified approach [13].

Our primary goal is to investigate the generalized nonlinear Schrödinger model in order to develop soliton solutions. Several alternative methods have been used to analyze the discussed model in the literature. A few of them are: Raza et al. [14] extracted soliton solutions and analyze stability by using three different methodologies on the model under observation, Seadawy et al. [15] applied generalized Kudryashov method on considered model with higher order to obtain exact travelling wave solutions, Qian et al. [16] extracted new results of discussed model by implementing modified Khater method and Khater et al. [17] obtained a variety of solutions to the considered equation by the implementation of the three schemes.

This study aims to employ the $G^{\mathrm{\prime }}/(b{G}^{\mathrm{\prime }}+G+a)$-expansion technique to attain precise solutions for the presented model. To our understanding, this specific equation has not previously been investigated using this approach [18]. The approach currently employed offers distinct benefits in contrast to methods previously explored, as it produces solutions in the shape of trigonometric and hyperbolic functions. The results are graphically displayed in three dimensional, contour plots and two dimensional plots for specific values of the parameters.

The bifurcation theory to explore differential equations have gained significant attention in recent research [19–22]. Many researchers have explored this concept in both perturbed and unperturbed frameworks. In a recent research endeavour, the authors delved into the bifurcation characteristics, chaotic phenomena and sensitivity analysis of the elliptic Schrödinger equation coupled with an external perturbation. They employed bifurcation and chaos theories as their analytical tools. Bifurcation theory, as referenced in literature [23–25], identifies critical points, often equilibrium points, where the behaviour of a dynamical system is contingent upon the underlying physical parameters. This analysis encompasses the examination of phase portraits associated with equilibrium points in dynamic system influenced by these physical parameters.

In contrast, chaos theory, as cited in references [26–29], aids in ascertaining whether the examined model exhibits chaos, signifying that when an external influence is exerted on the nonlinear processes within a medium, the outcomes may either sustain stability or devolve into a state of unpredictability and disarray. The investigation of these principles has yielded valuable perspectives into the dynamic characteristics of diverse nonlinear models across various domains of science, encompassing disciplines such as engineering and fibre optics. In [30], Sambas et al. proposed a novel mathematical model of a three-dimensional chaotic system that incorporates a closed butterfly-like curve of equilibrium points. Through a comprehensive bifurcation analysis of this novel system, we investigate its inherent characteristics such as chaotic behaviour, multiple stable states and transient chaos. In [31], an innovative three-dimensional chaotic system is introduced, characterized by line equilibrium, and its dynamic characteristics, including Lyapunov exponents, phase portraits, equilibrium points, bifurcation diagram, multistability and coexisting attractors, were thoroughly examined. In [32], a novel chaotic system is presented that consists of two circles containing equilibrium points. The dynamics of the suggested dynamical system are explored by means of assessing Lyapunov exponents, examining bifurcation diagrams and analyzing multistability.

Sensitivity [33] is also discussed for checking the sensitive behaviour the considered model.

In summary, this investigation examines a variety of optical solitons, analyzes the phase portraits of bifurcation analysis, identifies chaotic behaviour using different chaos detection tools, and conducts sensitivity analysis of the generalized higher-order nonlinear Schrödinger equation. These inquiries assist readers in understanding the dynamic observations of the nonlinear physical model from multiple perspectives. The nonlinear Schrödinger equation is a well-known equation that describes the propagation of optical solitons in nonlinear media. It is widely used in describing solitons in a variety of nonlinear optical systems, including optical fibres and waveguides whereas the extended higher-order nonlinear Schrödinger equation involves a higher-order nonlinearity, which means it involves higher-order derivatives or nonlinear terms compared to the nonlinear Schrödinger equation. It may be specifically designed for certain scenarios or optical systems where the nonlinear Schrödinger equation falls short. The choice between the nonlinear Schrödinger equation and a higher-order version would ultimately depend on the specific requirements and characteristics of the optical system under consideration and whether the added complexity of a higher-order equation is justified by improved accuracy.

The remaining sections are organized in the following order: Section 2 provides the algorithm of the proposed methodology. Section 3 introduces the mathematical analysis of governing equation. Section 4 is used to provide the diversity of optical solitons to governing equation. The bifurcation analysis of planer dynamical system is discussed in Section 5. Section 6 identifies the chaotic behaviour of the perturbed dynamical system using different chaos detecting tools. Section 7 studies the sensitivity analysis of the governing equation. Last section provides the concluding remarks of the study.

2. Algorithm of ${\mathit{G}}^{\mathrm{\prime }}\mathit{/}\mathit{(}\mathit{b}{\mathit{G}}^{\mathrm{\prime }}\mathit{+}\mathit{G}\mathit{+}\mathit{a}\mathit{)}$-expansion method

In the section, we will elucidate the approach known as the $G^{\mathrm{\prime }}/(b{G}^{\mathrm{\prime }}+G+a)$-expansion method by taking into consideration the subsequent partial differential equation: (1) $\mathrm{\wp }(w,w_{xx},w_{tt},\dots .)=0.$(1) We converted Equation (1(1) $\mathrm{\wp }(w,w_{xx},w_{tt},\dots .)=0.$(1) ) to an ODE by using the travelling waves hypothesis as follows: (2) $\mathrm{\Re }(\mathrm{{\rm Y}},{\mathrm{{\rm Y}}}^{\mathrm{\prime }\mathrm{\prime }},\dots .)=0.$(2) Suppose that the altered ordinary differential equation holds the following solution: (3) $\mathrm{{\rm Y}}\left(\varphi \right)=\sum _{i=0}^m{c}_i{H}^i,$(3) The given ODE can be characterized as $H=H\left(\varphi \right)=\frac{G^{\mathrm{\prime }}}{b{G}^{\mathrm{\prime }}+G+a}$, where a, b are constants. The value of m can be determined through the application of the homogeneous balancing principle, while the constants $c_i$ can be resolved. Here, $H=H\left(\varphi \right)$ adheres to the given ordinary differential equation: (4) $H^{\mathrm{\prime }}=\frac{\mathrm{d}H\left(\varphi \right)}{\mathrm{d}\varphi }=(j-l-1)H^2+\frac{2l-j}{b}-\frac{l}{b^2}.$(4) $G=G\left(\varphi \right)$ represents the solution to the subsequent ordinary differential equation (ODE): $G^{\mathrm{\prime }\mathrm{\prime }}=-\frac{j}{b}{G}^{\mathrm{\prime }}-\frac{j}{b^2}G-\frac{l}{b^2}a,$here j and l are real constants. Equation (4(4) $H^{\mathrm{\prime }}=\frac{\mathrm{d}H\left(\varphi \right)}{\mathrm{d}\varphi }=(j-l-1)H^2+\frac{2l-j}{b}-\frac{l}{b^2}.$(4) ) has the following forms of solutions,

Type 1: If $\delta =j^2-4l>0$, then $G=-a+m_1{e}^{\frac{1}{2b}(-j-\sqrt{\delta })\varphi }+m_2{e}^{\frac{1}{2b}(-j+\sqrt{\delta })\varphi },$in this case, $m_1$ and $m_2$ are arbitrary constants that satisfy $a^2+m_1^2+m_2^2\ne 0$.

The function $H=H\left(\varphi \right)$ can be written as below $H=\frac{m_1(j+\sqrt{\delta })+m_2(j-\sqrt{\delta })j^{\frac{\sqrt{\delta }}{b}\varphi }}{b{m}_1(j-2+\sqrt{\delta })+b{m}_2(j-2-\sqrt{\delta })j^{\frac{\sqrt{\delta }}{b}\varphi }}.$Furthermore, it is possible to express H as a function $H=\frac{\begin{array}{c}\left(j\right(m_2-m_1)-\sqrt{\delta }(m_2+m_1\left)\right)\sinh \left(\frac{\sqrt{\delta }}{2b}\varphi \right)\\ +\left(j\right(m_2+m_1)-\sqrt{\delta }(m_2-m_1\left)\right)\cosh \left(\frac{\sqrt{\delta }}{2b}\varphi \right)\end{array}}{\begin{array}{c}b\left(\right(j-2\left)\right(m_2-m_1)-\sqrt{\delta }(m_2+m_1\left)\right)\\ \sinh \left(\frac{\sqrt{\delta }}{2b}\varphi \right)\\ +b\left(\right(j-2\left)\right(m_2+m_1)-\sqrt{\delta }(m_2-m_1\left)\right)\\ \cosh \left(\frac{\sqrt{\delta }}{2b}\varphi \right)\end{array}}.$If $(j-2)(m_2-m_1)-\sqrt{\delta }(m_2+m_1)=0$, then $H=\frac{j-2l}{2b(j-l-1)}-\frac{\sqrt{\delta }}{2b(j-l-1)}\tanh \left(\frac{\sqrt{\delta }}{2b}\varphi \right).$If $(j-2)(m_2+m_1)-\sqrt{\delta }(m_2-m_1)=0$, then $H=\frac{j-2l}{2b(j-l-1)}-\frac{\sqrt{\delta }}{2b(j-l-1)}\coth \left(\frac{\sqrt{\delta }}{2b}\varphi \right).$Type 2: If $\delta =j^2-4l<0$, then $G=-a+e^{\frac{-j}{2b}\varphi }(m_1\cos (\frac{\sqrt{-\delta }}{2b}\varphi )+m_2\sin (\frac{\sqrt{-\delta }}{2b}\varphi \left)\right)$.

In this type, $H=H\left(\varphi \right)$ has the following representation $H=\frac{\begin{array}{c}(j{m}_1-\sqrt{-\delta }{m}_2)\cos \left(\frac{\sqrt{-\delta }}{2b}\varphi \right)\\ +(j{m}_2-\sqrt{-\delta }{m}_1)\sin \left(\frac{\sqrt{-\delta }}{2b}\varphi \right)\end{array}}{\begin{array}{c}b\left(\right(j-2)m_1-\sqrt{-\delta }{m}_2)\cos \left(\frac{\sqrt{-\delta }}{2b}\varphi \right)\\ +b\left(\right(j-2)m_2+\sqrt{-\delta }{m}_1)\sin \left(\frac{\sqrt{-\delta }}{2b}\varphi \right)\end{array}}.$If $(j-2)m_2+\sqrt{-\delta }{m}_1=0$, then $H=\frac{j-2l}{2b(j-l-1)}+\frac{\sqrt{-\delta }}{2b(j-l-1)}\tan \left(\frac{\sqrt{-\delta }}{2b}\varphi \right).$If $(j-2)m_1-\sqrt{-\delta }{m}_2=0$, then $H=\frac{j-2l}{2b(j-l-1)}-\frac{\sqrt{-\delta }}{2b(j-l-1)}\cot \left(\frac{\sqrt{-\delta }}{2b}\varphi \right).$By substituting Equations (3(3) $\mathrm{{\rm Y}}\left(\varphi \right)=\sum _{i=0}^m{c}_i{H}^i,$(3) ) and (4(4) $H^{\mathrm{\prime }}=\frac{\mathrm{d}H\left(\varphi \right)}{\mathrm{d}\varphi }=(j-l-1)H^2+\frac{2l-j}{b}-\frac{l}{b^2}.$(4) ) into the transformed ordinary differential equation and equating the coefficients of $H^i$ to zero, a system of equations is derived. Solving this system yields the values of the unknowns.

The following segment will introduce succinct resolutions of the model being examined, employing the previously mentioned approach.

3. Mathematical analysis of governing equation

Consider the following form of generalized higher-order nonlinear Schrödinger equation (5) $\begin{array}{rl}& A_{\tau }+c_1{A}_{xx}+(c_2|A{|}^{2n}+c_3|A{|}^{4n})A\\ & +c_4\left(\frac{(|A|{)}_{xx}}{|A|}\right)A=0,\end{array}$(5) where $A(x,\tau )$ describe soliton wave profile. In Equation (5(5) $\begin{array}{rl}& A_{\tau }+c_1{A}_{xx}+(c_2|A{|}^{2n}+c_3|A{|}^{4n})A\\ & +c_4\left(\frac{(|A|{)}_{xx}}{|A|}\right)A=0,\end{array}$(5) ), the coefficient $c_1$ physically describes the second order dispersion of the wave, $c_2$ and $c_3$ explain the nonlinear effects due to intensity-dependent refractive index changes, while $c_4$ indicates the higher-order dispersion or nonlinear effects involving spatial variations of the wave amplitude.

In [35], Ankiewicz et al. worked on an extended nonlinear Schrödinger equation with higher-order odd (third order) and even (fourth order) terms with variable coefficients and acquired soliton solutions and approximate rogue wave solutions. In [34], Triki et al. investigated the exact analytical dark soliton solution of the extended nonlinear Schrödinger equation that includes terms accounting for higher-order odd (third order) and even (fourth order) terms in the presence of all the physical parameters. In [37], Arshed et al. studied the bright–dark solitary wave solutions of generalized higher-order nonlinear Schrödinger equation and its applications in optics. In Table 1, we provide a comparison table with some existing Schrodinger equations, from which the superiority of this work is determined.

Table 1. Comparison of previous literature studies and this study.

Study	Equation	Approach	Findings
Triki et al. [34]	Extended NLSE with higher-order odd and even terms	Analytical	Parametric stability and dark soliton
Ankiewicz et al. [35]	Extended NLSE with higher-order odd and even terms	Analytical	Rogue wave solutions
Ozisik et al. [36]	Higher-order NLSE having Kudryashov nonlinear refractive index	Analytical	Optical soliton solutions
Arshad et al. [37]	Generalized higher-order NLSE	Analytical	Solitary wave solutions
Younas et al. [38]	Generalized higher-order NLSE	Analytical	Optical solitons and other solutions in fibre optics
This study	Generalized higher-order NLSE	Analytical	Optical solitons, Dynamical behaviour of light wave promulgation, Chaotic and bifurcation analysis, Multistability response

To solve Equation (5(5) $\begin{array}{rl}& A_{\tau }+c_1{A}_{xx}+(c_2|A{|}^{2n}+c_3|A{|}^{4n})A\\ & +c_4\left(\frac{(|A|{)}_{xx}}{|A|}\right)A=0,\end{array}$(5) ) consider the transformation (6) $A(x,\tau )=u\left(\varphi \right)\exp \left(i\alpha (x,\tau )\right),$(6) here ϕ represents the travelling wave variable $\varphi =kx-w\tau .$Where, $u\left(\varphi \right)$ shows the amplitude of wave, while w represents the soliton's speed. Moreover, $\alpha (x,\tau )$ represents the phase element, delineated as follows: $\alpha (x,\tau )=sx+\delta \tau ,$here, s denote the wave frequency, while δ represents the wave number. Upon inserting Equation (6(6) $A(x,\tau )=u\left(\varphi \right)\exp \left(i\alpha (x,\tau )\right),$(6) ) into Equation (5(5) $\begin{array}{rl}& A_{\tau }+c_1{A}_{xx}+(c_2|A{|}^{2n}+c_3|A{|}^{4n})A\\ & +c_4\left(\frac{(|A|{)}_{xx}}{|A|}\right)A=0,\end{array}$(5) ) and combining the real and imaginary components, the velocity of the soliton can be determined. $w=2c_{1}sk,$and the following ODE (7) $k(c_1+c_4)u^{\mathrm{\prime }\mathrm{\prime }}-(\delta +c_1{s}^2)u+c_3{u}^{2n+1}+c_4{u}^{4n+1}=0.$(7) Balancing $u^{\mathrm{\prime }\mathrm{\prime }}$ with $u^{4n+1}$ in Equation (7(7) $k(c_1+c_4)u^{\mathrm{\prime }\mathrm{\prime }}-(\delta +c_1{s}^2)u+c_3{u}^{2n+1}+c_4{u}^{4n+1}=0.$(7) ) provides $N=\frac{1}{2n}$. To get closed-form solution, consider further transformation as (8) $u=u^{\frac{1}{2n}}.$(8) Equation (7(7) $k(c_1+c_4)u^{\mathrm{\prime }\mathrm{\prime }}-(\delta +c_1{s}^2)u+c_3{u}^{2n+1}+c_4{u}^{4n+1}=0.$(7) ) can be represented as (9) $\begin{array}{rl}& (\frac{a(1-n)}{4n^2}-\frac{a}{4n})(u^{\mathrm{\prime }}{)}^2+\frac{a}{2n}u{u}^{\mathrm{\prime }\mathrm{\prime }}-e{u}^2\\ & +c{u}^3+\mathrm{d}{u}^4=0,\end{array}$(9) where $a=k(c_1+c_4)$, $e=(\delta +c_1{s}^2)$, $c=c_3$ and $d=c_4$.

4. Extraction of solutions using the proposed technique

Subsequently, we proceed by using the $G^{\mathrm{\prime }}/(b{G}^{\mathrm{\prime }}+G+a)$-expansion method to address Equation (9(9) $\begin{array}{rl}& (\frac{a(1-n)}{4n^2}-\frac{a}{4n})(u^{\mathrm{\prime }}{)}^2+\frac{a}{2n}u{u}^{\mathrm{\prime }\mathrm{\prime }}-e{u}^2\\ & +c{u}^3+\mathrm{d}{u}^4=0,\end{array}$(9) ), aiming to unveil solutions reminiscent of solitons. Guided by the prescribed protocol tailored for the scenario where m = 1, the anticipated outcome for Equation (9(9) $\begin{array}{rl}& (\frac{a(1-n)}{4n^2}-\frac{a}{4n})(u^{\mathrm{\prime }}{)}^2+\frac{a}{2n}u{u}^{\mathrm{\prime }\mathrm{\prime }}-e{u}^2\\ & +c{u}^3+\mathrm{d}{u}^4=0,\end{array}$(9) ) takes shape as follows: (10) $u\left(\varphi \right)=c_0+c_{1}H\left(\varphi \right),$(10) Here, the variables $c_0$ and $c_1$ represent constants that are yet to be determined. Employing the $G^{\mathrm{\prime }}/(b{G}^{\mathrm{\prime }}+G+a)$-expansion strategy leads to the subsequent solutions. $c_0=\frac{c_1}{b}\sqrt{-\frac{-b^{2}j+2b^{2}l-l}{j-l-1}},c_1=c_1.$According to Type 1, for $\delta =j^2-4k>0$, we have (11) $u_1(x,\tau )=(\frac{c_1}{b}\sqrt{-\frac{-b^{2}j+2b^{2}l-l}{j-l-1}}+2)H,$(11) where $H=\frac{\begin{array}{c}\left(j\right(m_2-m_1)-\sqrt{\delta }(m_2+m_1\left)\right)\sinh \left(\frac{\sqrt{\delta }}{2b}\varphi \right)\\ +\left(j\right(m_2+m_1)-\sqrt{\delta }(m_2-m_1\left)\right)\cosh \left(\frac{\sqrt{\delta }}{2b}\varphi \right)\end{array}}{\begin{array}{c}b\left(\right(l-2\left)\right(m_2-m_1)-\sqrt{\delta }(m_2+m_1\left)\right)\\ \sinh \left(\frac{\sqrt{\delta }}{2b}\varphi \right)\\ +b\left(\right(l-2\left)\right(m_2+m_1)-\sqrt{\delta }(m_2-m_1\left)\right)\\ \cosh \left(\frac{\sqrt{\delta }}{2b}\varphi \right)\end{array}}.$Moreover, if $(l-2)(m_2-m_1)-\sqrt{\delta }(m_2+m_1)=0$, then Equation (11(11) $u_1(x,\tau )=(\frac{c_1}{b}\sqrt{-\frac{-b^{2}j+2b^{2}l-l}{j-l-1}}+2)H,$(11) ) becomes (12) $u_2(x,\tau )=(\frac{c_1}{b}\sqrt{-\frac{-b^{2}j+2b^{2}l-l}{j-l-1}}+2)H,$(12) where $H=\frac{j-2l}{2b(j-l-1)}-\frac{\sqrt{\delta }}{2b(j-l-1)}\tanh \left(\frac{\sqrt{\delta }}{2b}\varphi \right).$Moreover, if $(l-2)(m_2+m_1)-\sqrt{\delta }(m_2-m_1)=0$, then Equation (11(11) $u_1(x,\tau )=(\frac{c_1}{b}\sqrt{-\frac{-b^{2}j+2b^{2}l-l}{j-l-1}}+2)H,$(11) ) becomes (13) $u_3(x,\tau )=(\frac{c_1}{b}\sqrt{-\frac{-b^{2}j+2b^{2}l-l}{j-l-1}}+2)H,$(13) where $H=\frac{j-2l}{2b(j-l-1)}-\frac{\sqrt{\delta }}{2b(j-l-1)}\coth \left(\frac{\sqrt{\delta }}{2b}\varphi \right).$where $\varphi =kx-w\tau$ (Figures 1 and 2).

According to Type 2, for $\delta =j^2-4l<0$, we have (14) $u_4(x,\tau )=(\frac{c_1}{b}\sqrt{-\frac{-b^{2}j+2b^{2}l-l}{j-l-1}}+2)H,$(14) where $H=\frac{\begin{array}{c}(j{m}_1-\sqrt{-\delta }{m}_2)\cos \left(\frac{\sqrt{-\delta }}{2b}\varphi \right)\\ +(j{m}_2-\sqrt{-\delta }{m}_1)\sin \left(\frac{\sqrt{-\delta }}{2b}\varphi \right)\end{array}}{\begin{array}{c}b\left(\right(j-2)m_1-\sqrt{-\delta }{m}_2)\cos \left(\frac{\sqrt{-\delta }}{2b}\varphi \right)\\ +b\left(\right(j-2)m_2+\sqrt{-\delta }{m}_1)\sin \left(\frac{\sqrt{-\delta }}{2b}\varphi \right)\end{array}}.$Moreover, if $(j-2)m_2+\sqrt{-\delta }{m}_1=0$, then Equation (14(14) $u_4(x,\tau )=(\frac{c_1}{b}\sqrt{-\frac{-b^{2}j+2b^{2}l-l}{j-l-1}}+2)H,$(14) ) becomes (15) $u_5(x,\tau )=(\frac{c_1}{b}\sqrt{-\frac{-b^{2}j+2b^{2}l-l}{j-l-1}}+2)H,$(15) where $H=\frac{j-2l}{2b(j-l-1)}+\frac{\sqrt{-\delta }}{2b(j-l-1)}\tan \left(\frac{\sqrt{-\delta }}{2b}\varphi \right)).$Moreover, if $(j-2)m_1-\sqrt{-\delta }{m}_2=0$, then Equation (14(14) $u_4(x,\tau )=(\frac{c_1}{b}\sqrt{-\frac{-b^{2}j+2b^{2}l-l}{j-l-1}}+2)H,$(14) ) becomes (16) $u_6(x,\tau )=(\frac{c_1}{b}\sqrt{-\frac{-b^{2}j+2b^{2}l-l}{j-l-1}}+2)H,$(16) where (17) $H=\frac{j-2l}{2b(j-l-1)}-\frac{\sqrt{-\delta }}{2b(j-l-1)}\cot \left(\frac{\sqrt{-\delta }}{2b}\varphi \right).$(17) where $\varphi =kx-w\tau$ (Figure 3).

Figure 1. Visualization of $u_2(x,\tau )$ using a 3D, 2D and contour plots with parameters b = 1, j = 2 and l = 0.5.

Figure 2. Visualization of $u_3(x,\tau )$ using a 3D, 2D and contour plots with parameters b = 1, j = 2 and l = 0.5.

5. Bifurcation of the governing model

In this subsequent section, the model is studied qualitatively by utilizing the planar dynamical system of Equation (9(9) $\begin{array}{rl}& (\frac{a(1-n)}{4n^2}-\frac{a}{4n})(u^{\mathrm{\prime }}{)}^2+\frac{a}{2n}u{u}^{\mathrm{\prime }\mathrm{\prime }}-e{u}^2\\ & +c{u}^3+\mathrm{d}{u}^4=0,\end{array}$(9) ), that can be written as below (18) $\begin{array}{rl}& u^{\mathrm{\prime }}=v\\ & v^{\mathrm{\prime }}=\frac{{\delta }_1{v}^2+{\delta }_2{u}^2-{\delta }_3{u}^3-{\delta }_4{u}^4}{u},\end{array}$(18) here ${\delta }_1=\frac{-(1-n)}{2n}+\frac{1}{2},{\delta }_2=\frac{2en}{a},{\delta }_3=\frac{2cn}{a},{\delta }_4=\frac{2\mathrm{d}n}{a}.$However, the above system isn't Hamiltonian system. By using Equation (18(18) $\begin{array}{rl}& u^{\mathrm{\prime }}=v\\ & v^{\mathrm{\prime }}=\frac{{\delta }_1{v}^2+{\delta }_2{u}^2-{\delta }_3{u}^3-{\delta }_4{u}^4}{u},\end{array}$(18) ), we have (19) $\frac{\mathrm{d}{v}^2}{\mathrm{d}u}=\frac{2{\delta }_1{v}^2+2{\delta }_2{u}^2-2{\delta }_3{u}^3-2{\delta }_4{u}^4}{u},$(19) Clearly the point u = 0 shows singular point of Equation (19(19) $\frac{\mathrm{d}{v}^2}{\mathrm{d}u}=\frac{2{\delta }_1{v}^2+2{\delta }_2{u}^2-2{\delta }_3{u}^3-2{\delta }_4{u}^4}{u},$(19) ). Hence, this equation has the following solution under these conditions (20) $v^2=a_0{u}^4+a_1{u}^3+a_2{u}^2,$(20) where $a_0=\frac{2{\delta }_1^2{\delta }_4-5{\delta }_1{\delta }_4+3{\delta }_4}{({\delta }_1-2)({\delta }_1-1)(2{\delta }_1-3)}$, $a_1=\frac{2{\delta }_1^2{\delta }_3-6{\delta }_1{\delta }_3+4{\delta }_3}{({\delta }_1-2)({\delta }_1-1)(2{\delta }_1-3)}$,$a_2=\frac{-2{\delta }_1^2{\delta }_2+7{\delta }_1{\delta }_2-6{\delta }_2}{({\delta }_1-2)({\delta }_1-1)(2{\delta }_1-3)}$, we get (21) $(u^{\mathrm{\prime }}{)}^2-[a_0{u}^4+a_1{u}^3+a_2{u}^2]=0,$(21) from Equation (20(20) $v^2=a_0{u}^4+a_1{u}^3+a_2{u}^2,$(20) ), we have (22) $H(v,u)=v^2-[a_0{u}^4+a_1{u}^3+a_2{u}^2],$(22) Next, bifurcation analysis will be performed on the aforementioned model by utilizing the discrimination theory, which leads to Equation (22(22) $H(v,u)=v^2-[a_0{u}^4+a_1{u}^3+a_2{u}^2],$(22) ) along with potential energy as: (23) $U\left(u\right)=-[a_0{u}^4+a_1{u}^3+a_2{u}^2],$(23) further we get (24) $U^{\mathrm{\prime }}\left(u\right)=-4a_0{u}^3-3a_1{u}^2-2a_{2}u.$(24) We get (25) $U^{\mathrm{\prime }}\left(u\right)=-4a_0[u^3+{\rho }_0{u}^2+{\rho }_{1}u].$(25) Here ${\rho }_0=\frac{3a_1}{4a_0}$, ${\rho }_1=\frac{2a_2}{4a_0}$.

Figure 3. Visualization of $u_5(x,\tau )$ using a 3D, 2D and contour plots with parameters b = 1, j = 1 and l = 1.

Let $J(u,v)$ linearized matrix on critical points $(u,v)$, known as $\mathit{Jacobian}$. The determinant is: $\left|\begin{array}{cc}0& 1\\ -4a_0(3u^2+2{\rho }_{0}u+{\rho }_1)& 0\end{array}\right|,$So we get $J(u,v)=4a_0(3u^2+2{\rho }_{0}u+{\rho }_1),$therefore, $(u,v)$ would be saddle if $J(u,v)<0$, but if $J(u,v)>0$, centre will be shown, whereas it would be cusp if $J(u,v)=0$.

By introducing the discriminant for polynomial (26) $\mathrm{\Delta }={\rho }_0^2{\rho }_1^2-4{\rho }_1^3.$(26) Let's discuss the possible outcome.

Case a:

Let $\mathrm{\Delta }=0$ and ${\rho }_0>0$ and ${\rho }_1=0$, $U^{\mathrm{\prime }}\left(u\right)=-4a_0\left[\right(u-p{)}^2(u-q)].$Two critical points $(p,0)$ as well as $(q,0)$ exist. Where $(p,0)$ showed cusp point although $(q,0)$ is the centre for $a_0>0$ and saddle for $a_0<0$. The behaviour is showed in Figure 4 when ${\rho }_1=0$, $a_0=\frac{1}{6}$ and ${\rho }_0=1$. We get p = 0 and q = −1.

Figure 4. Global phase portrait for Case a.

Case b:

If $\mathrm{\Delta }=0$ and ${\rho }_1=0$, ${\rho }_0<0$, $U^{\mathrm{\prime }}\left(u\right)=-4a_0\left[\right(u-p{)}^2(u-q)].$Two critical points $(p,0)$ as well as $(q,0)$ exist in such situation. Where $(p,0)$ shows cuspidal and $(q,0)$ is a saddle if $a_0<0$ and centre for $a_0>0$. In Figure 5, phase portrait is depicted for ${\rho }_1=0$, $a_0=-\frac{1}{6}$ and ${\rho }_0=-1$. We get p = 0 and q = 1.

Figure 5. Global phase portrait for Case b.

Case c:

If $\mathrm{\Delta }>0$, ${\rho }_1>0$, ${\rho }_0>0$, $U^{\mathrm{\prime }}\left(u\right)=-4a_0\left[\right(u-p\left)\right(u-q\left)\right(u-n\left)\right].$In such assumption, $(p,0)$ will centre point, $(q,0)$ but $(n,0)$ are saddle points for $a_0<0$. Their behaviours are displayed in Figure 6 by a mathematical tool selecting $a_0=-\frac{1}{6}$, ${\rho }_0=3$, ${\rho }_1=1$. Where we get p = 0, q = −2.6, and n = −0.4.

Figure 6. Global phase portrait for Case c.

Case d:

If $\mathrm{\Delta }>0$, ${\rho }_1<0$, ${\rho }_0<0$, $U^{\mathrm{\prime }}\left(u\right)=-4a_0\left[\right(u-p\left)\right(u-q\left)\right(u-n\left)\right],$here three equilibrium points are attained. Here, $(p,0)$ will be saddle point, $(q,0)$ whereas $(n,0)$ are central points for $a_0>0$. In Figure 7, the behaviour is depicted by taking $a_0=\frac{1}{6}$, ${\rho }_0=-3$, ${\rho }_1=-1$. Here we get p = 0, q = 3.29 and n = −0.3.

Figure 7. Global phase portrait for Case d.

Case e:

If $\mathrm{\Delta }<0$, (27) $U^{\mathrm{\prime }}\left(u\right)=-4a_0\left[\right(u-p\left)\right[(u-q{)}^2+m^2]].$(27) Hence $(p,0)$ will be centre for $a_0>0$ and saddle for $a_0<0$. This will be unique critical point for this case.

6. Chaotic behaviour of the suggested model

In this segment, our objective is to examine the tumultuous behaviour of the system delineated by Equation (1(1) $\mathrm{\wp }(w,w_{xx},w_{tt},\dots .)=0.$(1) ). Consequently, through the inclusion of the periodic component, the structure of the two-dimensional system (18(18) $\begin{array}{rl}& u^{\mathrm{\prime }}=v\\ & v^{\mathrm{\prime }}=\frac{{\delta }_1{v}^2+{\delta }_2{u}^2-{\delta }_3{u}^3-{\delta }_4{u}^4}{u},\end{array}$(18) ) can be expressed as follows: (28) $\{\begin{array}{ll}& u^{\mathrm{\prime }}=v,\\ & v^{\mathrm{\prime }}=-4a_0[u^3+{\rho }_0{u}^2+{\rho }_{1}u]+g_0\cos \left(w\right),\\ & w^{\mathrm{\prime }}=\rho ,\end{array}$(28) where $w=\rho \xi$. The purpose of introducing the perturbation term by applying external force is to investigate the dynamical behaviour of the system of solution Equation (28(28) $\{\begin{array}{ll}& u^{\mathrm{\prime }}=v,\\ & v^{\mathrm{\prime }}=-4a_0[u^3+{\rho }_0{u}^2+{\rho }_{1}u]+g_0\cos \left(w\right),\\ & w^{\mathrm{\prime }}=\rho ,\end{array}$(28) ). We examine how the frequency and amplitude of the perturbation term influence the behaviour of the system of solution Equation (28(28) $\{\begin{array}{ll}& u^{\mathrm{\prime }}=v,\\ & v^{\mathrm{\prime }}=-4a_0[u^3+{\rho }_0{u}^2+{\rho }_{1}u]+g_0\cos \left(w\right),\\ & w^{\mathrm{\prime }}=\rho ,\end{array}$(28) ). In order to identify the chaotic dynamics in the system (28(28) $\{\begin{array}{ll}& u^{\mathrm{\prime }}=v,\\ & v^{\mathrm{\prime }}=-4a_0[u^3+{\rho }_0{u}^2+{\rho }_{1}u]+g_0\cos \left(w\right),\\ & w^{\mathrm{\prime }}=\rho ,\end{array}$(28) ), we utilize various methods for detecting chaos that are beneficial in determining chaotic behaviour in the given system. These methods encompass the use of a 2D-phase portrait, 3D-phase portrait, analysis of time series data, examination of Lyapunov exponents, and consideration of multistability.

• Figure 8 illustrates 3D visualization, time analysis graph, and 2D plots for the given parameters: $a_0=1$, ${\rho }_0=-\frac{1}{4}$, ${\rho }_1=\frac{1}{2}$, $g_0=2.5$, and $\rho =\pi$. It's observed that the perturbed system (28(28) $\{\begin{array}{ll}& u^{\mathrm{\prime }}=v,\\ & v^{\mathrm{\prime }}=-4a_0[u^3+{\rho }_0{u}^2+{\rho }_{1}u]+g_0\cos \left(w\right),\\ & w^{\mathrm{\prime }}=\rho ,\end{array}$(28) ) demonstrates quasi-periodic chaotic behaviour.

• Figure 9 encompasses the presentation of 3D, 2D, and time series plots corresponding to the parameters: $a_0=1$, ${\rho }_0=\frac{1}{4}$, ${\rho }_1=\frac{3}{4}$, $g_0=3.5$, and $\rho =2\pi$. Evidently, the system (28(28) $\{\begin{array}{ll}& u^{\mathrm{\prime }}=v,\\ & v^{\mathrm{\prime }}=-4a_0[u^3+{\rho }_0{u}^2+{\rho }_{1}u]+g_0\cos \left(w\right),\\ & w^{\mathrm{\prime }}=\rho ,\end{array}$(28) ) exhibits quasi-periodic behaviour.

• Figure 10 showcases the visualization of 3D, 2D phase plots, and a time series graph with the following parameter values: $a_0=1$, ${\rho }_0=\frac{5}{4}$, ${\rho }_1=\frac{3}{4}$, $g_0=4$, and $\rho =2\pi$. The noteworthy observation is that the altered dynamical system demonstrates a quasi-periodic pattern.

• In Figure 11, we determine the existence of chaos in system (28(28) $\{\begin{array}{ll}& u^{\mathrm{\prime }}=v,\\ & v^{\mathrm{\prime }}=-4a_0[u^3+{\rho }_0{u}^2+{\rho }_{1}u]+g_0\cos \left(w\right),\\ & w^{\mathrm{\prime }}=\rho ,\end{array}$(28) ) through Lyapunov exponents. Three Lyapunov exponent values ${\lambda }_1=0.224958,{\lambda }_2=-0.224958$ and ${\lambda }_1=0$ are to be found and shown graphically, see Figure 11. From Figure 11, we can easily observe that system (28(28) $\{\begin{array}{ll}& u^{\mathrm{\prime }}=v,\\ & v^{\mathrm{\prime }}=-4a_0[u^3+{\rho }_0{u}^2+{\rho }_{1}u]+g_0\cos \left(w\right),\\ & w^{\mathrm{\prime }}=\rho ,\end{array}$(28) ) has chaotic behaviour. Similarly, Figure 12 is presented for Lyapunov exponents using different parameter values rather than Figure 11.

• In Figure 13, we determine the existence of chaos in system (28(28) $\{\begin{array}{ll}& u^{\mathrm{\prime }}=v,\\ & v^{\mathrm{\prime }}=-4a_0[u^3+{\rho }_0{u}^2+{\rho }_{1}u]+g_0\cos \left(w\right),\\ & w^{\mathrm{\prime }}=\rho ,\end{array}$(28) ) through multistability at different initial conditions. Multistability guarantees that for different initial value problem the considered equation has stable solutions. Normally we take these initial value problems closer to each other so that in a single frame stability can be shown. It also shows that for a minor change in the initial value problems how the considered equation reacts. In Figure 13, we detect the existence of chaotic behaviour of the system (28(28) $\{\begin{array}{ll}& u^{\mathrm{\prime }}=v,\\ & v^{\mathrm{\prime }}=-4a_0[u^3+{\rho }_0{u}^2+{\rho }_{1}u]+g_0\cos \left(w\right),\\ & w^{\mathrm{\prime }}=\rho ,\end{array}$(28) ) at three different initial conditions such as $(0.1,0.1,0),(0.2,0.1,0),(0.3,0.1,0)$ through multistability analysis.

• In Figure 14, we identify the chaotic behaviour through the chaotic attractor using $a_0=1,{\rho }_0=0.25,{\rho }_1=0.75,g_0=2.5,\rho =3.5$ at (0.3,0,0.3). The system (28(28) $\{\begin{array}{ll}& u^{\mathrm{\prime }}=v,\\ & v^{\mathrm{\prime }}=-4a_0[u^3+{\rho }_0{u}^2+{\rho }_{1}u]+g_0\cos \left(w\right),\\ & w^{\mathrm{\prime }}=\rho ,\end{array}$(28) ) shows the stable chaotic attractor.

• In Figure 15, we identify the chaotic behaviour through the basin of attractor using $a_0=1,{\rho }_0=0.25,{\rho }_1=0.75,g_0=2.5,\rho =3.5$ at initial condition (0.3,0,0.3). In figure, red line represents the chaotic attractor and other two lines are stable at (0.3,0,0.3).

Figure 8. Depicting the system with parameters $a_0=1$, ${\rho }_0=-\frac{1}{4}$, ${\rho }_1=\frac{1}{2}$, $g_0=2.5$ and $\xi =\pi$, at initial condition (0.3, 0.3). (a) 2D-Phase profile with perturbed term, (b) 3D-Phase profile with perturbed term, (c) Poincare section and (d) Time Analysis.

Figure 9. Depicting the system with parameters $a_0=1$, ${\rho }_0=\frac{1}{4}$, ${\rho }_1=\frac{3}{4}$, $g_0=3.5$ and $\xi =2\pi$, at initial condition (0.3, 0.3). (a) 2D-Phase profile with perturbed term, (b) 3D-Phase profile with perturbed term, (c) Poincare section and (d) Time Analysis.

Figure 10. Depicting the system with parameters $a_0=1$, ${\rho }_0=\frac{5}{4}$, ${\rho }_1=\frac{3}{4}$, $g_0=4$ and $\xi =2\pi$, at initial condition (0.3, 0.3). (a) 2D-Phase profile with perturbed term, (b) 3D-Phase profile with perturbed term, (c) Poincare section and (d) Time Analysis.

Figure 11. Detection of chaos via Lyapunov exponents at $a_0=1,{\rho }_0=-0.25,{\rho }_1=0.5,g_0=2.5,\rho =\pi$.

Figure 12. Detection of chaos via Lyapunov exponents at $a_0=1,{\rho }_0=0.25,{\rho }_1=0.75,g_0=2.5,\rho =\pi$.

Figure 13. Detection of chaos via multistability at $a_0=1,{\rho }_0=-0.25,{\rho }_1=0.5,g_0=2.5,\rho =\pi$.

Figure 14. Detection of chaos via attractor using $a_0=1,{\rho }_0=0.25,{\rho }_1=0.75,g_0=2.5,\rho =3.5$ at (0.3,0,0.3).

Figure 15. Detection of chaos via basin of attractor using $a_0=1,{\rho }_0=0.25,{\rho }_1=0.75,g_0=2.5,\rho =3.5$ at (0.3,0,0.3).

7. Sensitivity analysis

This part elucidates the delicate nature of the proposed system. Sensitivity evaluation involves assessing the system's responsiveness. When minor changes to initial conditions only slightly influence the system, sensitivity is low. However, if even small adjustments to the baseline conditions lead to substantial system changes, sensitivity is high. Figure 15 depicts the sensitivity of the wave profile in the dynamic structural system under two distinct initial conditions. The system is sensitive, as can be seen in Figures 15(a) and (b), where a change in the amplitude structure of the waves can be seen as the interval progresses from 0 to 20. This indicates that these curves are not overlapping in Figure  15(a), but that the system is sensitive because when we alter a small value in the initial condition, a significant difference can be seen (Figure 16).

Figure 16. (a) Shows the analysis of the system letting the initial conditions $(1,0.10)$ and $(1,0.15)$ with n = 1.5, $c_4=0.4$, $\delta =0.5$, $c_1=c_3=k=1$. (b) Here taking initial conditions $(1,0.15)$ and $(1,0.15)$ which shows the overlapping between the curves.

8. Conclusion

In this study, we have investigated the generalized higher-order nonlinear Schrödinger equation behaviour in predicting the time-based evolution of disturbances in marginally stable or unstable media, specifically in the context of light wave propagation. The $G^{\mathrm{\prime }}/(b{G}^{\mathrm{\prime }}+G+a)$-expansion method is utilized to derive precise solutions for this model. The analysis takes a multifaceted approach, encompassing bifurcation behaviour, quasi-periodic oscillations, chaotic motion and sensitivity assessment of the equation. Bifurcation analysis involves scrutinizing the planar dynamical system at equilibrium points. Moreover, by introducing an external periodic perturbation, a perturbed dynamic system is formed, whose quasi-periodic and chaotic attributes are identified through chaos detection tools like 3D phase portraits, 2D phase plots, time series analysis, Lyapunov exponents and multistability. These techniques reveal the perturbed system's unpredictable patterns and chaotic behaviour. Lastly, sensitivity analysis across three initial conditions underscores the equation's high sensitivity to even minor changes. Overall, these findings are both innovative and significant in understanding the evolution of disruptions in less stable media. They also indicate the method's potential to unravel dynamics in various nonlinear models within fields like engineering and fibre optics, offering new avenues for research and insight.

Future investigation into the generalized higher-order nonlinear Schröodinger equation encompasses numerous facets. One particular focus area pertains to the formulation of higher order smooth positon and breather positon solutions, which possess a broader scope compared to the solutions of the constituent equations. Another aspect involves the examination of the dynamical depiction of the higher-order modulation instability for fourth-order nonlinear Schröodinger equations, which mirrors a comparable geometric interpretation as that of the fourth-order equation. These avenues of research will contribute to a more profound comprehension of the behaviour and characteristics of the generalized higher-order nonlinear Schröodinger equation.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Availability of data and materials

The data that support the findings of this study are available from the corresponding author upon reasonable request.