Introduction

Propagation of nonlinear waves in astrophysical plasmas has generated a lot of interest in the plasma community: a large number of investigations are on-going on nonlinear waves, such as solitons, shocks, double layers, and so on, which are observed in space, astrophysical, and laboratory plasmas.

Most of these studies have focussed on the nonlinearity of the ion acoustic wave, which is an important wave in plasma. The first nonlinear analysis of the ion acoustic wave was by Sagdeev [1]; its first experimental observation was by Ikezi et al. [2]. Outside of the laboratory, the observations by Viking and Freja spacecrafts have identified solitary structures in the magnetosphere as density depressions.

A nonlinear and dispersive medium always supports a solitary wave. In a soliton structure, the nonlinearity is balanced by the dispersion. The Korteweg–deVries (KdV) equation, having negative nonlinear and positive dispersive effects, will give rise to rarefactive solitons. When both nonlinearity and dispersive effects are positive, the KdV equation may result in a compressive soliton. Such compressive solitons have been studied, for example, in a moving electron–positron plasma containing positively and negatively charged dust [3]. However, a medium having significant dissipative effect and dispersion supports the formation of shock waves instead of solitons which is best described by the Korteweg–deVries–Burgers (KdVB) equation. The dissipative Burger term in the nonlinear KdVB equation arises due to dissipative mechanisms, such as wave–particle interactions, turbulence, dust charge fluctuations in a dusty plasma, multi-ion streaming, Landau damping, anomalous viscosity, etc. [4–6]. When wave breaking due to nonlinearity is balanced by the combined effect of dispersion and dissipation, a monotonic or oscillatory dispersive shock wave is generated in a plasma [7]. If the dissipative effect is negligible the KdVB equation transforms into the KdV equation which permits a soliton solution. The KdVB equation has thus been extensively used to study the properties of solitons and shock waves in dusty plasmas [8–12].

The presence of high energy particles in the tail of a distribution causes plasmas observed in different space environments to deviate significantly from the well-known Maxwellian distribution. Using solar wind data, Vasyliunas first predicted a non-Maxwellian distribution [13]; this distribution, which later came to be known as the “kappa distribution,” has been found in many magnetospheric and astrophysical environments. Ion-acoustic shock waves in plasmas consisting of superthermal electrons and positrons have been studied recently by many authors [14–20]; the effects of electrons and positrons on ion acoustic solitary waves have also been studied [21].

A cometary plasma is composed of hydrogen, and new born heavier ions and electrons with relative densities depending on their distances from the nucleus. Initially, positively charged oxygen ions were treated as the main heavier ion [22, 23]. However, the discovery of negatively charged oxygen ions [24] enables one to consider the plasma environment around a comet as a pair-ion plasma (O+, O−) with other ions (both lighter and heavier) constituting the other components of the plasma.

Thus, a cometary plasma is a true multi-ion plasma consisting of both lighter and heavier ions and electrons with different temperatures. We thus model our plasma as consisting of a pair-ion plasma of oxygen ions, lighter hydrogen ions, and two components of electrons with different temperatures. The lighter hydrogen ions and electrons are modeled by kappa distributions.

A very complex structure of multiple sub-shocks and interplanetary structures of comet Halley was reported by Giotto. In addition, unambiguous observations of bow shock crossings were provided by Vega-1 [25]. These shock structures are seriously affected by heavier ions due to mass loading of the solar wind and pickup driven ion instabilities [25]. And these heavy ions can be both positively and negatively charged as discussed above. In addition, in a recent study, Voelzke and Izaguirre, analyzing 886 images of comet Halley, identified 41 solitary structures [26].

The stability of ion acoustic waves in a four component dusty plasma has been investigated recently by several authors. Early studies on nonlinear wave propagation in four component dusty plasmas were by Sakanaka and Spassovska [27] and Verheest [28]. Various aspects of dust acoustic waves have been considered in the studies by several authors [29–36]. Other nonlinear electrostatic waves have been considered by Ghosh et al. [37, 38] and Dutta et al. [39].

