2. Dynamical equations

The nonlinear interaction of linearly polarized finite amplitude intense laser pulse with longitudinal electron density perturbations that are driven by the laser ponderomotive force in a relativistic degenerate plasma can be described by the following set of dimensionless equations which are, the EM wave equation, the electron continuity and momentum equations and the Poisson equation in the Coulomb gauge (Gratton et al. Reference Gratton, Gnavi, Galvao and Gomberoff1997; Roy et al. Reference Roy, Chatterjee and Misra2020; Misra & Chatterjee Reference Misra and Chatterjee2018):

(2.1)\begin{gather} \left(\frac{\partial^2}{\partial t^2}- c^2\frac{\partial^2}{\partial z^2}\right) {A}_x= \frac{\omega_{pe}^2m_ec^2}{n_0H_e}n_e^2A_x, \end{gather}

(2.2)\begin{gather}\frac{\partial \left(\gamma_en_e\right)}{\partial t} +\frac{\partial}{\partial z} (\gamma_en_e {v}_{ez})=0, \end{gather}

(2.3)\begin{gather}\frac{\textrm{d}}{\textrm{d}t}\left(\frac{\gamma_eH_ev_{ez}}{n_ec^2}\right)=e\left(\frac{\partial\phi}{\partial z}-\frac{1}{2}\frac{en_e}{\gamma_eH_e}\frac{\partial A_x^2}{\partial z}\right)-\frac{1}{\gamma_en_e}\frac{\partial P_e}{\partial z}, \end{gather}

(2.4)\begin{gather}\frac{\partial^2 \phi}{\partial z^2}=4{\rm \pi} e(\gamma_en_e-n_0), \end{gather}

where $\textrm {d}/\textrm {d}t\equiv \partial /\partial t+ {\boldsymbol {v}}_j\boldsymbol {\cdot } \boldsymbol {\nabla }$, and $e$, $n_e$ (with its equilibrium value $n_0$ in laboratory frame), $m_e$ and $v_{ez}$ are, respectively, the charge, the number density, the mass and $z$-component of the velocity of electrons. Since ions form the neutralizing background, $\gamma _i=1$ and so, $\gamma _i n_i=n_0$, the equilibrium number density of electrons and ions. Thus, the electron number density in the laboratory frame with its equilibrium and perturbation parts may be written as $n_L\equiv \gamma _en_e=n_0+n'_L$. Furthermore, $c$ is the speed of light in vacuum, $\omega _{pe}$ is the electron plasma frequency, $A_x$ is the $x$-component of the vector potential, $\phi$ is the electrostatic potential, $P_e$ is the electron degeneracy pressure at zero temperature, given by (Chandrasekhar Reference Chandrasekhar1935)

(2.5)\begin{equation} \left. \begin{gathered} \left(P_e,\mathcal{E}_e\right)=\frac{m_e^4c^5}{3{\rm \pi}^2\hbar^3}[f(R),R^3(1+R^2)^{1/2}-f(R)],\\ f(R)=\tfrac18[ R(2R^2-3)(1+R^2)^{1/2}+3\sinh^{{-}1}R], \end{gathered}\right\}\end{equation}

where $\hbar =h/2{\rm \pi}$ is the reduced Planck's constant, $R=p_{Fe}/m_ec=\hbar (3{\rm \pi} ^2n_e)^{1/3}/m_ec$ is the dimensionless degeneracy parameter and $H_e\equiv P_e+\mathcal {E}_e=n_em_ec^2\sqrt {1+R^2}$ is the enthalpy per unit volume measured in the rest frame of each element of the fluid with $\mathcal {E}_e$ denoting the total energy density, i.e. $\mathcal {E}_e=m_en_ec^2+\bar {\epsilon }_e$ and $\bar {\epsilon }_e$ the internal energy of the fluid. Also, $\gamma _e$ is the relativistic factor, given by (Misra & Chatterjee Reference Misra and Chatterjee2018)

(2.6)\begin{equation} \gamma_e=\sqrt{\frac{1+ e^2n_e^2 A^2_x/H^2_e}{1-v_{ez}^2/c^2}}. \end{equation}

For the slow-motion approximation of relativistic dynamics of electrons, i.e. $en_eA_x/H_e \sim o(\epsilon )$; $\phi$, $n_{e1}$, $v_{ez}\sim o(\epsilon ^2)$, where $0<\epsilon <1$ is a scaling parameter, (2.1) to (2.4) can be reduced to the following coupled equations (Misra & Chatterjee Reference Misra and Chatterjee2018):

