The propagation of nonlinear structures in electron–positron (EP) plasmas is quite different from that in ordinary electron–ion (EI) plasmas because of the same mass of both the species. In a number of interesting physical situations, EP plasmas have been observed and investigated. The examples include pulsar environments, supernovae, active galactic nuclei, and cluster explosions [1–4]. Many astrophysical plasmas are an admixture of mainly electrons and positrons having a small fraction of ions, whereas the converse is true for the laboratory plasmas where the positron concentration is very low.

Many approaches have been proposed to create laboratory EP plasmas over the past few decades [5–7]. In initial experiments, relativistic positrons from radioactive neon gas were trapped directly into a magnetic mirror (MM) [5, 6]. Tsytovich and Wharton [7] suggested trapping positrons in an MM from a LINAC source. Later on, Boehmer [8] employed cyclotron heating to trap moderated positrons from a radioactive source and to heat the trapped positrons to relativistic energies. The most successful experimental approach, to obtaining positron plasmas, however, has turned out to be the scattering from a buffer gas into a Penning trap. In this case, the positron density is large enough to study the collective modes in plasmas [9]. It is important to note that pair annihilation can take place in EP plasmas. It is, therefore, a prerequisite that, for the study of collective modes in such plasma, the annihilation time scale ought to be much longer than the time scale for plasma effects, which is typically the inverse of the plasma frequency [9]. It is quite common in astrophysical and laboratory plasmas to have ions in addition to electrons and positrons, and, therefore, it is important to understand the linear and nonlinear behavior in electron–positron–ion (EPI) plasmas. Many studies have been carried out in the past few decades to investigate the linear and nonlinear wave propagation in EPI plasmas that highlight the importance of the inclusion of positrons in the system [10–18].

The observations of geophysical and laboratory plasmas have indicated the presence of two populations of electrons, namely, cold and hot electrons. Examples include hot cathode discharge plasmas, solar wind at about 1 AU, turbulent plasma for thermonuclear interest, and strongly interacting beam plasma systems [19]. In collisionless plasmas of space, Earth’s bow shock [20], and in the up/downstream of interplanetary shocks [21], two-temperature plasmas have been observed. Moreover, this type of plasmas has been observed in heliospheric termination shock [22], planetary magnetospheres [23], and space plasmas [24]. Additionally, different missions such as Polar [25], Geotail [26], FAST [27], and S3 − 3 [28] observed and discovered the two populations of electrons in the auroral zone and magnetosphere. The investigation of the linear ion acoustic wave (IAW) in plasma of two population electrons has a significant impact on its characteristics [29]. Many authors studied the ion acoustic solitons in plasma of two population electrons and reported many interesting facts about them [30–32].

No well-established methods have been developed, so far, that can solve all nonlinear differential equations arising in mathematics, physics, engineering, and biological sciences. Concerted efforts have thus been made to develop various analytical techniques to study nonlinear phenomena. Some of these methods include the homogeneous balance method [33–35], hyperbolic tangent expansion method [36–38], trial function method [39, 40], nonlinear transformation method [41, 42], and sine–cosine method [43]. Although, soliton and shock solutions have been obtained, periodic solutions could not be found by these methods. Porubov et al. [44, 45] obtained some exact periodic solutions for some nonlinear wave equations using the Weierstrass elliptic function which involved complicated deductions. Fu et al. [46] proposed the Jacobi elliptic function expansion method (JEFEM) and applied it to some nonlinear wave problems. In addition to shock and solitary wave solutions, various periodic solutions, based on the Jacobi elliptic sine function, were also found [46]. To discuss all three possible waveforms (for traveling wave solution), we use JEFEM to present all possible analytic solutions of the mKdVE. These include both bounded and unbounded periodic solutions. In recent years, theory of the dynamical system has generated a lot of interest in the study of nonlinear equations arising in plasma physics and fluid mechanics [47–56].

Saha and Chatterjee [57] used bifurcation of the phase portrait to analyze the qualitative behavior of the dust ion acoustic traveling wave solution of the modified Kadomtsev–Petviashvili equation (mKPE). The authors showed the existence of periodic, solitary, and kink wave solutions using the theory of the dynamical system by investigating the nature of the critical points. Analytical traveling wave solution of the mKPE is also obtained, and the two approaches are shown to be consistent. This work [57] is a generalization of the work presented by Samanta et al [52] who identified two solutions instead of three, using bifurcation analysis. Tamang et al; [53] used the planar dynamical system approach to present nonlinear homoclinic and nonlinear periodic trajectories of the modified mKdVE and generalized Gardner equation (GGE). A few articles, in which Asit Saha [52, 53, 57] is one of the authors, bifurcation has been used to identify the wave profiles and their properties depending on the parameter values. It may be noted that the bifurcation of the wave structure, as mentioned in these articles, should not be confused with the bifurcation diagram. We know that the qualitative behavior of the phase portraits such as periodic, solitary, or kink depends on the characterizing parameters of the problem. The phase portraits are drawn for some discrete values of the parameters without any prior knowledge of their behaviors. The qualitative change occurring in the phase diagrams is observed visually depending on the closed, homoclinic, or heteroclinic orbits.