For reasons given at the middle of this section, a reasonably accurate modeling of a cometary plasma requires at least five components. We, therefore, investigate the existence of ion-acoustic shock waves in a five component cometary plasma consisting of positively and negatively charged oxygen ions, kappa described hydrogen ions, hot solar electrons and slightly colder cometary electrons. Related studies of the effects of two types of electrons on ion-acoustic solitary waves are those by Khaled [40] and Shahmansouri et al. [41]. Similarly, the dust acoustic shock wave has also been studied in a dusty plasma with kappa distributed ions [11]. The KdVB equation is derived for this system. It is found that the amplitude of shock wave decreases with increasing kappa values. Its strength decreases with increasing temperatures of positively charged oxygen ions and densities of negatively charged oxygen ions.

Basic equations

We consider the existence of ion-acoustic shock waves in a five component plasma consisting of negatively and positively charged oxygen ions (represented, respectively, by subscripts ‘1’ and ‘2’), kappa described hydrogen ions, hot electrons of solar origin and colder electrons of cometary origin. At equilibrium, charge neutrality requires that

$$n_{{{\text{ce}}0}} + n_{{{\text{se}}0}} + Z_{1} n_{10} = n_{{{\text{H}}0}} + Z_{2} n_{20} .$$

In the above relation, \(n_{\text{ce0}}\) and \(n_{\text{he0}}\) represent the equilibrium densities of cometary electrons and solar electrons, respectively. In addition, \(n_{10}\), \(n_{20}\), and \(n_{{{\text{H}}0}}\) are, respectively, the equilibrium densities of negatively charged oxygen (O−) ions, positively charged oxygen (O+) ions, and hydrogen ions. Z1 and Z2 denote the charge numbers of O− and O+ ions, respectively.

The kappa distribution of species ‘s’ is given by:

$$n_{\text{s}} = n_{\text{s0}} \left[ {1 + \frac{{e_{\text{s}} \varphi }}{{k_{\text{B}} T_{\text{s}} (\kappa_{s} - {\raise0.7ex\hbox{$3$} \!\mathord{\left/ {\vphantom {3 2}}\right.\kern-0pt} \!\lower0.7ex\hbox{$2$}})}}} \right]^{{ - \kappa_{\text{s}} + {\raise0.7ex\hbox{$1$} \!\mathord{\left/ {\vphantom {1 2}}\right.\kern-0pt} \!\lower0.7ex\hbox{$2$}}}} .$$
(1)

In (1), s = H for hydrogen, \(= {\text{se}}\) for solar electrons, and \(= {\text{ce}}\) for cometary photo-electrons.\(n_{\text{s}}\) denotes the density (with the subscript ‘0’ denoting the equilibrium value), \(e_{\text{s}}\) is the charge, \(T_{\text{s}}\) is the temperature, \(\kappa_{\text{s}}\) is the spectral index for species ‘s’, \(k_{\text{B}}\) is Boltzmann’s constant, and \(\varphi\) is the potential.

The dynamics of the heavier ions can be described by the following hydrodynamic equations:

$$\frac{{\partial n_{j} }}{\partial t} + \frac{{\partial (n_{j} v_{j} )}}{\partial x} = 0,$$
(2)
$$\left( {\frac{\partial }{\partial t} + v_{j} \frac{\partial }{\partial x}} \right)v_{j} = \mp \frac{{Z_{j} e}}{{m_{j} }}\frac{\partial \varphi }{\partial x} - \frac{1}{{m_{j} n_{j} }}\frac{{\partial P_{j} }}{\partial x} + \mu_{j} \frac{{\partial^{2} v_{j} }}{{\partial x^{2} }},$$
(3)

where ‘−’ sign refers to positively charged oxygen ions (and vice versa) and \(v_{j}\) and \(m_{j}\), respectively, denote the fluid velocity and mass of the j-species of ions (j = O−, O+). In (3), the adiabatic equation of state for ions is \(\frac{{P_{j} }}{{P_{j0} }} = \left( {\frac{{n_{j} }}{{n_{j0} }}} \right)^{\gamma } = N_{j}^{\gamma }\), where \(P_{j0} = n_{j0} k_{\text{B}} T_{j}\), γ = (N + 2)/N for an N dimensional system, and \(\mu_{j}\) is the ion kinematic viscosity. Here, we are considering a one-dimensional system, and hence, γ = 3.

Poisson’s equation is given by

$$\frac{{\partial^{2} \varphi }}{{\partial x^{2} }} = - 4\pi e\left( {n_{\text{H}} + Z_{2} n_{2} - Z_{1} n_{1} - n_{\text{ce}} - n_{\text{se}} } \right).$$
(4)

