Our early universe is composed of completely ionized gases, which comprised of opposite charge particles of same masses and charge densities [1]. Negatively charged particles (electrons) and positively charged particles (positrons) are considered to be the fundamental composites of ionized gases. Electron–positron particles are main constituents of many astrophysical surroundings, such as, pulsar magnetospheres [2], active galactic nuclei [3], solar flares [4], fireballs producing γ -ray bursts [5], etc. Plasmas having electron and positron are detected in laboratory examinations [6 –8]. Electron–positron pair in plasmas may be generated during intense short laser pulse transmission [9]. Many astrophysical [4, 5] and laboratory [10 –13] observation shows that plasmas are composed of electrons, ions and positrons. Surko and Murphy [10] studied that for a large scale of parameters, desolation of electrons and positrons occur when 1 × 10¹² cm⁻³ is considered as density of electrons with temperature of 1 eV and the desolation of positron takes more than 1 s. Hence, the examination of electron–positron–ion (EPI) plasmas is beneficial to perceive features of astrophysical and laboratory environments. Lately, transmissions of different waves in plasmas bearing composition of EPI have gained immense attention of researchers [14 –19].

Recently, investigations of supernonlinear waves (SNWs) [20, 21] were initiated in plasmas. Dubinov et al [22] examined the properties of SNWs in multi-component plasma environments. The authors reported that SNWs can exist in plasma environments when there are at least three different plasma components. The multi-component plasma environment shows complexity which prompts existence of SNW features. In the form of plasma waves, SNWs exist in different physical environments such as, electrostatic and MHD waves [23]. All SNWs are greatly anharmonic waves with a substantially large amplitude [24]. The studies of solitary solution of SNWs known as supersolitons were reported earlier [25 –28], but were not specified as supersolitons. Verheest et al [29 –31] studied the existence of supersolitons in different three-component plasma systems. Singh et al [32] studied the presence of positive potential supersolitons with nonlinear potential solitons. Rufai et al [33] investigated a plasma limitation suggesting that the plasma system should have at least three-component for the transmission of SNWs. Lakhina et al [34] examined plasma system with three-components namely, two kinds of ions and Maxwellian electrons, observed absence of SNWs. Therefore, the necessary condition for SNWs to restrain in a plasma system is composition of three or more plasma components. Recently, El-Wakil et al [35] reported transmission of SNWs in four-constituent nonextensive dusty plasmas.

The kinetic motions of a charged particle are governed by Maxwellian for a plasma system exhibiting low range of interactions and collisions. However, due to existence of stationary states in some system, charged particles of the system tend to diverge from the standard Maxwellian distribution and converge to non-Maxwellian distributions. Some of non-Maxwellian distribution are q -nonextensive, kappa, nonthermal distributions, etc. Tsallis [36] proposed the Boltzmann–Gibbs–Shannon q -nonextensive entropy. In plasmas, q- nonextensive distribution has huge implication on plasma constituents that are widely used in astrophysical and cosmological environments such as, collisionless thermal plasma, hadronic matter, stellar polytropes, quark-gluon plasma, etc. The strength of nonextensivity q varies as −1 < q  < 1 and q  > 1. It is not normalizable for q  < −1 and tends to Maxwellian distribution for $q\to 1$ . Tsallis [37, 38] established the nonextensive distribution for systems that show extensive interactivity and collisions, and are not definable under the traditional Boltzmann–Gibbs statistics. The nonextensitivity of combined (A  + B) is projected by $S_q^{(A+B)}=S_q^{(A)}+S_q^{(B)}+(1-q)S_q^{(A)}{S}_q^{(B)}$ , where A and B are two distinct systems. This formulation classify the basis of system dynamics generalized by Tsallis. The entropy generalized by Tsallis is associated to plasma dynamics [39, 40], Hamiltonian systems [41] with extensive interaction, and nonlinear gravitational model [42]. Moreover, Liu et al [43] and Lima et al [39] studied plasma systems that follow velocity distribution other than Maxwellian. Many works [44 –48] were examined to study the Tsallis q -entropy and to enhance the theory of generalized statistics.

