1 Introduction

Nonlinear evolution equations, such as the Kadomtsev–Petviashvili equation (KPE), play a vital role in modeling wave dynamics in complex physical systems where linear models are inadequate [1,2,3,4,5]. These equations help describe solitons, shocks, and turbulence in two-dimensional media, making them essential in plasma physics, fluid dynamics, and nonlinear optics. In plasma environments, especially those involving electron-positron-ion (EPI) compositions, nonlinear wave phenomena are of great interest due to their presence in astrophysical, space, and laboratory conditions [6,7,8,9]. Positrons are commonly generated in the interstellar medium via cosmic-ray interactions with atoms, leading to EPI plasmas. Electrons and positrons have identical masses but opposite electric charges, enabling them to exhibit similar dynamical timescales. Accurate modeling requires an understanding of key plasma parameters. Typical ion and electron densities range from \(10^5\) to \(10^{20} \, \text {cm}^{-3}\), while positron densities are generally lower, often approximated as \(n_p \sim 10^{-3} n_e\). Electron and positron temperatures vary between \(1\) and \(100\, \text {eV}\), depending on the environment. Further details can be found in [10, 11].

Nonlinear waves, especially solitary structures, are fundamental to the study of plasma dynamics. Solitons, which retain their shape after interactions, arise due to the balance between nonlinearity and dispersion. The dynamics of not only single soliton propagation but also soliton collisions—whether head-on (\( \pi \) angle) or overtaking (0 angle)—are influenced by various plasma parameters [12,13,14,15,16,17,18,19,20,21]. Traditionally, particle distributions are assumed to be Maxwellian; however, satellite and laboratory observations confirm that many plasmas exhibit non-Maxwellian behaviors [23,24,25]. These include superthermal particles in space plasmas (e.g., solar wind, magnetotail), motivating the use of generalized distributions like the \((\alpha , q)\)-distribution.

Several studies have employed the KPE to investigate ion-acoustic wave dynamics in two-dimensional plasma systems [1, 13, 20, 21, 27,28,29,30]. For instance, Hafez et al. [1] conducted a detailed analysis of ion-acoustic solitons using the KPE framework. Mahmood and his collaborators [27, 28] and Dutta [29] derived the KPE to study two-dimensional electrostatic solitons in various plasma conditions, offering insights into the parametric regimes favorable for stable solitary structures. Hafeez-Ur-Rehman et al. [30] focused on weakly relativistic plasmas with nonthermal particles and showed that enhanced ion velocity significantly influences soliton propagation. Malik [31] investigated the role of electron inertia in relativistic plasmas and its effect on KP soliton properties. Furthermore, Zaman et al. [20] explored bifurcations, periodic ion-acoustic waves, and chaotic behavior in three-component unmagnetized plasmas using phase portraits in conjunction with the KPE. Akther and Hafez [21] examined the propagation of lump solitons in a (2 + 1)-dimensional relativistic plasma comprising q-distributed electrons and positrons, using the reductive perturbation and Hirota bilinear methods to derive and analyze soliton structures under various interaction conditions, including collisional effects. Kaur and Saini [32] further extended this work to dust ion acoustic solitons in non-Maxwellian multispecies plasmas, utilizing KP, MKP, and CKP equations to highlight how different plasma parameters affect solitary wave behavior. Recently, attention has shifted toward multi-soliton dynamics and overtaking collisions, which cannot be fully captured by the standard KPEs. To address these gaps, this study investigates:

  1. (i)

    a two-dimensional fluid model for an unmagnetized multi-species plasma with \((\alpha , q)\)-distributed [33] inertialess particles,

  2. (ii)

    derivation of the KPE, modified KPE (MKPE), and Gardner-KPE (G-KPE) using the reductive perturbation technique (RPT),

  3. (iii)

    construction of soliton and overtaking collisional multi-soliton solutions via the Hirota bilinear method (HBM) [34], and

  4. (iv)

    analysis of the influence of plasma parameters on the propagation and interaction of solitons.

The paper is organized as follows: Sect. 2 outlines the mathematical plasma model. Section 3 detail the derivations and solutions for the KPE, MKPE, and G-KPE, respectively. Section 4 discusses the results, and Sect. 5 concludes the study.

2 Two-dimensional mathematical model equations

This study considers a two-dimensional, nonlinear propagation in a completely ionized, unmagnetized, multi-component plasma system. The plasma consists of hot ions, nonextensive positrons, and nonextensive electrons. The system is governed by the condition: \(n_{i0} + n_{p0} = n_{e0}\), where \(n_{i0}\), \(n_{e0}\), and \(n_{p0}\) represent the number densities of ions, electrons, and positrons, respectively, in the unperturbed state (as described by Hafez et al. [1]). In thermal equilibrium, the Maxwellian velocity distribution is typically used to represent the most probable particle distribution. However, in space plasmas, the particle velocity distributions often show non-Maxwellian features, such as superthermal tails, which follow a power-law decay with velocity. Given the extensive nonthermal nature of inertialess charged particles, their thermal pressure is often modeled using the \((\alpha ,q)\)-distribution. Using this distribution function, the electron and positron density functions can be expressed accordingly as

$$\begin{aligned} \left. \begin{array}{cccc} \mathscr {N}_e=n_{e0} \left[ 1+(q-1)\left( \frac{e\psi }{\mathscr {T}_e}\right) \right] ^\frac{q+1}{2(q-1)} \times \left[ 1-\mathscr {K}_1\left( \frac{e\psi }{\mathscr {T}_e}\right) +\mathscr {K}_2\left( \frac{e\psi }{\mathscr {T}_e}\right) ^2\right] \\ \mathscr {N}_p=n_{p0} \left[ 1-(q-1) \left( \frac{e\sigma \psi }{\mathscr {T}_p}\right) \right] ^\frac{q+1}{2(q-1)} \times \left[ 1+\mathscr {K}_1 \left( \frac{e\sigma \psi }{\mathscr {T}_p}\right) +\mathscr {K}_2 \left( \frac{e\sigma \psi }{\mathscr {T}_p}\right) ^2\right] \\ \end{array}\right\} , \end{aligned}$$
(1)

where \( \mathscr {K}_1=\frac{16\alpha q}{3-14q+15q^2+12\alpha }\) and \(\mathscr {K}_2=\frac{16\alpha q(2q-1)}{3-14q+15q^2+12\alpha }\), q is the nonextensivity and \(\alpha \) is the population of nonthermal particles index. Based on this assumption, the following normalized fundamental equations that govern the nonlinear dynamics of IAWs are derived.

$$\begin{aligned} & \frac{\partial \mathscr {N}_i}{\partial t} + \frac{\partial }{\partial x} \left( \mathscr {N}_i \mathscr {U}_i\right) + \frac{\partial }{\partial y}\left( \mathscr {N}_i \mathscr {V}_i\right) = 0, \end{aligned}$$
(2)
$$\begin{aligned} & \quad \left( \frac{\partial }{\partial t} + \mathscr {U}_i \frac{\partial }{\partial x} + \mathscr {V}_i \frac{\partial }{\partial y}\right) \mathscr {U}_i + \frac{\partial \psi }{\partial x} + \frac{\delta }{\mathscr {N}_i} \frac{\partial \mathscr {P}_i}{\partial x} = 0, \end{aligned}$$
(3)
$$\begin{aligned} & \quad \left( \frac{\partial }{\partial t} + \mathscr {U}_i \frac{\partial }{\partial x} + \mathscr {V}_i \frac{\partial }{\partial y}\right) V_i + \frac{\partial \psi }{\partial y} + \frac{\delta }{\mathscr {N}_i} \frac{\partial \mathscr {P}_i}{\partial y} = 0, \end{aligned}$$
(4)
$$\begin{aligned} & \quad \left( \frac{\partial }{\partial t} + \mathscr {U}_i \frac{\partial }{\partial x} + \mathscr {V}_i \frac{\partial }{\partial y}\right) \mathscr {P}_i + 3\mathscr {P}_i \frac{\partial }{\partial x} \mathscr {U}_i + 3\mathscr {P}_i \frac{\partial \mathscr {V}_i}{\partial y} = 0, \end{aligned}$$
(5)
$$\begin{aligned} & \quad \frac{\partial ^2 \psi }{\partial x^2} + \frac{\partial ^2 \psi }{\partial y^2} = 1 + \mathscr {Q}_1 \psi + \mathscr {Q}_2 \psi ^2 + \mathscr {Q}_3 \psi ^3 + \dots - \mathscr {N}_i, \end{aligned}$$
(6)

where

