Bifurcation of low-frequency ion-acoustic nonlinear structures and chaos with fractional effect in Non-Maxwellian magnetoplasmas

Abstract

The dynamics of ion-acoustic waves (IAWs) in nonextensive magnetoplasmas with the two-temperature electrons are studied. The plasma model equations are transformed into the modified Korteweg–de Vries equation employing the reductive perturbation method. The Grunwald–Letnikov definition and the Newton–Leipnik algorithm are addressed to study the forced dynamical system. Through the fourth-order Runge–Kutta scheme, the fractional forced system is solved numerically. It is seen that the quasi-period route to chaos is impacted by the fractional parameter α in a forced system. It is also found that both regular and irregular oscillations appear due to the variations in the fractional-order term. Phase portrait plots, time series, and bifurcation diagrams corroborate these results. The relevance of this investigation can help to understand the problem of anomalous transport due to the drift wave turbulence state of the plasma.

Introduction

Nowadays, the low-frequency wave, also known as ion-acoustic wave (IAW), is widely studied in plasma, and its propagation is connected to plasma confinement.^1,2 Usually caused by a concentration gradient, IAWs are acoustic-type waves in which lighter species, such as electrons, provide the restoring force, while heavy ionic species supply inertia.^3–5 Ion-acoustic solitons (IAS), also known as nonlinear IAWs, are a typical example of solitons in plasmas and exhibit solitonic behavior. Space missions and observatories, including S3-3, GEOTAIL, Freja, POLAR, FAST, and Viking, have reported the presence of nonlinear solitary waves in space.^6–8

The existence of two-electron species (hot and cold) in the plasma has been shown to alter how an IAW forms and spreads.^9–11 It was observed that the amplitude of electron-acoustic soliton and IAS in two-electron temperature (TET) increased in both cases as the temperature of the hot electron increased.^12–14 Many intriguing physical problems have recently been investigated in non-Maxwellian TET plasmas through the analysis of various nonlinear equations.^15–17 Researchers investigated Langmuir and IAS propagation properties in a non-magnetized plasma with TET. It was observed that these solitons traveled faster than sound and were characterized by density rarefactions. ¹⁸

Thermal plasma systems are those in which the particles are in thermal equilibrium and are typically described by the Maxwellian velocity distribution. Electrons in different regions of space plasmas have been found by multiple satellite missions to behave differently from Maxwellian electrons.^19,20 One of the distributions widely used today to describe a plasma out of thermodynamic equilibrium is the nonextensive distribution. This distribution, initially recommended by Renyi, ²¹ then later confirmed by Tsalis, ²² made it possible to describe the statistics of particles in space successfully. The degree of nonextensivity in this distribution is quantified by the entropic index q. It is important to note that super-extensivity (sub-extensivity) is denoted by q < 1 (q > 1) and Maxwellian behavior by q = 1. As previously noted, energetic electrons with non-Maxwellian components are seen in a range of astrophysical plasma domains. ²³ Recently, the super positron-acoustic rogue waves using the nonextensive distribution of electrons and positrons have been investigated in a magnetized plasma. ²⁴ It has been shown that increasing the nonextensive term leads to decreasing the amplitude of the super positron-acoustic rogue waves. Furthermore, the nonextensive parameter influences the stability domain in the collisional electronegative plasma, as demonstrated by Alhejaili et al. ²⁵

