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.

Keywords

chaos fractional effect low-frequency ion-acoustic nonlinear structures magnetoplasma Newton–Leipnik algorithm

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 ${\hat{B}}_{0}$ is assumed to act along the z-axis, that is, ${\hat{B}}_{0}$ = $B_{0} \hat{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 (n_{i} U_{i}) = 0,

(1)

\partial_{t} U_{i} + (U_{i} \nabla) U_{i} = - \nabla ϕ - U_{i} \times B_{0} \hat{z},

(2)

\nabla^{2} ϕ = β n_{c} + (1 - β) 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_{j 0} {[1 + e \frac{ϕ (q_{j} - 1)}{k_{B} T_{j}}]}^{\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} ϕ}{\partial x^{2}} + \frac{\partial^{2} ϕ}{\partial y^{2}} + \frac{\partial^{2} ϕ}{\partial z^{2}} = μ_{2} ϕ^{2} + μ_{1} ϕ - n_{i} + 1,

(5)

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

\begin{matrix} μ_{1} = \frac{β (q_{c} + 1) + (1 - β) (q_{h} + 1) σ}{2}, \\ μ_{2} = \frac{β (q_{c} + 1) (3 - q_{c})}{8} + \frac{(1 - β) (q_{h} + 1) (3 - q_{h}) σ^{2}}{8}, \end{matrix}

