2.1. Fluid model

The fluid equations are governed by the normalized continuity equation,

(2.1)\begin{equation} \frac{\partial n}{\partial t}+\frac{\partial}{\partial x}\left(nu\right)=0,\end{equation}

the normalized momentum equation,

(2.2)\begin{equation} \frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+\frac{\partial\phi}{\partial x}=0,\end{equation}

and the normalized Poisson equation,

(2.3)\begin{equation} \frac{\partial^{2}\phi}{\partial x^{2}}+n-{\rm e}^{\phi}=0.\end{equation}

In the normalized fluid equations, $n$ denotes the ion number density normalized with respect to the equilibrium ion density $n_{i0}$. The ion fluid velocity is represented by $u$ and is normalized with respect to the ion-acoustic speed for Boltzmann electrons $C_{{\rm IA}}=(k_{B}T_{e}/m_{i})$, where $m_{i}$ denotes the ion mass, and the electrostatic potential $\phi$ is normalized with respect to $k_{B}T_{e}/e$, where $k_{B}$ is the Boltzmann constant, $T_{e}$ is the electron temperature and $e$ is the electron charge. In addition, the time variable $t$ is normalized with respect to the inverse ion plasma frequency $\omega _{{\rm pi}}^{-1}=(m_{i}\varepsilon _{0}/n_{i0}{\rm e}^{2})^{1/2}$, where $\varepsilon _{0}$ is the permittivity of free space. The spatial variable $x$ is normalized with respect to the Debye length $\lambda _{D}=(\varepsilon _{0}k_{B}T_{e}/n_{i0}e^{2})^{1/2}$.

2.2. Reductive perturbation analysis

In the following, we use the results obtained by Washimi & Taniuti (Reference Washimi and Taniuti1966). Rather than producing a full derivation, we report the important definitions and solutions used in this study.

For reductive perturbation analysis, we introduce the following perturbation expansions:

(2.4)\begin{gather} n=1+\varepsilon n_{1}+\varepsilon^{2}n_{2}+\cdots, \end{gather}

(2.5)\begin{gather}u=\varepsilon u_{1}+\varepsilon^{2}u_{2}+\cdots, \end{gather}

(2.6)\begin{gather}\phi=\varepsilon\phi_{1}+\varepsilon^{2}\phi_{2}+\cdots, \end{gather}

along with the moving frame coordinates

(2.7a,b)\begin{equation} \xi=\varepsilon^{1/2}\left(x-t\right),\quad\tau=\varepsilon^{3/2}t.\end{equation}

By substituting these expansions, as well as a Taylor series expansion of the function ${\rm e}^{\phi }$ in Poisson's equation, it follows that the electrostatic potential $\phi _{1}$ satisfies the following KdV equation:

(2.8)\begin{equation} \frac{\partial\phi_{1}}{\partial\tau}+\frac{1}{2}\frac{\partial^{3}\phi_{1}}{\partial\xi^{3}}+\phi_{1}\frac{\partial\phi_{1}}{\partial\xi}=0,\end{equation}

while $n_{1}=u_{1}=\phi _{1}$.

For the purpose of this paper, we are interested in two solutions. The first of these is the single soliton solution, given by

(2.9)\begin{equation} \phi_{1}(\xi,\tau)=3c\,\text{sech}^{2}\left[\sqrt{\frac{c}{2}}\left(\xi-c\tau\right)\right]. \end{equation}

To express the solution in the original coordinates, we note that $\phi \approx \varepsilon \phi _{1}$. By setting ${c=1}$ and using the original coordinates as expressed in (2.7a,b), one obtains the following solution for the electrostatic potential:

(2.10)\begin{equation} \phi(x,t)=3\varepsilon\,\text{sech}^{2}\left[\sqrt{\frac{\varepsilon}{2}}\left(x-(1+\varepsilon)t\right)\right]. \end{equation}

It follows that the amplitude of the soliton is given by $A=3\varepsilon$, while the speed is given by $M=1+\varepsilon$. Therefore, for any choice of $M>1$, one can determine the value of $\varepsilon$ to determine the solution. However, it should be noted that the approximation works best for $\varepsilon \ll 1$. In addition, since $n_{1}=\phi _{1}$ and $u_{1}=\phi _{1}$, one obtains the following solutions for the number density and ion fluid velocity:

(2.11)\begin{gather} n=1+3\varepsilon\,\text{sech}^{2}\left[\sqrt{\frac{\varepsilon}{2}}\left(x-(1+\varepsilon)t\right)\right], \end{gather}

(2.12)\begin{gather}u=3\varepsilon\,\text{sech}^{2}\left[\sqrt{\frac{\varepsilon}{2}}\left(x-(1+\varepsilon)t\right)\right]. \end{gather}

The second solution we consider is the two-soliton solution. In particular, we consider the two-soliton solution where the excess velocity of the fast soliton is twice as large as the excess velocity of the slow soliton. Here, the excess velocity is the speed in excess of the ion acoustic speed. More specifically, if the speed of the soliton is given by $M$, the excess speed is given by $M-1$. The solution is given by

(2.13)\begin{equation} \phi_{1}(\xi,\tau)=12\frac{2\cosh^{2}\dfrac{\eta_{1}}{\sqrt{2}}+\text{sinh}^{2}\eta_{2}}{\left[a_{-} \text{cosh}\left(\eta_{2}+\dfrac{\eta_{1}}{\sqrt{2}}\right)+a_{+}\text{cosh}\left(\eta_{2}-\dfrac{\eta_{1}}{\sqrt{2}} \right)\right]^{2}}, \end{equation}

where $a_{\pm }=\sqrt {2}\pm 1$ and $\eta _{j}=\xi -jt$ for $j=1,2$. We can express the solution in terms of the original coordinates, so that

(2.14)\begin{equation} \phi(x,t)=12\varepsilon\frac{2\cosh^{2}\zeta_{1}+\text{sinh}^{2}\zeta_{2}}{\left[a_{-}\text{cosh} \left(\zeta_{2}+\zeta_{1}\right)+a_{+}\text{cosh}\left(\zeta_{2}-\zeta_{1}\right)\right]^{2}},\end{equation}