Let us first recall the procedure adopted in these studies. A model nonlinear partial differential equation (PDE) is initially transformed into a nonlinear ordinary differential equation (ODE) using a traveling wave solution. The ODE is then formulated in terms of the corresponding dynamical system, and the phase plane diagram is drawn near the critical points using linear stability analysis. The solution behavior is based on the nature of the critical points, and the wave behavior is predicted from the phase diagram. For example, the center in the phase space corresponds to the closed orbit representing the periodic wave, the homoclinic orbit representing the solitary wave, and heteroclinic orbits corresponding to the kink (or shock) wave front. The closed-domain traveling wave solution generally corresponds to these waveforms. The stability of the wave solution is determined from the critical point and the Jacobian matrix. The phase diagram and the corresponding wave nature are determined for specifically choosing the parameters.

The aforementioned analysis raises a few questions for further thinking: a) Can we choose a set of premeditated parameters which gives rise to the wave structure of our choice? b) Can we divide the whole range of a parameter into sub-intervals with different topological behaviors (if so)? Can we know the behavior of the solution for some choice of the parameters without actually experimenting it through the phase plane diagram? Can we guess the behavior for a continuous spectrum of the parameter of interest instead of its discrete values? The answer to these questions lies in stability and bifurcation theory, which, to our understanding, has not been given due importance. The issues raised previously are addressed in this work besides finding the exact analytical solution. After taking up the analytical solution and the solution behavior of mKdVE, we show that the two approaches lead to same conclusions and corroborate each other. We are lucky to find the exact solution here; however, in more complex problems involving higher-order non-linearities, where analytical solution is not possible, the qualitative solution gives us all the physics we need to know.

In the present paper, no heteroclinic orbit appears, and, therefore, no shock solution exists for all possible parameters of the mKdVE. The existence of periodic (bounded and unbounded) and solitary wave solutions is found both mathematically and qualitatively. The bifurcation points are identified at which the topological behavior of the wave form changes, that is, the switch of the dynamical behavior from periodic to solitary waves or the other way round. A complete range of the important parameters like α, β, δ, γ, μ, and ν have been specified from the bifurcation diagram predicting the behavior in advance without any fear of wrong conclusion or incomplete conclusion (leaving some values of the parameters unattended). In passing, we note that the homoclinic behavior mentioned in [57] is actually from 0.8 < q < 1.0, instead of, q > − 1 . The important conclusions and observations highlighting important physical features are discussed in the concluding section. The analysis presented here will be useful to explore the physical behaviors in the problems with greater complexities in some future works where the analytical solution is not a possibility.

The paper is organized as follows: in Section 2 and Section 3, we present a set of model equations and derive the mKdVE using the reductive perturbation technique (RPT) [58–60]. In Section 4, we present the analytical solutions of the mKdVE using the JEFEM. All possible analytic solutions including bounded periodic solutions, solitary solutions, and unbounded periodic solutions are found. Section 5 explores the qualitative behavior of the mKdVE using the planar theory of the dynamical system. In Section 5, we present the results and discussion of both the quantitative and qualitative analyses using the plasma parameters that are consistent with the satellite observations of space plasmas. The main findings of the paper are recapitulated in Section 7.

We consider multicomponent plasma comprising hot and cold electrons, warm adiabatic ions, and isothermal positrons. In order to study the nonlinear behavior of the IAWs, we employ the following set of normalized fluid equations [32]: ∂ t n i + ∂ x n i u i = 0 , (1) ∂ t u i + u i ∂ x u i = − ∂ x ϕ − 5 3 δ n i − 1 / 3 ∂ x n i , (2) ∂ x 2 ϕ = − ρ = n e − 1 − α n i − n p , (3) where n _e = n _h + n _c .