We normalize (2)–(4) using the parameters of O− ions according to \(\varphi = \frac{e\phi }{{k_{\text{B}} T_{1} }}\), \(V_{j} = \frac{{v_{j} }}{{c_{\text{s}} }}\), and \(\tau = \frac{t}{{\omega_{{{\text{p}}1}}^{ - 1} }}\), where \(c_{\text{s}} = \left( {\frac{{Z_{1} k_{\text{B}} T_{1} }}{{m_{1} }}} \right)^{1/2}\) and \(\omega_{{{\text{p}}1}} = \left( {\frac{{4\pi Z_{1}^{2} e^{2} n_{10} }}{{m_{1} }}} \right)^{1/2}\). The variable x is normalized using \(\lambda_{{{\text{D}}1}} = \left( {\frac{{Z_{1} k_{\text{B}} T_{1} }}{{4\pi Z_{1}^{2} e^{2} n_{10} }}} \right)^{1/2}\), while \(N_{j} = \frac{{n_{j} }}{{n_{j0} }}\).

Thus, Eqs. (2)–(4) can be rewritten as:

$$\frac{{\partial N_{1} }}{\partial \tau } + \frac{{\partial (N_{1} V_{1} )}}{\partial x} = 0,$$
(5)
$$\frac{{\partial N_{2} }}{\partial \tau } + \frac{{\partial (N_{2} V_{2} )}}{\partial x} = 0,$$
(6)
$$\frac{{\partial V_{1} }}{\partial \tau } + V_{1} \frac{{\partial V_{1} }}{\partial x} = \frac{\partial \phi }{\partial x} - \frac{{3N_{1}^{{}} }}{{Z_{1} }}\frac{{\partial N_{1} }}{\partial x} + \rho_{1} \frac{{\partial^{2} V_{1} }}{{\partial x^{2} }},$$
(7)
$$\frac{{\partial V_{2} }}{\partial \tau } + V_{2} \frac{{\partial V_{2} }}{\partial x} = \frac{{ - Z_{2} m}}{{Z_{1} }}\frac{\partial \phi }{\partial x} - \frac{{3m\beta N_{2}^{{}} }}{{Z_{1} }}\frac{{\partial N_{2} }}{\partial x} + \rho_{2} \frac{{\partial^{2} V_{2} }}{{\partial x^{2} }},$$
(8)

where \(m = \frac{{m_{1} }}{{m_{2} }}\), \(\beta = \frac{{T_{2} }}{{T_{1} }}\), \(\rho_{1} = \frac{{\mu_{1} }}{{\omega_{\text{p1}} \lambda_{\text{D}}^{ 2} }}\), and \(\rho_{2} = \frac{{\mu_{2} }}{{\omega_{{{\text{p}}1}} \lambda_{\text{D}}^{2} }}.\)

\(\rho_{1}\) and \(\rho_{2}\) now represent the normalized kinematic viscosities of the pair ions.

The normalized Poisson’s equation after substitution of (1) is:

$$\begin{aligned} \frac{{\partial^{2} \phi }}{{\partial x^{2} }} & = N_{1} - N_{2} \left( {1 + \mu_{\text{ce}} + \mu_{\text{se}} - \mu_{\text{H}} } \right) + \mu_{\text{ce}} \left( {1 - \frac{\phi }{{\sigma_{ce} \left( {\kappa_{ce} - 3/2} \right)}}} \right)^{{ - (\kappa_{\text{ce}} - 1/2)}} \\ & \quad + \mu_{\text{se}} \left( {1 - \frac{\phi }{{\sigma_{\text{se}} \left( {\kappa_{\text{se}} - 3/2} \right)}}} \right)^{{ - (\kappa_{\text{se}} - 1/2)}} - \mu_{\text{H}} \left( {1 + \frac{\phi }{{\sigma_{\text{H}} \left( {\kappa_{\text{H}} - 3/2} \right)}}} \right)^{{ - (\kappa_{\text{H}} - 1/2)}} , \\ \end{aligned}$$
(9)

where \(\mu_{\text{ce}} = \frac{{n_{{{\text{ce}}0}} }}{{Z_{1} n_{10} }}\), \(\mu_{\text{se}} = \frac{{n_{\text{se0}} }}{{Z_{1} n_{10} }}\), \(\mu_{\text{H}} = \frac{{n_{{{\text{H}}0}} }}{{Z_{1} n_{10} }}\), \(\sigma_{\text{ce}} = \frac{{T_{\text{ce}} }}{{T_{1} }}\), \(\sigma_{\text{se}} = \frac{{T_{\text{se}} }}{{T_{1} }}\), and \(\sigma_{\text{H}} = \frac{{T_{\text{H}} }}{{T_{1} }}.\)