The influence of the magnetic field (MF) on the plasma environment has remained of significant interest in observational, theoretical, and numerical research since most astrophysical plasmas are magnetized. It is well established that the MF alters the plasma’s characteristics, causing non-identical propagation of electrostatic waves in directions parallel and perpendicular to the MF while also permitting oblique propagation.^26–28 The nonlinear dynamics of dense plasmas composed of electrons, positrons, ions, and charged dust particles under the combined influence of magnetic and gravitational fields have been investigated using a quantum hydrodynamic model. ²⁹ Further research has examined electron-acoustic waves (EAWs) in magnetized plasmas featuring nonthermal Kappa-distributed ions. ³⁰ The derivation of the Korteweg–de Vries (KdV) equation from this model aids in identifying distinct EAW propagation modes in space and laboratory plasmas. The properties of magnetosonic shock waves and solitons have also been explored in Fermi plasmas under quasi-periodic perturbations. ³¹ This analysis revealed how nonlinear wave structures evolve in such environments. More recently, a study employing reductive perturbation method (RPM) investigated relativistic EAWs in a thermal magnetoplasma with Kappa-distributed ions. ³² Properties of dust magnetosonic solitons in magnetized electron–ion–dust plasmas were systematically investigated by Singla et al.. ³³ Through derivation of the KdV equation, the authors conducted comprehensive numerical analyses examining how key plasma parameters influence soliton amplitude, width, and propagation characteristics. A subsequent study ³⁴ developed a relativistic degenerate magneto-rotating quantum plasma model. This work revealed that nucleus-acoustic envelope solitons exhibit strong dependence on three critical factors which includes finite-temperature effects, rotational frequency, and mass density ratios between nuclear species. Chandra et al. ³⁵ employed RPM to investigate IAW dynamics in magnetized multicomponent plasmas. Their methodology involved derivation of a nonlinear evolutionary equation, numerical solution via fourth-order Runge–Kutta integration, and stability analysis through phase-space examination of the resulting ordinary differential equations. This approach provided quantitative criteria for structural stability of stationary wave solutions. In astrophysical contexts, recent work ³⁶ has examined magneto-acoustic wave phenomena in dense stellar plasmas containing electron–positron–ion triplets under strong MF. Sarkar et al. ³⁷ conducted a comprehensive investigation of Rayleigh–Taylor instability in ultra-dense astrophysical environments, particularly neutron stars. Employing a quantum hydrodynamic framework, the study systematically analyzed the instability growth rate and its characteristic properties. In a related study, a novel quantum two-stream instability mode was discovered in relativistically degenerate, magnetized dense plasmas. ³⁸ This research provided detailed analysis of several critical factors influencing wave stability, including: quantum Bohm potential effects, MF strength and orientation, relative streaming velocities, and particle number density distributions. Furthermore, when strong enough, MF considerably influences the dynamics of electrostatic waves, making them complex and having a complex behavior in the phase space. As a result, magnetized plasmas can sustain various nonlinear structures, including periodic, quasi-periodic, and chaotic oscillations, in addition to solitary waves.^39–41

Numerous disciplines, including biology, mathematics, physics, chemistry, engineering, hydrology, finance, and social sciences, have shown great interest in fractional calculus.^42–51 For example, Yusry et al. ⁵² investigated the fractal Duffing oscillator under periodic forcing using a non-perturbative analytical approach. Their methodology successfully reduced the order of the governing differential equation and derived periodic solutions via Galerkin’s technique. Ali et al. ⁵³ developed an innovative mesh-free computational approach utilizing moving kriging (MK) interpolation to solve the magnetohydrodynamic (MHD) flow of fractional Maxwell fluids, employing finite difference methods for temporal discretization and MK shape functions for spatial resolution. In a complementary engineering application addressing elevator safety, recent research ⁵⁴ has designed a bio-inspired fractal buffer system featuring hierarchical vibration-absorbing metamaterials that mimic gecko pad structures, demonstrating significant potential for impact energy dissipation and damage mitigation in elevator accidents while offering broader implications for the development of advanced safety buffers in mechanical systems. A novel methodology for characterizing fractal networks in chaotic systems has been developed for microstructure analysis of robot-laser-hardened steels, incorporating multiple laser beam angles. ⁵⁵ Models with fractional derivatives better describe the dynamic behaviors of complex physical systems. Indeed, fractional derivative models are more precise and better adapted to experimental data than integer-order models. Fractional derivatives also help maintain the memory effect in real-world systems, in which the flow of a particle at a given time is influenced by its prior history. ⁵⁶ Different definitions of fractional derivatives have been introduced and utilized in literature in previous years, for instance, truncated M-fractional derivative, ⁵⁷ conformable fractional, ⁵⁸ Reimann–Liouville, ⁵⁹ Caputo fractional, ⁶⁰ Grunwald–Lentnikov definition, ⁶¹ β time-fractional, ⁶² etc. To approximate the relevant derivative operator and solve a fractional differential equation, previous states of the system must be taken into account. Dynamical systems with fractional derivatives present challenges for nonlinear analysis and involve an innovative approach to chaos detection. In the plasma environment, models having fractional derivatives in time, space, or space-time have attracted remarkable attention from many researchers. The KdV system having fractional time derivative was derived from studying IAWs in a collisionless plasma, and it was observed that fractional time operators considerably influence the amplitudes of the IAS. ⁶³ Additionally, the analytical solution of the modified Korteweg–de Vries (mKdV) system having a space-time fractional derivative was recently evaluated by Jamilu et al. ⁶⁴ using a sine-cosine approach. The Kadomstev–Petviashvili (KP) equation having space-time fractional derivative was derived to investigate the impact of the fractional parameter on the IAWs. ⁶⁵ Nazari–Golshan derived space fractional Korteweg–de Vries–Burgers (KdVB) equation to study shock wave excitation in nonextensive plasma with relativistic ions. ⁶⁶ Uddin et al. ⁶⁷ also used the KdVB equation and investigated the effects of β fractional derivative and other plasma parameters on nonlinear ion-acoustic shock waves.