(2.7)\begin{gather} \left(\frac{\partial^2}{\partial t^2}-\frac{\partial^2}{\partial z^2}+1\right)A+(1-\delta_{e})(N-\alpha A^2)A=0, \end{gather}

(2.8)\begin{gather}\left(\frac{\partial^2}{\partial t^2}-\delta_{e}\frac{\partial^2}{\partial z^2}+1\right)N=\frac{1}{2}(1-\delta_{e})\frac{\partial^2 A^2}{\partial z^2}, \end{gather}

where $N\equiv N_{e1}/n_0=\gamma _en_e/n_0-1$ is the dimensionless longitudinal electron density perturbation and $A\equiv \eta _e eA_x/m_ec^2$ is the dimensionless vector potential along the $x$-axis. The space and time coordinates are normalized according to $t\rightarrow t\sqrt {\eta _e}\omega _p$, $z\rightarrow z\sqrt {\eta _e}\omega _p/c$. Also, $\alpha =\eta _e^2/2$, $\delta _e=(1-\eta _e^2)/3$ and $\eta _e=1/\sqrt {1+R_0^2}$ with $R_0=\hbar (3{\rm \pi} ^2n_0)^{1/3}/mc$ denoting a measure of the strength of plasma degeneracy, i.e. $R_0\ll 1$, corresponds to the weakly relativistic degenerate plasma, whereas $R_0\gg 1$ is referred to as ultrarelativistic degenerate plasmas.

For the evolution and stability of EM solitons in relativistic degenerate plasmas, it appears much more difficult to solve the fully relativistic one-dimensional fluid model due to the generation of higher harmonics of linearly polarized incident laser pulse. However, a more convenient approach is to study the coupled equations (2.7) and (2.8) instead. So, we introduce a slowly varying complex wave envelope in that form. Here, we note that in contrast to circular polarization, the linearly polarized EM waves have odd harmonics for the vector potential $A$ and even harmonics for the electron density perturbation $N$. Thus, we have

(2.9)\begin{equation} \left.\begin{gathered} A=\tfrac{1}{2}(a\,\textrm{e}^{-\textrm{i}t}+a^{*}\,\textrm{e}^{\textrm{i}t}),\\ N=N_{0}+\tfrac{1}{2}(N_{2}\,\textrm{e}^{-\textrm{i}2t}+N_{2}^{*}\,\textrm{e}^{\textrm{i}2t}), \end{gathered}\right\} \end{equation}

where the asterisk denotes the complex conjugate of the corresponding physical quantity. Substituting the expansions in (2.9) for $A$ and $N$ into (2.8), and collecting the zeroth and second harmonic terms $({\sim }\textrm {e}^{-\textrm {i}2t})$, we obtain the following envelopes for $N_{0}$ and $N_{2}$:

(2.10)\begin{equation} \left.\begin{gathered} N_{0}=\tfrac{1}{4}(1-\delta_{e})(|a|^{2})_{zz},\\ N_{2}={-}\tfrac{1}{12}(1-\delta_{e})(a^{2})_{zz}. \end{gathered}\right\} \end{equation}

Next, substituting (2.9) into (2.7), and collecting the first harmonic terms $({\sim }\textrm {e}^{-\textrm {i}t})$, we obtain the following equation for the EM wave amplitude $a$ (for details see appendix A):

(2.11)\begin{align} &\textrm{i}\frac{\partial a}{\partial t}+\frac{1}{2}(a)_{zz}+\frac{3}{8}\alpha(1-\delta_{e})|a|^{2}a-\frac{1}{8}(1-\delta_{e})^{2}(|a|^{2})_{zz}a +\frac{1}{48}(1-\delta_{e})^{2}(a^{2})_{zz}a^{*}=0. \end{align}