Derivation of KdVB equation

We use the reductive perturbation method to derive the KdVB equation from (5) to (9) by introducing the transformations

$$\xi = \varepsilon^{1/2} (x - \lambda t),\,\,\quad \tau = \varepsilon^{3/2} t,\quad \rho_{j} = \varepsilon^{1/2} \rho_{j0}$$

where ε is a smallness parameter and λ is the wave phase speed.

To apply the reductive perturbation technique, the various parameters are expanded as:

$$N_{1,2} = 1 + \varepsilon N_{1,2}^{(1)} + \varepsilon^{2} N_{1,2}^{(2)} + \cdots$$
(10)
$$V_{(1,2)} = \varepsilon V_{(1,2)}^{(1)} + \varepsilon^{2} V_{(1,2)}^{(2)} + \cdots$$
(11)
$$\phi = \varepsilon \phi^{(1)} + \varepsilon^{2} \phi^{(2)} + \cdots$$
(12)

We substitute (10)–(12) in (5)–(9) and equate the coefficients of different powers of ε. From the coefficients of order \(\varepsilon^{3/2}\), we get the first-order terms as:

$$N_{1}^{1} = \frac{{\phi^{1} }}{{\left( {\frac{3}{{Z_{1} }} - \lambda^{2} } \right)}}$$
(13)

and

$$N_{2}^{1} = \frac{{(Z_{2} /Z_{1} )m\phi^{1} }}{{\left( {\lambda^{2} - \frac{3m\beta }{{Z_{1} }}} \right)}}.$$
(14)

Expressions for \(V_{1}^{1}\) and \(V_{2}^{1}\) can be obtained by multiplying (13) and (14) by \(\lambda\).

In addition, the linear dispersion relation is:

$$\lambda^{2} = \frac{{S \pm \sqrt {S^{2} - 12mZ_{1}^{2} T\left[ {3m\beta T + mZ_{2} \left( {1 + \mu_{ce} + \mu_{se} - \mu_{H} } \right) + m\beta Z_{1} } \right]} }}{{2TZ_{1}^{2} }},$$
(15)

where

$$S = Z_{1}^{2} + mZ_{1} Z_{2} \left( {1 + \mu_{ce} + \mu_{se} - \mu_{H} } \right) + 3Z_{1} T(1 + m\beta ),$$
(16)

and

$$T = \frac{{\mu_{\text{ce}} \left( {\kappa_{\text{ce}} - 1/2} \right)}}{{\left( {\kappa_{\text{ce}} - 3/2} \right)\sigma_{\text{ce}} }} + \frac{{\mu_{\text{se}} \left( {\kappa_{\text{se}} - 1/2} \right)}}{{\left( {\kappa_{\text{se}} - 3/2} \right)\sigma_{\text{se}} }} + \frac{{\mu_{\text{H}} \left( {\kappa_{\text{H}} - 1/2} \right)}}{{\left( {\kappa_{\text{H}} - 3/2} \right)\sigma_{\text{H}} }}.$$
(17)

Equating the coefficients of \(\varepsilon^{5/2}\) in (5) and (6), we get

$$\frac{{\partial N_{i}^{1} }}{\partial \tau } - \lambda \frac{{\partial N_{i}^{2} }}{\partial \xi } + \frac{{\partial V_{i}^{2} }}{\partial \xi } + \frac{{\partial (N_{i}^{1} V_{i}^{1} )}}{\partial \xi } = 0\quad i = 1,\,2.$$
(18)

And equating the coefficient of order \(\varepsilon^{5/2}\) in (7) and (8) results in:

$$\frac{{\partial V_{1}^{1} }}{\partial \tau } - \lambda \frac{{\partial V_{1}^{2} }}{\partial \xi } + V_{1}^{1} \frac{{\partial V_{1}^{1} }}{\partial \xi } = \frac{{\partial \phi^{2} }}{\partial \xi } - \frac{{3N_{1}^{1} }}{{Z_{1} }}\frac{{\partial N_{1}^{1} }}{\partial \xi } - \frac{3}{{Z_{1} }}\frac{{\partial N_{1}^{2} }}{\partial \xi } + \rho_{10} \frac{{\partial^{2} V_{1}^{1} }}{{\partial \xi^{2} }},$$
(19)
$$\frac{{\partial V_{2}^{1} }}{\partial \tau } - \lambda \frac{{\partial V_{2}^{2} }}{\partial \xi } + V_{2}^{1} \frac{{\partial V_{2}^{1} }}{\partial \xi } = - \frac{{Z_{2} m}}{{Z_{1} }}\frac{{\partial \phi^{2} }}{\partial \xi } - \frac{{3m\beta N_{2}^{1} }}{{Z_{1} }}\frac{{\partial N_{2}^{1} }}{\partial \xi } - \frac{3m\beta }{{Z_{1} }}\frac{{\partial N_{2}^{2} }}{\partial \xi } + \rho_{20} \frac{{\partial^{2} V_{2}^{1} }}{{\partial \xi^{2} }}.$$
(20)

Finally, equating the coefficients of terms of order ε 2 from Poisson’s equation (9) gives:

$$\begin{aligned} \frac{{\partial^{2} \phi^{1} }}{{\partial \xi^{2} }} & = N_{1}^{2} - N_{2}^{2} \left( {1 + \mu_{\text{ce}} + \mu_{\text{se}} - \mu_{\text{H}} } \right) + T\phi^{2} + \frac{{\mu_{\text{ce}} }}{2}\frac{{\left( {\kappa_{\text{ce}}^{2} - 1/4} \right)}}{{\left( {\kappa_{\text{ce}} - 3/2} \right)^{2} \sigma_{\text{ce}}^{2} }}(\phi^{1} )^{2} \\ & \quad + \frac{{\mu_{\text{se}} }}{2}\frac{{\left( {\kappa_{\text{se}}^{2} - 1/4} \right)}}{{\left( {\kappa_{\text{se}} - 3/2} \right)^{2} \sigma_{\text{se}}^{2} }}(\phi^{1} )^{2} - \frac{{\mu_{\text{H}} }}{2}\frac{{\left( {\kappa_{\text{H}}^{2} - 1/4} \right)}}{{\left( {\kappa_{\text{H}} - 3/2} \right)^{2} \sigma_{\text{H}}^{2} }}(\phi^{1} )^{2} . \\ \end{aligned}$$
(21)

Substituting the values from (13) and (14) into (18) to (21) and eliminating the second-order terms, we obtain the KdVB equation as:

$$A\frac{{\partial \phi^{1} }}{\partial \tau } + \phi^{1} \frac{{\partial \phi^{1} }}{\partial \xi } + B\frac{{\partial^{3} \phi^{1} }}{{\partial \xi^{3} }} - C\frac{{\partial^{2} \phi^{1} }}{{\partial \xi^{2} }} = 0.$$
(22)

In (22), the coefficients A, B, and C are given by:

$$A = \frac{{ - 2Z_{1}^{2} \lambda \left( {Z_{1} \lambda^{2} - 3} \right)\left( {Z_{1} \lambda^{2} - 3m\beta } \right)\left[ {\left( {Z_{1} \lambda^{2} - 3m\beta } \right)^{2} + m\left( {\frac{{Z_{2} }}{{Z_{1} }}} \right)\left( {1 + \mu_{\text{ce}} + \mu_{\text{se}} - \mu_{\text{H}} } \right)\left( {Z_{1} \lambda^{2} - 3} \right)^{2} } \right]}}{D}$$
$$B = \frac{{ - \left( {Z_{1} \lambda^{2} - 3} \right)^{3} \left( {Z_{1} \lambda^{2} - 3m\beta } \right)^{3} }}{D}$$
$$C = \frac{{ - Z_{1}^{2} \lambda \left( {Z_{1} \lambda^{2} - 3} \right)\left( {Z_{1} \lambda^{2} - 3m\beta } \right)\left[ {\rho_{10} \left( {Z_{1} \lambda^{2} - 3m\beta } \right)^{2} - m\rho_{20} \left( {\frac{{Z_{2} }}{{Z_{1} }}} \right)\left( {1 + \mu_{\text{ce}} + \mu_{\text{se}} - \mu_{\text{H}} } \right)\left( {Z_{1} \lambda^{2} - 3} \right)^{2} } \right]}}{D}$$

where