Since genuine physical systems frequently include a variety of perturbations, nonlinear dynamical systems under external perturbations have recently received much interest recently.^68–70 Research indicates that chaotic dynamics, a fascinating nonlinear phenomenon, might result by adding a perturbation to an integrable solution.^71,72 Chaotic systems refer to many physical systems which, although governed by perfectly deterministic laws of evolution, become unpredictable over a certain period of time. ⁷³ Chaotic systems are susceptible to initial conditions, which means that considering a system with two very close initial conditions leads to a considerable trajectory divergence over time. From the point of view of the dynamics of nonlinear systems and the theory of solitary waves, chaotic movements are of particular interest.^74,75 Bifurcation and chaotic oscillations in a forced dynamical system are observed in complex plasma through the quasi-periodic route to chaos.^76,77 Saha et al. ⁷⁸ carried out a qualitative study on positron acoustic waves in a multicomponent plasma system and showed that cold positron dynamics involve the excitation of solitons of opposite polarities (compressive and rarefactive). Furthermore, considering external periodic perturbation on dynamic systems allowed them to exhibit quasi-periodic and chaotic orbits. Based on the theory of planar dynamical systems, Tamang et al. ⁷⁹ studied the nonlinear characteristic of positron acoustic waves and showed the appearance of quasi-periodic and chaotic states in this model subjected to an external periodic disturbance. Saha and coworkers ⁸⁰ also investigated multiperiodic wave and chaotic features of IAWs under the KP equation in a plasma model with nonextensive distribution. Recently, Pradhan and his collaborators ⁸¹ analyzed the chaotic scenario of IAWs in a quantum pair ion plasma. Moreover, ion-acoustic periodic, super-periodic, kink, and anti-kink solitary wave solutions were examined by Punam et al. ⁸² in nonextensive solar wind plasma. The quasi-period route to chaos was also demonstrated in their analysis due to the external periodic excitation. Furthermore, Bhowmick et al. ⁸³ showed the presence of chaotic structures in a dusty plasma with nonthermal species and external MF.

A survey of the available literature shows that the transition from quasi-periodic state to chaotic motion of IAWs in magnetized electron-ion plasmas with fractional-order derivative and nonextensive hot and cold electrons has not been investigated so far. Therefore, our objective in this work is to investigate, on the one hand, the bifurcation analysis of IAWs by using a planar Hamiltonian system. On the other hand, the quasi-periodicity and chaotic profiles of the IAWs under the influence of fractional-order derivatives in magnetized plasmas with TET on the framework of mKdV equations are examined.

Presentation of the plasma model

To investigate the dynamics of IAWs, we consider a collisionless magnetized plasma model composed of cold and hot electrons modeled by a nonextensive distribution and positively charged inertial ions. In the absence of any perturbation, n_h0/n_i0 = 1 − n_c0/n_i0 represents the overall charge neutrality condition, which is expressed as a ratio of electron density at equilibrium to that of the ions. Here, n_c0(n_h0) and n_i0 represent the number densities of cold (hot) electrons and ions, respectively. The ambient MF ${\widehat{B}}_0$ is assumed to act along the z-axis, that is, ${\widehat{B}}_0$ = $B_0\widehat{z}$ . In this fluid model, the behavior of inertial ions is governed by the equations of continuity and motion and coupled to the Poisson equation. In normalized form, these three equations are written:

${\partial }_t{n}_i+\nabla \left(n_i{U}_i\right)=0,$

(1)

${\partial }_t{U}_i+(U_i\nabla )U_i=-\nabla \varphi -U_i\times B_0\widehat{z},$

(2)

${\nabla }^2\varphi =\beta n_c+(1-\beta )n_h-n_i,$

(3)where β = n_c0/n_i0.

The number densities n _c (cold electrons) and n _h (hot electrons) can be obtained by using the q-nonextensive distributions ²² :

$n_j=n_{j0}{\left[1+e\frac{\varphi (q_j-1)}{k_B{T}_j}\right]}^{\frac{q_j+1}{2(q_j-1)}},j=c,h,$

