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 Reference Kalman, Rosenberg and DeWitt2000; Shukla Reference Shukla2001; Xie & He Reference Xie and He2001; Wang & Wu Reference Wang and Wu2002; Fortov et al. Reference Fortov, Ivlev, Khrapak, Khrapak and Morfill2005; Shukla & Eliasson Reference Shukla and Eliasson2009; Teng et al. Reference Teng, Chang, Tseng and Lin2009; Couëdel et al. Reference Couëdel, Nosenko, Ivlev, Zhdanov, Thomas and Morfill2010; Ghorui, Chatterjee & Wong Reference Ghorui, Chatterjee and Wong2013; Ghorui et al. Reference Ghorui, Samanta, Maji and Chatterjee2014; Merlino Reference Merlino2014; Han et al. Reference Han, Gao, Zhang, Wang and Duan2015; Seadawy Reference Seadawy2015). The remarkable feature of a dusty plasma is that there are low-frequency waves (Rao, Shukla & Yu Reference Rao, Shukla and Yu1990; Shukla Reference Shukla1992; Barkan, Merlino & D'Angelo Reference Barkan, Merlino and D'Angelo1995; Pieper & Goree Reference Pieper and Goree1996; Praburam & Goree Reference Praburam and Goree1996; Merlino et al. Reference Merlino, Barkan, Thompson and D'Angelo1998; Bandyopadhyay et al. Reference Bandyopadhyay, Prasad, Sen and Kaw2008; Dubinov & Sazonkin Reference Dubinov and Sazonkin2013; Kumar Tiwari & Sen Reference Kumar Tiwari and Sen2016; Verheest & Hereman Reference Verheest and Hereman2019).
Many researchers focused on the nonlinear solitary waves in a dusty plasma (El-Labany et al. Reference El-Labany, El-Taibany, Mamun and Moslem2004; Bandyopadhyay et al. Reference Bandyopadhyay, Prasad, Sen and Kaw2008; Emamuddin, Yasmin & Mamun Reference Emamuddin, Yasmin and Mamun2013). A solitary wave was first discovered in 1834. Later, a Korteweg–de Vries (KdV) equation was proposed in 1895 (Korteweg & De Vries Reference Korteweg and De Vries1895). However, until Kruskal and Zabusky discovered a solitary wave solution from a KdV equation in 1965 (Zabusky & Kruskal Reference Zabusky and Kruskal1965; Olivier, Verheest & Maharaj Reference Olivier, Verheest and Maharaj2016), 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 Reference Washimi and Taniuti1966), fluid dynamics (Ono Reference Ono1991; Duan, Wang & Wei Reference Duan, Wang and Wei1996), nonlinear lattice (Zabusky & Kruskal Reference Zabusky and Kruskal1965; Wadati Reference Wadati1990), Bose–Einstein condensates (Huang, Velarde & Makarov Reference Huang, Velarde and Makarov2001), etc.
Recently, the application scope of the perturbation method (Lee Reference Lee2012; Olivier et al. Reference Olivier, Verheest and Maharaj2016) is studied by using the one-dimensional particle-in-cell (PIC) code (Qi et al. Reference Qi, Xu, Duan and Yang2014; Zhang et al. Reference Zhang, Yang, Xu, Yang, Qi and Duan2014, Reference Zhang, Qi, Duan and Yang2015; Gao et al. Reference Gao, Zhang, Zhang, Li and Duan2017; Zhang et al. Reference Zhang, Yang, Hong, Qi, Duan and Yang2017; Wang et al. Reference Wang, Han, Zhang, Gao, Li, Duan and Zhang2018). 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. Reference Qi, Xu, Duan and Yang2014; Zhang et al. Reference Zhang, Yang, Xu, Yang, Qi and Duan2014).
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:
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:
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
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
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:
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
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:
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.
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 Reference Duan and Shi2003): ${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
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 Reference Duan and Shi2003). We then obtain a CKdV equation as follows:
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
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
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:
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. Reference Qi, Xu, Duan and Yang2014; Zhang et al. Reference Zhang, Qi, Duan and Yang2015, Reference Zhang, Yang, Hong, Qi, Duan and Yang2017):
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. Reference Qi, Xu, Duan and Yang2014; Zhang et al. Reference Zhang, Qi, Duan and Yang2015). A summary of a computational cycle of the PIC method is shown in figure 2.
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:
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.
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$.
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:
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.
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 Reference Lin and Duan2005): $\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.
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.
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.
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$.
Acknowledgements
Editor E. Thomas Jr. thanks the referees for their advice in evaluating this article.
Funding
This work was supported by the National Natural Science Foundation of China (Nos. 11965019, 42004131) and the Foundation of Gansu Educational Committee (No. 2022QB-178).
Declaration of interests
The authors report no conflict of interest.