$$\mathscr {Q}_1 = \frac{(1 + \rho \,\sigma )}{1 - \rho } \left[ \frac{q+1}{2} - \mathscr {K}_1\right] ,$$
$$\mathscr {Q}_2 = \frac{(1 - \rho \sigma ^2)}{1 - \rho } \left[ \mathscr {K}_2 - \frac{q+1}{2} \mathscr {K}_1 + \frac{(q+1)(3-q)}{8}\right] ,$$
$$\mathscr {Q}_3= \frac{(1 + \rho \sigma ^3)}{1 - \rho } \left[ \frac{q+1}{2} \mathscr {K}_2 - \frac{(q+1)(3-q)}{8}\mathscr {K}_1 + \frac{(q+1)(3-q)(5-3q)}{48} \right] .$$

In this system, \(\mathscr {U}_i\) and \(\mathscr {V}_i\) represent the ion fluid velocities along the \(x\)- and \(y\)-axes, respectively, while \(\mathscr {P}_i\) is the ion pressure. \(\mathscr {T}_e\) and \(\mathscr {T}_p\) represent the temperatures of the electrons and positrons, respectively. The symbol \(e\) denotes the electric charge, and \(\psi \) is the electrostatic potential. The physical variables in the system are normalized as follows:(i) the ion number density \(\mathscr {N}_i\) is normalized by the unperturbed electron density \(n_{e0}\), (ii) \(\mathscr {U}_i\) and \(\mathscr {V}_i\), are normalized by \(\sqrt{\mathscr {T}_i e / \mathscr {M}_i}\), where \(\mathscr {M}_i\) is the ion mass, (iii) \(\psi \) is normalized by \(\mathscr {T}_e / e\), (iv) \(t\) is normalized by the plasma frequency \(\sqrt{\mathscr {M}_i / (4\pi n_{e0} e^2)}\) and (v) \(x\) and \(y\) are normalized by the Debye length \(\lambda _{Dc} = \sqrt{\mathscr {T}_e / (4\pi n_{e0} e^2)}\). Due the normalization, the ion-to-electron temperature ratio \(\delta = \mathscr {T}_i / \mathscr {T}_e\) and the positron-electron number density ratio \(\rho =\frac{\mathscr {N}_{p0}}{\mathscr {N}_{e0}}\) and electron to positron temperature ratio \(\sigma =\frac{\mathscr {T}_e}{\mathscr {T}_p}\) are obtained. Additionally, nextensivity index \(q\) is used to describe the distribution of the particles: (i) for \(q < 1\), the particles are superthermal, (ii) for \(q > 1\), the particles are subthermal, and (iii) \(q \rightarrow 1\), the particles follow an isothermal distribution with \(\alpha =0\). Otherwise, one can use \(0<\alpha <1\) and \(q=1\) for the nonthermality case. The geometrical configuration of the considered plasma assumption is displayed in Fig. 1.

Fig. 1
figure 1

2D representation of nonlinear wave propagation in a plasma, showing wavefronts, plasma, ion, electron, and positron regions, with associated velocity components and electrostatic potential gradient

Full size image

The considered plasma model consists of hot adiabatic ions, nonextensive electrons, and nonextensive positrons. This configuration is well-motivated by observations in both astrophysical and laboratory plasma environments. In particular, multi-component plasmas with significant positron populations and nonthermal particle distributions are characteristic of pulsar magnetospheres [35], active galactic nuclei [6], and the early universe [36], where high-energy processes produce abundant electron–positron pairs. Space missions such as Voyager, Cluster, and Wind [37] have reported non-Maxwellian velocity distributions in plasmas like the solar wind and Earth’s magnetosphere, which are effectively described by the nonextensive distribution. Similarly, laboratory experiments involving laser–plasma interactions [38] have observed superthermal tails and energetic particle populations, reinforcing the presence of nonthermal and non-Maxwellian features. In such environments, the temperature and density conditions typically satisfy \(\mathscr {T}_i < \mathscr {T}_e \ge \mathscr {T}_p\) and \(\mathscr {N}_{p0} < \mathscr {N}_{e0}\), which justify the choice of parameter ranges:

$$\begin{aligned} 0< \delta = \frac{\mathscr {T}_i}{\mathscr {T}_e}< 1, \quad \sigma = \frac{\mathscr {T}_e}{\mathscr {T}_p} \ge 1, \quad 0< \rho = \frac{\mathscr {N}_{p0}}{\mathscr {N}_{e0}} < 1. \end{aligned}$$

3 Mathematical analysis

3.1 Derivation of KPE

The derivation of evolution equations or the construction of their solutions often relies on perturbative methods. Among these, the RPT has proven to be a powerful tool, particularly effective for extracting nonlinear evolution equations—such as the KPEs—from complex plasma models [39,40,41]. This technique is especially well-suited for systems where wave dynamics and background fields evolve on different temporal and spatial scales. The nature of the perturbation variables in such models naturally lends itself to RPT, which simplifies the original equations while retaining the essential nonlinear and dispersive features. Unlike conventional perturbation methods, RPT is specifically designed to manage multi-scale behaviors, making it ideal for analyzing finite-amplitude wave propagation in plasma systems, where both slow and fast dynamics coexist. To capture the slow evolution of wave envelopes and account for different scales in space and time, we introduce stretched coordinates defined as follows [1, 20, 21]:

$$\begin{aligned} \xi = \epsilon (x - \lambda _0 t), \quad \eta = \epsilon ^{2} y, \quad \tau = \epsilon ^{3} t, \end{aligned}$$
(7)

where \(\lambda _0\) represents the linear phase velocity and \(\epsilon \) is a small expansion parameter that characterizes the weak dispersion and nonlinearity. The purpose of introducing these stretched coordinates is to rescale the space-time variables, enabling a clearer analysis of the system’s evolution, especially when multiple spatial and temporal scales are involved. We can expand physical quantities in terms of \(\epsilon \) as [1, 20, 21]:

$$\begin{aligned} \mathscr {N}_i & = 1 + \epsilon ^2 \mathscr {N}_{i1} + \epsilon ^4 \mathscr {N}_{i2} + \cdots , \nonumber \\ \mathscr {P}_i & = 1 + \epsilon ^2 \mathscr {P}_{i1} + \epsilon ^4 \mathscr {P}_{i2} + \cdots , \nonumber \\ \mathscr {U}_i & = \mathscr {U}_{i0} + \epsilon ^2 \mathscr {U}_{i1} + \epsilon ^4 \mathscr {U}_{i2} + \cdots , \nonumber \\ \psi & = \epsilon ^2 \psi _{i1} + \epsilon ^4 \psi _{i2} + \cdots , \nonumber \\ \mathscr {V}_i & = \epsilon ^3 \mathscr {V}_{i1} + \epsilon ^5 \mathscr {V}_{i2} + \cdots , \end{aligned}$$
(8)

where is ion streaming velocity. Next, we substitute Eqs. (7) and (8) into Eqs. (2)–(6) to express the model equations in terms of different powers of \(\epsilon \). The lowest power of \(\epsilon \) gives the following relations:

$$\begin{aligned} \mathscr {N}_{i1} = \frac{1}{\varrho _1} \psi _{i1} = \mathscr {Q}_1 \psi _{i1}, \quad \mathscr {U}_{i1} = \frac{\Lambda }{\varrho _1} \psi _{i1}, \quad \mathscr {P}_{i1} = \frac{3}{\varrho _1} \psi _{i1}, \end{aligned}$$
(9)

and

$$\begin{aligned} \lambda _0 = \mathscr {U}_{i0} + \sqrt{3\delta + \frac{1}{\mathscr {Q}_1}}, \varrho _1 = \left( \lambda _0 - \mathscr {U}_{i0} \right) ^2 - 3\delta , \quad \Lambda = \lambda _0 - \mathscr {U}_{i0}. \end{aligned}$$
(10)

The next power of \(\epsilon \) gives a system of partial differential equations involving nonlinearity (for simplicity, the nonlinearity is omitted here):