(4)where cold (hot) electrons are subjected to q _c (q _h ), a nonextensive parameter that indicates the deviation from thermodynamic equilibrium (q → 1). Moreover, as noted in the literature, ⁸⁴ using the nonextensive distribution with q < 1 is strongly advised to explain the characteristics of solitary wave structures in a plasma model. Therefore, we assume q_c,h < 1 in the rest of our study.Using the power series expansion of the exponential functions around zero, equation (3) becomes

$\frac{{\partial }^2\varphi }{\partial x^2}+\frac{{\partial }^2\varphi }{\partial y^2}+\frac{{\partial }^2\varphi }{\partial z^2}={\mu }_2{\varphi }^2+{\mu }_1\varphi -n_i+1,$

(5)where μ₁ and μ₂ depend on the two-temperature electron and the nonextensive parameter

$\begin{array}{c}{\mu }_1=\frac{\beta (q_c+1)+(1-\beta )(q_h+1)\sigma }{2},\hfill \\ {\mu }_2=\frac{\beta (q_c+1)(3-q_c)}{8}+\frac{(1-\beta )(q_h+1)(3-q_h){\sigma }^2}{8},\hfill \end{array}$

(6)and σ = T _c /T _h . For the entire computations, we will use the following defined parameter values, which are consistent with those used in astrophysics and space plasma: β = 0.1 − 1, σ = 0.1 − 1, and q = 0.3 − 1. ⁸⁵ The parameter β, defined as the ratio of cold electron density to ion density, plays a crucial role in plasma dynamics, where β ≃ 0.1 leads to increased turbulence due to insufficient cold electrons, while β ≃ 0.5 creates a stable regime characteristic of laboratory plasmas through balanced electron-ion densities, and β ≃ 1 supports IAW propagation while maintaining quasi-neutrality. Similarly, the temperature ratio σ between cold and hot electrons significantly influences plasma behavior, with σ ≃ 0.1 causing strong population separation and electrostatic instabilities, σ ≃ 0.5 improving stability through moderated thermal differences, and σ ≃ 1 minimizing thermal separation effects. Furthermore, the nonextensivity parameter q in Tsallis statistics governs plasma characteristics, where q = 0.3 indicates localized particle interactions near equilibrium, q = 0.5 introduces moderate turbulence through nonlocal effects, and q ≃ 1 reverts to Boltzmann–Gibbs statistics with standard exponential distributions.

Reduction to the mKdV equation

By applying the RPM, we study the weakly nonlinear propagation of IAWs in the magnetized plasma model considered. This method extends the independent variables in the following ways:

$\xi =ϵ\left(k_{1}x+k_{2}y+k_{3}z-v_{p}t\right),\tau ={ϵ}^{3}t,$

(7)where v _p is the IAWs’ phase speed and ϵ is a small parameter, whereas k₁, k₂, and k₃ are the direction cosines along x, y, and z axes, respectively. Then $k_1^2+k_2^2+k_3^2=1$ . The dynamic variables n(x, t), u(x, t), v(x, t), w(x, t), and ϕ(x, t) are put in the form:

$\begin{array}{c}n=1+{ϵ}^{1/2}{n}_1+ϵn_2+{ϵ}^{3/2}{n}_3,\hfill \\ u=ϵu_1+{ϵ}^{3/2}{u}_2+{ϵ}^2{u}_3,\hfill \\ v=ϵv_1+{ϵ}^{3/2}{v}_2+{ϵ}^2{v}_3,\hfill \\ w={ϵ}^{1/2}{w}_1+ϵw_2+{ϵ}^{3/2}{w}_3,\hfill \\ \varphi ={ϵ}^{1/2}{\varphi }_1+ϵ{\varphi }_2+{ϵ}^{3/2}{\varphi }_3.\hfill \end{array}$

(8)Substituting the previous perturbed quantities defined by equation (8) into equations (1)–(3), and applying the RPM, we obtain in the lowest order of ϵ:

$n_1=\frac{k_3^2{\varphi }_1}{v_p^2},u_1=\frac{-k_2}{{\omega }_{ci}}\frac{\partial {\varphi }_1}{\partial \xi },v_1=\frac{k_1}{{\omega }_{ci}}\frac{\partial {\varphi }_1}{\partial \xi },w_1=\frac{k_3{\varphi }_1}{v_p},$

(9)where the phase speed of the IAWs is given by