where $\zeta _{1}=\sqrt {{\varepsilon }/{2}}(x-(1+\varepsilon t)$ and $\zeta _{2}=\sqrt {\varepsilon }(x-(1+2\varepsilon )t$. During periods of the solution where $|t|\gg 1$, the fast soliton propagates with a speed of $M_{f}=1+2\varepsilon$, while the slow soliton propagates with speed $M_{s}=1+\varepsilon$. In addition, the corresponding solutions for the ion density and ion fluid velocity are given by

(2.15)\begin{equation} n(x,t)=1+12\varepsilon\frac{2\cosh^{2}\zeta_{1}+\text{sinh}^{2}\zeta_{2}}{\left[a_{-}\text{cosh}\left(\zeta_{2}+\zeta_{1}\right)+a_{+}\text{cosh}\left(\zeta_{2}-\zeta_{1}\right)\right]^{2}} \end{equation}

and

(2.16)\begin{equation} u(x,t)=12\varepsilon\frac{2\cosh^{2}\zeta_{1}+\text{sinh}^{2}\zeta_{2}}{\left[a_{-}\text{cosh}\left(\zeta_{2}+\zeta_{1}\right)+a_{+}\text{cosh}\left(\zeta_{2}-\zeta_{1}\right)\right]^{2}}, \end{equation}

respectively.

As an example, consider the two-soliton solution (2.14) with $\varepsilon =0.1$, corresponding to the collision of two solitons, where the slow and fast solitons propagate with speeds $M_{s}=1.1$ and $M_{f}=1.2$, respectively. In figure 1, we show different representations of the solution. In figure 1(a), the solution is shown in terms of the original coordinates by using the perturbation expansions (2.5) and transformations (2.7a,b). Unfortunately, the width of the solitons is small relative to the total spatial domain so that the resulting solutions only produce thin lines that provide little detail. To gain a better perspective, we plot the solutions in terms of moving frames coordinate defined in terms of the slow and fast soliton speeds, defined as

(2.17a,b)\begin{equation} \xi_{s}=x-M_{s}t,\quad \xi_{f}=x-M_{f}t. \end{equation}

Figure 1(b) shows the solution plotted with respect to the slow soliton time frame $\xi _{s}$. Notice that prior to the collision, the slow soliton (light green line) is vertical, indicating that its speed coincides with the moving frame. After the collision, we see that this line shifts to the left, indicative of the phase shift that results from the collision. Following the collision, the shifted curve remains vertical, indicating that the speed after the collision is unaffected. Similarly, figure 1(c) shows the propagation of the fast soliton (yellow curve) to remain unchanged after the collision, except for a phase shift to the right. Due to the full recovery of both solitons after the collisions, these collisions are referred to as elastic collisions.

Figure 1. Different graphical representations of the time evolution of the electrostatic potential $\phi$ corresponding to the two-soliton solution with $\varepsilon =0.1$. (a) Solution plotted in terms of the original coordinates of $x$ and $t$. (b) and (c) Solutions plotted in the moving frames $\xi _{s}$ and $\xi _{f}$, respectively.

2.3. Sagdeev pseudopotential analysis

An important aspect of reductive perturbation analysis is the fact that it is limited to the small-amplitude regime. For larger amplitudes, higher order nonlinear effects become significant and even dominant for large amplitudes. For such solutions, Sagdeev pseudopotential analysis provides an alternative approach that retains the full nonlinearity of the system. However, the resulting analysis is only relevant to the construction of single-soliton solutions.

Sagdeev pseudopotential analysis introduces a moving frame of the form

(2.18)\begin{equation} \xi=x-Mt, \end{equation}

where $M$ is the propagation speed. The idea is then to look for solutions that remain constant in this frame of reference. A necessary condition for obtaining soliton solutions is that one must have a propagation speed that exceeds the acoustic speed, that is, $M>1$. In addition, we impose asymptotic boundary conditions that ensure that the plasma returns to the equilibrium far away from the soliton, given by the following set of boundary conditions: $n\rightarrow 1$, $u\rightarrow 0$ and $\phi \rightarrow 0$ when $|\xi |\rightarrow \infty$. By substituting the moving frame variable into the continuity equation and performing a straightforward integration yields

(2.19)\begin{equation} u=M-\frac{M}{n}.\end{equation}

Similarly, the momentum equation can be integrated in a straightforward manner. By using (2.19) to eliminate $u$, one obtains

(2.20)\begin{equation} n=\frac{1}{\sqrt{1-\dfrac{2\phi}{M^{2}}}}.\end{equation}

Substitution of (2.20) into (2.19) allows us to also express the fluid velocity $u$ in terms of $\phi$, given by

(2.21)\begin{equation} u=M-\sqrt{M^{2}-2\phi}.\end{equation}

The final step is to use Poisson's equation to obtain the electrostatic potential $\phi$. By substituting the moving frame variable into (2.3), Poisson's equation becomes

(2.22)\begin{equation} \frac{{\rm d}^{2}\phi}{{\rm d}\xi^{2}}={\rm e}^{\phi}-\frac{1}{\sqrt{1-\dfrac{2\phi}{M^{2}}}}.\end{equation}

By multiplying this equation by ${{\rm d}\phi }/{{\rm d}\xi }$, integrating and applying the boundary conditions, one obtains the following first-order ordinary differential equation (ODE) in the form of an energy integral:

(2.23)\begin{equation} \frac{1}{2}\left(\frac{{\rm d}\phi}{{\rm d}\xi}\right)^{2}+V(\phi)=0,\end{equation}

where the Sagdeev potential $V(\phi )$ is given by

(2.24)\begin{equation} V(\phi)=M^{2}\left[1-\sqrt{1-\frac{2\phi}{M^{2}}}\right]+1-{\rm e}^{\phi}.\end{equation}

Unfortunately, (2.23) and (2.24) cannot be integrated analytically. However, one can integrate the equation numerically to construct the soliton solutions. Once the electrostatic potential is constructed numerically, one can substitute the solutions into (2.20) and (2.21) to obtain the ion number density and ion fluid velocity, respectively The unstable nature of the solutions of (2.23) leads to some inaccuracies for the region where $|\phi |\ll 1$. To deal with this, one can apply an asymptotic analysis, as discussed in § 4. Throughout this paper, this novel approach is used to construct initial conditions for solitons.

2.4. Higher-order effects on single-soliton solutions

For solitons with small amplitudes and speeds only marginally faster than the acoustic speed $M_{a}=1$, the difference between the single-soliton solutions obtained from the KdV equation and the Sagdeev pseudopotential are very small. Gradually, as the soliton amplitude increases, the differences become increasingly more apparent, thus indicating that higher order nonlinear effects become significant. To illustrate this, let us consider the number densities obtained for solutions with small amplitude that arises when $M=1.01$ and a relatively large amplitude that arises when $M=1.3$. The results are shown in figure 2. In panel (a), we see that the solution obtained from the KdV equation agrees well with the Sagdeev potential. The two sets of solutions are almost identical except at the peak, where we see that the KdV solution slightly underestimates the maximum number density. However, in panel (b), we see a dramatic difference between the KdV approximation and the exact solution obtained from the Sagdeev approach. We see that the solution obtained from the RPT overestimates the width, but significantly underestimates the peak density.

Figure 2. Comparison of solutions for the ion number density obtained from Sagdeev pseudopotential analysis (blue curves) and RPT (red curves): (a) solutions for $M=1.01$; (b) solutions for $M=1.3$.

This leads to the question: how does the increase in peak density affect overtaking collisions of solitons? In particular, two-soliton solutions of the KdV equations show that overtaking soliton collisions are elastic, with collisions only resulting in a phase shift. For large-amplitude soliton collisions, are these collision properties still valid or do higher order nonlinear effects, such as larger peak ion densities, lead to inelastic collisions? Since there is no analytical approach available to study this question, we investigate this problem by means of numerical simulation.

4.1. Constructing soliton initial conditions on arbitrary interval lengths

Despite the analytical form for the Sagdeev potential, (2.23) has no closed form solutions. As such, one must integrate (2.23) numerically to obtain the solution of $\phi$. To do so, we consider the initial value problem given by (2.22) with initial conditions given by $\phi (0)=\phi _{r}$ and ${{\rm d}\phi }/{{\rm d}\xi }=0$, where $\phi _{r}$ is the positive root of the Sagdeev potential, that is, $V(\phi _{r})=0$ with $\phi _{r}>0$. Due to the symmetry of soliton solutions, it then follows that $\phi (\xi )=\phi (-\xi )$, so that the solution only has to be constructed for $\xi >0$ and can then be reflected about the $\phi$ axis to produce the remainder of the solution.

One issue that arises from this approach is that the solution of the ODE initial value problem is unstable. As a result, numerical and round-off errors lead to inaccuracies in the limit when $\phi \rightarrow 0$. To understand this problem, it is useful to refer to the potential well analogy of (2.23), where the Sagdeev potential acts as a frictionless potential well for a fictitious particle. For the particle to finish on the top of the potential well, the solution must be solved exactly, otherwise the particle will either fall back into the well (underestimation) or fall off on the other side of the well (overestimation).

To deal with the tails, we briefly review the results of Olivier et al. (Reference Olivier, Verheest and Maharaj2017), where an asymptotic analysis was used to compare the tails of regular solitons with those of acoustic speed solitons. The idea is to use a Taylor series analysis to fit a parabola for $|\phi |\ll 1$. Since $V(0)=V^{\prime }(0)=0$, it follows that the Taylor series expansion of $V$ about $\phi =0$ is given by

(4.1)\begin{equation} V\left(\phi\right)=\tfrac{1}{2}V^{\prime\prime}\left(0\right)\phi^{2}+\mathcal{O}\left(\phi^{3}\right), \end{equation}

so that we can approximate the Sagdeev potential as

(4.2)\begin{equation} V\left(\phi\right)\approx{-}\frac{M^{2}-1}{2M^{2}}\phi^{2}, \end{equation}

provided that $|\phi |\ll 1$. By substituting this approximation into (2.23), one can easily integrate the equation analytically to obtain

(4.3)\begin{equation} \phi\left(\xi\right)=C\exp\left({\pm}\sqrt{\frac{M^{2}-1}{M^{2}}}\xi\right),\end{equation}

where $C$ is an integration constant. To construct the initial condition, we solve the initial value problem (2.22) numerically until one observes a sufficiently small solution of $\phi$. One then substitutes the given $\xi$ value, say $\xi _{0}$ into (4.3), to determine the constant of integration $C$. Clearly for $\xi \rightarrow +\infty$, one must use the lower minus sign in the exponent, so that

(4.4)\begin{equation} C=\phi_{0}\exp\left(\sqrt{\frac{M^{2}-1}{M^{2}}}\xi_{0}\right). \end{equation}

As an example, we show the Sagdeev potential for $M=1.1$ in figure 3(a). Here the red dot represents a pseudoparticle for illustrative purposes. Now, if the particle is released from this position, the absence of friction will mean that it will approach the local maximum of the potential well at the origin, and the position $\phi$ will approach zero as $\xi$ approaches infinity, as $\xi$ plays the role of time in the analogy. However, the solution obtained from numerically integrating (2.22) is shown in figure 3(b) by the blue curve. Here we see that the solution behaves appropriately for $0\leqslant \xi \leqslant 40$. However, at some point, the value of $\phi$ starts to increase. This corresponds to the particle returning to its original position as if it does not have enough energy to reach the top. However, as stated earlier, this is due to numerical errors and the unstable nature of the solution. Indeed, numerically, one always observes either that the particle returns to $\phi \approx \phi _{r}$ or that the particle overshoots the origin, resulting in $\phi \rightarrow -\infty$. The red curve shows the solution obtained by fitting the asymptotic tail at $\phi \approx 10^{-4}$. Here, the curve decreases exponentially as one would expect. The resulting solution can be constructed for the necessary interval length without any stability issues.

Figure 3. Construction of the tail of the soliton initial condition. In panel (a), the Sagdeev pseudopotential well is shown with the blue curve, while the red dot indicates the fictitious particle. The red part of the potential near the origin shows the part of the curve where the asymptotic expansion is applied. In panel (b), the blue curve shows the solution obtained from numerical integration only, while the red curve shows the soliton solution obtained by fitting the asymptotic tail.

5.1. Small-amplitude soliton collision $M_{s}=1.01$

To start off, we consider the simulation results from a collision between two solitons, where the slow soliton propagates with a speed of $M_{s}=1.01$ and the fast soliton has a speed of $M_{f}=1.02$. The amplitudes of the two solitons are given by $\phi _{s}\approx 0.029777$ and $\phi _{f}\approx 0.059117$, respectively. These values agree well with the corresponding values associated with reductive perturbation analysis, given by $\phi _{s}=0.03$ and $\phi _{f}=0.06$, respectively. This is to be expected, as both solitons are in the small-amplitude regime.

The large widths of the solitons require a sufficiently large choice of truncated interval length $L$ to avoid soliton overlap at the initial condition. Here we used a truncated interval length $L=1000$. Conversely, the small gradients associated with the large widths allow us to use a relatively sparse spatial grid. For $M_{s}=1.01$, the simulations provided accurate results for a spatial step size of $\Delta x=1$ and a temporal step size of $\Delta t=0.1$. Since the difference between the speeds is small, the time integration must be performed over a substantial time period. Here, we solved the solution for $0\leqslant t\leqslant 60\,000$.

The result of the simulation is shown in figure 5. In figure 5(a), the solution is shown in terms of the slow moving frame $\xi _{s}=x-M_{s}t$. After the collision at $t\approx 30\,000$, we see that the slow soliton (green line) remains vertical, indicative that the propagation of the slow soliton after the collision is unchanged. In addition, we see that the soliton re-emerges slightly to the left of its original position, indicating a phase shift to the left.

Figure 5. Simulation results for a soliton collision with slow and fast soliton speeds of $M_{s}=1.01$ and $M_{f}=1.02$, respectively. The moving frame references $\xi _{s}=x-M_{s}t$ and $\xi _{f}=x-M_{f}t$ are used in panels (a) and (b), respectively. In panel (c), the solid blue line shows the simulation results at $t=60\,000$, while the associated two-soliton solution from RPT (2.14) is shown with the black dots.

To consider the effect of the collision on the fast soliton, we plot the solution in terms of the fast moving frame $\xi _{f}=x-M_{f}t$ in figure 5(b). The fact that the fast soliton (yellow curve) remains vertical after the collision shows that the speed of the fast soliton is also unaffected by the collision. In addition, we see that the fast soliton re-emerges to the right of its original position, indicating a phase shift towards the right.

Figure 5(c) shows a comparison between the two-soliton solution from RPT (2.14) and the simulation result at the termination time $t=60\,000.$ The solid blue line shows the result obtained from the simulation, while the black dots show the two-soliton solution obtained from RPT. This panel shows good agreement between the simulation result and RPT. This is to be expected, as the solitons fall well within the small-amplitude regime.

This leads to the main question of the paper: how do higher order nonlinear effects affect soliton collisions in terms of elasticity and phase shift? To investigate this question, we next consider a collision of two solitons in the large-amplitude regime.

5.2. Large-amplitude soliton collision $M_{s}=1.2$

We now consider the simulation results for a collision between two large-amplitude solitons, where the slow soliton propagates with a speed of $M_{s}=1.2$ and the fast soliton has a speed of $M_{f}=1.4$. The amplitudes of the two solitons are given by $\phi _{s}\approx 0.52439$ and $\phi _{f}\approx 0.93827$, respectively. For solitons with such large amplitudes, the effects of higher-order nonlinearities can no longer be ignored, so that one would expect deviations from the small-amplitude regime.

For these simulations, the widths of the solitons are fairly small, so that a significantly smaller interval length of $L=200$ can be used. However, the large gradients associated with these small widths require us to choose a much finer spatial grid than for the previous example, namely a spatial step size of $\Delta x=0.01$ and a temporal step size of $\Delta t=0.001$. Fortunately, the large difference between the speeds means that the time integration can be performed over a relatively short time span. Here, we solved the solution for $0\leqslant t\leqslant 500.$

The results of the simulation are shown in figure 6. In figure 6(a), the solution is shown in terms of the slow moving frame $\xi _{s}=x-M_{s}t$. After the collision at $t\approx 250$, we see that the slow soliton (green line) remains vertical. Remarkably, this shows that the propagation of the slow soliton after the collision is unchanged. In addition, the leftward phase shift is clearly visible. Similarly, figure 6(b) shows that the fast soliton also recovers its original speed after the collision. In this panel, the rightward phase shift of the fast soliton is clearly visible.

Figure 6. Simulation results for a soliton collision with slow and fast soliton speeds of $M_{s}=1.2$ and $M_{f}=1.4$, respectively. The moving frame references $\xi _{s}=x-M_{s}t$ and $\xi _{f}=x-M_{f}t$ are used in panels (a) and (b), respectively. In panel (c), the solid blue line shows the simulation results at $t=500$, while the associated two-soliton solution from RPT is shown with the black dashed line.

To compare the simulation results with the associated result from RPT, we plot the result of the simulation along with the two-soliton solution (2.14) in figure 6(c) at the termination time $t=500.$ The solid blue line shows the result of the solution obtained from the simulation, while the black dashed line show the two-soliton solution of RPT. As mentioned in § 2.4, the difference in shapes is due to higher order nonlinear effects. The most important aspect that is shown here is that there is a significant difference in phase shifts. In particular, the leftward phase shift of the slow soliton is less than that predicted by the two-soliton solution (2.14) from RPT. Similarly, the rightward phase shift of the fast soliton is smaller than that predicted by RPT.

To summarize, the simulation shows that the elastic nature of overtaking collisions is conserved in the large-amplitude regime. The only effect of higher order nonlinearities is a reduction in the magnitude of the phase shifts of both slow and fast solitons. To investigate this further, we now take a closer look at the effect of higher order nonlinearities on phase shifts.

5.3. Higher order nonlinear effects on phase shift

The results from the large-amplitude simulation revealed that the only effect of higher order nonlinear effects is to reduce the phase shift of the solitons after the collision. In the following, we make a comparison between phase shifts obtained from simulations with those obtained from the two-soliton solution (2.14) via RPT.

In table 1, we show the two sets of phase shifts corresponding to different speeds. For both the simulation and the two-soliton solution, the phase shifts were determined by comparing the solutions before and after the collision. For the slow solitons, we see good agreement for speeds $M_{s}\leqslant 1.1$. Beyond this, we see that the differences between the KdV and simulation results grow increasingly fast. For large speeds, we see that the phase shift of the slow soliton is closer to the simulation of the fast soliton (in magnitude) than to the phase shift predicted by RPT. Similarly, the difference of the phase shifts between the simulated and RPT for the fast soliton becomes larger as the speed increases. Note that increasing speed also leads to an increase in amplitude, thereby amplifying the effects of higher order nonlinearities.

Table 1. Comparison between phase shifts predicted by RPT and obtained from simulations for different speeds.