$$\begin{aligned} & \frac{\partial \mathscr {N}_{i1}}{\partial \tau } - (\lambda _0 - \mathscr {U}_{i0}) \frac{\partial \mathscr {N}_{i2}}{\partial \xi } + \frac{\partial \mathscr {U}_{i2}}{\partial \xi } + \frac{\partial (\mathscr {N}_{i1} \mathscr {U}_{i1})}{\partial \xi } + \frac{\partial \mathscr {V}_{i1}}{\partial \eta } = 0, \end{aligned}$$
(11)
$$\begin{aligned} & \quad \frac{\partial \mathscr {U}_{i1}}{\partial \tau } - (\lambda _0 - \mathscr {U}_{i0}) \frac{\partial \mathscr {U}_{i2}}{\partial \xi } + \mathscr {U}_{i1} \frac{\partial \mathscr {U}_{i1}}{\partial \xi } + \frac{\partial \psi _{i2}}{\partial \xi } + \delta \frac{\partial \mathscr {P}_{i2}}{\partial \xi } - \delta N_{i1} \frac{\partial P_{i1}}{\partial \xi } = 0, \end{aligned}$$
(12)
$$\begin{aligned} & \quad -(\lambda _0 - \mathscr {U}_{i0}) \frac{\partial \mathscr {V}_{i1}}{\partial \xi } + \frac{\partial \psi _{i1}}{\partial \eta } + \delta \frac{\partial \mathscr {P}_{i1}}{\partial \eta } = 0, \end{aligned}$$
(13)
$$\begin{aligned} & \quad \frac{\partial \mathscr {P}_{i1}}{\partial \tau } - (\lambda _0 - \mathscr {U}_{i0}) \frac{\partial \mathscr {P}_{i2}}{\partial \xi } + \mathscr {U}_{i1} \frac{\partial \mathscr {P}_{i1}}{\partial \xi } + 3 \frac{\partial \mathscr {U}_{i2}}{\partial \xi } + 3 \mathscr {P}_{i1} \frac{\partial \mathscr {U}_{i1}}{\partial \xi } + 3 \frac{\partial \mathscr {V}_{i1}}{\partial \eta } = 0, \end{aligned}$$
(14)
$$\begin{aligned} & \quad \frac{\partial ^2 \psi _{i1}}{\partial \xi ^2} = \mathscr {Q}_1 \psi _{i2} + \mathscr {Q}_2 \psi _{i1}^2 - \mathscr {N}_{i2}. \end{aligned}$$
(15)

By eliminating \(\mathscr {N}_{i2}\), \(\mathscr {U}_{i2}\), \(\mathscr {P}_{i2}\), and \(\psi _{i2}\) from the above equations, we derive the KPE as below.

$$\begin{aligned} \frac{\partial }{\partial \xi } \left( \frac{\partial \psi _{i1}}{\partial \tau } + \mathscr {R}_1 \psi _{i1} \frac{\partial \psi _{i1}}{\partial \xi } + \mathscr {R}_2 \frac{\partial ^3 \psi _{i1}}{\partial \xi ^3} \right) + \mathscr {R}_3 \frac{\partial ^2 \psi _{i1}}{\partial \eta ^2} = 0, \end{aligned}$$
(16)

where

$$\begin{aligned} \mathscr {R}_1 = \frac{3\Lambda }{2 \varrho _1} - Q_2 \frac{\varrho _1^2}{\Lambda } + \frac{3\delta }{2 \varrho _1 \Lambda }, \quad \mathscr {R}_2 = \frac{\varrho _1^2}{2 \Lambda }, \quad \mathscr {R}_3 = \frac{\Lambda }{2}. \end{aligned}$$
(17)

3.2 Soliton and overtaking collisional soliton solutions of KPE

To analyze the solutions of nonlinear evolution equations characterizing the collisionless propagation of wave dynamics, various mathematical techniques [2,3,4,5, 34, 42,43,44] can be employed. Among these, the HBM stands out as a powerful and systematic approach for deriving exact soliton solutions. It effectively captures both the dynamics of solitary wave propagation in collisionless regimes and the complex interactions involved in multi-soliton collisions. It involves rewriting the nonlinear equation into a bilinear form using the Hirota differential operators \(D_x, D_t\), defined by

$$ D_x^m D_t^n \{f \cdot g\} = \left( \frac{\partial }{\partial x} - \frac{\partial }{\partial x'}\right) ^m \left( \frac{\partial }{\partial t} - \frac{\partial }{\partial t'}\right) ^n f(x,t) g(x',t') \bigg |_{x'=x,\, t'=t}. $$

The dependent variable is then expressed through a transformation (e.g., \(\psi = 2 \frac{\partial ^2}{\partial x^2} \ln f\)) that converts the original nonlinear equation into bilinear form. Solutions are constructed by expanding \(f\) in a perturbation series and solving order-by-order to obtain multi-soliton solutions systematically.

To determine the soliton and multi-soliton solutions of Eq. (16) using the HBM, we start by defining the auxiliary function:

$$\begin{aligned} \psi _{i1}(\xi , \eta , \tau ) = \frac{12 \mathscr {R}_2}{\mathscr {R}_1} (\ln \mathscr {F})_{\xi \xi } = \frac{12 \mathscr {R}_2}{\mathscr {R}_1} \left[ \frac{\mathscr {F} \mathscr {F}_{\xi \xi } - \mathscr {F}_\xi ^2}{\mathscr {F}^2} \right] . \end{aligned}$$
(18)

Equation (16) can be rewritten using Hirota operators as:

$$\begin{aligned} (D_{\xi } D_{\tau } + \mathscr {R}_2 D_{\xi }^4 + \mathscr {R}_3 D_{\eta }^2)(\mathscr {F} \cdot \mathscr {F}) = 0. \end{aligned}$$
(19)

For a single soliton solution, we consider:

$$\begin{aligned} \mathscr {F} = 1 + e^{\theta _1}, \quad \theta _1 = k_1 \xi + l_1 \eta - \left( \mathscr {R}_2 k_1^3 + \frac{\mathscr {R}_3 l_1^2}{k_1} \right) \tau . \end{aligned}$$
(20)

Thus, the soliton solution is:

$$\begin{aligned} \psi _{i1}(\xi , \eta , \tau ) = \frac{12 \mathscr {R}_2}{\mathscr {R}_1} \left[ \ln \left( 1 + e^{\theta _1} \right) \right] _{\xi \xi }. \end{aligned}$$
(21)

For overtaking collision of two solitons, we assume \(\mathscr {F} = 1 + e^{\theta _1} + e^{\theta _2} + a(1,2) e^{\theta _1 + \theta _2}\), leading to

$$\begin{aligned} \psi _{i1}(\xi , \eta , \tau ) = \frac{12 \mathscr {R}_2}{\mathscr {R}_1} \left[ \ln \left( 1 +\sum _{i=1}^{2} e^{\theta _i} + a(1,2) e^{\sum _{i=1}^{2}\theta _i} \right) \right] _{\xi \xi }. \end{aligned}$$
(22)

For overtaking collision of three solitons, we take \(\mathscr {F} = 1 + e^{\theta _1} + e^{\theta _2} + e^{\theta _3} + a(1,2) e^{\theta _1 + \theta _2}\) \(+ a(1,3) e^{\theta _1 + \theta _3} + a(2,3) e^{\theta _2 + \theta _3} + a(1,2,3) e^{\theta _1 + \theta _2 + \theta _3}\), leading to

$$\begin{aligned} \psi _{i1}(\xi , \eta , \tau ) = \frac{12 \mathscr {R}_2}{\mathscr {R}_1} \left[ \ln \left( 1 +\sum _{i=1}^{3} e^{\theta _i} +\sum _{1 \le i < j \le 3} a(i,j) e^{\theta _i + \theta _j} + a(1,2,3) e^{\sum _{i=1}^{3}\theta _3} \right) \right] _{\xi \xi }. \end{aligned}$$
(23)

In the above equations, the parameters \(a(i,j)\) and \(a(1,2,3)\) follow the following expressions:

$$\begin{aligned} & a(i,j) = \frac{k_i \omega _j + k_j \omega _i + \mathscr {R}_2 \left( 4 k_i^3 k_j + 4 k_j^3 k_i - 6 k_i^2 k_j^2\right) + 2 \mathscr {R}_3 l_i l_j}{k_i \omega _j + k_j \omega _i + \mathscr {R}_2 \left( 4 k_i^3 k_j + 4 k_j^3 k_i + 6 k_i^2 k_j^2\right) + 2 \mathscr {R}_3 l_i l_j}, \quad i,j = 1,2,3.,\\ & \quad a(1,2,3) = a(1,2) \cdot a(1,3) \cdot a(2,3). \end{aligned}$$

The solutions presented above are valid when \(\mathscr {R}_1 \ne 0\), but they break down as \(\mathscr {R}_1 \rightarrow 0\). This highlights a significant issue: the current approach fails to address the case when \(\mathscr {R}_1\) approaches zero, necessitating a higher-order correction to derive a more robust nonlinear evolution equation. This oversight weakens the overall validity of the solution, and the need for further refinement should have been anticipated in the analysis.

Fig. 2
figure 2

Plot of nonlinear coefficient \(\mathscr {R}_1\) of KPE with regards to (a) \(\alpha \) (\(\rho =0.1\), \(\sigma =1\), \(q =0.7214188862\) and \(\delta =0.5\)), b q (\(\rho =0.1\), \(\sigma =1\), \(\alpha =0.1\) and \(\delta =0.5\)), c \(\rho \) (\(q=0.65\), \(\sigma =1\), \(\alpha =0.1\) and \(\delta =0.5\)) and (b) \(\delta \) (\(\rho =0.5\), \(\sigma =1\), \(\alpha =0\) and \(q=-0.87\))