$v_p^2=\frac{k_3^2}{{\mu }_1}.$

(10)At higher-order ϵ^3/2, we obtain

$\begin{array}{c}{n}_2={\mu }_2{\varphi }_1^2-\frac{k_3^2{\varphi }_1^2{\mu }_1}{2v_p^2},\hfill \\ u_2=\frac{k_2{k}_3^2{\varphi }_1\left(\frac{\partial {\varphi }_1}{\partial \xi }\right){\omega }_{ci}+v_p^3{k}_1\frac{{\partial }^2{\varphi }_1}{\partial {\xi }^2}}{{\omega }_{ci}^2{v}_p^2},\hfill \\ v_2=\frac{-k_1{k}_3^2{\varphi }_1\left(\frac{\partial {\varphi }_1}{\partial \xi }\right){\omega }_{ci}+v_p^3{k}_2\frac{{\partial }^2{\varphi }_1}{\partial {\xi }^2}}{{\omega }_{ci}^2{v}_p^2},\hfill \\ w_2=0,\hfill \\ {\varphi }_2=-\frac{k_3^2{\varphi }_1^2}{2v_p^2}.\hfill \end{array}$

(11)

Considering the order of power ϵ², we obtain

$\begin{array}{ccc}\hfill k_2\frac{\partial }{\partial \xi }(n_1{v}_1)+k_3\frac{\partial }{\partial \xi }(n_1{w}_2)+k_1\frac{\partial }{\partial \xi }(n_1{u}_1)+k_3\frac{\partial }{\partial \xi }(n_2{w}_1)+k_1{n}_1\frac{\partial u_1}{\partial \xi }+k_2{n}_1\frac{\partial v_1}{\partial \xi }+& \hfill & \hfill \\ \hfill k_3{n}_2\frac{\partial w_1}{\partial \xi }-v_g\frac{\partial n_3}{\partial \xi }+k_1\frac{\partial u_2}{\partial \xi }+k_2\frac{\partial v_2}{\partial \xi }+k_3\frac{\partial w_3}{\partial \xi }+\frac{\partial n_1}{\partial \tau }=0.& \hfill & \hfill \end{array}$

(12)

Using equations (9)–(11) and setting ϕ = ϕ₁, equation (12) becomes

$\frac{\partial \varphi }{\partial \tau }+S_1{\varphi }^2\frac{\partial \varphi }{\partial \xi }+S_2\frac{{\partial }^3\varphi }{\partial {\xi }^3}=0,$

(13)Equation (13) represents the mKdV equation, which combines the nonlinear effects (via S₁) and the dispersive effects (via S₂) given by the following expressions:

$S_1=\frac{2{\mu }_2{v}_p^4+3k_3^4}{2v_p^3},S_2=\frac{v_p^3\left(k_1^2+k_2^2\right)}{k_3^2{\omega }_{ci}^2}.$

(14)These expressions show that S₁ does not depend on the MF parameter ω _ci . In contrast, both coefficients S₁ and S₂ are functions of the nonextensive velocity distribution and the TET.

Dynamical system and chaotic pattern with fractional-order derivative

To study the behavior of IAWs on the framework of equation (13), we introduce into a stationary frame the transformation χ = ξ − v₀τ; here, v₀ is the traveling waves’ velocity, and then equation (13) becomes

$-v_0\frac{\mathrm{d}\varphi }{\mathrm{d}\mathrm{\chi }}+S_1{\varphi }^2\frac{\mathrm{d}\varphi }{\mathrm{d}\mathrm{\chi }}+S_2\frac{{\mathrm{d}\varphi }^3}{\mathrm{d}{\mathrm{\chi }}^{\mathrm{3}}}=0.$

(15)Integrating equation (15), and using the appropriate boundary conditions (ϕ → 0 and dϕ/dχ → 0 at χ → ±∞), we end up with the equation

$\frac{{\mathrm{d}}^2\varphi }{\mathrm{d}{\chi }^2}=\frac{v_0\varphi -\frac{S_1{\varphi }^3}{3}}{S_2}.$

(16)The nonlinear differential equation (16) can be put in the form of a second-order autonomous nonlinear system as follows:

$\begin{array}{c}\frac{\mathrm{d}\varphi }{\mathrm{d}\chi }=y,\hfill \\ \frac{\mathrm{d}\mathrm{z}}{\mathrm{d}\chi }=\frac{v_0\varphi -\frac{S_1{\varphi }^3}{3}}{S_2}.\hfill \end{array}$