Nonlinear dynamics and chaotic phenomenon have been widely applied in various fields of engineering and science such as, biomedical engineering [49, 50], communication and cryptography [51 –55]. Of late, it has been investigated that dynamical systems can produce a series of different stable features for a particular set of parameters. The corresponding phenomena is defined as the existence of coexisting attractors in the system. Coexisting attractors have been observed in several real world phenomenon [56 –61]. However, the study of coexisting attractors in plasmas, as a new area of research in a chaotic system, is still in its beginning. The multistability is not reported much in plasma systems. In the work [62], it has been investigated that under the evidence of external forcing, the butterfly pattern can be distorted in a plasma perturbation model. In this article, we have proposed that the nonautonomous system can produce chaos and various coexisting attractors, which are quantified using phase space and Lyapunov exponents.

Lately, the application of bifurcation of planar dynamical systems [63] has enriched the study of underlying dynamics of plasma systems. Many physical systems are constantly in dynamic motions, application of bifurcation concept helps to analyze transitions and instabilities occurring in the system. The best possible way to describe the bifurcation analysis is through phase plane analysis which can be controlled by system parameters [63]. The bifurcation theory was employed in plasma for the first to describe nonlinear waves by Samanta et al [64]. Employing bifurcation theory, many researchers [65, 66] investigated various acoustic waves in plasmas. The features of supersolitons [67 –70] were examined using bifurcation concept in multi-component plasmas. Very recently, few works [71 –73] were published on supernonlinear periodic waves in plasmas. Tamang and Saha [73] reported positron-acoustic superperiodic waves and its dynamic motions in plasmas. Recently, Saha et al [74] reported bifurcation and coexisting attractors under the dynamics of Schrödinger equation in quantum plasmas. However, coexisting orbits of ion-acoustic waves (IAWs) using the planar dynamical systems are not reported in classical EPI plasma as far as we have surveyed. Therefore, the main objective of this study is to present superperiodicity, chaos and coexisting orbits of IAWs in a multi-component nonextensive plasma containing cold fluid ions, cold electrons, hot electrons and positrons.

The article is systematized as follows: in section 2, plasma model is considered. In section 3, conservative dynamical system (CDS) and superperiodic wave feature are obtained. In section 4, quasiperiodic and chaotic orbits under external perturbation are examined. In section 5, coexisting orbits are shown. Section 6 provides conclusions of the work.

We inspect the propagation of IAWs in a multi-component plasma system containing fluid ions, Maxwellian hot positrons and q -nonextensive cold and hot electrons. The following normalized model equations [75] represent the dynamics of IAWs: where $n_p={\mathrm{e}}^{-\gamma \varphi }$ is number density of hot positrons. Here, n and n _h (n _c ) are number densities of ions and hot electrons (cold electrons), respectively. Here, ion velocity and electrostatic potential are denoted by v and φ, respectively. γ is given by $T_{\mathrm{e}\mathrm{f}\mathrm{f}}/T_p$ with effective temperature $T_{\mathrm{e}\mathrm{f}\mathrm{f}}=T_c/({\mu }_c+{\mu }_h\beta )$ , where T _h , T _c and T _p denote temperatures of hot electrons, cold electrons and hot positrons, respectively. An equilibrium density ratio α is specified by $\left(\frac{n_{\mathrm{p}0}}{n_{\mathrm{e}0}}\right)$ , where ‘0’ is used for unperturbed quantity. Here, n _h , n _c and n _p are normalized by n _e0 . v is normalized by $C_s={(T_{\mathrm{e}\mathrm{f}\mathrm{f}}/m_i)}^{1/2}$ . Here, x and t are normalized by ${\lambda }_D={(T_{\mathrm{e}\mathrm{f}\mathrm{f}}/4\pi n_{\mathrm{e}0}{e}^2)}^{1/2}$ and ${\omega }_{\mathrm{p}\mathrm{i}}^{-1}={(4\pi n_{\mathrm{e}0}{e}^2/m_i)}^{-1/2}$ , respectively, where m _i is mass of ion and electronic charge e .

The nonextensive velocity distribution of electrons is achieved by considering the following function [76]:

where the distribution f _e (v) satisfies the laws of thermodynamics and prosper the Tsallis entropy (q). Here, k _B is the Boltzmann constant, T _e , m _e and v are temperature, mass and velocity of electrons, respectively and C _q is the normalizing constant expressed as