Full size image

3.3 Derivation of MKPE

To study electrostatic soliton propagation near the critical values (CVs) of any specific parameters for which \(\mathscr {R}_1 \rightarrow 0\) with the fixed values of the remaining parameters, the expansion relations provided by Eq. (8) no longer hold. The existence CVs (highlighted by green color) are displayed in Fig. 2a, b. Figure 2 obviously showed that the KPE is no longer validated to describe the electrostatic soliton propagation at these CVs. Because the amplitude of solitons becomes infinity at these CVs. This limitation highlights the need for a different approach, where the state variables are expanded in the following form to enable further mathematical analysis:

$$\begin{aligned} \mathscr {N}_i &= 1 + \epsilon \mathscr {N}_{i1} + \epsilon ^2 \mathscr {N}_{i2}+ \epsilon ^3 \mathscr {N}_{i3} + \cdots \nonumber \\ \mathscr {P}_i & = 1 + \epsilon \mathscr {P}_{i1} + \epsilon ^2 \mathscr {P}_{i2}+ \epsilon ^3 \mathscr {P}_{i3} + \cdots \nonumber \\ \mathscr {U}_i & = \mathscr {U}_{i0} + \epsilon \mathscr {U}_{i1} + \epsilon ^2 \mathscr {U}_{i2}+ \epsilon ^3 \mathscr {U}_{i3} + \cdots \nonumber \\ \psi & = \epsilon \psi _{i1} + \epsilon ^2 \psi _{i2}+ \epsilon ^3 \psi _{i3} + \cdots \nonumber \\ \mathscr {V}_i & = \epsilon ^2 \mathscr {V}_{i1} + \epsilon ^3 \mathscr {V}_{i2} + \epsilon ^4 {\mathscr {U}}_{i3}+ \cdots , \end{aligned}$$
(24)

Now, one needs to compose the model equations in terms of various powers \(\epsilon \) by setting Eqs. (7) and (24) into Eqs. (2)–(6). The smallest power of \(O(\epsilon ^3)\), one obtains the similar equation as in Eq. (9). Hence, from the next order of \(\epsilon \), it gives

$$\begin{aligned} & \mathscr {N}_{i2} = \left( \frac{3\,\Lambda ^2}{2 \,\varrho _1^3} + \frac{3\delta }{2 \varrho _1^3}\right) \, \psi _{i1}^2 + \frac{1}{ \varrho _1} \psi _{12}=\mathscr {Q}_1 \psi _{i2} + \mathscr {Q}_2 \psi _{i1}^2, \end{aligned}$$
(25)
$$\begin{aligned} & \quad \mathscr {U}_{i2} = \left( \frac{\Lambda ^3}{2 \varrho _1^3} + \frac{9\delta \Lambda }{2 \varrho _1^3}\right) \psi _{i1}^2 + \frac{\Lambda }{ \varrho _1} \psi _{12},\, \mathscr {P}_{i2} = \left( \frac{6 \Lambda ^2}{ \varrho _1^3} + \frac{3}{2 \varrho _1^2}\right) \psi _{i1}^2 + \frac{3}{ \varrho _1} \psi _{12}, \end{aligned}$$
(26)

and

$$\begin{aligned} A_{cv}=\frac{3\,\Lambda ^2}{2 \,\varrho _1^3} + \frac{3\delta }{2 \varrho _1^3}-\mathscr {Q}_2\sim \mathscr {R}_1=0. \end{aligned}$$
(27)

Equation (27) is clearly shown that one can easily determine the CVs of any one parameter with the constant values of the remaining parameter for which the above equation holds. However, one can collect the next power of \(\epsilon \) equations, which yields the following equations:

$$\begin{aligned} & \frac{\partial \mathscr {N}_{i1}}{\partial \tau } - (\lambda _0 - \mathscr {U}_{i0}) \frac{\partial \mathscr {N}_{i3}}{\partial \xi } + \frac{\partial \mathscr {U}_{i3}}{\partial \xi } + \frac{\partial (\mathscr {N}_{i1} \mathscr {U}_{i2})}{\partial \xi } + \frac{\partial (\mathscr {N}_{i2} \mathscr {U}_{i1})}{\partial \xi }+ \frac{\partial \mathscr {V}_{i1}}{\partial \eta }= 0, \end{aligned}$$
(28)
$$\begin{aligned} & \quad \frac{\partial \mathscr {U}_{i1}}{\partial \tau } - (\lambda _0 - \mathscr {U}_{i0}) \frac{\partial \mathscr {U}_{i3}}{\partial \xi } + \frac{\partial (\mathscr {U}_{i1} \mathscr {U}_{i2})}{\partial \xi } + \frac{\partial \psi _{i3}}{\partial \xi } + \delta \frac{\partial \mathscr {P}_{i3}}{\partial \xi } \nonumber \\ & \quad - \delta \mathscr {N}_{i1} \frac{\partial \mathscr {P}_{i2}}{\partial \xi } - \delta \mathscr {N}_{i2} \frac{\partial \mathscr {P}_{i1}}{\partial \xi } + \delta \mathscr {N}_{i1}^2 \frac{\partial \mathscr {P}_{i1}}{\partial \xi } = 0, \end{aligned}$$
(29)
$$\begin{aligned} & \quad -(\lambda _0 - \mathscr {U}_{i0}) \frac{\partial \mathscr {V}_{i2}}{\partial \xi }+ \mathscr {V}_{i1} \frac{\partial \mathscr {V}_{i1}}{\partial \xi }- \delta \mathscr {N}_{i1} \frac{\partial \mathscr {P}_{i1}}{\partial \eta } + \frac{\partial \psi _{i2}}{\partial \eta } + \delta \frac{\partial \mathscr {P}_{i2}}{\partial \eta } = 0, \end{aligned}$$
(30)
$$\begin{aligned} & \quad \frac{\partial \mathscr {P}_{i1}}{\partial \tau } - (\lambda _0 - \mathscr {U}_{i0}) \frac{\partial \mathscr {P}_{i3}}{\partial \xi } + \mathscr {U}_{i2} \frac{\partial \mathscr {P}_{i1}}{\partial \xi }+ \mathscr {P}_{i1} \frac{\partial \mathscr {P}_{i2}}{\partial \xi } + 3 \frac{\partial \mathscr{U}_{i3}}{\partial \xi } \nonumber \\ & \quad + 3 \mathscr {P}_{i1} \frac{\partial \mathscr {U}_{i2}}{\partial \xi }+ 3 \mathscr {P}_{i2} \frac{\partial \mathscr {U}_{i1}}{\partial \xi } + 3 \frac{\partial \mathscr {V}_{i1}}{\partial \eta } = 0, \end{aligned}$$
(31)
$$\begin{aligned} & \quad \frac{\partial ^2 \psi _{i1}}{\partial \xi ^2} = \mathscr {Q}_1\, \psi _{i3} + \mathscr {Q}_3\, \psi _{i1}^3 + 2 \,\mathscr {Q}_2\,\psi _{i1} \,\psi _{i2} - \mathscr {N}_{i3}. \end{aligned}$$
(32)

Eliminating \(\mathscr{N}_{i3}\), \(\mathscr{U}_{i3}\), \(\mathscr{P}_{i3}\) and \(\psi _{i3}\) from the above equations along with the relations and conditions as in (9), and (25)–(27), the following MKPE is derived:

$$\begin{aligned} \frac{\partial }{\partial \xi } \left( \frac{\partial \psi _{i1}}{\partial \tau } + \mathscr {R}_{12}\, \psi _{i1}^2 \,\frac{\partial \psi _{i1}}{\partial \xi } + \mathscr {R}_2 \frac{\partial ^3 \psi _{i1}}{\partial \xi ^3} \right) + \mathscr {R}_3 \frac{\partial ^2 \psi _{i1}}{\partial \eta ^2}=0, \end{aligned}$$
(33)

where

$$\begin{aligned} \mathscr {R}_{12} = \frac{81\Lambda \delta }{4\,\varrho _1^3} + \frac{27\delta ^2}{4\Lambda \varrho _1^3} + \frac{3}{2\varrho \Lambda } + \frac{3\Lambda }{4\,\varrho } - \frac{3\varrho _1^2}{2\Lambda } \mathscr {Q}_3. \end{aligned}$$

        

3.4 Soliton and overtaking collisional soliton solutions of MKPE

To determine the soliton and multi-soliton solutions of Eq. (33) using the HBM, we start by defining the auxiliary function:

$$\begin{aligned} & \psi _{i1}(\xi , \eta , \tau ) = \sqrt{\frac{24 \mathscr {R}_2}{\mathscr {R}_{12}}} \,\frac{\partial }{\partial \xi } \left( \tan ^{-1} \frac{\mathscr {F}}{\mathscr {G}} \right) = \sqrt{\frac{24 \mathscr {R}_2}{\mathscr {R}_{12}}}\,\left( \frac{\mathscr {G} \mathscr {F}_x - \mathscr {F} \mathscr {G}_x}{\mathscr {F}^2 + \mathscr {G}^2}\right) \end{aligned}$$
(34)

Equation (33) can then be converted to the following form involving Hirota operators:

$$\left( D_\xi D_\tau + \mathscr {R}_2 D_\xi ^4 + \mathscr {R}_3 D_\eta ^2 \right) (\mathscr {F} \cdot \mathscr {G}) = 0,$$
$$\left( D_\xi D_\tau + \mathscr {R}_2 D_\xi ^4 + \mathscr {R}_3 D_\eta ^2 \right) \big (\mathscr {F} \cdot \mathscr {F} - \mathscr {G} \cdot \mathscr {G}\big ) = 0,$$
$$D_\xi ^2 \big (\mathscr {F} \cdot \mathscr {F} + \mathscr {G} \cdot \mathscr {G}\big )=0.$$

For a single soliton solution, we consider:

$$\begin{aligned} \mathscr {F} = e^{\theta _1}, \mathscr {G}=1, \quad \theta _1 = k_1\xi + l_1\eta -\left( \mathscr {R}_2\, k_1^3 + \frac{\mathscr {R}_3\, l_1^2}{k_1}\right) \tau . \end{aligned}$$

As a result, the single soliton propagation solution is obtained as

$$\begin{aligned} \psi _{i1}(\xi , \eta , \tau ) = \sqrt{\frac{24 \mathscr {R}_2}{\mathscr {R}_{12}}} \,\frac{\partial }{\partial \xi } \left( \tan ^{-1} \frac{e^{k_1\xi + l_1\eta -\left( \mathscr {R}_2\, k_1^3 + \frac{\mathscr {R}_3\, l_1^2}{k_1}\right) \tau }}{1} \right) . \end{aligned}$$
(35)

For overtaking collision of two solitons, we assume \(\mathscr {F} = e^{\theta _1} + e^{\theta _2},\,\)\( \mathscr {G}\\ = 1- a(1,2)\, e^{\theta _1 + \theta _2}\). As a result, the overtaking collision of two soliton propagation solution of MKPE is determined as

$$\begin{aligned} \psi _{i1}(\xi , \eta , \tau ) = \sqrt{\frac{24 \mathscr {R}_2}{\mathscr {R}_{12}}} \,\frac{\partial }{\partial \xi } \left[ \tan ^{-1}\left( \frac{e^{\theta _1} + e^{\theta _2}}{1-a(1,2)\, e^{\theta _1 + \theta _2}}\right) \right] , \end{aligned}$$
(36)

where

$$\theta _i = k_i\, \xi + l_i\, \eta - \left( \mathscr {R}_2\, k_i^3 + \frac{\mathscr {R}_3\, l_i^2}{k_i}\right) \tau ,\quad \omega _i = -\left( \mathscr {R}_2\, k_i^3 + \frac{\mathscr {R}_3 \,l_i^2}{k_i}\right) , \quad i = 1,2,$$
$$a(1,2) = \frac{k_1 \omega _2 + k_2 \omega _1 + \mathscr {R}_2 \left( 4 k_1^3 k_2 + 4 k_2^3 k_1 - 6 k_1^2 k_2^2\right) + 2 \mathscr {R}_3 l_1 l_2}{k_1 \omega _2 + k_2 \omega _1 + \mathscr {R}_2 \left( 4 k_1^3 k_2 + 4 k_2^3 k_1 + 6 k_1^2 k_2^2\right) + 2 \mathscr {R}_3 l_1 l_2}.$$

For overtaking collision of three solitons, we assume \(\mathscr {F} = \sum _{i=1}^{3}e^{\theta _i} - a(1,2,3)\, e^{\sum _{i=1}^{3}\theta _i }\), \( \mathscr {G}= 1- a(1,2)\, e^{\theta _1 - \theta _2}- a(1,3)\, e^{\theta _1 + \theta _3}- a(2,3)\, e^{\theta _2 + \theta _3}\). As a result, the overtaking collision of three soliton propagation solution of MKPE is determined as

$$\begin{aligned} \psi _{i1}(\xi , \eta , \tau ) = \sqrt{\frac{24 \mathscr {R}_2}{\mathscr {R}_{12}}} \,\frac{\partial }{\partial \xi } \left[ \tan ^{-1}\left( \frac{ \sum _{i=1}^{3}e^{\theta _i} - a(1,2,3)\, e^{\sum _{i=1}^{3}\theta _i }}{1- a(1,2)\, e^{\theta _1 - \theta _2}- a(1,3)\, e^{\theta _1 + \theta _3}- a(2,3)\, e^{\theta _2 + \theta _3}}\right) \right] , \end{aligned}$$
(37)

where

$$\theta _i = k_i\, \xi + l_i\, \eta - \left( \mathscr {R}_2\, k_i^3 + \frac{\mathscr {R}_3\, l_i^2}{k_i}\right) \tau ,\quad \omega _i = -\left( \mathscr {R}_2\, k_i^3 + \frac{\mathscr {R}_3 \,l_i^2}{k_i}\right) ,$$
$$a(i,j) = \frac{k_i \omega _j + k_j \omega _i + \mathscr {R}_2 \left( 4 k_i^3 k_j + 4 k_j^3 k_i - 6 k_i^2 k_j^2\right) + 2 \mathscr {R}_3 l_i l_j}{k_i \omega _j + k_j \omega _i + \mathscr {R}_2 \left( 4 k_i^3 k_j + 4 k_j^3 k_i + 6 k_i^2 k_j^2\right) + 2 \mathscr {R}_3 l_i l_j}, \quad i,j = 1,2,3$$
$$a(1,2,3) = a(1,2)\,\, a(1,3)\,\, a(2,3).$$

The solutions presented above are valid when \(\mathscr {R}_{12} \ne 0\), but they break down as \(\mathscr {R}_{12} \rightarrow 0\). This highlights a significant issue: the current approach fails to address the case when \(\mathscr {R}_{12}\) approaches zero, necessitating a higher-order correction to derive a more robust nonlinear evolution equation. This oversight weakens the overall validity of the solution, and the need for further refinement should have been anticipated in the analysis.

3.5 Derivation of G-KPE

A new form of KPE arises when \(\mathscr {R}_1\) is close to zero, but not exactly zero–specifically, when it is of the same order as \(\epsilon \) (\(A \approx O(\epsilon )\)). By applying the stretching coordinates from Eq. (7) and the expansion from Eq. (23), which were used in the derivation of MKPE, we analyze the equations from Eqs. (2) to (6) and equate the coefficients of different powers of \(\epsilon \). This leads to a series of evolution equations for the first, second, and third-order terms. At the second order of \(\epsilon \), we get the same equations as for KPE (see Eq. (9)), with the same phase velocity. At the third order of \(\epsilon \), the values of \(\mathscr {N}_{i2}\), \(\mathscr {U}_{i2}\), and \(\mathscr {P}_{i2}\) match those obtained in the MKPE derivation (see Eqs. (25), (26) and (27)). In our mathematical analysis, the nonlinear coefficient \(\mathscr {R}_1\) is expressed as: \( \mathscr {R}_1 = \frac{3\Lambda }{2\varrho _1} - \mathscr {Q}_2 \frac{\varrho _1^2}{\Lambda } + \frac{3\delta }{2\varrho _1\Lambda }, \) while the dispersion coefficient \(\mathscr {R}_2\) is: \( \mathscr {R}_2 = \frac{\varrho _1^2}{2\Lambda }, \) and the coefficient \(\mathscr {R}_3\) is: \( \mathscr {R}_3 = \frac{\Lambda }{2}. \) We can rewrite these as: \( \mathscr {R}_1 = \mathscr {A} \mathscr {B}_1, \quad \mathscr {R}_2 = \frac{\mathscr {A}}{2}, \quad \mathscr {R}_3 = \frac{\Lambda \mathscr {A}}{2}, \) where \(\mathscr {A} = \frac{\varrho _1^2}{\Lambda }\) and \(\mathscr {B}_1 = \frac{3\Lambda ^2}{2\varrho _1^3} - \mathscr {Q}_2 + \frac{3\delta }{2\varrho _1^3}\). It’s important to note that \(\mathscr {A}\) appears in all these coefficients and can never be zero. However, \(\mathscr {B}_1\) could potentially become zero. Now, inserting the expression of \(\mathscr {N}_{i2}\) as in Eq. (25), we can easily verify that the Poisson equation (6) at O(\(\epsilon ^2\)) is automatically satisfied. This is because the only nonzero term \(-\mathscr {B}_1 \psi _{i1}^2\) is the O(\(\epsilon ^3\)). Therefore, this term must be included in the Poisson equation at the next order, yields