(17)

We assume that $\overrightarrow{\mathrm{\Gamma }}=(y,v_0\varphi -\frac{S_1{\varphi }^3}{3}/S_2)$ , which corresponds to system (17). Then $\overrightarrow{\mathrm{\Gamma }}$ satisfies the condition $\overrightarrow{\nabla }\overrightarrow{\mathrm{\Gamma }}=0$ ; thus, system (17) is conservative. Moreover, this system is a planar Hamiltonian system with Hamiltonian functions:

$H\left(\varphi ,y\right)=\frac{y^2}{2}+\frac{S_1{\varphi }^4}{12S_2}-\frac{v_0{\varphi }^2}{2S_2}.$

(18)The potential energy of the system, which depends on the relevant physical plasma parameters, reads

$E(\varphi )=\frac{S_1{\varphi }^4}{12S_2}-\frac{v_0{\varphi }^2}{2S_2}.$

(19)

However, when the nonlinear dynamic system interacts with external shocks, the resulting system is perturbed and can lead to irregular or chaotic behavior. To investigate the nonlinear dynamics of the perturbed system, we introduce a periodic disturbance as follows:

$\begin{array}{c}\frac{\mathrm{d}}{\mathrm{d}\chi }\varphi =y,\hfill \\ \frac{\mathrm{d}}{\mathrm{d}\chi }y=\frac{v_0\varphi -\frac{S_1{\varphi }^3}{3}}{S_2}+f_0\cos (\omega \chi ),\hfill \end{array}$

(20)where f₀ and ω represent the magnitude and frequency of the periodic disturbance, respectively, and this forced dynamical system can reveal interesting periodic and chaotic structures. These different structures will be analyzed using time series and phase portrait plots. To investigate the effect of fractional-order derivatives on the framework of system (20), we consider the Grunwald–Letnikov fractional derivative of order α of the function g(t) given by ⁸⁶

$a{D_t}^{\alpha }g(t)=\underset{h\to 0}{\lim }\frac{1}{h^{\alpha }}\sum _{j=0}^{\frac{t-a}{h}}{(-1)}^j\frac{\alpha !}{j!(\alpha -j)!}g(t-jh).$

(21)Here, α is a non-integer positive real number (0 < α < 1). According to this definition, the dynamical forced mKdV system (20) becomes

$\begin{array}{c}a{D_t}^{\alpha }x\left(t\right)=y\left(t\right),\hfill \\ a{D_t}^{\alpha }y\left(t\right)={\alpha }_{1}x\left(t\right)-{\alpha }_2{x}^3\left(t\right)+f_{0}cos\left(\omega t\right),\hfill \end{array}$

(22)where α₁ = v₀/S₂ and α₂ = S₁/3S₂. By assuming that t _k = kh, (k = 1, 2, 3, …), and a = (t _k − L _m /h), where L _m is the length of the memory and h is the time step, the explicit numerical approximation of fractional derivation at time t _k point is

$\begin{array}{c}a{D_{t_k}}^{\alpha }g(t)=h^{-\alpha }\sum _{j=0}^k{(-1)}^j\frac{\alpha !}{j!(\alpha -j)!}g((k-j)h),\hfill \\ a{D_{t_k}}^{\alpha }g(t)=h^{-\alpha }\sum _{j=0}^k{(-1)}^j\frac{\alpha !}{j!(\alpha -j)!}g(t_{k-j}),\hfill \end{array}$

(23)which leads to the formula

$a{D_{t_k}}^{\alpha }g(t)=h^{-\alpha }\sum _{j=0}^k{c}_j^{\alpha }g(t_{k-j}).$

(24)Now, we apply this formula to equation (22), that is, $a{D_t}^{\alpha }x(t)=g(x(t),t)$ , so the numerical solution is obtained as

$\begin{array}{ccc}\hfill h^{-\alpha }\sum _{j=0}^k{c}_j^{\alpha }x(t_{k-j})=g(x(t_k),t_k),& \hfill & \hfill \\ \hfill c_0^{\alpha }x(t_k)+\sum _{j=1}^k{c}_j^{\alpha }x(t_{k-j})=h^{\alpha }g(x(t_k),t_k),& \hfill & \hfill \\ \hfill x(t_k)=h^{\alpha }g(x(t_k),t_k)-\sum _{j=1}^k{c}_j^{\alpha }x(t_{k-j}).& \hfill & \hfill \end{array}$