Equation (2.11) has the form of a generalized nonlinear Schrödinger (GNLS) equation with local (cubic) as well as non-local (derivative) nonlinearities. It is noticed that both the cubic and nonlocal nonlinear coefficients are significantly modified by the effects of the relativistic electron degeneracy pressure. It is to be noted that in the limit of $R_0\rightarrow 0$, (2.10) and (2.11) assume the same form as (7) and (8) in Hadzievski et al. (Reference Hadzievski, Jovanovic, Skoric and Mima2002) after replacing $A$ by $a$. However, there are some disagreements with the factor $1/2$ in the expression for $N_2$ (2.10) and in the last term of (2.11). This may be due to some sign mismatch in the coefficient of $({A^2}/{2})A$ (2.7) with the similar term in Hadzievski et al. (Reference Hadzievski, Jovanovic, Skoric and Mima2002). The cubic nonlinearity should be correctly written as $({A^2}/{2})A$ or $({a^2}/{2})a$ instead of $-({a^2}/{2})a$, as in Hadzievski et al. (Reference Hadzievski, Jovanovic, Skoric and Mima2002). Similar nonlinear terms with the correct positive sign can be found in some other previous investigations, e.g. (29) of Gratton et al. (Reference Gratton, Gnavi, Galvao and Gomberoff1997).

We look for a localized stationary solution of (2.11) in the form of a moving soliton $a=\rho (\xi )\exp [\textrm {i}\theta (\xi )+\textrm {i}\lambda ^{2}t]$ where $\xi =z-v_0t$ and $v_0$ is the soliton velocity in the moving frame of reference. Thus, substituting this solution into (2.11), we obtain the following equations for the soliton phase and the amplitude:

(2.12)\begin{gather} \theta_{\xi\xi}\rho\left[1+\frac{1}{12}(1-\delta_{e})^{2}\rho^2\right] +\theta_{\xi}\rho_{\xi}\left[2+\frac{1}{3}(1-\delta_{e})^{2}\rho^{2}\right]-2v_0\rho_{\xi}=0, \end{gather}

(2.13)\begin{gather}\rho_{\xi\xi}-\frac{5}{12}(1-\delta_{e})^{2}\frac{\rho}{\zeta}\rho_{\xi}^{2}=\frac{\rho}{\zeta}\left(2\lambda^{2}-2v_0\theta_{\xi}+\theta_{\xi}^2\right) +\frac{\rho^3}{\zeta}\left[\frac{1}{6}(1-\delta_{e})^{2}\theta_{\xi}^{2}-\frac{3}{4}\alpha(1-\delta_{e})\right]\! =0, \end{gather}

where $\zeta (\rho )=1-(5/12)(1-\delta _{e})^{2}\rho ^{2}$. Using the same boundary conditions, namely $\rho$, $\rho _{\xi }$, $\rho _{\xi \xi }\rightarrow 0$ as $\xi \rightarrow \pm \infty$, the integration of (2.12) approximately yields $\theta (\xi )_{\xi }=v_0$, while that of (2.13) gives

(2.14)\begin{equation} \rho_{\xi}^{2}=\frac{\rho^2}{\zeta}\left(2\lambda^2-v_0^2-\left[\frac{3}{8}\alpha(1-\delta_{e})-\frac{1}{12}(1-\delta_{e})^{2}v_0^{2}\right]\rho^2\right). \end{equation}

Further integration of (2.14) yields a soliton solution in the following implicit form:

(2.15)\begin{align} \pm \xi &= \sqrt{\frac{10(1-\delta_e)}{\varDelta}}\ln\frac{\sqrt{\zeta}+\sqrt{\zeta-\zeta_0}}{\sqrt{\left\vert\zeta_0\right\vert}} -\frac{1}{2\sqrt{2\lambda^2-v_0^2}}\ln\frac{\sqrt{1-\rho^2/\rho_0^2}+\sqrt{\zeta}}{\left\vert\sqrt{1-\rho^2/\rho_0^2}-\sqrt{\zeta}\right\vert}, \end{align}

where $\zeta _0(\rho _0)=1-(5/12)(1-\delta _{e})^{2}\rho _0^{2}$, $\varDelta =9\alpha -2(1-\delta _{e})v_0^{2}$ and

(2.16)\begin{equation} \rho_{0}^{2}=\frac{24(2\lambda^{2}-v_0^{2})}{\varDelta(1-\delta_{e})} \end{equation}