$$\begin{aligned} \frac{\partial ^2 \psi _{i1}}{\partial \xi ^2} = \mathscr {Q}_1\, \psi _{i3} + \mathscr {Q}_3\, \psi _{i1}^3 + 2 \,\mathscr {Q}_2\,\psi _{i1} \,\psi _{i2} - \mathscr {N}_{i3}-\mathscr {B}_1\psi _{i1}^2. \end{aligned}$$
(38)

Simplifying Eqs. (28)–(31) and (38), we obtain the G-KPE in the following form:

$$\begin{aligned} \frac{\partial }{\partial \xi } \left( \frac{\partial \psi _{i1}}{\partial \tau } + \mathscr {R}_1 \psi _{i1} \frac{\partial \psi _{i1}}{\partial \xi } + \mathscr {R}_{12} \psi _{i1}^2 \frac{\partial \psi _{i1}}{\partial \xi } + \mathscr {R}_2 \frac{\partial ^3 \psi _{i1}}{\partial \xi ^3} \right) + \mathscr {R}_3 \frac{\partial ^2 \psi _{i1}}{\partial \eta ^2} = 0, \end{aligned}$$
(39)

where

$$\begin{aligned} \begin{aligned} \mathscr {R}_1&= \frac{3\,\Lambda }{{2\,\varrho _1}} - \mathscr {Q}_2\frac{\varrho _1^2}{\Lambda } + \frac{{3\delta }}{{2\,\varrho _1\Lambda }}, \quad \mathscr {R}_2 = \frac{\varrho _1^2}{2\Lambda }, \quad \mathscr {R}_3 = \frac{\Lambda }{2}, \\ \mathscr {R}_{12}&= \frac{81\Lambda \delta }{4\,\varrho _1^3} + \frac{27\delta ^2}{4\Lambda \varrho _1^3} + \frac{3}{2\varrho \Lambda } + \frac{3\Lambda }{4\,\varrho } - \frac{3\varrho _1^2}{2\Lambda } \mathscr {Q}_3. \end{aligned} \end{aligned}$$
(40)

3.6 Soliton and overtaking collisional soliton solutions of G-KPE

To determine the soliton and multi-soliton solutions of G-KPE, one can apply a suitable transformation that simplifies the equation and makes it more easily integrable. Thus, Eq. (39) is converted to the following form by considering \(\psi _{i1} = \phi - \frac{\mathscr {R}_1}{2 \mathscr {R}_{12}}\):

$$\begin{aligned} \frac{\partial }{\partial \xi } \left( \frac{\partial \phi }{\partial \tau } + \mathscr {R}_0 \frac{\partial \phi }{\partial \xi } + \mathscr {R}_{12} \phi ^2 \frac{\partial \phi }{\partial \xi } + \mathscr {R}_2 \frac{\partial ^3 \phi }{\partial \xi ^3} \right) + \mathscr {R}_3 \frac{\partial ^2 \phi }{\partial \eta ^2} = 0, \end{aligned}$$
(41)

where \(\mathscr {R}_0 = -\frac{\mathscr {R}_1^2}{4 \mathscr {R}_{12}}\). According to the HBM, one needs to consider firstly the following auxiliary function:

$$\begin{aligned} \phi (\xi , \eta , \tau ) = \sqrt{\frac{24 \mathscr {R}_2}{\mathscr {R}_{12}}} \,\frac{\partial }{\partial \xi } \left( \tan ^{-1} \frac{\mathscr {F}}{\mathscr {G}} \right) = \sqrt{\frac{24 \mathscr {R}_2}{\mathscr {R}_{12}}}\,\left( \frac{\mathscr {G} \mathscr {F}_x - \mathscr {F} \mathscr {G}_x}{\mathscr {F}^2 + \mathscr {G}^2}\right) . \end{aligned}$$
(42)

Equation (41) can then be converted to the following form involving Hirota operators:

$$\left( D_\xi D_\tau + \mathscr {R}_2 D_\xi ^4 + \mathscr {R}_0 D_\xi ^2 + \mathscr {R}_3 D_\eta ^2 \right) (\mathscr {F} \cdot \mathscr {G}) = 0,$$
$$\left( D_\xi D_\tau + \mathscr {R}_2 D_\xi ^4 + \mathscr {R}_0 D_\xi ^2 + \mathscr {R}_3 D_\eta ^2 \right) \big (\mathscr {F} \cdot \mathscr {F} - \mathscr {G} \cdot \mathscr {G}\big ) = 0,$$
$$D_\xi ^2 \big (\mathscr {F} \cdot \mathscr {F} + \mathscr {G} \cdot \mathscr {G}\big )=0.$$

For a single soliton solution, one needs to consider

$$\begin{aligned} \mathscr {F} = e^{\theta _1}, \mathscr {G}=1, \quad \theta _1 = k_1\xi + l_1\eta -\left( \mathscr {R}_0\, k_1+\mathscr {R}_2\, k_1^3 + \frac{\mathscr {R}_3\, l_1^2}{k_1}\right) \tau \end{aligned}$$

Hence, the single-soliton solution of G-KPE is determined as

$$\begin{aligned} \psi _{i1} = \sqrt{\frac{24 \mathscr {R}_2}{\mathscr {R}_{12}}} \,\frac{\partial }{\partial \xi } \left( \tan ^{-1} \frac{e^{k_1\xi + l_1\eta -\left( \mathscr {R}_0\, k_1+\mathscr {R}_2\, k_1^3 + \frac{\mathscr {R}_3\, l_1^2}{k_1}\right) \tau }}{1}\right) - \frac{\mathscr {R}_1}{2 \mathscr {R}_{12}}, \end{aligned}$$
(43)

For overtaking collision of two solitons, we assume \(\mathscr {F} = e^{\theta _1} + e^{\theta _2},\) \(\mathscr {G}= 1- a(1,2)\, e^{\theta _1 + \theta _2}\). Hence, the overtaking collision of two soliton propagation solution of G-KPE is determined as

$$\begin{aligned} \psi _{i1}(\xi , \eta , \tau ) = \sqrt{\frac{24 \mathscr {R}_2}{\mathscr {R}_{12}}} \,\frac{\partial }{\partial \xi } \left[ \tan ^{-1}\left( \frac{e^{\theta _1} + e^{\theta _2}}{1-a(1,2)\, e^{\theta _1 + \theta _2}}\right) \right] - \frac{\mathscr {R}_1}{2 \mathscr {R}_{12}}, \end{aligned}$$
(44)

where \(\begin{gathered}\theta _i = k_i\, \xi + l_i\, \eta - \left( \mathscr {R}_0\, k_i+\mathscr {R}_2\, k_i^3 + \frac{\mathscr {R}_3\, l_i^2}{k_i}\right) \tau, \hfill \\ \omega _i = -\left( \mathscr {R}_0\, k_i+\mathscr {R}_2\, k_i^3 + \frac{\mathscr {R}_3 \,l_i^2}{k_i}\right) , i = 1,2, \end{gathered}\)

$$\begin{aligned} a(1,2) = \frac{k_1 \omega _2 + k_2 \omega _1 + \mathscr {R}_2 \left( 4 k_1^3 k_2 + 4 k_2^3 k_1 - 6 k_1^2 k_2^2\right) +2 \mathscr {R}_0 k_1 k_2+ 2 \mathscr {R}_3 l_1 l_2}{k_1 \omega _2 + k_2 \omega _1 + \mathscr {R}_2 \left( 4 k_1^3 k_2 + 4 k_2^3 k_1 + 6 k_1^2 k_2^2\right) +2 \mathscr {R}_0 k_1 k_2 + 2 \mathscr {R}_3 l_1 l_2} \end{aligned}$$

For overtaking collision of three solitons, we assume \(\mathscr {F} = \sum _{i=1}^{3}e^{\theta _i} - a(1,2,3)\, e^{\sum _{i=1}^{3}\theta _i }\), \( \mathscr {G}= 1- a(1,2)\, e^{\theta _1 - \theta _2}- a(1,3)\, e^{\theta _1 + \theta _3}- a(2,3)\, e^{\theta _2 + \theta _3}\). Hence, the overtaking collision of three soliton propagation solution of G-KPE is determined as