The cold and hot populations of electrons and positrons are assumed to be inertialess on the ion time scale and assumed to follow the Maxwell–Boltzmann distribution function and are given as [32] n p = exp − γ ϕ , (4) n c = μ ⁡ exp ϕ μ + ν β , (5) n h = ν ⁡ exp β ϕ μ + ν β , (6) and n e = 1 + ϕ + μ + ν β 2 2 μ + ν β 2 ϕ 2 + μ + ν β 3 6 μ + ν β 3 ϕ 3 + ⋯ . (7) In the aforementioned equations, n _i is the density of ions normalized by n _i0; u _i is the fluid velocity of the ion species normalized by ion acoustic speed C _i ; n _c0, n _h0, and n _i0 are the equilibrium densities of two electron populations and the ion component, respectively; ϕ is the electrostatic potential normalized by thermal potential T _eff /e; α is the equilibrium density ratio of positron to electron species; ν is the equilibrium density ratio of hot electrons to the total electron number density; μ is the equilibrium density ratio of cold electrons to the total electron number density; β is the ratio of cold to hot electron temperatures; and δ is the ratio of ion to electron effective temperature. In the aforementioned equations, time and space coordinates have been normalized by the inverse of the ion–plasma frequency ω p i − 1 and Debye length λ _D , respectively, whereas electron densities n _h and n _c are normalized by n _e0. Here, β = T _c /T _h , α = n _po /n _eo , μ = n _co /n _eo , ν = n _ho /n _eo , δ = T _i /T _eff , γ = T _eff /T _p , T _eff = T _h T _h /(μT _h + νT _c ), c i = ( T eff / m i ) 1 / 2 , ω p i − 1 = ( 4 π n e o e 2 / m i ) − 1 / 2 , and λ D = ( T eff / 4 π n e o e 2 ) 1 / 2 .

Before deriving the mKdV equation, we remark that in some cases, it so happens that the nonlinear term in the KdV equation becomes identically zero because of the value of plasma parameters. As a consequence, no solitary wave is possible. This necessitates to consider the next-order (higher-order) non-linearity to balance dispersion in the governing equation leading to the mKdV equation. Moreover, it is worthwhile to explore the effect of cubic non-linearity on the formation of nonlinear structures, which becomes possible with the mKdV equation. For that, the following stretched coordinates ξ and τ are introduced [32] ξ = ε x − λ t , τ = ε 3 t , (8) where λ indicates the phase velocity and ɛ is a small and real perturbation parameter. We next expand the physical variables in Eqs. (4)–(7) as follows: n i = 1 + ϵ n i 1 + ϵ 2 n i 2 + ϵ 3 n i 3 + ⋯ , u i = ϵ u i 1 + ϵ 2 u i 2 + ϵ 3 u i 3 + ⋯ , ϕ = ϵ ϕ 1 + ϵ 2 ϕ 2 + ϵ 3 ϕ 3 + ⋯ , ρ = ϵ ρ 1 + ϵ 2 ρ 2 + ϵ 3 ρ 3 + ⋯ . (9) Substituting Eqs (8)–(9) into Eqs (1)–(3) and equating the terms with same power of ϵ will yield the following equations in the lowest order, that is, O (ϵ ²) u i 1 = λ n i 1 = 3 λ 3 λ 2 − 5 σ ϕ 1 , λ = 1 − α 1 + α γ + 5 δ 3 . (10) For the next higher-order equations, that is, O (ϵ ³), and after substituting the values of n _i1, u _i1, and λ, we obtain u i 2 = 81 λ 3 − 18 λ 2 λ ̃ 2 ϕ 1 2 − 15 δ 2 λ ̃ 3 ϕ 1 2 + 3 λ ̃ ϕ 2 , (11) n i 2 = 81 λ 2 − 15 δ 2 λ ̃ 3 ϕ 1 2 + 3 λ ̃ ϕ 2 , (12) α γ 2 − μ + ν β 2 μ + ν β 2 + 81 1 − α λ 2 λ ̃ 3 − 15 δ 1 − α λ ̃ 3 = 0 , (13) where λ ̃ = ( 3 λ 2 − 5 δ ) .

To the next higher-order equation of ϵ, that is, O (ϵ ⁴), we obtain ∂ τ n i 1 − λ ∂ ξ n i 3 + ∂ ξ u i 2 n i 1 + ∂ ξ n i 2 u i 1 + ∂ ξ u i 3 = 0 , (14) ∂ τ u i 1 − λ ∂ ξ u i 3 + u i 1 ∂ ξ u i 2 + u i 2 ∂ ξ u i 1 + ∂ ξ ϕ 3 − 5 9 δ n i 2 ∂ ξ n i 1 + 5 3 δ ∂ ξ n i 3 − 5 9 δ n i 1 ∂ ξ n i 2 + 10 27 δ n i 1 2 ∂ ξ n i 1 = 0 , (15) and ∂ ξ 2 ϕ 1 + 1 − α n i 3 − α γ ϕ 3 − 1 6 α γ 3 ϕ 1 3 + α γ 2 ϕ 1 ϕ 2 − ϕ 3 − μ + ν β 2 ϕ 1 ϕ 2 μ + ν β 2 − 1 6 μ + ν β 3 ϕ 1 3 μ + ν β 3 = 0 . (16) Differentiating Eqs (11), (12) and using Eq. (13), the following mKdVE is obtained: ∂ τ ϕ 1 + A ϕ 1 2 ∂ ξ ϕ 1 + B ∂ ξ 3 ϕ 1 = 0 , (17) with A = 3 λ 2 − 5 δ 2 18 1 − α λ 6561 1 − α λ 4 2 λ ̃ 5 − 2565 1 − α δ λ 2 2 λ ̃ 5 − 486 1 − α λ 2 λ ̃ 4 + 225 δ 2 1 − α 2 λ ̃ 5 + 30 1 − α δ λ ̃ 4 − α γ 3 2 − 1 2 μ + ν β 3 μ + ν β 3 ] (18) and B = 3 λ 2 − 5 δ 2 18 1 − α λ . (19)