is the squared amplitude (maximum) of the linearly polarized moving EM soliton with the reduced eigenfrequency $\varLambda \equiv \lambda ^2-v_0^2$. In particular, for $v_0=0$, one can obtain the standing soliton solution of the GNLS (2.11). In general, (2.15) describes a two-parameter family of solutions of (2.11) which are symmetric about the origin, and while combined together generates the form of a soliton. As an illustration, we plot the soliton solution (2.15) as shown in figure 1 for different values of the parameters $v_0$, $\lambda$ and $R_0$. It is seen that for a fixed value of $R_0$ and $\lambda$ but with an increasing value of the soliton velocity, the amplitude decreases while the soliton profile (width) broadens. An opposite feature is observed with increasing value of $\lambda$, i.e. the soliton amplitude increases but the width decreases. However, when the relativistic degeneracy effect is more pronounced, an enhancement of both the amplitude and width of the soliton is seen to occur. In order to justify the physics behind it we consider the small amplitude limit of solitons in which the non-local (ponderomotive) nonlinear terms may be neglected so that one recovers the cubic nonlinear Schrödinger equation from (2.11) as

(2.17)\begin{equation} \textrm{i}\frac{\partial a}{\partial t}+P\frac {\partial^2 a} {\partial z^2}+Q|a|^2a=0, \end{equation}

where $P=1/2$ and $Q=(3/8)\alpha (1-\delta _e)$. Equation (2.17) has a travelling wave solution of the form $a\sim \sqrt {a_0}\,\text {sech} [(z-v_0t)/w]$, where $wa_0=\sqrt {2|P/Q|}$. Physically, the relativistic degeneracy pressure of electrons provides the wave dispersion (which increases with increasing values of $R_0$) quite distinct from the dispersion due to separation of charged particles. So, depending on the dispersion and nonlinear effects, the soliton amplitude and width can increase or decrease. Since $P$ appears as a constant due to its normalization, the wave amplitude/width can increase or decrease according to when the value of the nonlinear coefficient $Q$ decreases or increases. It is noted that the values of $Q$ decrease with increasing values of $R_0$. Thus, it may be likely that the soliton amplitude and width can increase with increasing values of $R_0$. However, when the wave amplitude is not small, the ponderomotive nonlinear effects can intervene with the dynamics and compete with the cubic nonlinearity which may result in an increase of the wave amplitude and/or width, or some other nonlinear phenomena including collapse.

Figure 1.

A profile of the EM soliton, given by (2.15), is shown for different values of the relativistic degeneracy parameter R ₀, the soliton eigenfrequency λ and the soliton velocity v ₀, as in the legend.

3. Stability analysis

In order to examine the stability of the moving EM soliton, we follow the Vakhitov–Kolokolov stability criterion (Vakhitov & Kolokolov Reference Vakhitov and Kolokolov1973). According to the criterion, solitons are stable under the longitudinal perturbations if

(3.1)\begin{equation} \frac{\textrm{d}P_{0}}{\textrm{d}\lambda^{2}}>0, \end{equation}

where $P_{0}$ is the soliton photon number defined by

(3.2)\begin{equation} P=\int |a|^{2} \,\textrm{d}z. \end{equation}

An expression for $P_{0}(\lambda )$ can be obtained for the soliton solution (2.15) as

(3.3)\begin{equation} P_{0}(\lambda)=\frac{\sqrt{6}}{\sqrt{\varDelta(1-\delta_{e}) }} \left( \rho_{0}+\sqrt{\frac{3}{5}}\frac{\zeta_0}{1-\delta_{e}} \ln\frac{1+ \sqrt{1-\zeta_0}}{\left\vert 1- \sqrt{1-\zeta_0}\right\vert}\right). \end{equation}

So, according to the condition (3.1), the moving EM soliton (2.15) turns out to be stable in the region $\lambda <\lambda _{s}$, where $\lambda _s$ is some instability threshold value of $\lambda$ at which $P_0$ reaches a local maximum and above which $\textrm {d}P_0/\textrm {d}\lambda ^2<0$. An analytic expression of $\lambda _{s}$ cannot be determined in its explicit form. However, we try to find its values numerically for different values of the degeneracy parameter $R_0$ and the soliton velocity $v_0$.