and

Performing integration on f _e (v), electron number density with nonextensive distribution is acquired as:

Thus, the normalized cold and hot electron number densities can be expressed as where q  > −1, and it measures strength of nonextensivity. Here, $\beta =\frac{T_c}{T_h}$ , ${\mu }_c=\frac{n_{\mathrm{c}0}}{n_{\mathrm{e}0}}$ and ${\mu }_h=\frac{n_{\mathrm{h}0}}{n_{\mathrm{e}0}}$ . The charge-neutrality condition is given by μ _c  + μ _h  = 1.

It is effective to obtain all possible nonlinear waves of the considered plasma system. To study this, the model equations are changed into a CDS [77 –79] employing the wave alteration $\xi =x-Vt$ , where the speed of wave is denoted by V . Substituting this alteration to equation (1 ) and performing integration subjected to conditions n  = 1, v  = 0, as $\xi \to \pm \mathrm{\infty }$ , one can deduce Similarly, from equation (2 ) we obtain the following expression subjected to conditions v  = 0, φ  = 0, as $\xi \to \pm \mathrm{\infty }$ , From equations (6 ) and (7 ), we get Substituting equations (4 ), (5 ) and (8 ) in equation (3 ), we deduce We rewrite the above equation considering up to 4th degree term as where



Equation (10 ) can be framed into a dynamical system as: The Hamiltonian function H (φ, z) of the system (11 ) is given by If we consider $\overrightarrow{F}=(z,a\varphi +b{\varphi }^2+c{\varphi }^3+d{\varphi }^4),$ then Equation (13 ) shows that the system (11 ) is a CDS.

To achieve different qualitative phase portrait profiles, we first obtain all equilibrium points of the dynamical system (11 ) by $\frac{\mathrm{d}\varphi }{\mathrm{d}\xi }=0$ and $\frac{\mathrm{d}z}{\mathrm{d}\xi }=0$ , which give while $m=\frac{c}{d},n=\frac{b}{d}$ and $o=\frac{a}{d}$ . The CDS (11 ) has four equilibrium points at E ₀ (φ ₀ ,0), E ₁ (φ ₁ ,0), E ₂ (φ ₂ ,0), and E ₃ (φ ₃ ,0), while φ ₀  = 0, and φ ₁, φ ₂, φ ₃ are real roots of the equation ${\varphi }^3+m{\varphi }^2+n\varphi +o=0$ . The determinant of the Jacobian matrix at E _i (φ _i , 0), i  = 0, 1, 2, 3, of the CDS (11 ) is obtained as Using the concept of planar dynamical systems [78, 79], we attain saddle point and center about an equilibrium point of the CDS (11 ) for D  < 0 and D  > 0, respectively. In figure 1, we present all possible phase profiles of the CDS (11 ) depending on different possible sets of values of system parameters q α, γ, μ _c , β and  V . Figures 1 (a) and (b) show phase profiles of the CDS (11 ) for $q=-0.17$ and q  = −0.26, respectively, with $\alpha =0.9,$ γ  = 0.01, μ _c  = 0.06, β  = 0.05 and V  = 0.5. The phase profile in figure 1 (a) shows saddle point at the origin and as we increase the value of q in the range −1 < q  < 0 keeping other parameters fixed, we observe that the nature of stability at equilibrium points of the system (11 ) changes. Therefore, we observe a center at the origin in figure 1 (b). In this case, we find that there are equilibrium points featuring pairs of centers and saddle points, and four kind of orbits namely, nonlinear homoclinic and periodic orbits, supernonlinear periodic and homoclinic orbits of the CDS (11 ). Any particular orbit in phase space is associated to a particular wave. Therefore, IAWs in the EPI plasma system has nonlinear solitary and periodic IAW solutions corresponding to nonlinear homoclinic and periodic orbits. These types of nonlinear IAW solutions are already reported [80]. Also, we observe existence of supernonlinear solitary and periodic IAWs solutions corresponding to supernonlinear homoclinic and periodic orbit, respectively [71]. Figure 1 (c) show phase profiles of the CDS (11 ) for q  = 0.77 with α  = 0.62, γ  = 0.06, μ _c  = 0.03, β  = 0.07 and V  = 0.705. The phase portrait manifested in figure 1 (c) for the range 0 < q  < 1 presents a saddle at the origin. Also, we perceive the existence of nonlinear solitary and periodic IAW solutions corresponding to nonlinear homoclinic and periodic orbits, and supernonlinear solitary and periodic IAWs solutions corresponding to supernonlinear homoclinic and periodic orbit, respectively. Figures 1 (d) and (e) display phase profiles of the CDS (11 ) for q  = 1.7 and q  = 1.8, respectively, with $\alpha =0.9,\gamma =0.01,$ μ _c  = 0.06, β  = 0.01 and V  = 0.776 in the subextensive range of parameter q . It is evident from figures 1 (d) and (e) that the system (11 ) in the range of q  > 1 also shows the existence of nonlinear and supernonlinear periodic and solitary waves.

