1 Introduction
Reductive perturbation analysis was introduced into plasma physics by Washimi & Taniuti (Reference Washimi and Taniuti1966), who illustrated its ability to study small-amplitude solitons, based on the Korteweg–de Vries (KdV) equation for acoustic-type nonlinear modes. Such solitons have been known for more than a century from the study of nonlinear waves on the surface of shallow water by Korteweg & De Vries (Reference Korteweg and De Vries1895). The introduction of these ideas into plasma physics quickly gave rise to a never-ending stream of papers, expanding into generalizations for other nonlinear evolution equations and extensions to various multi-species plasma compositions. From this extensive literature only a relevant selection can be cited.
The general properties of nonlinear solitary waves are twofold. First, there is an intricate relation between speed, amplitude and width, with zero boundary conditions far from the maximum of the structure. Amplitude and width are inversely related in the sense that the amplitude decreases when the width increases, a property that precludes their having a linear counterpart. The other typical characteristic is that larger and faster solitons overtake smaller and slower ones, and the solitons come out of the overtaking nonlinear interaction process seemingly unchanged. Again, such solitary waves have no linear counterpart.
Noting that periodic modes are possible in many different media, we are quite familiar with their description when the nonlinearities are omitted, yielding typical harmonic profiles. Extensions to nonlinear periodic modes, in particular in plasma physics and for not too large nonlinearities, have introduced cnoidal waves, with properties that are quite different from the soliton properties. For cnoidal waves particular care has to be taken with respect to the boundary or initial conditions, since for periodic modes the structure pattern is repeated and one cannot rely on conditions at infinity, periodic functions having no well-defined limits there.
In early studies of cnoidal wave solutions in plasmas, great care was taken in the treatment of the boundary conditions, see, for example Konno, Mitsuhashi & Ichikawa (Reference Konno, Mitsuhashi and Ichikawa1979), Ichikawa (Reference Ichikawa1979), Tiwari, Jain & Chawla (Reference Tiwari, Jain and Chawla2007) and Prudskikh (Reference Prudskikh2009). An important characteristic of these solutions is that they were constructed in a way that ensures conservation of number density. The conservation of number density is related to the conservation of mass, as is discussed in detail in Appendix A. Another important result is the fact that solutions propagating at both superacoustic and subacoustic solutions were reported, see for example figure 1 in Konno et al. (Reference Konno, Mitsuhashi and Ichikawa1979).
Unfortunately, more recent studies of cnoidal waves in plasmas have ignored these facts, where conservation of number density has been completely ignored. Examples of this can be found in Saini & Sethi (Reference Saini and Sethi2016), Ur-Rehman, Mahmood & Hussain (Reference Ur-Rehman, Mahmood and Hussain2017), Tolba et al. (Reference Tolba, Moslem, Elsadany, El-Bedwehy and El-Labany2017), Singh et al. (Reference Singh, Ghai, Kaur and Saini2018), Liu et al. (Reference Liu, Goree, Flanagan, Sen, Tiwari, Ganguli and Crabtree2018), Farhadkiyaei & Dorranian (Reference Farhadkiyaei and Dorranian2018), Kaur et al. (Reference Kaur, Singh, Kohli and Saini2018), Ur-Rehman, Mahmood & Kaladze (Reference Ur-Rehman, Mahmood and Kaladze2019), Ur-Rehman & Mahmood (Reference Ur-Rehman and Mahmood2019), Ali et al. (Reference Ali, Masood, Rizvi and Mirza2020), Tamang & Saha (Reference Tamang and Saha2020) and Mehdipoor & Asri (Reference Mehdipoor and Asri2022). In addition, the authors used similar boundary conditions to those satisfied by solitons. For the purpose of this paper, we refer to this formalism as the soliton boundary approach to deriving cnoidal wave solutions. The resulting solutions clearly violate number density conservation, as is quickly verified by comparing the solutions with those illustrated in § 3 and in Appendix A.
In this paper, we deconstruct the soliton boundary approach, and show that the erroneous boundary conditions can instead be expressed as a set of initial conditions that are consistent with the solutions. Once these details are well understood, a straightforward modification of the initial conditions allows one to obtain cnoidal wave solutions that conserve number density.
It is worth pointing out that, in the original approach of Konno et al. (Reference Konno, Mitsuhashi and Ichikawa1979), a power series expansion is used in order to obtain the final expression, that is expressed in terms of the Jacobi elliptic parameter $k^2$ (note that, in this paper, we use $m=k^2$ as the Jacobi elliptic parameter). One drawback of this approach is that the solutions break down in the small-amplitude limit when $k^2\ll 1$. In our formalism, no approximations are made. For the subacoustic cnoidal wave solutions considered here, one is able to construct solutions in this limit. We use these solutions to show consistency with linear theory. A detailed study of superacoustic cnoidal wave solutions that includes the soliton limit $k^2\rightarrow 1$ is beyond the scope of the present paper, and is a subject for future investigation.
The paper is structured as follows: in § 2, the model is introduced, along with the derivation of linear periodic solutions. In § 3 we provide a full derivation using the soliton boundary approach, where all the mathematical inconsistencies of this approach are highlighted. In addition, we deconstruct the method in the form of an initial value problem rather than a boundary value problem. Section 4 deals with the new formalism in the form of a new set of initial conditions. In § 5, the cnoidal wave solutions are derived, where we limit ourselves to the subacoustic case. The resulting solutions are discussed in detail in § 6. Finally, conclusions and future work are discussed in § 7.
2 Fluid equations and linear study
For the purpose of this paper, we consider one of the simplest fluid models that supports cnoidal wave solutions, namely a plasma consisting of cold singly charged ions and inertialess Boltzmann-distributed electrons. This is done in order to simplify the underlying algebra as much as possible, thus allowing us to focus on the main principles behind the approach. To this end, we consider the normalized fluid equations, given by the continuity, momentum and Poisson equations,
respectively. Here, the usual normalizations apply as follows: the ion number density $n$ is normalized with respect to the equilibrium number density $n_{0}$, and the electrostatic potential $\phi$ is normalized with respect to $\kappa _{B}T_{e}/e$, where $\kappa _{B}$ is the Boltzmann constant, $T_{e}$ is the electron temperature and $e$ is the electron charge. The space coordinate $x$ is normalized with respect to the Debye length $\lambda _{D}=({\varepsilon _{0}\kappa _{B}T_{e}}/{n_{0}e^{2}})^{{1}/{2}}$, where $\varepsilon _{0}$ is the permittivity of free space, and $T_{e}$ is the electron temperature. The time variable $t$ is normalized with respect to the inverse plasma frequency $({\varepsilon _{0}m_i}/{n_{0}e^{2}})^{{1}/{2}},$ where $m_i$ denotes the ion mass. As a result, the ion fluid velocity $u$ is normalized with respect to the ion acoustic speed $c=({\kappa _{B}T_{e}}/{\varepsilon _{0}m_i})$.
The equilibrium state of the plasma is given by $n=1$, $u=0$ and $\phi =0$. For linear theory, we consider a small perturbation of the form
where $n_{1},u_{1},\phi _{1}\sim O(\varepsilon )$ for some $\varepsilon \ll 1$. Substituting this into (2.1), (2.2) and (2.3), and retaining only terms of order $\varepsilon$, one obtains the well-known dispersion relation
From the dispersion relation, one can now construct sinusoidal solutions by means of the superposition principle, giving
A fundamental property of the linear periodic solutions is that it conserves number density. This means that
for all $t$, where $n_{0}=1$ is the equilibrium density, the wavelength is given by $L={2{\rm \pi} }/{k}$, and $x_{0}$ and $t_{0}$ are arbitrary points in space and time, respectively. In Appendix A we show that this conservation is an essential ingredient to obtain physically significant periodic solutions. As such, we incorporate this condition into our construction of the periodic solutions, as detailed in § 4.
3 Soliton boundary approach to cnoidal wave derivation
In order to justify the need for our new approach, we start off with the so-called soliton boundary approach. To this end, we apply this approach to the model under consideration. In this derivation, we point out the questionable treatment and inconsistencies of the integration constants. These issues arise due to an ill-posed set of the boundary conditions. Despite these mathematical concerns, one does obtain cnoidal wave solutions. In § 3.2 we deconstruct this method to show that the treatment of the boundary conditions instead corresponds to a specific choice of initial conditions. In addition, we illustrate that the resulting solutions fail to conserve the number density.
3.1 Methodology
In the soliton boundary approach, one looks for periodic solutions that satisfy the boundary conditions
Now, right at the onset, it should be noted that non-trivial (i.e. non-constant) periodic functions have no asymptotic limit. As such, the very notion of introducing a limit when $x$ or $\xi$ approaches infinity is deeply flawed. Despite this obvious mathematical error, we proceed by introducing the usual reductive perturbation expansions
along with the stretched coordinates
Substitution of the expansions (3.2), (3.3), (3.4), along with the stretched coordinates (3.5a,b) into the continuity equation (2.1) leads to an equation that yields terms of differing orders of $\varepsilon$. At the lowest order $O(\varepsilon ^{{3}/{2}})$, one obtains
while the higher order $O(\varepsilon ^{{5}/{2}})$ yields
Similarly, the momentum equation (2.2) gives
at order $O(\varepsilon ^{{3}/{2}})$, and
at order $O(\varepsilon ^{{5}/{2}})$. Finally, Poisson's equation gives
at order $O(\varepsilon )$, and
at order $O(\varepsilon ^{2}).$ Following the usual reductive perturbation approach, we start by solving the lowest-order equations (3.6), (3.8) and (3.10), and use the resulting solutions to solve the higher-order equations (3.7), (3.9) and (3.11).
One can solve (3.6) by integrating with respect to $\xi,$ resulting in the equation
where $C_{1}$ is an integration constant. To find $C_{1}$, we use the ill-posed boundary conditions, namely that $n\rightarrow 1$ and $u\rightarrow 0$ when $x\rightarrow \infty$. It then follows that $n_{1}\rightarrow 0$ and $u_{1}\rightarrow 0$ when $\xi \rightarrow \infty$, so that $C_{1}=0$ and
We now proceed to solve the lowest-order equation associated with the momentum equation (3.8). Once again, a straightforward integration leads to
where $C_{2}$ is an integration constant. As before, the ill-posed boundary conditions $u\rightarrow 0$ and $\phi \rightarrow 0$ yield that $C_{2}=0$. From (3.13) it follows that
so that
Substituting this expression of $n_{1}$ into the lowest-order equation derived from Poisson's equation (3.10), one obtains the compatibility condition
Hence, it follows that $M=\pm 1$. Since we are interested in waves propagating in the direction of the magnetic field, we choose $M=1$. From (3.13) and (3.16) it then follows that $n_{1}=u_{1}=\phi _{1}$.
We now proceed to the higher-order equations. Substituting $M=1$, $n_{1}=\phi _{1}$ and $u_{1}=\phi _{1}$ into (3.7) gives the following partial differential equation:
Similarly, (3.9) gives
Substituting (3.18) into (3.19) then leads to
Finally, the KdV equation is obtained by differentiating (3.11) with respect to $\xi$, and substituting (3.20) into the resulting equation, giving
In order to derive the cnoidal wave solutions for the KdV equation (3.21), it is customary to look for solutions that remain constant in the moving frame that propagate slightly faster than the acoustic speed. We therefore introduce the variable
In terms of the original coordinates, this is given by
By substituting the variable (3.22) into the KdV equation (3.21), one reduces the KdV equation to the following third-order ordinary differential equation (ODE):
Integrating this equation once with respect to $\eta$ gives
where $C_{3}$ is a constant of integration. The usual approach is to choose $C_{3}=0$, based on the argument that if $\phi \rightarrow 0$ in the limit when $\eta \rightarrow \infty$, then ${{\rm d}^{2}\phi _{1}}/{{\rm d}\eta ^{2}}\rightarrow 0$ in the same limit. The resulting second-order ODE is then given by
In order to perform another integration, (3.25) is then multiplied with ${{\rm d}\phi _{1}}/{{\rm d}\eta }$ and integrated to get
where $C_{4}$ is a constant of integration. If one were to stay mathematically consistent, one would note that, since $\phi _{1}\rightarrow 0$ when $\eta \rightarrow \infty$, one must then have that ${{\rm d}\phi _{1}}/{{\rm d}\eta }\rightarrow 0$ when $\eta \rightarrow \infty$, similarly to the argument used to obtain $C_3=0$. In this case, the constant $C_{4}=0$, so that the equilibrium is recovered, and no periodic solutions are found. However, here, the soliton boundary approach contradicts itself by choosing $C_{4}$ to be an arbitrary constant. This implies that ${{\rm d}\phi _{1}}/{{\rm d}\eta }$ approaches a non-zero value in the limit $\eta \rightarrow \infty$. Not only does this contradict the use of the previous limits, but it clearly also implies that the solutions are unbounded rather than periodic. Nevertheless, by following this reasoning one obtains the following equation:
The differential equation (3.28) has cnoidal wave solutions of the form
where $\mu _{1}<\mu _{2}<\mu _{3}$ are the real roots of the ‘Sagdeev potential’, i.e. the cubic polynomial
and
The solution (3.29) is periodic with period $L=12K(m)/(\mu _{3}-\mu _{1})$, where $K(m)$ is the complete elliptic integral of the first kind.
Clearly, a necessary condition for the existence of cnoidal wave solutions is that the equation (3.30) must have three distinct real roots. This depends on the choices of $v$ and $C_{4}$. Indeed, for a fixed choice of $v$, one can show that the polynomial (3.30) has three real roots provided that
As an example, let us consider cnoidal wave solutions associated with the soliton boundary approach corresponding to periodic waves propagating with a velocity of $M=1.01$. To achieve this, we set $\varepsilon =0.01$, $v=1$, and use different choices of ${C_{4}\in (0,4/3)}$. To this end, we let $C_{4}=\frac {4}{3}a,$ where $0< a<1$. The solutions are shown in figure 1. Here, panel (a) shows the phase portraits for six different choices of $a$. It is worth pointing out that all these solutions are inconsistent with the original set of boundary conditions. Clearly all choices of $a$ result in solutions that remain positive for all values of $\eta$. Indeed, if $\phi _{1}$ remains positive, it can never approach $0$ in any limit, and therefore violates the original set of boundary conditions. Moreover, these results are inconsistent with the conservation of number density. To show this, consider figure 1(b), where the number density is plotted in the original coordinates at $t=0$. The blue curve corresponds to $a=0.2$, and the red curve corresponds to $a=0.99$. In addition, the black dotted line shows the equilibrium number density. It is clear that neither of these solutions conserves the number density due to the fact that the perturbed periodic wave oscillates above the equilibrium $n=1$ level. As such, it is clear that
where $L$ is the spatial period of the solution, so that the conservation of number density is violated.
3.2 Deconstruction of the method
As we mentioned above, there are a number of mathematical flaws and contradictions in the derivation of the cnoidal wave solutions. In particular, the boundary conditions for the number density $n$, fluid velocity $u$ and electrostatic potential $\phi$ are all inconsistent with periodic solutions. In addition, the integration constant $C_{4}$ is constructed in a way that implies that the derivative ${{\rm d}\phi _{1}}{{\rm d}\eta }$ approaches a non-zero value in the limit $\eta \rightarrow \infty$, clearly contradicting the boundary condition stating that $\phi _{1}\rightarrow 0$ when $\eta \rightarrow \infty$. Despite these fundamental flaws, one obtains cnoidal wave solutions.
As it turns out, the calculation of the integration constants can instead be obtained through the application of a set of initial conditions. The key idea here is to notice that the integration constants are determined at a specific choice of $\xi$ or $\eta$, rather than at the boundaries.
To simplify the deconstruction, we look to express the conditions that allow us to determine the integration constants as an initial condition. To do so, we set $t=0$, so that $\xi =\varepsilon ^{{1}/{2}}x$ and $\eta =\varepsilon ^{{1}/{2}}x$. Let us start with the integration constant $C_{1}$. By setting $C_{1}=0$, it follows that there exists a $\xi _{0}$ such that $n_{1}(\xi _{0})=0$ and $u_{1}(\xi _{0})=0$. By setting $t=0$, it follows that $n(x_{0},0)=1$ and $u(x_{0},0)=0$, where $x_{0}=\varepsilon ^{-{1}/{2}}\xi _{0}$. Similarly, the second integration constant can be obtained by setting $\phi (x_{0},0)=0$. This implies that the function $\phi (x,0)$ has a position $x=x_0$ where the solution satisfies the equilibrium conditions.
The next integration constant was chosen in a way that $C_{3}=0$. At $x=x_{0}$ and $t=0$, this implies that
In other words, this boundary condition shows that the function $\phi _{1}(x,0)$ has a spatial inflection point at $x=x_{0}$ and $t=0$.
For the final integration constant $C_{4}$, it follows that $({{\rm d}\phi _{1}}/{{\rm d}\eta })(\eta _{0})=\pm \sqrt {2C_{4}}$. By choosing the point where the slope reaches its maximum, we can use the positive $+$ sign. If we rewrite this in the original coordinates at $t=0$, it follows that
In other words, at $t=0$ and $x_{0}=\varepsilon ^{-{1}/{2}}\xi _{0}$, the function $\phi (x,0)$ satisfies three conditions, namely (i) the function reaches its equilibrium value, (ii) the function has a non-zero slope and (iii) the function has an inflection point. The latter implies that the function reaches its maximum slope at $x_{0}$.
When expressed in the form of an initial condition problem, the soliton limit approach is mathematically sound. However, the problem of number density conservation still remains. This is due to the fact that the solutions are forced to reach their equilibrium state at the inflection points. In the next section, we derive a new formalism by relaxing this condition.
4 New reductive perturbation formalism for nonlinear periodic waves
In order to address the problem of number density conservation, we now introduce a new formalism that is a modified variation of the deconstructed initial value problem discussed above. To this end, we define the initial conditions at an inflection point $x_{0}$ satisfying
The deconstructed method assumes that one recovers the equilibrium values of the density, fluid velocity and electrostatic potential at the inflection point. However, in our formalism, we relax this requirement by introducing the following initial conditions at the inflection point:
Based on the expression shown in (3.35), we introduce the following condition on the first derivative at the inflection point:
The role of $\alpha$ in the final solution is to produce solutions of different amplitudes, where larger choices of $\alpha$ correspond to solutions with larger amplitudes. Finally, we look for solutions that satisfy the conservation of number density (2.9).
In this formalism we apply the same perturbation expansions (3.2)–(3.4), and the same stretching coordinates (3.5a,b). As such, the continuity equation yields the same lower- and higher-order equations (3.6) and (3.7), respectively. Similarly, the momentum equation yields equations (3.8) and (3.9), respectively, while Poisson's equation yields equations (3.10) and (3.11), respectively.
In order to derive the KdV equation, one must integrate each of the lower-order equations. For the continuity equation, integration of the resulting lower-order equation (3.6) gives
where $C_{1}$ is an integration constant. To find $C_{1}$, we use the initial conditions. Since $n(x_{0},0)=1+\varepsilon n_{10}$ and $u(x_{0},0)=\varepsilon u_{10}$, it follows that $n_{1}(\xi _{0},0)=n_{10}$ and $u_{1}(\xi _{0},0)=u_{10}$, where $\xi _{0}=\varepsilon ^{{1}/{2}}x_{0}$. It follows that $C_{1}=-Mn_{10}+u_{10}$, so that
We now proceed to solve the lowest-order equation associated with the momentum equation, given by (3.8). Once again, a straightforward integration leads to
where $C_{2}$ is an integration constant. Once more, the initial conditions $u(x_{0},0)=\varepsilon u_{10}$ and $\phi (x_{0},0)=\varepsilon \phi _{10}$ yield that $C_{2}=-Mu_{10}+\phi _{10}$. From (4.5) it follows that
so that
Substituting this expression of $n_{1}$ into the lowest-order equation derived from Poisson's equation (3.10), one obtains the compatibility condition
Since $\phi _{1}(x_{0},0)=\phi _{10}$, it follows that $M=\pm 1$, where once more we choose $M=1$ to study waves that propagate in the direction of the magnetic field. Assuming that $u_1=\phi_1$, it follows from (4.5) and (4.8) that $n_{1}=u_{1}=\phi _{1}$ and $n_{10}=u_{10}=\phi _{10}$. Importantly, since $\phi _{1}=n_{1}$, it follows that the conservation of number density requires that
for all $t,$ where $L$ is the spatial period of the periodic wave.
We now proceed to the higher-order equations. Substituting $M=1$, $n_{1}=\phi _{1}$ and $u_{1}=\phi _{1}$ into the higher-order continuity equation (3.7), one recovers the partial differential equation (3.18). From a similar treatment of the higher-order momentum equation (3.9), one recovers (3.19). Substituting (3.18) into (3.19) then leads to (3.20). Finally, the KdV equation is obtained by differentiating (3.11) with respect to $\xi$, and substituting (3.20) into the resulting equation, giving
Notice that, despite the different choice of boundary conditions, one obtains a KdV equation with identical coefficients as the one produced by the soliton limit method (3.21).
To study subacoustic cnoidal wave solutions, we consider solutions for the KdV equation (4.11) that propagate slightly below the acoustic speed moving frame $\xi =\varepsilon ^{{1}/{2}}(x-t)$. To this end, we introduce the moving frame
In terms of the original coordinates, this is given by
That is to say, solutions that are stationary in the frame with velocity $M=1-\varepsilon v$.
By substituting the variable (4.12) into the KdV equation (4.11), one obtains
Integrating this equation once with respect to $\eta$ gives
where $C_{3}$ is a constant of integration. From the initial conditions $\phi _{1}(x_{0},0)=\phi _{10}$ and $({\partial ^{2}\phi }/{\partial x^{2}})(x_{0},0)=0$ it follows that $C_{3}=v\phi _{10}+\frac {1}{2}\phi _{10}^{2}$, so that
To perform another integration, we multiply (4.15) with ${{\rm d}\phi _{1}}/{{\rm d}\eta }$ and integrate to get
where $C_{4}$ is a constant of integration. At $x=x_{0}$ and $t=0$, it follows from the initial conditions that $\phi _{1}=\phi _{10}$, and
From the initial condition $({\partial \phi }/{\partial x})(x_{0},0)=\sqrt {2}\varepsilon ^{{3}/{2}}\alpha$, it therefore follows that
where $\eta _{0}=\varepsilon ^{{1}/{2}}x_{0}$. Consequently, the coefficient in (4.17) is given by $C_{4}=\alpha ^{2}-\frac {2}{3}\phi _{10}^{3}-v\phi _{10}^{2}$, so that
Here, $\phi _{10}\neq 0$ to prevent the reduction to the deconstructed soliton boundary approach, $v$ is a fixed velocity and $\alpha ^{2}$ can be varied to produce different amplitudes of periodic waves travelling at the same speed, similar to the role of the amplitude $a$ in the linear periodic solutions (2.8).
Finally, the differential equation (4.20) has cnoidal wave solutions of the form
where $\mu _{1}<\mu _{2}<\mu _{3}$ are the real roots of the cubic polynomial satisfying
and
The solution (4.21) is periodic with period $L=12K(m)/(\mu _{3}-\mu _{1})$, where $K(m)$ is the complete elliptic integral of the first kind.
5 Analysis of cnoidal wave solutions
One of the main advantages of our derivation is the fact that it allows the roots of the cubic polynomial $\mu _{1}$, $\mu _{2}$ and $\mu _{3}$ to be expressed analytically. This property allows us to construct cnoidal wave solutions that conserve the number density of the plasma in equilibrium. In the following, details of this approach are provided.
5.1 Necessary condition for existence of solutions
The analysis of the cnoidal wave solution is complicated by the fact that the solution has three free variables, namely $v$, $\alpha ^{2}$ and $\phi _{10}$. Here, $v$ is directly related to the propagation speed of the periodic wave. On the other hand, $\alpha ^{2}$ is related to the slope of the periodic solution at the inflection point $\xi =0$. In addition, $\phi _{10}$ is chosen in a way that ensures that the number density is conserved. For the cnoidal wave solution, larger values of $\alpha ^{2}$ correspond to larger fluctuations.
It is clear from the solution (4.21) that finding three real roots to the polynomial (4.22) is crucial in order to analyse the periodic waves. For our analysis, we treat $v$ as a fixed parameter, while $\alpha ^{2}$ is varied within this choice of $v$, and $\phi _{10}$ depends on the choice of $\alpha$.
It is well known that the nature of the roots of a generalized cubic polynomial of the form
depends on the discriminant
see for example Abramowitz & Stegun (Reference Abramowitz and Stegun1965). In particular, if $\varDelta >0$, the polynomial has three distinct real roots. This is a necessary condition for the existence of cnoidal wave solutions. For the polynomial (4.22) we have $a=\frac {1}{3}$, $b=v$, $c=-\left(2v\phi_{10}+\phi_{10}^2\right)$ and $d=-(\alpha ^{2}-\frac {2}{3}\phi _{10}^{3}-v\phi _{10}^{2})$.
In the following, we assume that $v$, $\alpha$ and $\phi _{10}$ are chosen in a way that ensures that the discriminant $\varDelta$ is strictly positive. The relationship between $\alpha ^{2}$ and $\phi _{10}$ will be considered in more detail after the final solution has been derived.
5.2 Analytical expressions for the roots $\mu _{1}$, $\mu _{2}$ and $\mu _{3}$
The three real roots of the cubic polynomial (4.22) satisfying $\varDelta >0$ can be expressed in terms of trigonometric functions as follows:
for $j=0,1,2$, where
and
For the polynomial (4.22), one then obtains the following roots:
and
where
Since $0\leq \theta \leq {\rm \pi}$, one can easily confirm by plotting the functions $\cos ({\theta }/{3})$, $\cos ({(\theta +2{\rm \pi} )}/{3})$ and $\cos ({(\theta +4{\rm \pi} )}/{3})$ that
whenever the polynomial has three real roots. As such, we have that,
By substituting these expressions into (4.21), it follows that the solution can be expressed as
where
and
5.3 Choice of $\phi _{10}$
The general solution (5.13) depends on three free parameters, namely $\alpha$, $v$ and $\phi _{10}$. In order to construct solutions, we consider a fixed choice of the velocity $v$. For the resulting propagation speed we then construct periodic solutions with different amplitudes. These solutions are obtained by varying $\alpha$. Finally, for each choice of $\alpha$ we look for the value of $\phi _{10}$ such that the conservation of number density (4.10) is satisfied. That is, we look for a choice of $\phi _{10}$ such that the following identity holds:
for all $\tau$. Here, $L$ is the wavelength, given by
where $K(\beta _{3})$ is the complete elliptic integral of the first kind.
The integral (5.18) can be expressed analytically, thanks to the following indefinite integral:
Here, the functions $\text {sn}(z\,|\,m)$, $\text {cn}(z\,|\,m)$ and $\text {dn}(z\,|\,m)$ are Jacobi elliptic functions, $\text {am}(z\,|\,m)$ is the Jacobi elliptic amplitude function and $E(z\,|\,m)$ is the incomplete elliptic integral of the second kind.
By using the incomplete integral (5.20) along with the identities
and
it follows that
Notice that, for a fixed choice of $v$ and $\alpha$, this integral depends on $\phi _{10}$, as is obvious from (5.10) and (5.14)–(5.17). As a result, finding $\phi _{10}$ becomes a root-finding problem for the function $I(\phi _{10})$, given by the right-hand side of (5.23).
In order to show the relationship between $\alpha$ and $\phi _{10}$, consider the special case where $v=1$. The dependence of $\phi _{10}$ on $\alpha$ is shown in figure 2. The figure clearly shows that $\phi _{10}$ increases with $\alpha$. In terms of the solution, this means that the vertical position of $\phi _{1}$, the point where the spatial slope reaches its maximum, must increase in order to compensate for the increased asymmetry associated with larger amplitudes. These aspects are explored in more detail in § 6.
6 Results
After obtaining an analytical expression for the solution, and obtaining a compatibility condition resulting in a root-finding problem, we now consider solutions. To investigate the solutions numerically, we note that we can set $v=1$ without loss of generality. This is due to the fact that, for the same set of results for $\phi _{1}$, one can vary the velocity by choosing different values of $\varepsilon$. The only restriction here is that $\varepsilon$ must be sufficiently small to ensure that higher-order nonlinear effects remain negligible. With this in mind, we start by investigating the solutions of $\phi _{1}$ with $v=1$ in detail in terms of the stretched coordinate $\eta$. After this, we consider the effect of variations in $\varepsilon$ on the solutions in terms of the original normalized variables $\phi$, $x$ and $t$.
6.1 Effect of $\alpha$ on $\phi _{1}$ solutions
As discussed earlier, the parameter $\alpha$ is a free parameter that produces different solutions, all propagating with the same speed. The choice of $\alpha$ is related to the slope at the inflection points, and in general, larger choices of $\alpha$ will produce solutions with larger amplitudes. Here, we consider the effect of $\alpha$ on the $\phi _{1}$ solutions for $0<\alpha \leq 5$. It should be pointed out that, when $\alpha$ is chosen too large, the solutions may violate the small-amplitude assumptions associated with the reductive perturbation theory. As such, we did not consider solutions beyond $\alpha =5$.
To start off our conversation, we consider the solution associated with $\alpha =1$. We use this solution to illustrate some important aspects of the relationship between the choice of initial conditions and the resulting solutions.
The solution for $\alpha =1$ is shown in figure 3 for $0\leq \eta \leq 2L$, where $L\approx 4.364096$ is the wavelength. For $\alpha =1$, one obtains $\phi _{10}\approx 0.206277$. Here, we want to point out two important aspects about these solutions. Firstly, $\phi _{10}$ is chosen in a way that ensures that
This means that the shaded area below the line (shown in the light blue shading) is equal to the shaded area above the line (shown with the grey shading). For each choice of $\alpha$, there is a unique value of $\phi _{10}$ that ensures this conserved quantity. Secondly, the value of $\phi _{10}$ corresponds to the value that $\phi _{1}$ obtains at the inflection points. These points correspond to the maximum (positive) and minimum (negative) slopes ${{\rm d}\phi _{1}}/{{\rm d}\eta }$. To illustrate this, we show the tangent line of the inflection point at $\eta \approx 7.837211$ and $\phi _{1}=\phi _{10}\approx 0.206277,$ indicated with the black dot. The slope of the tangent line is given by $\sqrt {2}\alpha$, with $\alpha =1$ in this instance. Notice that this is only one of four inflection points visible in the figure, each occurring where the solution intersects the horizontal $\phi _{1}=\phi _{10}$ line.
Having established some general properties of the cnoidal wave solutions, we now turn our attention to the effect that $\alpha$ has on the solutions. To do this, we plot the solutions for various choices of $\alpha$. Figure 4(a) shows the solutions for $\alpha =0.1$, $\alpha =0.5$ and $\alpha =1$. Here, it is obvious that the amplitude of the wave increases with $\alpha$, as is expected. For $\alpha =0.1$, the solution closely resembles a sinusoidal wave. On the other hand, for $\alpha =1$ we see that the maximum value is greater than $\phi _{1}=1$, whereas the minimum value does not reach $\phi _{1}=-1$. This is indicative of the asymmetry associated with nonlinear effects.
In figure 4(b) we show the solutions for $\alpha =1,2,3,4$ and $\alpha =5$. Here, we see that the asymmetry increases with $\alpha$. This is most obvious for the $\alpha =5$ solution, where the maximum exceeds $\eta =4$, whereas the minimum does not even reach the $\eta =-3$ level. In addition, the asymmetry of this solution is clear when comparing the part of the solution below the $\phi _{1}$ axis with that of the solution above the $\phi _{1}$ axis for the $\alpha =5$ solution. Interestingly, the panel also shows that the wavelength decreases with $\alpha$.
Many papers in the literature represent the solutions in terms of phase space portraits. These solutions are plotted in figure 4(c) for various choices of $0.5\leq \alpha \leq 5$. Here, we see that all trajectories intersect the vertical $d\phi_1/d\eta$ axis. This is in contrast to those shown in many recent papers, where small-amplitude trajectories never intersect this axis. That is, they remain either to the left or to the right of the $d\phi_1/d\eta$ axis, in clear violation of the conservation of number density. Our results therefore provide an important departure from these results. In addition to this, the phase portrait also clearly shows the increased asymmetry resulting from larger choices of $\alpha$. For $\alpha =0.5,$ the trajectory is nearly elliptic, corresponding to a sinusoidal solution. However, as $\alpha$ is increased, the trajectories are further and further deformed to a teardrop shape, with the difference in magnitudes between the minimum and maximum $\phi _{1}$ values increasingly widening.
The solutions were also plotted in Sagdeev potential form in figure 4(d), as is commonly done in the literature. Here, we clearly see that the two relevant roots are on opposite sides of the $\phi _{1}$ axis, again indicating that the solutions oscillate around the equilibrium value of $\phi _{1}=0$. Interestingly, this representation clearly shows that the values of $\phi _{10}$, represented by the $\phi _{1}$ value at the local minimum, increase with $\alpha$.
Based on these results, we take a closer look at the effect of $\alpha$ on two critical characteristics of the solutions, namely the wavelength and extreme values of the solutions. To investigate the former, we plotted the wavelength of $\phi _{1}$ on the interval $0<\alpha \leq 5$ in figure 5(a). Note that the wavelength is simply given by
Here, we can clearly see that the wavelength decreases with $\alpha$, in agreement with the results of figure 4(b). This result is somewhat surprising due to the fact that $\beta _{3}$ increases with $\alpha$, combined with the fact that $K(\beta _{3})$ is unbounded in the limit when $\beta _{3}\rightarrow 1$. However, it turns out that the increase of $\beta _{2}$ is faster than the increase of $K(\beta _{3})$ associated with increasing values of $\alpha$. It should be noted that this is a direct result of including $\phi _{10}$ in the initial conditions. As such, this is a novel result.
In addition, we calculated the maximum and minimum values of $\phi _{1}$ for different values of $\alpha$. These quantities are simply given by
The results are shown in figure 5(b) for $0<\alpha \leq 5$. Here, we can clearly see that the maximum value of $\phi _{1}$ increases much faster than the decrease in the minimum of $\phi _{1}$ associated with an increase in $\alpha$. This also reflects the asymmetry that is associated with large values of $\alpha$.
6.2 Solutions of $\phi$ in original coordinates
Until now, we only considered the solution $\phi _{1}$ in the stretched coordinates. A crucial aspect of the derivation of the KdV equation (4.11) from reductive perturbation analysis is the fact that $\varepsilon$ must be sufficiently small to ensure that higher-order nonlinear effects are negligible. In the following, we consider how this limitation affects the solution in the original frame of reference, that is, in terms of the original coordinates $x$ and $t$.
To start the analysis, we remind the reader that
so that smaller choices of $\varepsilon$ reduce the amplitude of the wave. Secondly, notice that
so that the velocity associated with the moving frame is given by
For our choice of $\varepsilon$ and $v$, it therefore follows that the choice of a small $\varepsilon$ corresponds to a propagation speed that is marginally smaller than the acoustic speed $M_{a}=1$. Thirdly, for a fixed time $t$ one obtains that
where $\gamma =-(1-\varepsilon v)t$ is a constant. This shows that a small choice of $\varepsilon$ corresponds to a large stretching in the $x$-coordinate (relative to the $\eta$ coordinate).
To illustrate these effects, we consider solutions of $\phi (x,t)$ in the original coordinates for $v=1$, $\varepsilon =0.001$ and $\alpha =1$. Here, the propagation speed is given by
and the wavelength is given by
The period, in other words the smallest $T>0$ satisfying $\phi (x,t+T)=\phi (x,t)$, is given by
In figure 6 the solution is shown across five wavelengths $0\leq x\leq 5L$ and over two periods $0\leq t\leq 2T$. This figure clearly illustrates that the cnoidal wave solutions correspond to the small-amplitude long wavelength limit, and are associated with large periods. Indeed, the $x$ coordinate is normalized with respect to the Debye length $\lambda _{D}$, and the time coordinate $t$ is normalized with respect to the inverse plasma frequency $\omega _{p}^{-1}$. From the solution, it therefore follows that the wavelength stretches over more than 100 Debye lengths, while the period exceeds the inverse plasma frequency by a factor of more than 100. To summarize, the cnoidal wave solutions correspond to the low-amplitude, long wavelength and low frequency oscillations.
6.3 Comparison between linear and cnoidal wave solutions
Another advantage of transforming the cnoidal wave solutions back to the original coordinates is that it provides an easy comparison with the linear waves associated with the dispersion relation. Indeed, the linear waves are given by
where the free parameter $a\ll 1$ must be sufficiently small to avoid nonlinear effects, while $\omega$ satisfies the dispersion relation (2.7). In order to compare the linear and cnoidal waves, it should be noted that the propagation speed $M$ of the cnoidal wave must be equal to the phase speed, giving
In terms of amplitude, one can compare the two solution sets by setting the amplitude of the linear wave equal to the maximum value of the cnoidal wave, that is, by choosing
As a concrete example, we once again consider the cnoidal wave solutions associated with $v=1$ and $\varepsilon =0.001$. The solutions for $\alpha =0.01$, $\alpha =0.1$, $\alpha =1$ and $\alpha =5$ at $t=0$ are shown in figure 7(a–d), respectively. In these panels, the blue curves correspond to the cnoidal wave, whereas the red curve corresponds to the linear wave. In figure 7(a) we see that these two curves are indistinguishable. This shows that, for sufficiently small $\alpha$, associated with small amplitudes, the linear approximation provides a very accurate approximation of the cnoidal wave solution. Similarly, figure 7(b) shows that this strong agreement is maintained for $\alpha =0.1$. While there is a discrepancy around the minima, the remainder of the solutions are very similar in terms of shape and wavelength. This agreement, however, wanes for larger choices of $\alpha$, as is clear from figure 7(c), corresponding to $\alpha =1$. Here, we see that the linear solution overestimates the magnitude of the local minima. In addition, the wavelength of the linear solution is slightly larger than that of the cnoidal wave, as is clearly seen in the interval $250\leq t\leq 300$. For larger choices of $\alpha$, the contrast between the two sets of solutions grows even starker, as is clear from figure 7(d) showing the solutions for $\alpha =5$. Clearly, the magnitude of the local minima is greatly overestimated by the linear wave solution, along with a larger overestimation of the wavelength.
Based on these results, we can conclude that the cnoidal wave solutions are consistent with linear theory in the small-amplitude limit $\alpha \ll 1$. As the amplitudes increase, however, the nonlinear effects become more important. The result of this is that the two sets of solutions diverge more and more with increasing amplitude. In particular, the magnitude of the local minimum of the cnoidal waves is overestimated by the linear waves, while the wavelength of the cnoidal waves is also overestimated by the linear waves.
7 Conclusions
In recent years, a new approach to studying cnoidal wave solutions was introduced. In this approach, one uses boundary conditions similar to those used to obtain soliton solutions. In this paper, we demonstrate some of the flawed reasoning behind this implementation, and show that the resulting solutions violate the conservation of number density.
To address this, we propose a new formalism based on a set of initial conditions instead of boundary conditions. In this approach, solutions are constructed in a way that ensures the conservation of number density, as was done in the early works on this topic. While the early works used a power series expansion to simplify the solutions, our approach requires no approximations. One advantage is that one can study the small-amplitude limit where $m\rightarrow 0$. We use this to show that the solutions are consistent with linear theory in this limit.
It is important to note that the scope of this paper is limited to subacoustic cnoidal wave solutions. For this case, the cnoidal wave solutions do not approach the soliton limit $m\rightarrow 1$ in the large-amplitude limit. For this limit, one must consider superacoustic cnoidal waves, a topic for future investigation.
To summarize, this paper provides a new approach that allows the derivation of the cnoidal wave solution in a much more detailed form than previously reported. This improved level of detail makes it easier to construct the solutions, and also ensures that a fundamental quantity is conserved, namely the number density.
Acknowledgements
Editor Thierry Passot thanks the referees for their advice in evaluating this article.
Funding
This work is based on the research supported wholly/in part by the National Research Foundation of South Africa (Grant numbers 145712).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Conservation of number density
In order to see where the problems lie with many papers on cnoidal plasma waves, we start from the description of nonlinear waves on the surface of water. These are the oldest well-studied nonlinear waves, and provide the imagery that we often keep at the back of our mind when we think of waves and their properties. Totally leaving aside the nonlinear solitary waves, we now refer to linear or nonlinear periodic waves, often described by cnoidal functions. To this end, consider a narrow water canal where one may assume that the water height is uniform across the width of the canal (say, the $y$-coordinate). Therefore, the water volume depends only on the height ($z$-coordinate) of the water along the length ($x$-coordinate) of the canal. If the water in the canal is undisturbed, i.e. in equilibrium, then the volume of water along one unit in the $x$-direction is given by
where $w$ is the width of the canal, and $h$ is the equilibrium height (i.e. depth) of the water canal.
When we look for periodic wave solutions, we look for solutions that conserve the mass of the water in the equilibrium state. This conservation ensures that these waves can be generated through some perturbation (wind, for example) of the undisturbed water mass. Since the water mass is directly proportional to the water volume, one can ensure conservation of mass by ensuring conservation of volume.
Now, suppose that a periodic solution $z=f(x,t)$ exists with period $L$, that is, $f(x+L,t)=f(x,t)$ for all $x$ and $t$. The average mass per unit of $x$ is then proportional to the average volume, given by
To ensure conservation of mass, one must have that $V_e=\bar {V}_p(t)$ for all $t$. This condition can be reduced to the following condition:
that is, the normalized average height of the periodic solution must be equal to the equilibrium height.
Suppose that we would have solutions with graphs as in figure 8, we first address the positive solutions, in terms of the elevation above the equilibrium level surface. Here, the black line corresponds to the equilibrium height of the water, and the blue curve shows the solution. Since the wave remains above the equilibrium water height, it is clear that the average height will exceed the equilibrium height, so that the condition (A3) is not satisfied. This shows that these positive humps cannot be generated by a small or large disturbance of the water at rest, because it would need the addition of an infinite layer of water (local height times a length from minus to plus infinity). A similar reasoning shows that the negative solutions will have an average height below the equilibrium height. This would imply that some of the water must be removed before the solution can be generated, thus violating the conservation of mass.
Let us now turn our attention to electrostatic waves in plasmas. Unlike the water medium, the volume is unaffected by the wave. Instead, the number density fluctuates in a periodic fashion. For the plasma, the ion mass within a given volume is directly proportional to the total number of ions within this volume. The number of ions in a given volume can be determined by integrating the ion number density over that volume.
To derive a condition that ensures conservation of ion mass, consider a square of length $\ell$ that is perpendicular to the direction of the flow, i.e. the direction of the magnetic field. We assume that the number density remains constant within the square, and only varies along the direction of the magnetic field $x$. In normalized terms, the total number of ions within one unit of $x$ (i.e. one Debye length) is given by
since the (normalized) equilibrium number density is given by $n_0=1$.
Suppose now that $n(x,t)$ is a periodic solution with period $L$ that describes the number density along the magnetic field direction $x$. To ensure the conservation of ion mass, one must then have that
This can be reduced to the condition that
This condition ensures the conservation of ion mass.