$$L = \frac{{\mu_{\text{ce}} \left( {\kappa_{\text{ce}}^{2} - 1/4} \right)}}{{\left( {\kappa_{\text{ce}} - 3/2} \right)^{2} \sigma_{\text{ce}}^{2} }} + \frac{{\mu_{\text{se}} \left( {\kappa_{\text{se}}^{2} - 1/4} \right)}}{{\left( {\kappa_{\text{se}} - 3/2} \right)^{2} \sigma_{\text{se}}^{2} }} - \frac{{\mu_{\text{H}} \left( {\kappa_{\text{H}}^{2} - 1/4} \right)}}{{\left( {\kappa_{\text{H}} - 3/2} \right)^{2} \sigma_{\text{H}}^{ 2} }}$$

and

$$\begin{aligned} D & = L\left( {Z_{1} \lambda^{2} - 3} \right)^{3} \left( {Z_{1} \lambda^{2} - 3m\beta } \right)^{3} + 3Z_{1}^{2} \left( {1 + Z_{1} \lambda^{2} } \right)\left( {Z_{1} \lambda^{2} - 3m\beta } \right)^{3} - 3m^{2} Z_{2}^{2} \left( {1 + \mu_{\text{ce}} + \mu_{\text{se}} - \mu_{\text{H}} } \right)\left( {Z_{1} \lambda^{2} - 3} \right)\left( {Z_{1} \lambda^{2} + m\beta } \right). \\ \end{aligned}$$
(23)

Solution of KdVB equation

To find the solution of (22), we use the transformed coordinate \(\chi = \left( {\xi - V\tau } \right)\) of the comoving frame with speed V and use the boundary conditions:\(\phi^{1} \to 0\) and \(\frac{{\partial \phi^{1} }}{\partial \chi }\), \(\frac{{\partial^{2} \phi^{1} }}{{\partial \chi^{2} }}\), \(\frac{{\partial^{3} \phi^{1} }}{{\partial \chi^{3} }} \to 0\) as \(\chi \to \infty\) for a localized solution [42].

When the partial differential equation of a system is formed by the combined effect of dispersion and dissipation, a convenient method to solve it is “the tanh method” [43, 44]. Using the above transformation, (22) can be written as:

$$- AV\frac{{\partial \phi^{1} }}{\partial \chi } + \phi^{1} \frac{{\partial \phi^{1} }}{\partial \chi } + B\frac{{\partial^{3} \phi^{1} }}{{\partial \chi^{3} }} - C\frac{{\partial^{2} \phi^{1} }}{{\partial \chi^{2} }} = 0.$$
(24)

Again using the transformation \(\alpha = \tanh \chi\) and assuming a series solution of the form \(\phi^{1} (\alpha ) = \sum\nolimits_{i = 0}^{n} {a_{i} \alpha^{i} }\), we arrive at the solution of (22) as:

$$\phi^{1} = AV + 8Bk^{2} + \frac{{C^{2} }}{25B} - 12k^{2} \tanh^{2} [k(\xi - V\tau )] - \frac{12}{5}kC[1 - \tanh [k(\xi - V\tau )]].$$
(25)

The speed of comoving frame is related to the coefficients A, B, and C as \(V = \frac{{(100B^{2} k^{2} - C^{2} + 60kBC)}}{25AB}\) and \(k = \frac{ \pm C}{10B}\) which can be obtained using the boundary conditions.

Results

Solution (25) is applicable to any astrophysical plasma. However, in this study, we concentrate on parameters relevant to comet Halley: the observed value of the density of hydrogen ions was \(n_{\text{H}} = 4.95\text{cm}^{ - 3}\); their temperature was \(T_{\text{H}} = 8 \times 10^{4} K\). The temperature of the solar (or hot) electrons was \(T_{\text{se}} = 2 \times 10^{5} {\text{K}}\) [45]. The temperature of the second component of electrons, namely photoelectrons was set at \(T_{\text{ce}} = 2 \times 10^{4} {\text{K}}\). Negatively charged oxygen ions with an energy ~1 eV and densities ≤1 cm−3 was identified by Chaizy et al. [24]. We thus set the densities of positively charged oxygen ions at \(n_{20} = 0.5\text{cm}^{ - 3}\) and that of negatively charged oxygen ions at \(n_{10} = 0.05\text{cm}^{ - 3}\) [24, 45].