(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:

ξ = ϵ (k_{1} x + k_{2} y + k_{3} z - v_{p} t), τ = ϵ^{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{matrix} n = 1 + ϵ^{1 / 2} n_{1} + ϵ n_{2} + ϵ^{3 / 2} n_{3}, \\ u = ϵ u_{1} + ϵ^{3 / 2} u_{2} + ϵ^{2} u_{3}, \\ v = ϵ v_{1} + ϵ^{3 / 2} v_{2} + ϵ^{2} v_{3}, \\ w = ϵ^{1 / 2} w_{1} + ϵ w_{2} + ϵ^{3 / 2} w_{3}, \\ ϕ = ϵ^{1 / 2} ϕ_{1} + ϵ ϕ_{2} + ϵ^{3 / 2} ϕ_{3} . \end{matrix}

(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} ϕ_{1}}{v_{p}^{2}}, u_{1} = \frac{- k_{2}}{ω_{c i}} \frac{\partial ϕ_{1}}{\partial ξ}, v_{1} = \frac{k_{1}}{ω_{c i}} \frac{\partial ϕ_{1}}{\partial ξ}, w_{1} = \frac{k_{3} ϕ_{1}}{v_{p}},

(9)

where the phase speed of the IAWs is given by

v_{p}^{2} = \frac{k_{3}^{2}}{μ_{1}} .

(10)

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

\begin{matrix} n_{2} = μ_{2} ϕ_{1}^{2} - \frac{k_{3}^{2} ϕ_{1}^{2} μ_{1}}{2 v_{p}^{2}}, \\ u_{2} = \frac{k_{2} k_{3}^{2} ϕ_{1} (\frac{\partial ϕ_{1}}{\partial ξ}) ω_{c i} + v_{p}^{3} k_{1} \frac{\partial^{2} ϕ_{1}}{\partial ξ^{2}}}{ω_{c i}^{2} v_{p}^{2}}, \\ v_{2} = \frac{- k_{1} k_{3}^{2} ϕ_{1} (\frac{\partial ϕ_{1}}{\partial ξ}) ω_{c i} + v_{p}^{3} k_{2} \frac{\partial^{2} ϕ_{1}}{\partial ξ^{2}}}{ω_{c i}^{2} v_{p}^{2}}, \\ w_{2} = 0, \\ ϕ_{2} = - \frac{k_{3}^{2} ϕ_{1}^{2}}{2 v_{p}^{2}} . \end{matrix}

(11)

Considering the order of power ϵ², we obtain

\begin{matrix} k_{2} \frac{\partial}{\partial ξ} (n_{1} v_{1}) + k_{3} \frac{\partial}{\partial ξ} (n_{1} w_{2}) + k_{1} \frac{\partial}{\partial ξ} (n_{1} u_{1}) + k_{3} \frac{\partial}{\partial ξ} (n_{2} w_{1}) + k_{1} n_{1} \frac{\partial u_{1}}{\partial ξ} + k_{2} n_{1} \frac{\partial v_{1}}{\partial ξ} + \\ k_{3} n_{2} \frac{\partial w_{1}}{\partial ξ} - v_{g} \frac{\partial n_{3}}{\partial ξ} + k_{1} \frac{\partial u_{2}}{\partial ξ} + k_{2} \frac{\partial v_{2}}{\partial ξ} + k_{3} \frac{\partial w_{3}}{\partial ξ} + \frac{\partial n_{1}}{\partial τ} = 0 . \end{matrix}

(12)

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

\frac{\partial ϕ}{\partial τ} + S_{1} ϕ^{2} \frac{\partial ϕ}{\partial ξ} + S_{2} \frac{\partial^{3} ϕ}{\partial ξ^{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 μ_{2} v_{p}^{4} + 3 k_{3}^{4}}{2 v_{p}^{3}}, S_{2} = \frac{v_{p}^{3} (k_{1}^{2} + k_{2}^{2})}{k_{3}^{2} ω_{c i}^{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{d ϕ}{d χ} + S_{1} ϕ^{2} \frac{d ϕ}{d χ} + S_{2} \frac{{d ϕ}^{3}}{d χ^{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{d^{2} ϕ}{d χ^{2}} = \frac{v_{0} ϕ - \frac{S_{1} ϕ^{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{matrix} \frac{d ϕ}{d χ} = y, \\ \frac{d z}{d χ} = \frac{v_{0} ϕ - \frac{S_{1} ϕ^{3}}{3}}{S_{2}} . \end{matrix}

(17)

We assume that $\vec{Γ} = (y, v_{0} ϕ - \frac{S_{1} ϕ^{3}}{3} / S_{2})$ , which corresponds to system (17). Then $\vec{Γ}$ satisfies the condition $\vec{\nabla} \vec{Γ} = 0$ ; thus, system (17) is conservative. Moreover, this system is a planar Hamiltonian system with Hamiltonian functions:

H (ϕ, y) = \frac{y^{2}}{2} + \frac{S_{1} ϕ^{4}}{12 S_{2}} - \frac{v_{0} ϕ^{2}}{2 S_{2}} .

(18)

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

E (ϕ) = \frac{S_{1} ϕ^{4}}{12 S_{2}} - \frac{v_{0} ϕ^{2}}{2 S_{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{matrix} \frac{d}{d χ} ϕ = y, \\ \frac{d}{d χ} y = \frac{v_{0} ϕ - \frac{S_{1} ϕ^{3}}{3}}{S_{2}} + f_{0} \cos (ω χ), \end{matrix}

(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}}^{α} g (t) = \lim_{h \to 0} \frac{1}{h^{α}} \sum_{j = 0}^{\frac{t - a}{h}} {(- 1)}^{j} \frac{α!}{j! (α - j)!} g (t - j h) .

(21)

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

\begin{matrix} a {D_{t}}^{α} x (t) = y (t), \\ a {D_{t}}^{α} y (t) = α_{1} x (t) - α_{2} x^{3} (t) + f_{0} c o s (ω t), \end{matrix}

(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{matrix} a {D_{t_{k}}}^{α} g (t) = h^{- α} \sum_{j = 0}^{k} {(- 1)}^{j} \frac{α!}{j! (α - j)!} g ((k - j) h), \\ a {D_{t_{k}}}^{α} g (t) = h^{- α} \sum_{j = 0}^{k} {(- 1)}^{j} \frac{α!}{j! (α - j)!} g (t_{k - j}), \end{matrix}

(23)

which leads to the formula

a {D_{t_{k}}}^{α} g (t) = h^{- α} \sum_{j = 0}^{k} c_{j}^{α} g (t_{k - j}) .

(24)

Now, we apply this formula to equation (22), that is,

a {D_{t}}^{α} x (t) = g (x (t), t)

, so the numerical solution is obtained as

\begin{matrix} h^{- α} \sum_{j = 0}^{k} c_{j}^{α} x (t_{k - j}) = g (x (t_{k}), t_{k}), \\ c_{0}^{α} x (t_{k}) + \sum_{j = 1}^{k} c_{j}^{α} x (t_{k - j}) = h^{α} g (x (t_{k}), t_{k}), \\ x (t_{k}) = h^{α} g (x (t_{k}), t_{k}) - \sum_{j = 1}^{k} c_{j}^{α} x (t_{k - j}) . \end{matrix}

(25)

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

\begin{matrix} x (t_{k}) = y (t_{k - 1}) h^{α} - \sum_{j = 1}^{k} {c_{j}}^{α} x (t_{k - j}) \\ y (t_{k}) = (α_{1} x (t_{k}) - α_{2} x^{3} (t_{k}) + f_{0} c o s (ω t_{k})) h^{α} - \sum_{j = 1}^{k} {c_{j}}^{α} y (t_{k - j}), \end{matrix}

(26)

where the coefficients

c_{j}^{α}

meet the following recursive relations:

c_{j}^{α} = (1 - α + 1 / j) c_{j - 1}^{α}

, and

c_{0}^{α} = 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.

Footnotes

Acknowledgments

The authors express their gratitude to Princess Nourah bint Abdulrahman University Researchers Supporting Project Number (PNURSP2025R2), Princess Nourah bint Abdulrahman University, Riyadh, Saudi Arabia. This study is supported via funding from Prince Sattam bin Abdulaziz University project number (PSAU/2025/R/1446).

ORCID iDs

Samir A. El-Tantawy

Alim

B. B. Mouhammadoul

Camus G. L. Tiofack

Funding

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Use of Artificial Intelligence (AI) tools declaration

The authors state that no AI tools were employed in the preparation of this article.

Data Availability Statement

This article does not involve data sharing as no new data were produced in the study. If needed, the Mathematica codes for generating the figures can be acquired from the corresponding author.*

References

Truemann

Baumjohann

. Advanced space plasma physics. Imperial College, 1984, p. 243.

Nicholson

. Introduction to plasma theory. Wiley, 1983.

Rufai

Bharuthram

Singh

, et al. Low frequency solitons and double layers in a magnetized plasma with two temperature electrons. Phys Plasmas 2012; 19: 122308.

Baluku

Hellberg

. Ion acoustic solitons in a plasma with two-temperature kappa-distributed electrons. Phys Plasmas 2012; 19: 012106.

Baluku

Hellberg

Verheest

. New light on ion acoustic solitary waves in a plasma with two-temperature electrons. Europhys Lett 2010; 91: 15001.

Annou

. Ion-acoustic solitons in plasma: an application to Saturn’s magnetosphere. Astrophys Space Sci 2015; 357: 163.

Shahmansouri

Alinejad

. Electrostatic wave structures in a magnetized superthermal plasma with two-temperature electrons. Phys Plasmas 2013; 20: 082130.

Panwar

Ryu

Bains

. Oblique ion-acoustic cnoidal waves in two temperature superthermal electrons magnetized plasma. Phys Plasmas 2014; 21: 122105.

Saeed

Sarfraz

Yoon

, et al. Characteristics of heat flux and electromagnetic electron-cyclotron instabilities driven by solar wind electrons. Mon Not Roy Astron Soc 2017; 466: 4928.

10.

Nakamura

Sugai

. Experiments on ion-acoustic solitons in a plasma. Chaos Solitons Fractals 1996; 7: 1023–1031.

11.

Abbasi

Masood

Khan

, et al. Nonlinear ion acoustic waves in dissipative and dispersive magneto-rotating relativistic plasmas with two temperature superthermal electrons. Contrib Plasma Phys 2020; 60: e202000050.

12.

Nawaz

Jahangir

Masood

, et al. Cubic nonlinearity driven dust ion acoustic solitons with superthermal two-temperature electrons in the Saturn’s magnetosphere. Phys Scripta 2022; 97: 075601.

13.

Goswami

Buti

. Ion acoustic solitary waves in a two-electron-temperature plasma. Phys Lett 1976; 57: 149–150.

14.

Lakhina

Kakad

Singh

, et al. Ion-and electron-acoustic solitons in two-electron temperature space plasmas. Phys Plasmas 2008; 15: 062903.

15.

Nawaz

Masood

Jahangir

, et al. Interaction of gardner solitons in plasmas: applications in the Saturn’s magnetosphere. Phys Scripta 2021; 96: 045604.

16.

Noman

Islam

Hassan

, et al. Dust-ion-acoustic rogue waves in a dusty plasma having super-thermal electrons. Gas 2021; 1: 106–116.

17.

Sharma

Deka

Paul

, et al. Experimental study on controlled production of two-electron temperature plasma. Plasma Sources Sci Technol 2022; 31: 025013.

18.

Rao

Shukla

. Coupled langmuir and ion-acoustic waves in two-electron temperature plasmas. Phys Plasmas 1997; 4: 636–645.

19.

Liu

. Dust acoustic instability driven by drifting ions and electrons in the dust plasma with Lorentzian kappa distribution. Phys Plasmas 2009; 16: 123707.

20.

Boldyrev

Horaites

. Kinetic theory of the electron strahl in the solar wind. Mon Not Roy Astron Soc 2019; 489: 3412–3419.

21.

Renyi

. On a new axiomatic theory of probability. Acta Math Acad Sci Hungar 1055; 6: 285–335.

22.

Tsallis

. Dynamical scenario for nonextensive statistical mechanics. Physica A 2004; 340: 1–10.

23.

Futaana

Machida

Saito

, et al. Moon-related nonthermal ions observed by nozomi: Species, sources, and generation mechanisms. J Geophys Res 2003; 108: 1025.

24.

Mouhammadoul

Alim

Tiofack

CGL

, et al. On the super positron-acoustic rogue waves in q-nonextensive magnetoplasmas. Phys Fluids 2023; 35: 054109.

25.

Alhejaili

Mouhammadoul

Alim

TCGL

, et al. Modulational instability and associated breathers in collisional electronegative Non-Maxwellian plasmas. Phys Fluids 2023; 35: 103101.

26.

Mouhammadoul

Alim

TCGL

Tiofack

CGL

, et al. Modulated positron-acoustic waves and rogue waves in a magnetized plasma system with nonthermal electrons and positrons. Astrophys Space Sci 2020; 365: 90–98.

27.

Douanla

Tiofack

CGL

Alim

, et al. Three-dimensional rogue waves and dust acoustic dark soliton collisions in degenerate ultradense magnetoplasma in the presence of dust pressure anisotropy. Phys Fluids 2022; 34: 087105.

28.

Douanla

Tiofack

CGL

Alim

, et al. Dynamics and head-on collisions of multidimensional dust acoustic shock waves in a self-gravitating magnetized electron-depleted dusty plasma. Phys Fluids 2023; 35: 023103.

29.

Thakur

Chinmay

. Chandra. Stationary structures in a four component dense magnetoplasma with lateral perturbations. IEEE Trans Plasma Sci 2022; 50: 1545.

30.

Sarkar

Chandra

Dey

, et al. Forced KdV and envelope soliton in magnetoplasma with kappa distributed ions. IEEE Trans Plasma Sci 2022; 50: 1565–1578.

31.

Sharry

Dutta

Ghosh

, et al. Magnetosonic shocks and solitons in fermi plasma with quasiperiodic perturbation. IEEE Trans Plasma Sci 2022; 50: 1585–1597.

32.

Chandra

Banerjee

Sarkar

, et al. Multistability studies on electron acoustic wave in a magnetized plasma with supra-thermal ions. J Astrophys Astron 2022; 43: 71.

33.

Singla

Chandra

Saini

. Simulation study of dust magnetosonic excitations in a magnetized dusty plasma. Chin J Phys 2023; 85: 524–533.

34.

Chandra

Manna

Mahanta

. Nuclear acoustic envelope soliton in a relativistically degenerate magneto-rotating stellar plasma. Phys Scripta 2024; 99: 095613.

35.

Chandra

Mahanta

Batani

, et al. Nonlinear evolutionary stages of a dispersive kappa distributed magnetized plasma. Int Conf on Nonlinear Dynamics and Applications 2024: 318–384.

36.

Chandra

Kapoor

Chaudhuri

, et al. Interaction and decomposition of magnetoacoustic stationary structures in magnetospheric plasma. IEEE Trans Plasma Sci 2025; 53: 184–198.

37.

Sarkar

Chandra

Goswami

, et al. Growth of RT instability at the accreting magnetospheric boundary of neutron stars. AIP Conf Proc 2021; 2319: 30006.

38.

Goswami

Sarkar

Chandra

, et al. Amplitude-modulated electron-acoustic waves with bipolar ions and kappa-distributed positrons and warm electrons. Pramana 2021; 95(2): 54.

39.

Samanta

Saha

Chatterjee

. Bifurcations of nonlinear ion acoustic travelling waves in the frame of a zakharov-kuznetsov equation in magnetized plasma with a kappa distributed electron. Phys Plasmas 2013; 20: 052111.

40.

Selim

El-Depsy

El-Shamy

. Bifurcations of nonlinear ion-acoustic travelling waves in a multicomponent magnetoplasma with superthermal electrons. Astrophys Space Sci 2015; 360: 66.

41.

Abdikian

. Modulational instability of ion-acoustic waves in magnetoplasma with pressure of relativistic electrons. Phys Plasmas 2017; 24: 052123.

42.

Herrmann

. Fractional calculus: an introduction for physicists. World Scientific, 2011.

43.

Meerschaert

Sikorskii

. Stochastic models for fractional calculus. De Gruyter, 2012.

44.

Metzeler

Klafter

. The random walk’s guide to anomalous diffusion: a fractional dynamics approach. Phys Rep 2000; 339: 1–77.

45.

Metzler

Klafter

. The restaurant at the end of the random walk: recent developments in the description of anomalous transport by fractional dynamics. J Phys A 2004; 37: 161–R208.

46.

Almuqrin

Tiofack

CGL

Mohamadou

, et al. On the “Tantawy Technique” and other methods for analyzing the family of fractional Burgers’ equations: applications to plasma physics. J Low Freq Noise Vib Act Control 2025. https://doi.org/10.1177/14613484251314580

47.

Almuqrin

Tiofack

CGL

Douanla

, et al. On the “Tantawy Technique” and other methods for analyzing fractional fokker Plank-type equations. J Low Freq Noise Vib Act Control 2025. https://doi.org/10.1177/14613484251319893

48.

El-Tantawy

Al-Johani

Almuqrin

, et al. Novel approximations to the fourth-order fractional cahn–hillard equations: application to the tantawy technique and other two techniques with yang transform. J Low Freq Noise Vib Act Control 2025. https://doi.org/10.1177/1461348425132224

49.

El-Tantawy

Bacha

SIH

Khalid

, et al. Application of the tantawy technique for modeling fractional ion-acoustic waves in electronegative plasmas having Cairns distributed-electrons, part (I): fractional KdV solitary waves. Braz J Phys 2025; 55: 123.

50.

El-Tantawy

Alhejaili

Khalid

, et al. Application of the tantawy technique for modeling fractional ion-acoustic waves in electronegative nonthermal plasmas, part (II): fractional modified KdV-Solitary waves. Braz J Phys 2025; 55: 176.

51.

El-Tantawy

Khan

, et al. A novel approximation to the fractional KdV equation using the tantawy technique and modeling fractional electron-acoustic cnoidal waves in a nonthermal plasma. Braz J Phys 2025; 55: 163.

52.

El-Dib

Elgazery

Alyousef

. The up-grating rank approach to solve the forced fractal duffing oscillator by non-perturbative technique. FU Mech Eng 2024; 22: 199–216.

53.

Habibirad

Hesameddini

Shekari

, et al. A meshless method based on moving kriging interpolation for the numerical solution of the transient flow of magnetohydrodynamic fractional maxwell fluid equation. J Appl Comput Mech 2024; 10: 853–866.

54.

Zuo

Luo

Zeng

. FU Mech Eng 2024; 22: 633–642.

55.

Babic

Marinkovic

. A new approach to determining the network fractality with application to robot-laser-hardened surfaces of materials. Fractal and Fractional 2023; 7: 710.

56.

Prigogine

Resibois

. On the kinetics of the approach to equilibrium. Physica 1961; 27: 629.

57.

Saichev

Zaslavsky

. Fractional kinetic equations: solutions and applications. Chaos 1997; 7: 753–764.

58.

Khalil

Al Horani

Yousef

, et al. A new definition of fractional derivative. J Comput Appl Math 2014; 264: 65–70.

59.

Kilbas

Srivastava

Trujillo

. Theory and applications of fractional differential equations. Elsevier, 2006.

60.

Hilfer

. Experimental evidence for fractional time evolution in glass forming materials. Chem Phys 2002; 284: 399–408.

61.

Yang

. Advanced local fractional calculus and its applications. World Science Publisher, 2012.

62.

Oldham

Spanier

. The fractional calculus. Academic Press, 1974.

63.

El-Wakil

Abulwafa

El-Shewy

, et al. Ion-acoustic waves in plasma of warm ions and isothermal electrons using time-fractional KdV equation. Chin Phys B 2011; 20: 040508.

64.

Sabi’u

Jibril

Gadu

. New exact solution for the (3+ 1) conformable space–time fractional modified korteweg-de-vries equations via sine-cosine method. J Taibah Univ Sci 2019; 13: 91–95.

65.

El-Wakil

Abulwafa

Elgarayhi

, et al. Effect of space-time fractional on the ion acoustic waves in electron-positron-ion plasma. Astrophys Space Sci 2014; 350: 591–598.

66.

Nazari-Golshan

. Investigation of shock waves in nonextensive electron–positron–ion plasma with relativistic ions. Plasma Phys Rep 2020; 46: 943–949.

67.

Uddin

Hafez

Hwang

, et al. Effect of space fractional parameter on nonlinear ion acoustic shock wave excitation in an unmagnetized relativistic plasma. Front Physiol 2022; 9: 766.

68.

Cox

Mortell

. The evolution of resonant water-wave oscillations. J Fluid Mech 1986; 162: 99–116.

69.

. Generation of upstream advancing solitons by moving disturbances. J Fluid Mech 1987; 184: 75–99.

70.

Akylas

. On the excitation of long nonlinear water waves by a moving pressure distribution. J Fluid Mech 1984; 141: 455–466.

71.

Morales

Lee

. Ponderomotive-force effects in a nonuniform plasma. Phys Rev Lett 1974; 33: 1016–1019.

72.

Moon

. Homoclinic crossings and pattern selection. Phys Rev Lett 1990; 64: 412–414.

73.

Sprott

. A proposed standard for the publication of new chaotic systems. Int J Bifur Chaos 2011; 21: 2391–2394.

74.

Nozaki

Bekki

. Chaos in a perturbed nonlinear schrödinger equation. Phys Rev Lett 1983; 50: 1226–1229.

75.

Hirsch

Smale

Devancy

. Differential equations, dynamical systems, and an introduction to chaos. Elsevier, 2004.

76.

Abdikian

Saha

Alimirzaei

. Bifurcation analysis of ion-acoustic waves in an adiabatic trapped electron and warm ion plasma. J Taibah Univ Sci 2020; 14: 1051–1058.

77.

Mouhammadoul

Tiofack

CGL

Alim

, et al. Positron-acoustic traveling waves solutions and quasi-periodic route to chaos in magnetoplasmas featuring Cairns nonthermal distribution. Eur Phys J D 2021; 75: 61.

78.

Saha

Tamang

. Qualitative analysis of the positron-acoustic waves in electron-positron-ion plasmas with deformed kaniadakis distributed electrons and hot positrons. Phys Plasmas 2017; 24: 082101.

79.

Tamang

Saha

. Dynamical behavior of supernonlinear positron-acoustic periodic waves and chaos in nonextensive electron-positron-ion plasmas. Z Naturforch A 2019; 74: 499–511.

80.

Saha

Pradhan

Natiq

. Multiperiodic and chaotic wave phenomena of collective ion dynamics under KP-type equation in a magnetised nonextensive plasma. Phys Scripta 2022; 97: 095604.

81.

Pradhan

Mukherjee

Saha

, et al. Multistability and chaotic scenario in a quantum pair-ion plasma. Z Naturforch A 2021; 76: 109–119.

82.

Prasad

Saha

. Bifurcation analysis of ion-acoustic waves for schrödinger equation in nonextensive solar wind plasma. Adv Space Res 2021; 67: 9–19.

83.

Bhowmick

Sahu

. Chaos and nonlinear excitations of dust acoustic waves in presence of external magnetic field with nonthermal species. Eur Phys J D 2021; 75: 273.

84.

Verheest

. Ambiguities in the tsallis description of non-thermal plasma species. J Plasma Phys 2013; 79: 1031–1034.

85.

Shalini

SNS

Misra

. Modulation of ion-acoustic waves in a nonextensive plasma with two-temperature electrons. Phys Plasmas 2015; 22: 092124.

86.

Podlubny

. Fractional differential equations. Academic Press, 1999.

87.

Petras

. Fractional-order nonlinear systems: modeling, analysis and simulation. Springer, 2010.

88.

Yorke

. Period three implies chaos. Am Math Mon 1975; 82: 985–992.

89.

Syta

Litak

Lenci

, et al. Chaotic vibrations of the duffing system with fractional damping. Chaos 2014; 24: 013107.