Application scope of the reductive perturbation method to derive the KdV equation and CKdV equation in dusty plasma

Abstract

The application scopes of two different reductive perturbation methods to derive the Korteweg–de Vries (KdV) equation and coupled KdV (CKdV) equation in two-temperature-ion dusty plasma are given by using the particle-in-cell (PIC) numerical method in the present paper. It suggests that the reductive perturbation method (RPM) is valid if the amplitude of the CKdV solitary wave is small enough. However, for the KdV solitary wave, RPM is valid not only if the amplitude of the KdV solitary wave is small enough, but also if the nonlinear coefficient of the KdV equation is not tending to zero.

1. Introduction

Dusty plasmas usually contain highly negatively charged dust particles, electrons, ions as well as neutrals. Sometimes, there are different kinds of dust particles and different kinds of ions in a dusty plasma. Tremendous progress in dusty plasmas has been made to date (Kalman, Rosenberg & DeWitt 2000; Shukla 2001; Xie & He 2001; Wang & Wu 2002; Fortov et al. 2005; Shukla & Eliasson 2009; Teng et al. 2009; Couëdel et al. 2010; Ghorui, Chatterjee & Wong 2013; Ghorui et al. 2014; Merlino 2014; Han et al. 2015; Seadawy 2015). The remarkable feature of a dusty plasma is that there are low-frequency waves (Rao, Shukla & Yu 1990; Shukla 1992; Barkan, Merlino & D'Angelo 1995; Pieper & Goree 1996; Praburam & Goree 1996; Merlino et al. 1998; Bandyopadhyay et al. 2008; Dubinov & Sazonkin 2013; Kumar Tiwari & Sen 2016; Verheest & Hereman 2019).

Many researchers focused on the nonlinear solitary waves in a dusty plasma (El-Labany et al. 2004; Bandyopadhyay et al. 2008; Emamuddin, Yasmin & Mamun 2013). A solitary wave was first discovered in 1834. Later, a Korteweg–de Vries (KdV) equation was proposed in 1895 (Korteweg & De Vries 1895). However, until Kruskal and Zabusky discovered a solitary wave solution from a KdV equation in 1965 (Zabusky & Kruskal 1965; Olivier, Verheest & Maharaj 2016), different methods such as the perturbation method were used to study the nonlinear solitary waves. After that, solitary waves were found nearly in all branches of physics, especially in plasma physics (Washimi & Taniuti 1966), fluid dynamics (Ono 1991; Duan, Wang & Wei 1996), nonlinear lattice (Zabusky & Kruskal 1965; Wadati 1990), Bose–Einstein condensates (Huang, Velarde & Makarov 2001), etc.

Recently, the application scope of the perturbation method (Lee 2012; Olivier et al. 2016) is studied by using the one-dimensional particle-in-cell (PIC) code (Qi et al. 2014; Zhang et al. 2014, 2015; Gao et al. 2017; Zhang et al. 2017; Wang et al. 2018). A critical point of the reductive perturbation parameter $\epsilon ^*$ is found. If the perturbation parameter $\epsilon <\epsilon ^*$, the reductive perturbation method is valid, while if $\epsilon >\epsilon ^*$, the reductive perturbation method is no longer reasonable (Qi et al. 2014; Zhang et al. 2014).

As is well known, many nonlinear equations can be reduced to a more simple nonlinear equation, such as the KdV equation, by using the reductive perturbation method (RPM) under small amplitude and long wavelength approximation. However, in some cases, the reductive perturbation method is invalid to derive a KdV equation. In the present paper, we will show that the reductive perturbation method is invalid when the nonlinear coefficient of the KdV tends to zero. Due to this reason, we choose another different reductive perturbation method to derive a coupled KdV (CKdV) equation. Furthermore, we find that a larger amplitude of the solitary wave results in a greater deviation between the numerical results and the theoretical results, i.e. the RPM is valid if the amplitude of the CKdV solitary wave is small enough. However, for the KdV solitary wave, the RPM is valid not only if the amplitude of the KdV solitary wave is small enough, but also if the parameters have to be limited.

2. Theoretical model

We now focus on the dust acoustic waves in two-temperature-ion dusty plasma containing highly negatively charged dust particles, electrons and two different kinds of ions (high-temperature ions and low-temperature ions). Charge neutrality at equilibrium reads ${n_{il0}}+{n_{ih0}}={Z_{d0}}{n_{d0}}+{n_{e0}}$, where ${n_{il0}}$, ${n_{ih0}}$, ${n_{e0}}$ and ${n_{d0}}$ are the number densities of unperturbed low-temperature ion, high-temperature ion, electron and the dust particles, respectively. Here, ${Z_{d0}}$ is the unperturbed number of charges residing on the dust grain measured in the units of electron charge.

For convenience and generality, we assumed that the dusty plasma is unmagnetized and collisionless. The waves propagate in the $x$ direction. In this case, the dust fluid satisfies the following equations:

(2.1)\begin{gather} \frac{{\partial {n_d}}}{{\partial t}} + \frac{\partial ({n_d}{u_d})}{{\partial x}} = 0, \end{gather}

(2.2)\begin{gather} \frac{{\partial {u_d}}}{{\partial t}} + {u_d}\frac{{\partial {u_d}}}{{\partial x}} + \frac{1}{{{n_d}{m_d}}}\frac{{\partial {p_d}}}{{\partial x}} = \frac{{{Z_d}e}}{{{m_d}}}\frac{{\partial\phi }}{{\partial x}}, \end{gather}

(2.3)\begin{gather} \frac{{{\partial ^2}\phi }}{{\partial {x^2}}} = \frac{e}{{{\epsilon_0}}}\left({{Z_d}{n_d} + {n_e} - {n_{il}} - {n_{ih}}}\right), \end{gather}

where ${n_d}$, ${m_d}$, $u_d$ and ${p_d}$ refer to the number density, the mass, the velocity and the pressure of dust grain fluid, respectively. Here, ${Z_{d0}}$ is the number of charge measured in units of electron charge $e$ when a dusty plasma is in the equilibrium state. We assume ${Z_d} = {Z_{d0}}$. Additionally, $\phi$ is the electrostatic potential. We also assume that the electron number density ($n_e$), the low temperature ion number density ($n_{il}$) and the high temperature ion number density ($n_{ih}$) satisfy the Boltzmann distributions, i.e. ${n_e} = {n_{e0}}\exp ( {{{e\phi } /{{k}{T_e}}}} )$, ${n_{il}}= {n_{il0}}\exp ( { - {{e\phi } / {k{T_{il}}}}} )$ and ${n_{ih}}= {n_{ih0}}\exp ( { - {{e\phi }/ {k{T_{ih}}}}} )$, where ${T_e}$, ${T_{il}}$, ${T_{ih}}$ and $k$ refer to the temperatures of the electrons, the low temperature ions, the high temperature ions and the Boltzmann constant, respectively.

Before normalization, we define an effective temperature as follows:

(2.4)\begin{equation} \frac{1}{{{T_{{\rm eff}}}}} = \frac{1}{{{Z_{d0}}{n_{d0}}}}\left( {\frac{{{n_{e0}}}}{{{T_e}}} + \frac{{{n_{il0}}}}{{{T_{il}}}} + \frac{{{n_{ih0}}}}{{{T_{ih}}}}} \right). \end{equation}

All physical quantities are normalized as follows: the dust grain number density ${n_{d}}$ is normalized by ${n_{d0}}$, the electron number density and ion number density are normalized by ${{Z_{d0}}{n_{d0}}}$, the pressure of dust grain fluid ${p_d}$ is normalized by ${{Z_{d0}}{n_{d0}}{T_d}}$ and the charge of the dust grain is normalized by ${Z_{d0}}$. The space coordinates $x$, time $t$, velocity $u_d$ and electrostatic potential $\phi$ are normalized by the Debye length ${\lambda _D} = {({{k{\epsilon _0}{T_{{\rm eff}}}}/{{n_{d0}}{z_{d0}}{e^2}}})^{1/2}}$, the inverse of dust plasma frequency ${\omega ^{ - 1}} = {({{{\epsilon _0}{m_d}}/{{n_{d0}}{Z_{d0}}^2{e^2}}})^{1/2}}$, the dust-acoustic speed ${c_d} = {({{k{z_{d0}}{T_{{\rm eff}}}}/{{m_d}}})^{1/2}}$ and ${{k{T_{{\rm eff}}}}/e}$, respectively. Here, ${T_d}$ is the temperature of the dust grain fluid. We assumed that the equation of state of the dust grain fluid is ${p_d} = {n_d}k{T_d}$ and the process is approximately adiabatic, i.e. ${p_d} = cn_d^\gamma$. Then (2.1)–(2.3) become

(2.5)\begin{gather} \frac{{\partial {n_d}}}{{\partial t}} + \frac{\partial ({n_d}{u_d}) }{{\partial x}} = 0, \end{gather}

(2.6)\begin{gather} \frac{{\partial {u_d}}}{{\partial t}} + {u_d}\frac{{\partial {u_d}}}{{\partial x}} + \frac{\gamma'}{n_d}\frac{{\partial {n_d}}}{{\partial x}} = \frac{{\partial \phi }}{{\partial x}}, \end{gather}

(2.7)\begin{gather} \frac{{{\partial ^2}\phi }}{{\partial {x^2}}} ={n_d} + {n_e} - {n_{il}} - {n_{ih}}, \end{gather}

where $\gamma '= {{\gamma {T_d}} / {{Z_{d0}}{T_{{\rm eff}}}}}$. Here, ${n_e} = v{e^{s{\beta _1}\phi }}$, ${n_{il}} = {\mu _l}{e^{ - s\phi }}$ and ${n_{ih}} = {\mu _h}{e^{ - s{\beta _2}\phi }}$ refer to the dimensionless number densities of electrons, lower temperature ions and higher temperature ions, respectively, with $\nu = {{{n_{e0}}} / {( {{Z_{d0}}{n_{d0}}} )}}$, ${\mu _l} = {{{n_{il0}}} / {( {{Z_{d0}}{n_{d0}}} )}}$, ${\mu _h} = {{{n_{ih0}}} / {( {{Z_{d0}}{n_{d0}}} )}}$, ${\beta _1} = {{{T_{il}}} / {{T_e}}}$, ${\beta _2} = {{{T_{il}}} / {{T_{ih}}}}$ and $s = 1/( {v{\beta _1} + {\mu _1} + {\mu _h}{\beta _2}} )$.

3. Nonlinear waves

3.1. Derivation of the KdV equation

We now use the reductive perturbation method to study the dust acoustic solitary waves in the limited case where the wave amplitude is small but finite, while the wavelength is long enough.

We introduce new coordinates of $\xi = \epsilon ( {x - {v_0}t} )$ and $\tau = {\epsilon ^3}t$, where $\epsilon$ is a small parameter characterizing the order of the wavenumber and ${v_0}$ is the velocity of the dust acoustic solitary waves. The physical quantities in the system are expanded as follows: ${n_d} = 1 + {\epsilon ^2}{n_{d1}} + {\epsilon ^4}{n_{d2}} + \cdots$, ${u_d} = {\epsilon ^2}{u_{d1}} + {\epsilon ^4}{u_{d2}} + \cdots$ and $\phi = {\epsilon ^2}{\phi _1} + {\epsilon ^4}{\phi _2} + \cdots$. Substituting these expansions into normalized equations (2.5)–(2.7) and collecting the terms in different powers of $\epsilon$, we obtain the following equations: ${n_{d1}} = - {\phi _1}$, ${u_{d1}} = - {v_0}{\phi _1}$, ${v_0}^2 = 1 + \gamma '$ and the KdV equation

(3.1)\begin{equation} \frac{\partial {\phi _1} }{{\partial \tau }} + A{\phi _1}\frac{\partial {\phi _1} }{{\partial \xi }} + B\frac{{{\partial ^3 {\phi _1}}}}{{\partial {\xi ^3}}} = 0, \end{equation}

where $A = - ({1}/{{2{v_0}}})[ {3 + 2\gamma ' + ( {v{\beta _1}^2 - {\mu _l} - {\mu _h}{\beta _2}^2} ){s^2}} ]$ and $B = {1}/{{2{v_0}}}$.

One of the important solutions of the KdV equation (3.1) is a single solitary wave solution as follows:

(3.2)\begin{equation} {\phi _1} = {\phi _m}{\rm sech}^2\frac{{\xi - {u_0}\tau }}{W}, \end{equation}

where ${\phi _m} = {{3{u_0}}}/{A}$ and $W = \sqrt {{4B / {{u_0}}}}$. Rewriting (3.2) in experimental coordinates by $\xi = \epsilon ( {x - {v_0}t} )$ and $\tau = {\epsilon ^3}t$, we have

(3.3)\begin{equation} {\phi _1} = \frac{{3{u_0}}}{A}{\rm sech}^2\frac{{x - \left( {{v_0} + {u_0}{\epsilon ^2}} \right)t}}{{\sqrt {\dfrac{{4B}}{{{u_0}{\epsilon ^2}}}} }}. \end{equation}

We also know that ${n_{d1}} = - {\phi _1}$, ${u_{d1}} = - {v_0}{\phi _1}$ and the physical quantities in the system are expanded as follows: ${n_d} = 1 + {\epsilon ^2}{n_{d1}} + {\epsilon ^4}{n_{d2}} + \cdots$, ${u_d} = {\epsilon ^2}{u_{d1}} + {\epsilon ^4}{u_{d2}} + \cdots$ and $\phi = {\epsilon ^2}{\phi _1} + {\epsilon ^4}{\phi _2} + \cdots$. Then we consider the lower order terms, and we can obtain the following equations:

(3.4)\begin{gather} \phi = \frac{{3{u_0}{\epsilon ^2}}}{A}{\rm sech}^2\frac{{x - \left( {{v_0} + {u_0}{\epsilon ^2}} \right)t}}{{\sqrt {\dfrac{{4B}}{{{u_0}{\epsilon ^2}}}} }}, \end{gather}

(3.5)\begin{gather} {n_d} = 1 - \frac{{3{u_0}{\epsilon ^2}}}{A}{\rm sech}^2\frac{{x - \left( {{v_0} + {u_0}{\epsilon ^2}} \right)t}}{{\sqrt {\dfrac{{4B}}{{{u_0}{\epsilon ^2}}}} }}, \end{gather}

(3.6)\begin{gather} {u_d} ={-} {v_0}\frac{{3{u_0}{\epsilon ^2}}}{A}{\rm sech}^2\frac{{x - \left( {{v_0} + {u_0}{\epsilon ^2}} \right)t}}{{\sqrt {\dfrac{{4B}}{{{u_0}{\epsilon ^2}}}} }}, \end{gather}

where ${\phi _m} = {{3{u_0}{\epsilon ^2}} / A}$, $W = \sqrt {{4B / {{u_0}{\epsilon ^2}}}}$ are the amplitude and the width of the KdV solitary wave, respectively.

3.2. Derivation of the CKdV equation

It seems that the amplitude of the solitary wave solution of the KdV equation can be infinity when $A=0$. There is no physical meaning in this case. Therefore, the reductive perturbation is invalid in this case. Due to this reason, figure 1 shows the dependence of $A$ on ${\mu _l}$. Notice from figure 1 that $A$ can be negative, zero and positive.

Figure 1. Dependence of $A$ on ${\mu _l}$, where the parameters are ${T_e} = 5$ eV, ${T_{il}} = 0.3$ eV, ${T_{ih}} = 4$ eV, ${T_d} = 200$ eV, ${Z_{d0}} = 1000$, ${n_{d0}} = 1 \times {10^{12}}$ m$^{ - 3}$ and ${n_{ih0}} = 1 \times {10^{15}}$ m$^{ - 3}$.

Because $A$ can be zero, the KdV solitary wave solution may be invalid. Due to this reason, we now try to use the other method to find the nonlinear wave solution when $A=0$ or $A$ tends to zero. For this purpose, we use the following expansions (Duan & Shi 2003): ${n_d} = 1 + \epsilon {n_{d1}} + {\epsilon ^2}{n_{d2}} + {\epsilon ^3}{n_{d3}} + \cdots$, ${u_d} = \epsilon {u_{d1}} + {\epsilon ^2}{u_{d2}} + {\epsilon ^3}{u_{d3}} + \cdots$, $\phi = \epsilon {\phi _1} + {\epsilon ^2}{\phi _2} + {\epsilon ^3}{\phi _3} + \cdots$. Substituting these expansions into (2.5)–(2.7) and collecting the terms in different powers of $\epsilon$, we obtain the following equations at the different orders: ${n_{d1}} = - {\phi _1}$, ${u_{d1}} = - {v_0}{\phi _1}$, ${v_0}^2 = 1 + \gamma '$ and a modified KdV (MKdV) equation as

(3.7)\begin{equation} \frac{\partial {\phi _1} }{{\partial \tau }} + C{\phi _1}^2\frac{\partial {\phi _1} }{{\partial \xi }} + B\frac{{{\partial ^3 {\phi _1}}}}{{\partial {\xi ^3}}} = 0, \end{equation}

where $B = {1}/{{2{v_0}}}$, $C = ({1}/{{2{v_0}}})({15 + 26\gamma ' + 12\gamma {'^2}} )$.

If $A$ is small enough, but not zero, we use the same perturbation method as that to derive the MKdV equation (Duan & Shi 2003). We then obtain a CKdV equation as follows:

(3.8)\begin{equation} \frac{\partial {\phi _1} }{{\partial \tau }} + A{\phi _1}\frac{\partial {\phi _1}}{{\partial \xi }} + C{\phi _1}^2\frac{\partial {\phi _1} }{{\partial \xi }} + B\frac{{{\partial ^3 {\phi _1} }}}{{\partial {\xi ^3}}}= 0, \end{equation}

where $A = - ({1}/{{2{v_0}}})[ {3 + 2\gamma ' + ( {v{\beta _1}^2 - {\mu _l} - {\mu _h}{\beta _2}^2} ){s^2}} ]$, $B = {1}/{{2{v_0}}}$, $C = ({1}/{{2{v_0}}})( 15 + 26\gamma ' + 12\gamma {'^2})$ and $B > 0, C > 0$.

A single solitary wave solution of the CKdV equation is

(3.9)\begin{equation} {\phi _1} = \frac{1}{{{a_1} + {a_2}\cosh [{a_3}\left( {\xi - {u_o}\tau } \right)]}}, \end{equation}

where ${a_1} = {A / {6{u_0}}}$, ${a_2} = {{\sqrt {{A^2} + 6C{u_0}} } / {6{u_0}}}$ and ${a_3} = \sqrt {{{{u_0}} / B}}$.

Rewriting (3.9) in experimental coordinates, we have

(3.10)\begin{equation} {\phi _1} = \frac{1}{{\dfrac{A}{{6{u_0}}} + \dfrac{{\sqrt {{A^2} + 6C{u_0}} }}{{6{u_0}}}\cosh \left[ {\sqrt {\dfrac{{{u_0}{\epsilon ^2}}}{B}} \left( {x - \left( {{v_0} + {u_o}{\epsilon ^2}} \right)t} \right)} \right]}}. \end{equation}

We also know that ${n_{d1}} = - {\phi _1}$, ${u_{d1}} = - {v_0}{\phi _1}$ and the physical quantities in the system are expanded as follows: ${n_d} = 1 + \epsilon {n_{d1}} + {\epsilon ^2}{n_{d2}} + {\epsilon ^3}{n_{d3}} + \cdots$, ${u_d} = \epsilon {u_{d1}} + {\epsilon ^2}{u_{d2}} + {\epsilon ^3}{u_{d3}} + \cdots$, $\phi = \epsilon {\phi _1} + {\epsilon ^2}{\phi _2} + {\epsilon ^3}{\phi _3} + \cdots$. Then we consider the lower order terms, and we can obtain the following equations:

(3.11)\begin{gather} \phi = \frac{\epsilon }{{\dfrac{A}{{6{u_0}}} + \dfrac{{\sqrt {{A^2} + 6C{u_0}} }}{{6{u_0}}}\cosh \left[ {\sqrt {\dfrac{{{u_0}{\epsilon ^2}}}{B}} \left( {x - \left( {{v_0} + {u_o}{\epsilon ^2}} \right)t} \right)} \right]}}, \end{gather}

(3.12)\begin{gather} {n_d} = 1 - \frac{\epsilon }{{\dfrac{A}{{6{u_0}}} + \dfrac{{\sqrt {{A^2} + 6C{u_0}} }}{{6{u_0}}}\cosh \left[ {\sqrt {\dfrac{{{u_0}{\epsilon ^2}}}{B}} \left( {x - \left( {{v_0} + {u_o}{\epsilon ^2}} \right)t} \right)} \right]}}, \end{gather}

(3.13)\begin{gather} {u_d} ={-} \frac{{{v_0}\epsilon }}{{\dfrac{A}{{6{u_0}}} + \dfrac{{\sqrt {{A^2} + 6C{u_0}} }}{{6{u_0}}}\cosh \left[ {\sqrt {\dfrac{{{u_0}{\epsilon ^2}}}{B}} \left( {x - \left( {{v_0} + {u_o}{\epsilon ^2}} \right)t} \right)} \right]}}, \end{gather}

where ${\phi _m}= {\epsilon }/({{{A}/{{6{u_0}}} + ({{\sqrt {{A^2} + 6C{u_0}} }})/{{6{u_0}}}}})$ and $W = \sqrt {{B / {{u_0}{\epsilon ^2}}}}$ are the amplitude and the width of the CKdV solitary wave, respectively.

Notice that the CKdV equation (3.8) becomes the MKdV equation (3.7) in the limited case $A= 0$. Therefore, we mainly focus on the CKdV equation and its solutions in the following sections. Then we will give the numerical results of the KdV equation and CKdV equation by using the PIC numerical method in the following sections.

4. Numerical results

In this section, we will use the PIC numerical method to find the application scopes of the perturbation methods. During simulation, the dust grains are represented by a limited ensemble of super-particles (SPs), and both electrons and ions are treated as Boltzmann distributed fluids. The weight factor of SPs is $S$ which stands for the number of real particles. Initially, SPs are uniformly distributed in the space, its initial weight parameters $S$ and the velocities of each SP are obtained from the initial conditions. Therefore, the equation of motion of the system is Newton's equation as follows (Qi et al. 2014; Zhang et al. 2015, 2017):

(4.1)\begin{gather} {m_{sp}}\frac{{{\rm{d}}{v_{sp}}}}{{{\rm{d}}t}} = {q_{sp}}E, \end{gather}

(4.2)\begin{gather} \frac{{{\rm{d}}{x_{sp}}}}{{{\rm{d}}t}} = {v_{sp}}, \end{gather}

where ${m_{sp}}$, ${v_{sp}}$, ${q_{sp}}$ and ${x_{sp}}$ are the mass, velocity, charge and position of the SPs, respectively. In the PIC simulation, we divided the simulation region into several grid cells. The dust particles constantly exchange information with the background grid, when the dust particles move along their trajectories. At each time step, the positions and the velocities of SPs are weighted to all the grids, then we can calculate the charge density $\rho _{g}$ (or electric current density $J_{g}$ ). Once $\rho _{g}$ is obtained, the Maxwell's equations (electromagnetic model) or Poisson–Boltzmann equation (electrostatic model) will be solved numerically to derive the value of $E$ at each grid. In the electrostatic model, $B_{g} = 0$. Then the field imposed on each SP can be worked out and the electric field will drive each SP according to (4.1) and (4.2), which can be solved numerically by using the leap-frog algorithm. At last, the new positions and velocities are obtained, and the procedure will repeat until the simulation is completed (Qi et al. 2014; Zhang et al. 2015). A summary of a computational cycle of the PIC method is shown in figure 2.

Figure 2. Summary of a computational cycle of the PIC method.

4.1. Numerical results of the KdV equation

First we simulate the solitary wave solutions of the KdV equation. The initial conditions are given by (3.5) and (3.6) at $t=0$ as follows:

(4.3)\begin{gather} {n_d}|_{t=0} = 1 - \frac{{3{u_0}{\epsilon ^2}}}{A}{\rm sech}^2\left[ {\sqrt {\frac{{{u_0}{\epsilon ^2}}}{{4B}}} \left( {x - {x_0}} \right)} \right], \end{gather}

(4.4)\begin{gather} {u_d}|_{t=0} ={-} {v_0}\frac{{3{u_0}{\epsilon ^2}}}{A}{\rm sech}^2\left[ {\sqrt {\frac{{{u_0}{\epsilon ^2}}}{{4B}}} \left( {x - {x_0}} \right)} \right]. \end{gather}

The periodic boundary conditions are chosen. The simulation parameters are as follows: the spatial step is $\Delta x = 0.2$, the time step is $\Delta t = 0.004$, the number of grid cells is ${N_x} = 40\,000$, the number of super particles contained in each cell is $100$ and the total length of the $x$-axis is ${L_x} = \Delta x{N_x}$, ${x_0} = {{{L_x}} / 4}$.

Figure 3 shows the numerical results of the evolutions of the KdV solitary waves at different time $t$, where $A=0.37$. Notice that the solitary wave is a compressional one since $A>0$. The comparisons between PIC simulation results and the analytical ones at different times for $A=0.37$ are given in figure 4, which show good agreement between the two. Moreover, the numerical results are also undertaken in figure 5 for the case $A<0$, in which the rarefactive solitary wave is obtained. Good agreements between PIC simulation results and the analytical ones are also observed in figure 6 for the rarefactive solitary wave. It indicates that our analytical results are valid. Namely, the perturbation method in this case is valid.

Figure 3. PIC simulation results of the evolution of the KdV solitary waves at different times $t=19.9998/\omega ^{-1}, 1599.98/\omega ^{-1}, 3199.97/\omega ^{-1}$, where $A=0.37$ and the other parameters are $\epsilon = 0.05$, $\gamma = 3$, ${u_0} = 1$, ${x_0} = 2000$, ${T_e} = 5$ eV, ${T_{il}} = 0.3$ eV, ${T_{ih}} = 5$ eV, ${T_d} = 200$ eV, ${Z_{d0}} = 1000$, ${n_{d0}} = 1.0 \times {10^{12}}$ m$^{ - 3}$, ${n_{il0}} = 1.0 \times {10^{14}}$ m$^{ - 3}$ and ${n_{ih0}} = 1.0 \times {10^{15}}$ m$^{ - 3}$.

Figure 4. Comparisons between PIC simulation results and the analytical ones at different times $t=19.9998/\omega ^{-1}, 1599.98/\omega ^{-1}, 3199.97/\omega ^{-1}$, where $A=0.37$ and the other parameters are $\epsilon = 0.05$, $\gamma = 3$, ${u_0} = 1$, ${x_0} = 2000$, ${T_e} = 5$ eV, ${T_{il}} = 0.3$ eV, ${T_{ih}} = 5$ eV, ${T_d} = 200$ eV, ${Z_{d0}} = 1000$, ${n_{d0}} = 1.0 \times {10^{12}}$ m$^{ - 3}$, ${n_{il0}} = 1.0 \times {10^{14}}$ m$^{ - 3}$ and ${n_{ih0}} = 1.0 \times {10^{15}}$ m$^{ - 3}$.

Figure 5. PIC simulation results of the evolution of the KdV solitary waves at different times $t=19.9986/\omega ^{-1}, 999.93/\omega ^{-1}, 1999.86/\omega ^{-1}$, where $A=-1.23$ and the other parameters are $\epsilon = 0.1$, $\gamma = 3$, ${u_0} = 1$, ${x_0} = 1000$, ${T_e} = 3$ eV, ${T_{il}} = 0.5$ eV, ${T_{ih}} = 3$ eV, ${T_d} = 200$ eV, ${Z_{d0}} = 1000$, ${n_{d0}} = 1.0 \times {10^{12}}$ m$^{ - 3}$, ${n_{il0}} = 1.0 \times {10^{15}}$ m$^{ - 3}$ and ${n_{ih0}} = 1.1 \times {10^{15}}$ m$^{ - 3}$.

Figure 6. Comparisons between PIC simulation results and the analytical ones at different times $t=19.9986/\omega ^{-1}, 999.93/\omega ^{-1}, 1999.86/\omega ^{-1}$, where $A=-1.23$ and the other parameters are $\epsilon = 0.1$, $\gamma = 3$, ${u_0} = 1$, ${x_0} = 1000$, ${T_e} = 3$ eV, ${T_{il}} = 0.5$ eV, ${T_{ih}} = 3$ eV, ${T_d} = 200$ eV, ${Z_{d0}} = 1000$, ${n_{d0}} = 1.0 \times {10^{12}}$ m$^{ - 3}$, ${n_{il0}} = 1.0 \times {10^{15}}$ m$^{ - 3}$ and ${n_{ih0}} = 1.1 \times {10^{15}}$ m$^{ - 3}$.

We now focus on the question whether the KdV solitary wave expressed by (3.1) exists when $A=0$ or $A$ tends to zero. Figure 7(a–d) shows the numerical results of the evolutions of the waves at different time $t$ with the initial conditions of (4.3) and (4.4), where $A=0.001$, $A=-0.001$, $A=0.01$ and $A=-0.01$, respectively. It seems that the KdV solitary waves do not exist when $A$ is close to zero. More numerical results suggest that the KdV solitary waves exist when $|A|>0.06$, but do not exist if $|A|<0.06$.

Figure 7. PIC simulation results of the evolution of the nonlinear waves at different times, where (a) $A=0.001$; (b) $A=0.01$; (c) $A=-0.001$; (d) $A=-0.01$. The other parameters are $\epsilon = 0.01$, $\gamma = 3$, ${u_0} = 1$, ${x_0} = 2000$, ${T_e} = 5$ eV, ${T_{il}} = 0.3$ eV, ${T_{ih}} = 4$ eV, ${T_d} = 200$ eV, ${Z_{d0}} = 1000$, ${n_{d0}} = 1.0 \times {10^{12}}$ m$^{ - 3}$, ${n_{ih0}} = 1.0 \times {10^{15}}$ m$^{ - 3}$, and ${n_{il0}} = 1.277 \times {10^{14}}$ m$^{ - 3}$, $1.254 \times {10^{14}}$ m$^{ - 3}$, $1.283 \times {10^{14}}$ m$^{ - 3}$ and $1.31 \times {10^{14}}$ m$^{ - 3}$ in panels (a,b,c,d), respectively.

4.2. Numerical results of CKdV equation

If $A=0$ or $A$ tends to zero, the solution of the KdV equation is no longer valid, so we use the solution of the CKdV equation. Now we consider whether the CKdV solitary wave propagates when $A=0$ or $A$ tends to zero. The initial conditions are given by (3.12) and (3.13) at $t = 0$ as follows:

(4.5)\begin{gather} {n_d}|_{t=0} = 1 - \frac{\epsilon }{{\dfrac{A}{{6{u_0}}} + \dfrac{{\sqrt {{A^2} + 6C{u_0}} }}{{6{u_0}}}\cosh \left[ {\sqrt {\dfrac{{{u_0}{\epsilon ^2}}}{B}} \left( {x - {x_0}} \right)} \right]}}, \end{gather}

(4.6)\begin{gather} {u_d}|_{t=0} ={-} \frac{{{v_0}\epsilon }}{{\dfrac{A}{{6{u_0}}} + \dfrac{{\sqrt {{A^2} + 6C{u_0}} }}{{6{u_0}}}\cosh \left[ {\sqrt {\dfrac{{{u_0}{\epsilon ^2}}}{B}} \left( {x - {x_0}} \right)} \right]}}. \end{gather}

The periodic boundary conditions are chosen. The simulation parameters are as follows: the spatial step is $\Delta x = 0.5$, the time step is $\Delta t = 0.01$, the number of grid cells is ${N_x} = 20000$, the number of super particles contained in each cell is $100$, the total length of the $x$-axis is ${L_x} = \Delta x{N_x}$, $\gamma = 3$, ${x_0} = {{{L_x}} / 5}$, ${u_0} = 1$ and $\epsilon =0.01$.

Figure 8(a–d) shows the numerical results of the evolutions of the CKdV solitary waves at different time $t$, where $A=0.001$, $A=-0.001$, $A=0.01$ and $A=-0.01$, respectively. It shows that the CKdV solitary wave can propagate stably when $A=0$ or $A$ tends to zero.

Figure 8. PIC simulation results of the evolution of the CKdV solitary waves at different times, where (a) $A=0.001$; (b) $A=0.01$; (c) $A=-0.001$; (d) $A=-0.01$. The other parameters are $\epsilon = 0.01$, $\gamma = 3$, ${u_0} = 1$, ${x_0} = 2000$, ${T_e} = 5$ eV, ${T_{il}} = 0.3$ eV, ${T_{ih}} = 4$ eV, ${T_d} = 200$ eV, ${Z_{d0}} = 1000$, ${n_{d0}} = 1.0 \times {10^{12}}$ m$^{ - 3}$, ${n_{ih0}} = 1.0 \times {10^{15}}$ m$^{ - 3}$, and ${n_{il0}} = 1.277 \times {10^{14}}$ m$^{ - 3}$, $1.254 \times {10^{14}}$ m$^{ - 3}$, $1.283 \times {10^{14}}$ m$^{ - 3}$ and $1.31 \times {10^{14}}$ m$^{ - 3}$ in panels (a,b,c,d), respectively.

4.3. The application scope of the CKdV equation

To further understand the CKdV solitary wave propagation for different amplitudes, we change the value of $\epsilon$, and keep the other parameters constant (Lin & Duan 2005): $\gamma = 3$, ${u_0} = 1$, ${x_0} = 5000$, ${T_e} = 5$ eV, ${T_{il}} = 0.3$ eV, ${T_{ih}} = 4$ eV, ${T_d} = 200$ eV, ${Z_{d0}} = 1000$, ${n_{d0}} = 1.0 \times {10^{12}}$ m$^{ - 3}$, ${n_{il0}} = 1.28 \times {10^{14}}$ m$^{ - 3}$ and ${n_{ih0}} = 1.0 \times {10^{15}}$ m$^{ - 3}$. We will see how CKdV solitary waves vary concerning the values of the parameter $\epsilon$.

Figure 9(a) shows the numerical results of the evolutions of the CKdV solitary waves at different time $t$ when $\epsilon =0.01$, i.e. the wave amplitude is ${\phi _m}=0.009$. Figure 9(b) shows the comparisons between PIC simulation results and the analytical ones at different times.

Figure 9. (a) PIC simulation results of the evolution of the CKdV solitary waves at different times, when $\epsilon =0.01$; (b) comparisons between PIC simulation results and the analytical ones at different times, when $\epsilon =0.01$.

Notice that the amplitude of the solitary wave remains unchanged during its propagation. Good agreements between PIC simulation results and the analytical ones are observed in figure 9(b).

For further study, figures 10(a) and 11(a) show the numerical results of the evolutions of the CKdV solitary waves at different times when $\epsilon =0.02$ (${\phi _m}=0.018$) and $\epsilon =0.03$ (${\phi _m}=0.027$), respectively. The corresponding comparisons between numerical results and the analytical ones are shown in figures 10(b) and 11(b). Good agreement between PIC numerical results and the analytical ones is observed in figures 9(b) and 10(b). However, in figure 11(b), differences between PIC numerical results and the analytical ones are observed.

Figure 10. (a) PIC simulation results of the evolution of the CKdV solitary waves at different times, when $\epsilon =0.02$; (b) comparisons between PIC simulation results and the analytical ones at different times, when $\epsilon =0.02$.

Figure 11. (a) PIC simulation results of the evolution of the CKdV solitary waves at different times, when $\epsilon =0.03$; (b) comparisons between PIC simulation results and the analytical ones at different times, when $\epsilon =0.03$.

Though the solitary wave exists, it seems from figure 11 that there are some spatial asymmetry especially for large perturbation amplitude and long-time simulation. This is due to the following reason. The expression of the solitary wave from the RPM is only valid for the small-amplitude and long-wavelength limitation. In the simulation, we use the initial condition from the expression of the solitary wave from the RPM. If the limit condition is not satisfied, for example, a large amplitude solitary wave, the real solitary waveform is different from that of the expression of the solitary wave from the RPM. Therefore, as the solitary wave propagates, the uncorrected initial form of the solitary wave will evolve into a real solitary wave and some linear waves or radiations. Finally, spatial asymmetry in the solitary wave propagation in PIC simulations appears for large perturbation and long-time simulation.

To further understand how the differences between PIC numerical results and the analytical ones depend on the parameter $\epsilon$, or the wave amplitude, the variations of both the numerical results and the analytical ones of the propagation velocity of solitary waves with respect to different amplitudes of solitary waves are given in figure 12.

Figure 12. Numerical results and the analytical ones of the dependence of the propagation velocity (${V }/{C_d}$) of the CKdV solitary wave on the wave amplitude (${\phi _m }/({kT_{{\rm eff}}/e})$).

It is noted that when the amplitude of the solitary wave is small, i.e. the parameter $\epsilon$ is small, the PIC numerical results are in good agreement with the analytical ones. However, the differences between the two are observed when the wave amplitude is large enough. Moreover, the difference between the two increases as the wave amplitude increases.

5. Discussion and conclusion

In the present paper, two different reductive perturbation methods are used to derive the KdV equation and the CKdV equation, respectively, in two-temperature-ion dusty plasma. It is found by using the PIC numerical method that the reductive perturbation method to derive the KdV equation is invalid when the nonlinear coefficient $A$ of the KdV equation tends to zero. The application scope of the perturbation method to derive the KdV equation with respect to the nonlinear coefficient $A$ is given. It is noted that the KdV solitary wave does not exist when $A<0.06$, while it exists when $A>0.06$. However, from the numerical results, it is found that under the same parameters of the plasma, the CKdV solitary wave exists in the system even when the coefficient $A$ is close to zero. Furthermore, we find that a larger amplitude of the solitary wave results in a greater deviation between the numerical results and the theoretical results, i.e. when the amplitude of the solitary wave is small enough, the perturbation method to derive the CKdV equation is valid. However, for the KdV solitary wave, the reductive perturbation method is only valid not only if the amplitude of the KdV solitary wave is small enough, but also if $A>0.06$.