The profiles of the curves of $P_0(\lambda )$, as shown in figure 2, predict the existence of stable $(\lambda <\lambda _s)$ and unstable regions $(\lambda >\lambda _s)$. It is found that an increase of the soliton velocity $v_0$ shifts the instability threshold $\lambda _s$ towards larger values than that for the standing soliton with $v_0=0$ (see the solid and dashed curves). However, $\lambda _s$ can shift towards its lower values if the degeneracy parameter $R_0$ is increased (see the dotted and dash-dotted curves). It follows that the stability domain $\lambda <\lambda _{s}$ of EM solitons may be significantly reduced in the regime of strongly or ultrarelativistic degenerate plasmas $(R_0\gg 1)$. We mention that the Vakhitov–Kolokolov criterion predicts only a linear stability of EM solitons involving the exponentially growing or decaying modes. So, it may not make predictions about the subsequent nonlinear evolution of unstable EM envelopes or about the stability of localized structures with arbitrary profiles. In general, the GNLS equation (2.11) can admit, apart from the stationary soliton solution, the soliton collapse (Litvak & Sergeev Reference Litvak and Sergeev1978) and long-lived relaxation oscillations around the stable soliton amplitude due to cubic nonlinearity, and perhaps some other dynamical states due to the presence of nonlocal nonlinearities (Hadzievski et al. Reference Hadzievski, Jovanovic, Skoric and Mima2002; Pelinovsky, Afanasijev & Kivshar Reference Pelinovsky, Afanasijev and Kivshar1996). An alternative stability criterion for the solitons can also be formulated in terms of the Hamiltonian and photon number interrelation (see, e.g. Mancic et al. Reference Mancic, Hadzievski and Skoric2006). We, however, skip this analysis here and instead look for some other regimes for the existence and stability of solitons by revisiting the Vakhitov–Kolokolov stability criterion and the constraints on the parameters $\lambda$ and $v_0$ in the analytical soliton solution (2.15). The limits of $\lambda$ and $v_0$ are given by

(3.4a,b)\begin{equation} 2\lambda^2-v_0^2>0,\quad~ \zeta\equiv1-5\rho^{2}\left(1-\delta_{e}\right)^{2}/12>0. \end{equation}

One more limitation on the parameter $v_0$ is given by

(3.5)\begin{equation} \varDelta\equiv 9\alpha-2(1-\delta_e)v_0^2 > 0,\quad \textrm{i.e.},\ v_0 < v_s\equiv \sqrt{9\alpha/2(1-\delta_e)}. \end{equation}

An estimation reveals that for $0<(\alpha ,\delta _e,v_0) < 1$, the condition $\varDelta >0$ holds for $0\leq R_0\lesssim 1.36$ and $\varDelta < 0$ for $R > 1.36$.

Figure 2.

Photon number $P_{0}(\lambda )$, given by (3.3), is plotted against the eigenfrequency λ to show the stable $(\lambda < \lambda _s; \textrm {d}P_0/\textrm {d}\lambda ^2 > 0)$ and unstable $(\lambda > \lambda _s; \textrm {d}P_0/\textrm {d}\lambda ^2 < 0)$ regions for different values of the relativistic degeneracy parameter R ₀ and the soliton velocity $v_0$, as in the legend. The straight lines are drawn to indicate the corresponding threshold values $\lambda _s$ of $\lambda$ at which $P_0$ reaches a maximum value.

In what follows, considering all the stability conditions, namely (3.1), (3.4a,b) and (3.5) imposed on the parameters $\lambda$ and $v_0$, we define the regions of the soliton existence and stability in $(v_0,\lambda )$-plane as shown in figure 3. We find that, above the dashed curve of $\zeta \equiv 1-5\rho ^{2}(1-\delta _{e})^{2})/12=0$, no analytic soliton solution exists as the conditions are not satisfied there. However, one can obtain numerically a soliton solution of (2.11) there. Such a discrepancy in this particular region occurs due to an initial small phase $(\theta (\xi ))$ difference between these approximate analytical and numerical soliton solutions. Furthermore, below the solid curve of $2\lambda ^{2}-v^{2}=0$, no localized solution (analytical or numerical) can exist. The existence region of solitons is separated by a dotted curve. Between the dashed and the dotted curves is the region where $dP_{0}/d\lambda ^{2}<0$, but both $2\lambda ^{2}-v^{2}$ and $\zeta > 0$, and so in this region only an unstable soliton solution exists. However, the stable soliton region is in-between the solid and dotted curves where all the conditions for soliton stability are satisfied. From figure 3(a–c) it is noticed that, as we gradually enter from the regimes of weak, moderate to strong relativistic degenerate plasmas (by increasing the value of the degeneracy parameter $R_{0}$) both the stable and unstable regions get significantly compressed. Also, a stable region shifts to an unstable one as $R_0$ increases. Furthermore, figure 3(c) shows that a threshold value, $v_s$, of $v_0$ appears when $R_0>1.36$. In fact, for $R_0\gg 1$, no stable or unstable region may be found in the $(v_0,\lambda )$-plane. This may be the limitation of the linear stability analysis, which may not correctly predict the existence and the stability regions of moving solitons, especially in the regime of ultrarelativistic degeneracy.