(25)The fractional-order Newton–Leipnik algorithm below can be used to obtain the discrete form of the set of equations ⁸⁷ :

$\begin{array}{c}x\left(t_k\right)=y\left(t_{k-1}\right)h^{\alpha }-\sum _{j=1}^k{c_j}^{\alpha }x\left(t_{k-j}\right)\hfill \\ y\left(t_k\right)=\left({\alpha }_{1}x\left(t_k\right)-{\alpha }_2{x}^3\left(t_k\right)+f_{0}cos\left(\omega t_k\right)\right)h^{\alpha }-\sum _{j=1}^k{c_j}^{\alpha }y\left(t_{k-j}\right),\hfill \end{array}$

(26)where the coefficients $c_j^{\alpha }$ meet the following recursive relations: $c_j^{\alpha }=\left(1-\alpha +1/j\right)c_{j-1}^{\alpha }$ , and $c_0^{\alpha }=1$ .

Therefore, the Newton–Leipnik algorithm (26) is determined numerically by applying the fourth-order Runge–Kutta algorithm. The step-size h is chosen to balance accuracy and computational efficiency. For fractional systems, smaller h is often required due to the memory effect and typical values range from h = 0.01 to h = 0.001. Convergence criteria are typically ensured by step-size refinement. For chaotic systems, convergence is inferred from the consistency of Lyapunov exponents or bifurcation diagrams under step-size changes. According to the Runge–Kutta method that we used in this work, the integration step is h = 0.01, the simulation time is T _sim = 2000, and the local truncation error per step is O(h⁵). We will now analyze the impacts of fractional-order derivative terms and external periodic perturbation on the system. First, we consider that the system is unforced, that is, f₀ = 0.In this case, Figure 1 represents the evolution of the system dynamics for different values of the fractional-order derivative term. For α = 0.98, the phase portrait in panel 1(a) shows an aperiodic limit cycle confirmed by their corresponding time series (panel 1(b)). On the other hand, for α = 1, the behavior of the system changes, and the phase portrait in panel 1(c) displays a supernonlinear periodic orbit, whereas in panel 1(d) one can observe a supernonlinear periodic wave solution. So, the fractional effect alters the behavior of the system. The slight deviation from α = 1 introduces weak memory effects, leading to irregular but bounded oscillations.

Figure 1.

Limit cycles (a), aperiodic orbit (b), super-periodic orbit (c), and supernonlinear periodic wave (d), for ω _ci = 0.5, β = 0.2, q _c = 0.5, q _h = 0.3, σ = 0.2, v₀ = 0.2, f₀ = 0, and ω = 0.

Second, by fixing the magnitude of the periodic perturbation at f₀ = 0.8 and its frequency at ω = 0.05, we have obtained new phase portraits, which are illustrated in Figure 2. By setting α = 0.7, the phase portrait (panel 2(a)) and time series (panel 2(b)) display periodic oscillations. However, for α = 0.75, panel 2(c) presents a quasi-periodic state with their corresponding time series shown in panel 2(d). After that, for α = 0.85, the phase space (panel 2(e)) and time series (panel 2(f)) present the characteristics of chaotic oscillations due to the influence of the fractional-order term. Therefore, when the fractional parameter α increases from 0.7 to 0.85, the system passes successively from periodic oscillations to a quasi-periodic state and then to a chaotic state. On the other hand, the unpredictability of a chaotic system also manifests itself under the influence of the fractional effect. Indeed, we observed that when α increases from 0.85 to 0.9, the initially chaotic system becomes periodic again. This behavior is observed for α = 0.9 through panels 2(g) and 2(h), which exhibit periodic oscillations in the phase portrait and time series, respectively. This suggests that for strong memory effects (low values of α), the system is stable, while weaker memory (high values of α) destabilizes it. As α → 1, the system approaches classical (integer-order) dynamics, where dispersion balances nonlinearity, restoring stable periodic waves.

Figure 2.

Transition between periodic, quasi-periodic, and chaotic state with ω _ci = 0.8, β = 0.2, q _c = 0.5, σ = 0.8, q _h = 0.3, v₀ = 0.2, f₀ = 0.8, and ω = 0.05.