Now, we show the effects of system parameters on superperiodic wave solutions with respect to the distinct ranges of q . Hence, we consider parametric sets of figure 1 (a) for the range −1 < q  < 0, figure 1 (c) for the range 0 < q  < 1 and figure 1 (e) for the range q  > 1. In figures 2 (a)–(c), we display superperiodic wave solutions with effects of parameters q, β and γ, respectively, with parametric set of figure 1 (a) for the range −1 < q  < 0. Similarly, we also show effects of q, β and γ on superperiodic wave solutions in figures 2 (d)–(f) associated to figure 1 (c) for the range 0 < q  < 1, and in figures 2 (g)–(i) associated to figure 1 (e) for the range q  > 1. In this case, we observe that with increase in q, the amplitude of supernonlinear periodic IAWs decreases in the range −1 < q  < 0(see figure 2 (a)) and increases in the ranges 0 < q  < 1 and q  > 1 (see figures 2 (d) and (g)). With increase in β, the amplitude of supernonlinear periodic IAWs increases in the ranges −1 < q  < 0 and 0 < q  < 1 (see figures 2 (b) and (e)) and, decreases in the range q  > 1 (see figure 2 (h)). When γ increases, the amplitude of supernonlinear periodic IAWs increases in the range −1 < q  < 0 (see figure 2 (c)) and decreases smoothly in the ranges 0 < q  < 1 and q  > 1 (see figures 2 (f) and (i)).

In figures 3 (a)–(c), we display nonlinear periodic wave solutions with effects of parameters q, β and γ, respectively, with parametric set of figure 1 (a) for the range −1 < q  < 0. Similarly, we also show effects of q, β and γ on nonlinear periodic wave solutions in figures 3 (d)–(f) associated to figure 1 (c) for the range 0 < q  < 1, and in figures 3 (g)–(i) associated to figure 1 (e) for the range q  > 1. Here, we see that with increase in q, the nonlinear periodic IAWs become smooth in the ranges −1 < q  < 0 and q  > 1 (see figures 3 (a) and (g)) and spiky in the ranges 0 < q  < 1 (see figure 3 (d)). With increase in β, the nonlinear periodic IAWs become smooth in the ranges −1 < q  < 9 and 0 < q  < 1 (see figures 3 (b) and (e)) and spiky in the range q  > 1 (see figure 3 (h)). However, when γ increases, the nonlinear periodic IAWs become spiky in the ranges −1 < q  < 0 and 0 < q  < 1 (see figures 3 (c) and (f)) and smooth in the range q  > 1 (see figure 3 (i)).

Recently, the effect of the source perturbation on nonlinear waves under Gaussian shape is reported in the work [81]. But, the nonlinear perturbation as external forcing, can be of various types [82 –84]. In this study, we consider a nonlinear perturbation as f ₀ cos(ωξ). In presence of external periodic forcing f ₀ cos(ωξ), the system (11 ) can be expressed as: where f ₀ and ω denote strength and frequency of the external perturbation, respectively.