Figure 1 is a plot of the solution (25) of the KdVB Eq. (22), and shows the variation of the potential \(\varphi^{1}\) versus \(\chi\) as a function of the spectral indices; the parameters for the figure are \(n_{10} = 0.05\text{cm}^{ - 3}\), \(n_{20} = 0.5\text{cm}^{ - 3}\), \(n_{\text{H}} = 4.95\text{cm}^{ - 3}\), \(T_{\text{ce}} = 2 \times 10^{4} {\text{K}}\), \(T_{\text{se}} = 2 \times 10^{5} {\text{K}}\), \(T_{1} = T_{2} = 1.16 \times 10^{4} K\), and \(Z_{1} = Z_{2} = 1\). Curve (a) is for the spectral index \(\kappa\) = 3, curve (b) is for \(\kappa\) = 5, and curve (c) is for \(\kappa\) = 7. It is clear from the figure that the amplitude of shock wave increases with a decrease of kappa indices.

Fig. 1
figure 1

\(\phi^{1} \,{\text{vs}}\,\chi\) as a function of kappa indices

Full size image

Figure 2 is a plot of the potential \(\varphi^{1}\) versus \(\chi\) as a function of temperature ratio (\(\beta = \frac{{T_{2} }}{{T_{1} }}\)) of the pair ions. We keep the spectral indices at \(\kappa_{\text{ce}} = \kappa_{\text{se}} = \kappa_{\text{H}} = 3\); the other parameters are the same as in Fig. 1. Curve (a) is for β = 2, curve (b) is for β = 2.2, and curve (c) is for β = 2.4. We find that the strength of shock profile decreases with an increase of the temperature of positively charged oxygen ions.

Fig. 2
figure 2

\(\phi^{1} \,{\text{vs}}\,\chi\) as a function of temperature ratio \(\beta\)

Full size image

Figure 3 is again a plot of the potential \(\varphi^{1}\) versus \(\chi\) as a function of the negative oxygen ion densities n 10. The parameters used in this case are \(n_{20} = 0.5\text{cm}^{ - 3}\), \(n_{\text{H}} = 4.95\text{cm}^{ - 3}\), \(T_{\text{ce}} = 2 \times 10^{4} {\text{K}}\), \(T_{\text{se}} = 2 \times 10^{5} {\text{K}}\), \(T_{1} = T_{2} = 1.16 \times 10^{4} K\), \(Z_{1} = Z_{2} = 1\), and \(\kappa_{\text{ce}} = \kappa_{\text{se}} = \kappa_{\text{H}} = 5\). Curve (a) is for n 10 = 0.06 cm−3, curve (b) is for n 10 = 0.07 cm−3, and curve (c) is for n 10 = 0.08 cm−3. We find that the strength of shock wave decreases with an increase of the negatively charged oxygen ion densities.

Fig. 3
figure 3

\(\phi^{1} \,{\text{vs}}\,\chi\) as a function of O− density

Full size image

Figure 4 shows the variation of the amplitude of shock wave as a function of the positively charged oxygen ion densities. The parameters chosen are the same as in Fig. 3. Here, the negatively charged oxygen ion density is n 10 = 0.05 cm−3. Curve (a) is for n 20 = 0.3 cm−3, curve (b) is for n 20 = 0.5 cm−3, and curve (c) is for n 20 = 0.7 cm−3. We find that increasing positively charged oxygen ion density will increase the strength of the shock wave profile.

Fig. 4
figure 4

\(\phi^{1} \,{\text{vs}}\,\chi\) as a function of O+ density

Full size image

Figure 5 is a plot of the solution of the KdVB Eq. (22) and shows the variation of the potential \(\varphi^{1}\) versus \(\chi\) as a function of the normalized kinematic viscosity of positively charged oxygen ions (\(\rho_{2}\)); the parameters for the figure are \(n_{10} = 0.05\text{cm}^{ - 3}\), \(n_{20} = 0.5\text{cm}^{ - 3}\), \(n_{\text{H}} = 4.95\text{cm}^{ - 3}\), \(T_{\text{ce}} = 2 \times 10^{4} {\text{K}}\), \(T_{\text{se}} = 2 \times 10^{5} {\text{K}}\), \(T_{1} = T_{2} = 1.16 \times 10^{4} K\), and \(Z_{1} = Z_{2} = 1\). Curve (a) is for the kinematic viscosity \(\rho_{2}\) = 0.25, curve (b) is for \(\rho_{2}\) = 0.5, and curve (c) is for \(\rho_{2}\) = 0.75. It is clear from the figure that the strength of shock wave increases with an increase in kinematic viscosity of positively charged oxygen ions.