Figure 3.

The regions for the existence of EM solitons and their stability are shown in the $(v_0,\lambda )$-plane for different values of the relativistic degeneracy parameter $R_{0}$: (a) $R_{0}=0.02$; (b) $R_{0}=0.8$; and (c) $R_0=2$. The solid, dashed and dotted lines, respectively, represent the contour plots of $2\lambda ^{2}-v^{2}=0$, $\zeta \equiv 1-(5/12)\rho ^{2}(1-\delta _{e})^{2}=0$ and $\textrm {d}P_0/\textrm {d}\lambda ^2=0$. In panel (c), the dash-dotted line represents a threshold value of $v_0$: $v_s=0.78$ above which no soliton solution exists as $\varDelta \equiv 9\alpha -2(1-\delta _e)v_0^2 < 0$ there. Note that the value of $v_s$ decreases with increasing values of $R_0>1.36$. No such threshold exists in the region $0\leq v_0\leq 1$, $0\leq \lambda \leq 1$ when $R_0$ lies in $0\leq R_0\lesssim 1.36$.

4. Simulation results

So far we have obtained the possible regimes for the existence and stability of EM solitons based on an analytic stationary soliton solution of the GNLS equation (2.11). Next, we examine these regimes by a direct numerical simulation of (2.11). To this end, we use the Runge–Kutta scheme with a time step $\textrm {d}t=0.001$ and with an initial condition in the form of a soliton: $a(\xi ,0)\sim a_0\,\text {sech}^2(\xi /5)\exp (-\textrm {i} v_0\xi )$. The soliton evolution after time $t=300$ is shown in figures 4 and 5 for different values of $v_0$, $\lambda$ and $R_0$ relevant for stable and unstable regions (cf. figure 3) as predicted in § 3. Here, we note that the initial condition may not be a solution of the GNLS equation (2.11). However, as time goes on, the initial pulse radiates and the nonlinear and dispersion effects intervene or it to evolve as stable or unstable solitons. Numerical simulation reveals that initially launched soliton (2.15) with parameters $v_0$, $\lambda$ and $R_0$ in the stable region remains stable for a long time. However, the stable or unstable behaviours may be changed if one considers the parameter values slightly above or below the stable or unstable region.

Figure 4.

Spatio-temporal evolution of the GNLS equation (2.11) in different dynamical regimes for a fixed value of the degeneracy parameter $R_0=0.02$. Panels (a,c,e) are the surface plots of soliton solutions, while panels (b,d,f) represent the corresponding contour plots. Panels (a,b) show an evolution of standing soliton in the stable region with soliton velocity $v_0=0$, the soliton eigenfrequency $\lambda =0.207$, the amplitude $a_0=0.68$ and the soliton photon number $P_0=1.46$; panels (c,d) a moving soliton in the stable region with an increasing velocity $v_0=0.2$, the same $\lambda =0.207$ but different $a_0=0.67$ and $P_0=1.12$; panels (e,f) a moving soliton in the stable region with the same $v_0=0.2$ as in panels (c,d) but different $a_0=0.72$, $\lambda =0.26$ and $P_0=1.55$.

Figure 5.

Spatio-temporal evolution of the GNLS equation (2.11) for different values of the degeneracy parameter $R_0$: $R_0=0.02$ in the stable region, panels (a,b); $R_0=0.8$ in the stable region, panels (c,d); and $R_0=2$ in the unstable region panels (e,f). The other parameter values are fixed at $v_0=0.2$ and $\lambda =0.23$.