In figure 4, we display dynamical properties corresponding to the perturbed dynamical system (PDS) (17 ) by changing intensity of external periodic force (f ₀ ) and frequency (ω) considering different set of values from figures 1 (b)–(d). In figures 4 (a)–(c), we consider set of parameters same as figure 1 (b) with f ₀  = 0.005 at initial condition (0.17, 0) by varying (a) ω  = 0.001, (b) ω  = 0.005 and (c) ω  = 0.05. We observe that orbits in phase spaces figures 4 (a)–(c) transit from regular pattern covering surface to irregular patterns. Hence, it is observed that dynamical motion of the PDS (17 ) transits slowly from quasiperiodic behavior (figure 4 (a)) to chaotic behavior (figure 4 (c)) as we increase frequency of external periodic force (ω). In figures 4 (d) and (e), we show phase space for (d) ω  = 0.5 and (e) ω  = 2 of the PDS (17 ) considering f ₀  = 0.005 at initial condition (−0.004, 0) with set of parameters same as figure 1 (d). It is perceived that as we increase frequency (ω), orbits in phase spaces (figures 4 (d)–(e)) cover surface of space forming pattern of torus. Hence, we observe quasiperiodic behaviors in the subextensive range of parameter q . In figures 4 (f)–(g), we show phase space of the PDS (17 ) varying (f) f ₀  = 0.005 and (g) f ₀  = 0.01 for the range 0 < q  < 1 and considering ω  = 0.08 with set of parameters same as figure 1 (c). It is observed that the PDS (17 ) exhibits two types of quasiperiodic behaviors. From figure 4 (f), it is transparent that there exists torus-like band forming quasiperiodic patterns. On the other hand, figure 4 (g) shows that there exist triple connected torus exhibiting complex quasiperiodic behaviors. However, in figures 4 (h)–(l), we show phase space of the perturbed dynamical system (17 ) varying (h) ω  = 0.01, (i) ω  = 0.1, (j) ω  = 0.57, (k) ω  = 0.6, and (l) ω  = 0.67 at initial condition (0.3, 0) for f ₀  = 0.005 with the set of parameters same as figure 1 (c) in the range 0 < q  < 1. It is observed that the PDS (17 ) exhibits continuous change in dynamical behaviors for different values of ω . In figure 4 (h), we display the phase space in which orbits filling around like a toroid exhibit quasiperiodic behavior. The quasiperiodic behaviors shown in figures 4 (i) and (j) are patterns of single and double twisted toruses, respectively. As we increase value of ω, the distortion of torus continues to give more complex structures, such as chaos in phase space (figure 4 (k)). In phase space figure 4 (l), orbits spread out steadily around like a toroid covering the surface densely.

To set up the existence of chaos, the main effectual approach is the alteration of the Lyapunov exponents. The Lyapunov exponent plays as an excellent measure for determining the rate of exponential convergence and divergence of nearby orbits in phase space. For any orbit $\mathbf{X}(\mathbf{t})=(x_1(t),x_2(t),\dots ,x_n(t))$ of a continuous dynamical system with dimension n and consider in the phase space that two orbits beginning with initial conditions X ₀ and ${\mathbf{X}}_0^{\mathrm{\prime }}={\mathbf{X}}_0+\delta {\mathbf{X}}_0$ develop with time as the vectors X (t) and X (t) + δ X (t), respectively. Then, one can define the Lyapunov exponent [85] corresponding to that orbits as Here, d(X ₀, t) is distance between orbits X (t) and ${\mathbf{X}}^{\mathrm{\prime }}(t)$ . The main feature of Lyapunov exponent shows that for a chaotic dynamical system, at least one Lyapunov exponent is positive while, it is nonpositive for nonchaotic systems. Hence, the positive value of Lyapunov exponents indicates a chaotic motion. Using Wolf algorithm [86], we examine Lyapunov exponent for system (17 ). Here, we observe that the nonlinear dynamics of the dynamical system (17 ) varies effectively due to the impact of external force. Therefore, we investigate this change of system (17 ) by plotting the Lyapunov exponents with respect to frequency (ω) and the forcing parameter (f ₀ ) of an external periodic force in figure 5 corresponding to the phase space reported in figure 4 (c). The positive values of the Lyapunov exponents in figure 5 show the evidence of chaos for the different values of ω and f ₀ . However, the Lyapunov exponents at ω  = 0.05045, 0.05425 (in figure 5 (a)) and f ₀  = 0.0048, 0.00492 (in figure 5 (b)) tend to be zero describing the presence of quasiperiodic orbits of the dynamical system (17 ). It is important to notice that the system (17 ) still preserves the conservative behavior in presence of external force.