Fig. 5
figure 5

\(\phi^{1} \,{\text{vs}}\,\chi\) as a function of kinematic viscosity of O+ ions

Full size image

Figure 6 shows the variation of the potential \(\varphi^{1}\) versus \(\chi\) as a function of the normalized kinematic viscosity of negatively charged oxygen ions (\(\rho_{1}\)). The parameters used in this case are the same as in Fig. 5. Curve (a) is for the kinematic viscosity \(\rho_{1}\) = 0.1, curve (b) is for \(\rho_{1}\) = 0.3, and curve (c) is for \(\rho_{1}\) = 0.5. We find that the strength of shock wave decreases with an increase in the kinematic viscosity of negatively charged oxygen ions.

Fig. 6
figure 6

\(\phi^{1} \,{\text{vs}}\,\chi\) as a function of kinematic viscosity of O− ions

Full size image

In a study of electron acoustic solitary and shock waves in dissipative space plasmas with superthermal, hot electrons, Han et al. [46] found that the amplitude of the shock wave decreased with increasing spectral index \(\kappa\). They also found that the strength of the shock wave increased with increasing ratios of the cold to hot electron densities and electron kinematic viscosities. In another study on drift ion-acoustic shocks in an inhomogeneous 2-D quantum plasma, Masood et al. [47] found that the strength of the shock waves increased with increasing densities of ions.

The results of Fig. 1, which are in agreement with that of Han et al. [46], can be considered as extending their result to a multi-component plasma as we now have three components being described by kappa distributions. Similarly, the conclusions of Fig. 4 where the strength of the shock wave increases with increasing positively charged oxygen ion density is in agreement with the result of Masood et al. [47].

Finally, the other observation of Han et al. [46] that the strength of shock wave increased with increasing electron kinematic viscosity was interpreted as being due to increasing dissipation among the constituents. Increasing the density of positively charged oxygen ions would lead to an increase in the electron density; conversely, an increase in the negatively charged oxygen ion densities would lead to a decrease in the electron densities. The results of Figs. 5 and 6 are consistent with this observation of Han et al. [46].

More recently, Sabetkar and Dorranian [48] investigated the characteristics of dust acoustic solitary waves (DASWs) in a plasma composed of two populations of ions, electrons and negatively charged dust. The ions and electrons were modeled by three-dimensional non-extensive and kappa distributions, respectively. They derived the Zakharov–Kuznetsov (ZK) and modified Zakharov–Kuznetsov (mZK) equations to describe dust acoustic waves in the system. Both positive and negative polarity solutions were found to exist. They found that the amplitude of the ZK solitons decreased with increasing κ values, which is the conclusion from Fig. 1. In addition, the amplitude of these solitons increased with an increase in the temperature ratio of the cold to hot ions, similar to the conclusion from Fig. 2.

In yet another study, the same authors [49] studied the “fast” and “slow” dust acoustic waves in a plasma composed of Maxwellian electrons, kappa distributed positive ions, kappa-Schamel distributed negative ions and negatively charged dust particles. The 1 and 3D Schamel–Korteweg–deVries (S–KdV) equation that they derived admitted only compressive solitons. A point of convergence with their conclusion is the decrease in the amplitude of the soliton with the temperature ratio β p (=\(\frac{{T_{\text{p}} }}{{T_{\text{n}} }}\)). Further comparisons are not possible because of the different distributions used and also the different equations studied.

Conclusion

We have, in this paper, studied ion-acoustic shock waves in a five component plasma of positively and negatively charged oxygen ions, lighter hydrogen ions, and hot and cold electrons by deriving the KdVB equation. The impact of spectral indices kappa, temperature of positively charged oxygen ions, and density of oxygen ions (O+, O−) on the shock wave amplitude is studied. We find that in a cometary plasma with the above components, the nonlinear wave is in a transition state from shock to soliton [20]. A reduction in the shock wave amplitude is seen with increasing spectral indices and positively charged oxygen ion density. The strength of the shock profile also decreases with increasing positively charged oxygen ion temperatures and negatively charged oxygen ion densities. These results can be expected to contribute to an understanding of shocks in comets as we have two heavy ion components and heavy ions were surmised to affect shocks in cometary plasmas [25].