Figure 4 shows the evolution of EM solitons in the weakly relativistic regime $(R_0=0.02)$. It is found that in the stable region, the standing soliton $(v_0=0)$ with amplitude $a_0=0.68$, $\lambda =0.207$ and photon number $P=1.46$ oscillates around the stable equilibrium with a frequency close to the plasma oscillation frequency and so the soliton propagation remains stable for a long time (figure 4a,b). However, the moving soliton with the same $\lambda =0.207$, but with an increased velocity $v_0=0.2$ and reduced amplitude $a_0=0.67$ in the stable region, travels towards the downstream and exhibits a decay of its amplitude (figure 4c,d). Physically, for a fixed $\lambda$, as the soliton velocity increases, its frequency of oscillation decreases. This results in a diminution of the photon number from $P=1.46$ to $P=1.12$, and thereby leading to a decay of the soliton amplitude. However, as time goes on it relaxes towards a corresponding stable soliton. On the other hand, retaining the soliton speed at $v_0=0.2$, but increasing its amplitude to $a_0=0.72$ and frequency to $\lambda =0.26$, we find that both the soliton eigenfrequency and the photon number $(P_0=1.55)$ increase. As a result, though the soliton evolves with long-lived oscillating behaviours of the breather type, its amplitude grows and it may be prone to instability (figure 4e,f). A slight deviation from the stable equilibrium with an increase of $\lambda$ or the initial perturbation with an increased amplitude $a_0$ for a fixed soliton velocity can lead to an aperiodic growth of the amplitude and thereby the onset of collapse.

In order to examine the effects of the relativistic degeneracy on the soliton dynamics, we consider three different values of $R_0$, namely $R_0=0.02$, $R_0=0.8$ and $R_0=2$ to define, respectively, the weak, moderate and strong degeneracy of electrons. With reference to figure 3 (which predicts that the soliton stable region can shift to an unstable one by an increasing value of $R_0$) we see from figure 5 that for some fixed values of $\lambda =0.23$ and $v_0=0.2$, as $R_0$ increases the soliton amplitude grows, it loses its stability and eventually collapses. Figures 5(a) and 5(b) show that in the weakly relativistic regime of electron degeneracy, a moving soliton with an amplitude $a_0=0.6$ in the stable region travels towards downstream preserving its profile, i.e. the soliton remains stable for a longer time. The solitons also remain stable even in the limit of $R_0\rightarrow 0$, i.e. when there is no electron degeneracy. If the initial soliton profile is considered with an increased amplitude with $a_0=0.82$ in the regime of moderate degeneracy with $R_0=0.8$, the soliton amplitude grows but it evolves about the stable equilibrium until it remains in the stable region (figure 5c,d). However, as the parameter $R_0$ is further increased to $R_0=2$, the moving soliton with amplitude $a_0=1.6$ falls in the unstable region. In this situation, solitons cannot travel undistorted with a constant velocity. Its amplitude aperiodically grows and eventually collapses (figure 5e,f).

We note that when the electron degeneracy effect is small, i.e. $R_0\ll 1$, one finds $\eta _e\sim 1$ and $\delta _e\ll 1$ for which the coefficients of the GNLS equation (2.11) appear as constants. In this case, one recovers similar results as in Hadzievski et al. (Reference Hadzievski, Jovanovic, Skoric and Mima2002) with no degeneracy of electrons. However, for moderate or large values of $R_0$, as the parameter increases, the magnitudes of both the cubic and nonlocal nonlinearities tend to decrease, which results in an enhancement of the soliton amplitude and/or width. This is expected, since in absence of the non-local terms, the GNLS equation (2.11) admits a soliton solution with amplitude/width proportional to $\sqrt {|P/Q|}$, where $P$ corresponds to the group velocity dispersion and $Q$ the cubic nonlinear effects. However, in the presence of the non-local nonlinearities, not only do the solitons grow with higher values of $R_0$, a rapid aperiodic growth in amplitude creates a collapsing soliton until the amplitude reaches a critical value, and we can no longer observe the dynamic behaviours. Figure 5(e) shows that when the values of $v_0$ and $\lambda$ fall slightly outside the stable region (cf. figure 3c), the soliton amplitude grows to become unstable, and eventually collapses. From figure 5(f), it may be estimated that the soliton collapse starts to occur after time $t\sim 160$, and within the time interval $0 < t < 160$ the soliton amplitude remains of the order of unity. Thus, our analytical predictions for the existence and stability of EM solitons agree with the numerical results, that an increase of the soliton velocity shifts the instability threshold towards larger values, and the strong degeneracy effects can shift the steady sate dynamics of solitons to an unstable one which results into wave collapse.