In this section, we investigate existence of coexisting orbits for IAWs in presence of the nonlinear controller f ₀ cos(ωξ). We observe that the system (17 ) can show coexisting periodic, quasiperiodic and chaotic orbits for IAWs.

In figure 6, existence of multistability for IAWs in the perturbed dynamical system (17 ) is shown with the help of numerical simulation. Figure 6 (a) shows coexistence of quasiperiodic (blue curve) and chaotic orbits (purple curve) for IAWs by changing initial condition from (0.17, 0) for the blue curve to (0.001, 0) for the purple curve, when other parameters are kept same as figure 4 (b). This shows that the system (17 ) is sensitive to initial condition. In figure 6 (b), by varying the initial condition from (−0.004, 0) to (0.5, 0), we observe that the system (17 ) transits from quasiperiodic to chaotic orbits for IAWs showing sensitive dependency to the initial condition. Therefore, figure 6 (b) depicts coexistence of quasiperiodic orbits (green curve) for IAWs at initial condition (−0.004, 0) and chaotic orbits (blue curve) for IAWs at initial condition (0.5, 0), while other parameters are same as figure 4 (d). Figure 6 (c) displays coexistence of periodic orbit for IAWs at initial condition (0.09, 0) (green curve) and on changing the initial condition to (0.3, 0), orbits in phase space fill torus spreading over the surface i.e. quasiperiodic orbits for IAWs is achieved (red curve) for parametric values same as figure 4 (h). Figure 6 (d) is enlarged view of the quasiperiodic orbits for IAWs shown in figure 6 (c) with the red curve.

In figure 7, time series plots of IAWs are shown for the system (17 ) corresponding to coexisting orbits shown in phase spaces of figure 6 . Figure 7 (a) represents time series plot of coexisting quasiperiodic (dark blue color) and chaotic orbits (purple color) for IAWs shown in figure 6 (a). Figure 7 (b) shows time series plot of coexisting quasiperiodic (green color) and chaotic orbits (blue color) for IAWs presented in figure 6 (b). Time series plots in figures 7 (a) and (b) clearly show dynamical properties namely, quasiperiodic and chaotic orbits for IAWs existing in the perturbed system (17 ). However, figure 7 (c) displays time series plot of coexisting periodic (green curve) and quasiperiodic orbits (red color) for IAWs shown in figures 6 (c) and (d).

Superperiodicity, chaos and coexisting orbits of the nonlinear IAWs have been reported in a multi-component plasma containing fluid ions, Maxwellian hot positrons, cold and hot electrons with q -nonextensive distribution. The dynamics of IA superperiodic waves, which can be quantified by the phase portrait and the time series plots have been discussed. The significant impacts of system parameters like, nonextensive parameter (q), temperature ratio of cold and hot electrons (β) and the ratio of effective temperature and positron temperature (γ) have been observed on superperiodicity and nonlinear periodicity of IAWs. It has been encountered that superperiodic IAWs become smooth with increase in q and spiky with increase in β and γ for the range −1 < q  < 0. The superperiodic IAWs have been spiky with increase in q and β, and less spiky with increase in γ for the range 0 < q  < 1. However for the range q  > 1, the superperiodic IAWs have been spiky with increase with q and, smooth when β and γ increase. Furthermore, it has been observed that nonlinear periodic IAWs show similar features with increase in q, β and γ . The nonlinear periodic IAW has been smooth with increase in q and β while, has been sharp with increase in γ in the range −1 < q  < 0. However, in the range 0 < q  < 1, the nonlinear periodic IAW has been sharp and with increase in q and γ whereas, it has been smooth with increase in β . For the range q  > 1, it has been observed that the nonlinear periodic IAWs have been smooth for increase in q and γ, and sharp with increase in β . Effect of a nonlinear periodic perturbation has been presented on various types of orbits, like periodic, quasiperiodic and chaotic orbits for IAWs. The coexistence of such orbits for IAWs has been reported for suitable parametric region in their physical domain with different initial values.