Considering the mKdV Eq. (17), we obtain ∂ τ ϕ 1 + A ϕ 1 2 ∂ ξ ϕ 1 + B ∂ ξ 3 ϕ 1 = 0 , (20) where ϕ 1 ≡ ϕ 1 ξ , τ .

Substituting the traveling wave solution ϕ ₁ = ψ(η), where η = k (ξ − uτ), into Eq. (20), the following Duffing equation is obtained [61–66]: ψ ̈ + p ψ + q ψ 3 = 0 , (21) where ψ ̈ ≡ d 2 ψ / d η 2 , p = − u / B k 2 , and q = A / 3 B k 2 .

For brevity, we take BC = A here and rewrite Eq. (20) in the following form: ∂ τ ψ + A ψ 2 ∂ ξ ψ + B ∂ ξ 3 ψ = 0 . (37) According to the Lie symmetry method, Eq. (37) is invariant under the translational symmetry. For the invariant traveling wave solution, taking η = k ξ − u τ in Eq. (37) yields the following standard Duffing equation Eqs (61)–(66). ψ ̈ = u B ψ − C 3 ψ 3 , (38) where k = 1. Differential Eq. (38) is now expressed as an autonomous system of first-order equations (dynamical system) ψ ̇ = z = f ψ , z , z ̇ = u B ψ − C 3 ψ 3 = g ψ , z . (39) The Hamiltonian function (total energy) of the aforementioned system reads as follows: H ψ , z = z 2 2 − u 2 B ψ 2 + C 12 ψ 4 . (40) The dynamical system Eq. 39 reveals three equilibrium points (stationary solutions): ψ 1 , z 1 = 0,0 , ψ 2 , z 2 = 3 u B C , 0 , ψ 3 , z 3 = − 3 u B C , 0 . (41) The nature of the equilibrium can be determined from the eigenvalues of the Jacobian matrix at the location of the critical points. The eigenvalues of the Jacobian matrix of stationary solutions turn out to be λ 1 = ± u B , λ 2,3 = ± i 2 u B . (42) The first stationary solution is a saddle, while the second and third are centers. The local behavior of the solution at the saddle point (hyperbolic point) matches with the global behavior of the Hartman–Grobman theorem. For the centers (non-hyperbolic point), the Hartman–Grobman theorem does not apply [67], and we use the energy method. Since the Hamiltonian function Eq. (40) is conserved d H d ψ = 0 , closed orbits are represented by the center. We make a remark that for the traveling wave solution, we expect only three types of the waves, namely, periodic, solitary, or the shock waves. Furthermore, we discuss the qualitative behavior of the solution for a continuous range of the characteristic parameters with the help of bifurcation diagrams (branches). The results of bifurcation analysis are provided in Figures 8F, 9F, 10.

The nonlinear ion acoustic waves in unmagnetized collisionless plasma consisting of two populations of electrons (cold and hot), positrons, and warm adiabatic ions are investigated. The mKdVE is derived using multi-scale analysis and solved by the Jacobi elliptic function expansion method. The bounded periodic and solitary solutions are obtained. In addition to that, unbounded periodic solutions (of little physical significance) also appear in the solution. No shock wave is found for any value of plasma parameters.

The dynamical system theory and the bifurcation analysis are applied to determine the qualitative behavior of the mKdVE. The behaviors are determined from the nature of the critical points and the corresponding phase plane diagrams. The centers and homoclinic orbits ensure periodic and solitary wave behaviors, whereas the absence of heteroclinic orbit confirms the non-existence of the shock wave and corroborates the analytical results. A complete convergence between the analytical results and the qualitative behavior is established. In the end, we conclude with the remark that the qualitative behavior is useful when the closed-form solutions of the nonlinear differential equations are known and crucial and when analytic solutions are not possible. The results presented here provide a general framework and can be used to study a variety of plasma waves in laboratory, space, and astrophysical plasmas.