$$\begin{aligned} \psi _{i1} = \sqrt{\frac{24 \mathscr {R}_2}{\mathscr {R}_{12}}} \,\frac{\partial }{\partial \xi } \left[ \tan ^{-1}\left( \frac{ e^{\theta _1} + e^{\theta _2} +e^{\theta _3}- a(1,2,3)e^{\theta _1 + \theta _2++ \theta _3}}{1- a(1,2)\, e^{\theta _1 + \theta _2}- a(1,3)\, e^{\theta _1 + \theta _3}- a(2,3)\, e^{\theta _2 + \theta _3}}\right) \right] - \frac{\mathscr {R}_1}{2 \mathscr {R}_{12}}, \end{aligned}$$
(45)

where

$$a(i,j) = \frac{k_i \omega _j + k_j \omega _i + \mathscr {R}_2 \left( 4 k_i^3 k_j + 4 k_j^3 k_i - 6 k_i^2 k_j^2\right) +2 \mathscr {R}_0 k_i k_j+ 2 \mathscr {R}_3 l_i l_j}{k_i \omega _j + k_j \omega _i + R_2 \left( 4 k_i^3 k_j + 4 k_j^3 k_i + 6 k_i^2 k_j^2\right) +2 \mathscr {R}_0 k_i k_j + 2 R_3 l_i l_j}, \quad i,j = 1,2,3,$$
$$a(1,2,3) = a(1,2)\,\, a(1,3)\,\, a(2,3).$$

4 Results and discussion

In this section, we analyze the nonlinear propagation and interaction of ion-acoustic solitons in a (2+1)-dimensional multi-component plasma, modeled through the KPE, MKPE, and G-KPE. The RPM is applied to derive these nonlinear evolution equations, capturing the interplay between nonlinearity and dispersion. The nonlinear and dispersive coefficients in the derived equations are found to depend on key plasma parameters: the ion-to-electron temperature ratio (\(\delta \)), the electron-to-positron temperature ratio (\(\sigma \)), the positron-to-electron density ratio (\(\rho \)), the nonextensivity index (\(q\)), the nonthermality parameter (\(\alpha \)), and the ion streaming velocity (\(\mathscr {U}_{i0}\)). The chosen parameter ranges reflect realistic conditions in some space and astrophysical plasmas, such as pulsar magnetospheres and the solar wind, based on observational constraints [10, 11, 45, 46]. Specifically, we adopt: \(0< \delta \le 0.5,\quad \sigma \ge 1,\quad 0< \rho < 1,\quad \mathscr {U}_{i0} = 10^3\,\text {cm/s}\).

Fig. 3
figure 3

Electrostatic KP soliton propagation with (a) \(\rho =0.5\), \(\sigma =1\), \(q=0.1\), \(\alpha =0.3\) and \(\delta =0.1\); b \(\rho =0.5\), \(\sigma =1\), \(q=0.71\), \(\alpha =0.3\) and \(\delta =0.1\); c \(\rho =0.5\), \(\sigma =1\), \(q=0.3\), \(\alpha =0.3\) and \(\delta =0.1\); d \(\rho =0.5\), \(\sigma =1\), \(q=1\), \(\alpha =0\) and \(\delta =0.1\); e \(\rho =0.5\), \(\sigma =1\), \(q=1\), \(\alpha =0.2\) and \(\delta =0.1\); f \(\rho =0.5\), \(\sigma =1\), \(q=5\), \(\alpha =0.2\) and \(\delta =0.1\); The remaining parameters are \(\mathscr {U}_{i0}=10^3 cm/s\), \(k_1=l_1=1\) and \(\tau \)= 0

Full size image
Fig. 4
figure 4

Electrostatic KP soliton propagation with (a) \(\rho =0.5\), \(\sigma =1\), \(q=0.1\), \(\alpha =0.3\) and \(\delta =0.1\); b \(\rho =0.3\), \(\sigma =1\), \(q=0.1\), \(\alpha =0.3\) and \(\delta =0.1\); c \(\rho =0.5\), \(\sigma =10\), \(q=0.1\), \(\alpha =0.3\) and \(\delta =0.1\), d \(\rho =0.5\), \(\sigma =1\), \(q=0.1\), \(\alpha =0.3\) and \(\delta =0.3\). The remaining parameters are \(\mathscr {U}_{i0}=10^3 cm/s\), \(k_1=l_1=1\) and \(\tau \)=0

Full size image

Figures 3a–f and 4a–d: KP soliton propagation: These figures illustrate how KP solitons evolve in the \(\xi \)–\(\eta \) plane at initial time (\(\tau = 0\)) under varying plasma conditions. The physical insights gained are:

  • KP solitons are localized structures whose stability arises from a balance between nonlinear steepening and dispersive spreading.

  • As \(q < 1\) (superthermality), the system has more high-energy electrons/positrons, which enhances nonlinearity, resulting in taller and narrower solitons.

  • For \(q > 1\) (subthermality), dispersion dominates, and solitons become wider and shallower.

  • Increasing \(\delta \) (hotter ions) amplifies ion pressure, which increases soliton amplitude and sharpens its structure due to stronger ion response.

  • Increasing \(\rho \) (more positrons) neutralizes the charge imbalance, reducing the potential well, hence lowering soliton amplitude and width.

  • Higher \(\sigma \) (hotter electrons relative to positrons) reduces net nonlinearity, weakening the nonlinear steepening and broadening the soliton.

Figures 5a–d, 6a–d, and 7a–d: Overtaking and time evolution of KP soliton collisions

Fig. 5
figure 5

Electrostatic overtaking collision of KP (a) compressive double soliton along with (b) its contour plot with \(\rho =0.5\), \(\sigma =1\), \(q=0.1\), \(\alpha =0.1\) and \(\delta =0.5\), \(\tau \)=3; and (c) rarefactive double soliton along with (d) its contour plot with \(\rho =0.01\), \(\sigma =1\), \(q=0.71\), \(\alpha =0.1\) and \(\delta =0.5\), \(\tau \)= 0.5. The remaining parameters are \(\mathscr {U}_{i0}=10^3 cm/s\), \(k_1=1,l_1=2,k_2=3,l_2=4\)

Full size image
Fig. 6
figure 6

Electrostatic overtaking collision of KP (a) compressive triple soliton along with (b) its contour plot with \(\rho =0.5\), \(\sigma =1\), \(q=0.1\), \(\alpha =0.1\) and \(\delta =0.5\), \(\tau =0;\) and (c) rarefactive triple soliton along with (d) its contour plot with \(\rho =0.01\), \(\sigma =1\), \(q=0.71\), \(\alpha =0.2\) and \(\delta =0.5\), \(\tau =0.1\). The remaining parameters are \(\mathscr {U}_{i0}=10^3cm/s\), \(k_1=1,l_1=2,k_2=3,l_2=4,k_3=5,l_3=6\)

Full size image
Fig. 7
figure 7

Variation of KP overtaking collisional two and three solitons propagation for different values of \(\tau \), that is (a) \(\tau \le 0\), b \(\tau > 0\), c \(\tau \le 0\) and d \(\tau > 0\). The remaining parameters are \(\rho =0.5\), \(\sigma =1\), \(q=0.1\), \(\alpha =0.1\) and \(\delta =0.5\), \(\mathscr {U}_{i0}=10^3 cm/s\), \(k_1=1,l_1=2,k_2=3,l_2=4,k_3=5,l_3=6\)

Full size image

These figures show the interaction of multiple KP solitons via overtaking collisions and their evolution over time.

  • In overtaking collisions, a faster, high-amplitude soliton overtakes a slower, smaller one. Such interactions demonstrate the integrable nature of the KP system.

  • Post-collision, solitons preserve their identity and shape, verifying their solitonic character and the non-destructive nature of the interaction.

  • The coexistence of compressive (positive potential) and rarefactive (negative potential) solitons reveals that the plasma system can support bipolar structures under certain parameter regimes.

  • Contour plots in the \(\xi \)–\(\eta \) plane visualize energy localization and interaction dynamics clearly.

  • The results emphasize soliton robustness and their relevance to coherent structure propagation in the considered plasma environments.

Figures 8a–f and 9a–b: MKP soliton propagation

Fig. 8
figure 8

Electrostatic MKP soliton propagation with (a) \(\rho =0.5\), \(\sigma =1\), \(q=0.1\), \(\alpha =0.1\) and \(\delta =0.5;\) b \(\rho =0.1\), \(\sigma =1\), \(q=0.1\), \(\alpha =0.1\) and \(\delta =0.5;\) c \(\rho =0.5\), \(\sigma =10\), \(q=0.1\), \(\alpha =0.1\) and \(\delta =0.5;\) d \(\rho =0.5\), \(\sigma =1\), \(q=0.3\), \(\alpha =0.1\) and \(\delta =0.5;\) e \(\rho =0.5\), \(\sigma =1\), \(q=0.1\), \(\alpha =0.1\) and \(\delta =0.1;\) f \(\rho =0.5\), \(\sigma =1\), \(q=0.1\), \(\alpha =0.3\) and \(\delta =0.5\). The remaining parameters are \(\mathscr {U}_{i0}=10^3 cm/s\), \(k_1=l_1=1\) and \(\tau =1\)