To complete our investigation and confirm the results of the above analysis, we now use the bifurcation diagram, another tool to reveal chaos in a nonlinear dynamic system. The bifurcation diagram for the different values of α used previously is shown in Figure 3, where we took the magnitude of the periodic perturbation f₀ as a bifurcation parameter. We can observe the phenomena of period-doubling, which disappears when the fractional parameter value α is 1. This suggests that a periodic state gives way to a quasi-periodic or chaotic structure in the system. The chaotic interval begins for the parameter values f₀ = 0.9 (panel 3(a) with α = 0.7), f₀ = 0.5 (panel 3(b) with α = 0.85), f₀ = 0.25 (panel 3(c) with α = 0.9), and f₀ = 2.1 (panel 3(d) with α = 0.98). The results of Figure 3 also indicate that strong memory effects (low α) destabilize the system, while weaker memory (high α) stabilizes it. It is evident that when α tends toward 1, the interval where the system exhibits a chaotic condition reduces. It should be noted that Li and Yorke, ⁸⁸ through an article entitled “Period three implies Chaos,” were the first in 1975 to study chaotic motion in a mathematical model. Therefore, this reinforces our hypothesis that there will be chaos in the area where the period-doubling intersects.

Figure 3.

Bifurcation diagram of IAWs against f₀ when parameter α varies. The other parameters are identical to those used in Figure 2.

A comparable bifurcation diagram can be used to conduct a more methodical examination of the system’s behavior. Taking the fractional derivative α as a bifurcation parameter, the bifurcation diagram is given in Figure 4. The areas of α for which the system response shifts from chaotic (0 < α < 0.35) to periodic/quasi-periodic via a period-doubling cascade (0.4 < α < 0.8) and back to chaotic (0.8 < α < 0.9), then to periodic (0.9 < α < 1) are visible in Figure 4. Since the period-doubling phenomena is seen in the 0.4 < α < 0.8 range, it is quite likely that both periodic and quasi-periodic states can exist for the fractional parameter values selected in this zone. The fractional parameter values we employed in Figure 2(a)–(d), however are consistent with this range. On the other hand, the period-doubling intersects in 0.8 < α < 0.9, which is the time frame during which we were able to obtain a chaotic phase portrait Figure 2(e). The system is periodic for α = 0.9 and f₀ = 0.8, as demonstrated by the result shown in Figure 2(g)–(h), which is finally confirmed by the final interval 0.9 < α < 1. Remember that we obtain a periodic solution matching the solution for the classic forced mKdV equation when α approaches near 1. Thus, this analysis shows us the unpredictability of a chaotic system and its very high sensitivity to variations in the values of physical parameters. It is seen from these investigations that the fractional-order derivative term strongly influences the transition to chaos in our plasma model in a forced system. Moreover, in an unforced system, aperiodic, super-periodic, and periodic limit cycles are observed under the influence of the fractional-order derivative term. Furthermore, similar results were obtained in Ref. 89 where authors used the Duffing system with a fractional derivative term and different numerical tools for studying the dynamical system. The quasi-periodic route to chaos was found to depend on the fractional derivative term.

Figure 4.

Bifurcation diagram of IAWs against α, with the same values of Figure 2.

Conclusion

This research investigated the bifurcation and transition to chaos of the ion-acoustic waves in magnetoplasmas. The fluid model has inertial positive ions and cold and hot electrons, which follow the q-distribution function. The mKdV equation is derived, and their corresponding dynamical system is obtained. On the other hand, the perturbed dynamical system is used to study both regular and irregular states according to the memory effect. The Grunwald–Letnikov definition and the Newton–Leipnik algorithm are used for the dynamical forced system of the mKdV equation. The phase portrait, time series, and bifurcation diagram are plotted through the fourth-order Runge–Kutta scheme, which presents periodic, quasi-periodic, and chaotic oscillations as the fractional-order derivative term varies. These investigations show that, in a forced system, the fractional-order derivative term determines the transition to chaos in this plasma model. Furthermore, in an unforced system, the fractional-order derivative term impacts the occurrence of aperiodic, super-periodic, and periodic limit cycles. To extend MHD applicability to collisionless plasmas while accounting for kinetic effects, the Vlasov–Maxwell system can be enhanced through fractional calculus, providing a more sophisticated framework for studying high-dimensional plasma interactions. Incorporating collisional effects via Vlasov–Boltzmann or Vlasov–Fokker–Planck models would further improve accuracy and enable deeper understanding of complex plasma dynamics. Our future work will explore these systems through generalized fractional approaches (GFAs) to assess their viability for plasma wave studies, complemented by advanced numerical techniques including Galerkin methods and machine learning algorithms. This combined methodology promises to enhance solution accuracy and facilitate comprehensive analysis of IAW dynamics in drift wave turbulence environments, with results to be presented in an upcoming publication.