Full size image
Fig. 9
figure 9

Electrostatic MKP soliton propagation with (a) \(\rho =0.5\), \(\sigma =1\), \(q=1\), \(\alpha =0\) and \(\delta =0.5\); (b) \(\rho =0.1\), \(\sigma =1\), \(q=5\), \(\alpha =0\) and \(\delta =0.5\); (c) \(\rho =0.5\), \(\sigma =1\), \(q=1\), \(\alpha =0.1\) and \(\delta =0.5\), (d) \(\rho =0.5\), \(\sigma =1\), \(q=1\), \(\alpha =0.3\) and \(\delta =0.5\). The remaining parameters are \(\mathscr {U}_{i0}=10^3 m/s\), \(k_1=l_1=1\) and \(\tau \)= 1

Full size image

MKPE introduces higher-order nonlinearity and dispersion, which are critical near the CVs of plasma parameters.

  • Only compressive solitons are supported in the MKP regime, indicating stronger electron/positron response that sustains positive potential structures.

  • Higher-order nonlinear effects make MKP solitons sharper and more localized than their KP counterparts.

  • For \(q < 1\), solitons remain taller and narrower due to dominant superthermal particle effects near CVs.

  • Increasing \(\alpha \) increases soliton amplitude by raising the population of nonthermal high-energy particles.

  • Increasing \(\delta \) (hotter ions) decreases soliton amplitude in this case due to enhanced thermal spreading overpowering nonlinear compression.

  • Higher \(\rho \) and \(\sigma \) values reduce soliton amplitude and width, as excess positrons and hotter electrons weaken electrostatic potential formation.

Figures 10a–d and 11a–d: MKP overtaking soliton collisions and evolution

Fig. 10
figure 10

Electrostatic MKP overtaking collisional (a) double soliton and its (b) contour plot with \(\rho =0.5\), \(\sigma =1\), \(q=0.1\), \(\alpha =0.1\) and \(\delta =0.1\); and (c) triple and its (d) contour plot with \(\rho =0.5\), \(\sigma =1\), \(q=1\), \(\alpha =0.3\) and \(\delta =0.5\). The remaining parameters are \(\mathscr {U}_{i0}=10^3\,cm/s\), \(k_1=1,l_1=2,k_2=3,l_2=4,k_3=5,l_3=6\) and \(\tau = 0\)

Full size image
Fig. 11
figure 11

Variation of MKP overtaking collisional two and three solitons propagation for different values of \(\tau \), that is, (a) \(\tau \le 0\), (b) \(\tau > 0\), (c) \(\tau \le 0\) and (d) \(\tau > 0\). The remaining parameters are \(\rho =0.5\), \(\sigma =1\), \(q=0.1\), \(\alpha =0.1\) and \(\delta =0.5\), \(\mathscr {U}_{i0}=10^3 cm/s\), \(k_1=1,l_1=2,k_2=3,l_2=4,k_3=5,l_3=6\)

Full size image

These figures illustrate overtaking collisions of double and triple MKP solitons around CVs.

  • Collisions exhibit mild deformation due to higher-order effects but remain elastic in the long-term.

  • Solitons maintain structural identity and amplitude post-collision, indicating resilience in high-gradient plasma environments.

  • Stronger localization observed in MKP solitons reflects realistic sharp wavefronts.

Figures 12a–d and 13a–d: G-KP solitons and G-KP overtaking soliton collisions

Fig. 12
figure 12

Electrostatic (a) G-KP single and its (b) contour plot with \(\rho =0.1\), \(\sigma =1\), \(q=q_c=0.7214188862\), \(\alpha =0.1\) and \(\delta =0.5\), \(\tau \)= 1, (c) G-KP overtaking collisional double and its (d) contour plot with \(\rho =0.1\), \(\sigma =1\), \(q=q_c=0.7214188862\), \(\alpha =0.1\) and \(\delta =0.5\),\(\tau \)= 2; and (e) G-KP overtaking collisional triple and its (d) contour plot with \(\rho =0.1\), \(\sigma =1\), \(q=q_c=0.7214188862\), \(\alpha =0.1\) and \(\delta =0.5\),\(\tau \)= 1. The remaining parameters are \(\mathscr {U}_{i0}=10^3 cm/s\), \(k_1=1,l_1=2,k_2=3,l_2=4,k_3=5,l_3=6\)

Full size image
Fig. 13
figure 13

Variation of G-KP overtaking collisional two and three solitons propagation for different values of \(\tau \), that is, (a) \(\tau \le 0\), (b) \(\tau > 0\), (c) \(\tau \le 0\) and (d) \(\tau > 0\). The remaining parameters are \(\rho =0.1\), \(\sigma =1\), \(q=q_c=0.7214188862\), \(\alpha =0.1\) and \(\delta =0.5\), \(\mathscr {U}_{i0}=10^3 cm/s\), \(k_1=1,l_1=2,k_2=3,l_2=4,k_3=5,l_3=6\)

Full size image

The Gardner KP equation accounts for both KPE and MKPE features, providing a more general framework.

  • G-KP solitons encapsulate the effects of both quadratic and cubic nonlinearities, giving rise to a broader class of solitary structures.

  • These solitons exhibit stronger amplitude and sharper width than both KP and MKP solitons, especially near CVs.

  • Overtaking collisions of G-KP solitons reveal significant energy redistribution during interaction, but eventual recovery of soliton identity, a hallmark of soliton solutions in complex plasmas.

  • The system supports only compressive solitons under G-KPE, indicating a dominant unipolar potential field in high-gradient astrophysical regimes.

  • The amplitude and width of G-KP soliton and overtaking G-KP solitons—including double, and triple soliton collisions at the composition of CVs is vary higher than the amplitude and width of KP, and MKP soliton and their overtaking solitons collisions.

  • These results reinforce the applicability of G-KPE to extreme plasma scenarios with mixed thermal and nonthermal populations.

The propagation and interaction of ion-acoustic solitons in the present plasma model arise from a balance between the ion inertia (driving force) and the thermal pressure of electrons and positrons (restoring force). The plasma parameters significantly influence soliton characteristics and their overtaking collision dynamics. The nonextensivity parameter \(q\) modulates the particle distribution: for \(q < 1\), superthermal populations dominate, enhancing nonlinearity and producing sharper, taller solitons; for \(q > 1\), solitons become broader and shallower due to weaker nonlinearity. An increase in the nonthermality index \(\alpha \) intensifies the restoring force from energetic particles, thus amplifying soliton amplitude. Overall, the results reveal how plasma nonthermality, nonextensivity, temperature ratios, and density ratios fundamentally shape soliton characteristics, including amplitude, width, stability, and collision behavior. The physical interpretations align with observations in space and laboratory plasmas, where superthermal populations and nonlinear structures play key roles in energy transport and particle acceleration mechanisms.

5 Conclusion

This study examined the nonlinear propagation and interaction of ion-acoustic solitons in a two-dimensional plasma system consisting of hot ions, nonextensive electrons, and nonextensive positrons. By employing the reductive perturbation method, we derived three evolution equations—the KPE, MKPE, and G-KPE—to capture nonlinear wave dynamics across different regimes.

Our analysis reveals that soliton amplitude, width, and profile are strongly influenced by plasma parameters such as the temperature ratios (\(\delta \), \(\sigma \)), density ratio (\(\rho \)), nonextensivity (\(q\)), nonthermality (\(\alpha \)), and ion streaming velocity (\(\mathscr {U}_{i0}\)). Superthermal particles (\(q < 1\)) and high \(\alpha \) enhance nonlinearity, leading to taller, narrower solitons. In contrast, subthermal conditions (\(q > 1\)) and increased positron effects (\(\rho \), \(\sigma \)) reduce soliton strength. The interplay between ion inertia and electron-positron pressure governs the observed wave features. Overtaking collisions were investigated across all regimes. KP solitons showed ideal soliton behavior–retaining shape and speed post-collision–while MKP solitons, influenced by higher-order nonlinearities, exhibited slight deformation but preserved coherence. G-KP solitons, integrating features from both KPE and MKPE, demonstrated the most robust profiles and supported large-amplitude multi-soliton structures.

Overall, the results confirm that ion-acoustic solitons in such multi-component plasmas are not only stable but also highly responsive to physical conditions. The model and findings are relevant to astrophysical and laboratory plasma environments, offering insight into nonlinear wave propagation, collision behavior, and soliton stability in complex plasma systems.