# Stabilization and blow-up in the relativistic model of cold collisional
plasma

Olga S. Rozanova  Olga S. Rozanova
Mathematics and Mechanics Department
Lomonosov Moscow State University
Leninskie Gory Moscow 119991
Russian Federation<EMAIL_ADDRESS>and  Eugeniy V. Chizhonkov
Eugeniy V. Chizhonkov
Mathematics and Mechanics Department
Lomonosov Moscow State University
Leninskie Gory Moscow 119991
Russian Federation<EMAIL_ADDRESS>
###### Abstract.

We study the influence of the factor of electron-ion collisions on the
solution of the Cauchy problem in the one-dimensional relativistic model of
cold plasma and show that, depending on their intensity and initial data, two
scenarios are possible: either the solution remains smooth and stabilizes to a
stationary state, or during a finite time the oscillations blowup. In contrast
to the nonrelativistic model, when exact conditions can be obtained separating
the two behaviors, in a much more complicated relativistic situation, it turns
out to be possible to analytically estimate from below the time during which
the existence of a smooth solution and the guaranteed number of oscillations
during this time. In addition, we show that in contrast to the relativistic
case without taking into account collisions, when oscillations corresponding
to arbitrarily small deviations from the zero equilibrium position blow up,
the presence of electron collisions can suppress the blow-up of sufficiently
small oscillations. Further, based on the analysis of characteristics, a
numerical algorithm is constructed, the order of accuracy of which is
determined only by the smoothness of the initial data. Numerical experiments
are presented to illustrate the theoretical results. The initial conditions
are chosen as reasonably as possible from the point of view of full-scale
physical experiments.

###### Key words and phrases:

Quasilinear hyperbolic system, plasma oscillations, breaking effect, loss of
smoothness

###### 1991 Mathematics Subject Classification:

Primary 35Q60; Secondary 35L60, 35L67, 34M10

Partially supported by the Moscow Center for Fundamental and Applied
Mathematics.

##  1\.  Introduction

The hydrodynamic model of cold plasma is one of the most commonly used and has
long been used in physics [10], [6]. Depending on what effects we want to
study, the influence of ion motion can be neglected or taken into account to
one degree or another. One of the specific features of the cold plasma model
is that it has physical meaning only for smooth solutions. It is believed that
when the singularities of the solution are formed (here these are delta-shaped
singularities in the density component), energy is released, which leads to
the ”heating” of the plasma and the need to replace the model. Indeed, the
formation of a singularity within the framework of the Lagrangian approach
corresponds to the intersection of adjacent trajectories. When two different
particles occupy the same position in space and time, then further tracking
their motion requires the use of more complex models than classical
electrodynamics, since the infinite concentration of electric charge requires
a special interpretation.

For the nonrelativistic approximation, the problem of the criterion for the
formation of a singularity in the solution of the collisional cold plasma
equations in terms of the initial data is completely solved in [15] (the same
problem in terms of the Euler-Poisson equations is considered in [8]). It
turned out that for a constant collision coefficient, one can strictly divide
the initial data into those that in a finite time lead to the formation of a
singularity and those that correspond to a solution that is globally smooth in
time. In the latter case, the solution asymptotically tends to a stationary
state as $t\to\infty$. In addition, there is a threshold value of the
collision coefficient $\nu>0$, after which the oscillatory nature of the
damping is replaced by a monotonic one (this value is $\nu=2$, while the
physically natural values of $\nu$ are much smaller). However, for an
arbitrarily large fixed value of the collision coefficient, there are initial
data at which a singularity is formed. From the point of view of physical
applications, results of this kind are of limited value, since nonrelativistic
models are regarded as a strong simplification.

It is believed that relativistic plasma oscillations are always blow up in the
general case. As shown in [13], for the collisionless case, this is to a
certain extent true, since nonlinear resonance occurs in the solution (except
for solutions that correspond to special initial data). In particular, even an
arbitrarily small perturbation of the trivial equilibrium position generally
leads to blowing up.

A natural question arises: will the presence of collisions be able to extend
the lifetime of the solution, or, perhaps, completely suppress the blow-up? In
the present work, we obtain a positive answer to this question, extending the
results of [13]. We show that the lifetime of a smooth solution can be
estimated from below, including in terms of the number of oscillations, and
this time increases as $\nu$ grows. Moreover, for any fixed initial data,
choosing a sufficiently large $\nu>0$ allows one to obtain a global in time
smooth solution stabilizing to a trivial constant state. Further, we will show
that for each arbitrarily small $\nu>0$, there exists such a small
perturbation of the stationary state that the corresponding solution remains
globally smooth.

It should be noted that the influence of collisions on plasma oscillations has
already been analyzed by other authors (see, for example, [11],[16]). However,
in these works, the plasma resistance was taken into account simultaneously
with the viscosity. As a result, the blow-up effect fell outside the attention
of researchers. For the first time, it was possible to trace the influence of
electron collisions on the breaking of plane plasma oscillations in [9].
However, there the analysis was carried out by asymptotic and numerical
methods, therefore, it was not possible to obtain exact mathematical
formulations about the conditions of blowup (or, conversely, non-blowup).

The article is structured as follows. Section 2 contains a detailed statement
of the problem in Eulerian variables and a reduction to a quasilinear system
of hyperbolic equations. Section 3 contains the analytical results announced
above. Further, the method for proving analytical statements is transformed
into a numerical algorithm, the order of accuracy of which is determined only
by the smoothness of the solution. The paper also contains numerical
experiments that illustrate the theoretical results. The initial conditions
are chosen as reasonable as possible from the point of view of full-scale
physical experiments. For the sake of completeness, the physical
substantiation of the intensity of electron collisions is given. In the
conclusion, the results of the research are systematized.

The main goal of this work is to substantiate and demonstrate the fact that an
increase in the parameter $\nu$ leads to an expansion of the range of initial
data for which there is a smooth global solution in time.

##  2\.  Formulation of the problem

We will consider plasma as a compressible relativistic electron liquid,
neglecting recombination effects and ion motion. In vector form, the system of
hydrodynamic equations describing it, together with Maxwell’s equations, is

(1) $\begin{array}[]{l}\displaystyle{\partial n\over\partial
t}+\mbox{div}\,(n{\bf v})=0\,,\quad\displaystyle{\partial{\bf p}\over\partial
t}+\left({\bf v}\cdot\nabla\right){\bf p}=e\,\left({\bf
E}+\displaystyle{1\over c}\left[{\bf v}\times{\bf
B}\right]\right)-\nu_{ei}{\bf p},\vspace{0.5em}\\\
\gamma=\sqrt{1+\displaystyle{|{\bf p}|^{2}\over m^{2}c^{2}}}\,,\quad{\bf
v}=\displaystyle{{\bf p}\over m\gamma}\,,\vspace{0.5em}\\\
\displaystyle{1\over c}\frac{\partial{\bf E}}{\partial
t}=-\displaystyle{4\pi\over c}en{\bf v}+{\rm rot}\,{\bf
B}\,,\quad\displaystyle{1\over c}\frac{\partial{\bf B}}{\partial t}=-{\rm
rot}\,{\bf E}\,,\quad\mbox{div}\,{\bf B}=0\,,\end{array}$

where $e,m$ are the charge and mass of an electron (here the electron charge
has a negative sign:$e<0$), $c$ is the speed of light; $n,{\bf p},{\bf v}$ are
the density, momentum and speed of electrons; $\gamma$ is the Lorentz factor;
${\bf E},{\bf B}$ are vectors of electric and magnetic fields.

The system of equations (1) is one of the simplest models of plasma, which is
often called the equations of hydrodynamics of ”cold” plasma; it is well known
and described in sufficient detail in textbooks and monographs (see, for
example, [3] \- [2]).

The equation for the momentum contains the term $\nu_{ei}{\bf p}$, which
describes electron-ion collisions. Taking this effect into account can be
interpreted as the force of friction between particles; in the non-
relativistic case (see, for example, [1]) they often use a formula like

$-\nu_{\alpha\beta}\left({\bf v}_{\alpha}-{\bf v}_{\beta}\right),$

where $\nu_{\alpha\beta}$ is the effective frequency of collisions of charged
particles of the type $\alpha$ with particles of the type $\beta$ when
$\alpha\neq\beta$. For stationary ions (${\bf v}_{\beta}=0$) the formula is
simplified.

In order to analyze and construct a numerical solution of plane one-
dimensional relativistic plasma oscillations taking into account collisions
the basic equations (1) can be greatly simplified.

We will denote the independent variables in the Cartesian coordinate system in
the usual way — $(x,y,z)$, and assume that

  * •

the solution is determined only $x-$ components of vector functions ${{\bf
p}},{{\bf v}},{{\bf E}};$

  * •

there is no dependence in these functions on the variables $y$ and $z$.

Then from the system (1) we get

(2) $\begin{array}[]{c}\displaystyle{\partial n\over\partial
t}+\displaystyle{\partial\over\partial
x}\left(n\,v_{x}\right)=0,\quad\displaystyle{\partial p_{x}\over\partial
t}+v_{x}\displaystyle{\partial p_{x}\over\partial
x}=e\,E_{x}-\nu_{ei}\,p_{x},\vspace{1 ex}\\\
\gamma=\sqrt{1+\displaystyle{p_{x}^{2}\over
m^{2}c^{2}}}\,,\quad{v_{x}}=\displaystyle{p_{x}\over
m\,\gamma},\quad\displaystyle{\partial E_{x}\over\partial
t}=-4\,\pi\,e\,n\,v_{x}\,.\end{array}$

We introduce the dimensionless quantities

$\rho=k_{p}x,\quad\theta=\omega_{p}t,\quad V=\displaystyle{v_{x}\over c},\quad
P=\displaystyle{p_{x}\over m\,c},\quad E=-\,\displaystyle{e\,E_{x}\over
m\,c\,\omega_{p}},\quad N=\displaystyle{n\over
n_{0}},\quad\nu=\displaystyle{\nu_{ei}\over\omega_{p}},$

where $\omega_{p}=\left(4\pi e^{2}n_{0}/m\right)^{1/2}$ is the plasma
frequency, $n_{0}$ is the value of the unperturbed electron density,
$k_{p}=\omega_{p}/c$. In the new variables (2) takes the form

(3) $\begin{array}[]{c}\displaystyle{\partial
N\over\partial\theta}+\displaystyle{\partial\over\partial\rho}\left(N\,V\right)=0,\quad\displaystyle{\partial
P\over\partial\theta}+E+V\displaystyle{\partial P\over\partial\rho}+\nu
P=0,\vspace{1.5 ex}\\\ \gamma=\sqrt{1+P^{2}},\quad
V=\displaystyle{P\over\gamma},\quad\displaystyle{\partial
E\over\partial\theta}=N\,V\,.\end{array}$

From the first and last equations (3) it follows

$\displaystyle{\partial\over\partial\theta}\left[N+\displaystyle{\partial\over\partial\rho}E\right]=0.$

This relationship is valid both in the absence of plasma oscillations and in
their presence. Therefore, under the traditional assumption of a uniform
background charge density of stationary ions, this implies a simpler
expression for the electron density $N(\rho,\theta)$:

(4) $N(\rho,\theta)=1-\displaystyle{\partial
E(\rho,\theta)\over\partial\rho}.$

Formula (4) is a special case of the Gauss theorem [6] pp.33-53, which in
differential dimensional form has the form $\mbox{div}\,{\bf
E}=4\,\pi\,e(n-n_{0}).$ Using (4) in (3), we arrive at the equations
describing plane one-dimensional relativistic plasma oscillations taking into
account collisions:

(5) $\displaystyle{\partial P\over\partial\theta}+V\,\displaystyle{\partial
P\over\partial\rho}+E+\nu P=0,\quad\displaystyle{\partial
E\over\partial\theta}+V\,\displaystyle{\partial E\over\partial\rho}-V=0,\quad
V=\displaystyle{P\over\sqrt{1+P^{2}}}\,.$

Here $\rho$ and $\theta$ are dimensionless coordinates with respect to space
and time, respectively. The variable $P$ describes the momentum of the
electrons, $V$ is the velocity of the electrons, and $E$ is a function that
characterizes the electric field. To the system (5) it is necessary to add the
equation (4), which describes the behavior of the electron density, the most
important function in the cold plasma model. Below we study in the half-plane
$\\{(\rho,\theta)\,:\,\rho\in{\mathbb{R}},\;\theta>0\\}$ a solution of the
Cauchy problem for (4), (5) with initial conditions

(6) $P(\rho,0)=P_{0}(\rho),\quad
E(\rho,0)=E_{0}(\rho),\quad\rho\in{\mathbb{R}}.$

System (5) is of the hyperbolic type. It is well known that for such systems
there exists, locally in time, a unique solution to the Cauchy problem of the
same class as the initial data [5]. Below, it is enough for us to require the
smoothness of the initial data $C^{2}(\mathbb{R})$. Since the coefficient of
electron collisions in real problems is small, we mainly consider
$\nu\in[0,2)$ (in [15] a detailed explanation of how the upper threshold value
of $\nu$ arises).

System (5) for functions $P(\rho,\theta)$ $E(\rho,\theta)$ can be considered
together with an extended system for derivatives

$Q(\rho,\theta)=\displaystyle{\partial P(\rho,\theta)\over\partial\rho},\quad
D(\rho,\theta)=\displaystyle{\partial E(\rho,\theta)\over\partial\rho},$

namely,

(7) $\displaystyle{\partial Q\over\partial\theta}+V\,\displaystyle{\partial
Q\over\partial\rho}+D+\displaystyle{Q^{2}\over(1+P^{2})^{3/2}}+\nu
Q=0,\quad\displaystyle{\partial
D\over\partial\theta}+V\,\displaystyle{\partial
D\over\partial\rho}=(1-D)\displaystyle{Q\over(1+P^{2})^{3/2}},$

with initial conditions

(8) $Q(\rho,0)=\displaystyle{d\,P_{0}(\rho)\over d\,\rho},\quad
D(\rho,0)=\displaystyle{d\,E_{0}(\rho)\over
d\,\rho},\quad\rho\in{\mathbb{R}}.$

The extended system is useful for studying the formation of singularities of
the solution associated with infinite derivatives.

We consider the characteristics of the system (5), outgoing from a fixed point
$\rho_{0}$. The system of characteristics has the form

(9) $\displaystyle{dP\over d\theta}=-E-\nu\,P,\quad\displaystyle{dE\over
d\theta}=\displaystyle{P\over\sqrt{1+P^{2}}},\quad\displaystyle{d\rho\over
d\theta}=\displaystyle{P\over\sqrt{1+P^{2}}},$

whence it immediately follows that

(10) $2\sqrt{1+P^{2}}+E^{2}\leq
2\sqrt{1+P^{2}_{0}(\rho_{0})}+E^{2}_{0}(\rho_{0})={\mathcal{E}}_{0}(\rho_{0})=\mbox{\rm
const}.$

Thus, the solution itself is always bounded, which cannot be guaranteed for
its derivatives.

Derivatives along a characteristic starting from a point $\rho_{0}$, satisfy

(11) $\displaystyle{dQ\over d\theta}=-\nu Q-D-
KQ^{2},\quad\displaystyle{dD\over d\theta}=K(1-D){Q},\quad
K=(1+P^{2})^{-3/2}.$

In the nonrelativistic case, system (11) is split off from (9), but in the
relativistic case they are connected via $K(\theta)$, that is, $P(\theta)$.
This complicates the problem considerably.

Although we do not know the explicit representation of $K(\theta)$, we can
evaluate this function from two sides. Namely,

(12) $0<\displaystyle{8\over{\mathcal{E}}_{0}^{3}}=K_{-}\leq K(\theta)\leq
K_{+}=1.$

The upper bound for $K_{+}$ is obvious, the lower bound follows from (10).

Let us note that the system (5) possesses solutions in the form of simple
waves, i.e. with the constraint $E=E(P)$. For this particular case, system (5)
reduces to one quasilinear equation, which can be analyzed separately.
However, this analysis is not so simple as for $\nu=0$, when one can obtain an
explicit criterion of singularities formation (see [13]). The simple waves for
$\nu>0$ satisfy the equation

(13) $\frac{dE(P)}{dP}=-\frac{P}{\sqrt{1+P^{2}}(E(P)+\nu P)},$

which cannot be integrated explicitly.

##  3\.  Analytical results

The results on estimates of the lifetime of a smooth solution to the equations
of relativistic cold plasma, proved in this section, continue the method
developed in [13] (see also [14]).

Let us summarize it briefly. Since the formation of a singularity is
associated with infinite derivatives of the solution of (5), the solution
remains smooth as long as the projection of the characteristic curve described
by equations (9), (11) onto the phase plane $(D,Q)$ remains bounded. We denote
by $\alpha=D(\rho,0),\beta=Q(\rho,0)$ the starting point of this projection.
As shown in [13], in the nonrelativistic case, when system (11) is autonomous,
for $\nu=0$ condition

(14) $\beta^{2}+2\,\alpha-1<0$

means that the phase curve in the $(D,Q)$ plane is an ellipse, which ensures
that the derivatives are bounded in time. If the condition (14) is not met,
the curve on the phase plane is a parabola or hyperbola, and its unboundedness
corresponds to the fact that the derivatives become infinite in a finite time.
In this case, we are dealing with a criterion for the formation of a
singularity in terms of the initial data. For $\nu>0$, in the nonrelativistic
case one can also obtain a criterion of this kind [15], which is rather
cumbersome. However, its analysis shows that the condition (14) enough to keep
the phase trajectory bounded. For $\nu>0$, the equilibrium position $(0,0)$ is
asymptotically stable and, globally in time, the smooth solution $(V,E)$ tends
to trivial stationary state at $\theta\to\infty$.

In this case, the behavior of the solution for $\nu\in[0,2)$ and $\nu>2$ is
sharply different. At $\nu\in[0,2)$, the equilibrium position $(0,0)$ is the
focus, therefore the phase trajectory makes an infinite number of revolutions
around the origin, which corresponds to the oscillatory motion of the medium.
For $\nu>2$, the equilibrium position $(0,0)$ is a node, and the medium does
not oscillate. In this case, the phase trajectory always turns out to be in
the quadrant $D<0$, $Q<0$, it either goes to the origin of coordinates in
infinite time, or goes to infinity in a finite time. The latter corresponds to
the formation of a singularity of the solution, but by increasing $\nu$, any
fixed initial data can be translated into an area such that the trajectory
leaving it goes to the origin. Since the physically natural value of $\nu$ is
small, no suppression of oscillations is observed in experiments. However,
formally, for arbitrary initial data (6), one can find such $\nu$ that the
solution of the problem (5), (6) will be globally smooth in time.

In the relativistic case, such results are impossible due to the fact that the
system (11) ceases to be autonomous. Therefore, one has to be content with
estimates of the location of the projections of the phase curves on $(Q,D)$,
based on the fact that the solution of (11) can be estimated through its
analogue with the substitution $K_{\pm}=\rm const$ instead of $K=K(\theta)$.
Note that $K(\theta)$ can be found in quadratures from (9), but this does not
help the investigation, since $K(\theta)$ cannot be explicitly represented.
Immediately, we note that the smaller the difference $K_{+}-K_{-}=1-K_{-}$,
that is, the less the initial data differ from the zero equilibrium position,
the more accurate estimates can be obtained. For the nonrelativistic case,
this difference is zero.

Let us assume that the solution do not correspond to a simple wave given as
(13) and introduce new variables $u={D}/{Q}$, $\lambda=(1-D)/{Q}$, in which
the system (11) is written as

(15) $\displaystyle{du\over d\theta}=u^{2}+\nu
u+K,\quad\displaystyle{d\lambda\over d\theta}=\lambda u+\nu\lambda.$

From the first equation (15) it follows that $u$ can tend to infinity in a
finite time (or should do so depending on the choice of $\nu$). This can
happen for two reasons: $Q$ vanishes for finite $D$ or $Q$ and $D$ both go to
infinity. The latter situation can also occur if $u$ is bounded. Similarly, if
$u$ vanishes, then $Q$ tends to infinity or $D$ vanishes. Analysis of the
phase trajectories shows that for $D\geq 1$, the value of $Q$ always tends to
infinity in a finite time, but from the requirement that the density be
positive it follows that $D<1$. In the domain $D\in(0,1),\,Q>0$, the
trajectory reaches the boundary $Q=0$, and in the domain $D\in(0,1),\,Q<0$ and
$D<0,\,Q>0$ the trajectory reaches the boundary $D=0$. In the quadrant
$D<0,\,Q<0$, where $u>0$, the value of $u$ can go to infinity for various
reasons. If $Q$ turns to zero, then the trajectory makes a revolution around
the origin, otherwise $Q$ turns to infinity in a finite time. We need to
distinguish these two situations.

Let us introduce the following notation:

$\mbox{quadrant \, I}:D<0,\,Q>0,\quad\mbox{quadrant\, II}:D>0,\,Q>0,$
$\mbox{quadrant \, III}:D>0,\,Q<0,\quad\mbox{quadrant\, IV}:D<0,\,Q<0.$

It follows from the above that the derivatives can become unbounded only in
quadrant IV. We will say that a trajectory makes a revolution (the solution
makes one oscillation) if it returns to the same quadrant from which it left.

For a simple wave, $u=E^{\prime}(P)$ (see (13)) and the behavior of $u$ does
not correspond to the behavior of the derivatives. This case requires a
separate analysis.

###  3.1. Behavior on the first revolution

First, we obtain the simplest sufficient conditions for the initial data,
allowing us to conclude that the trajectory certainly makes one revolution
(returns from quadrant IV to quadrant I) or goes to infinity in quadrant IV.

The behavior of the solution for $\nu\in[0,2)$ and $\nu\geq 2$ is sharply
different. In the first case, we encounter nonlinear resonance, so the phase
trajectory can go to infinity both at the first and at subsequent revolutions.
In the second case, there is no oscillations and the behavior of the solution
is determined already on the first revolution.

######  Proposition 3.1.

Let $0<\nu<2$, and the initial data (6) are bounded on $\mathbb{R}$, belong to
the class $C^{2}(\mathbb{R})$ and do not correspond to a simple wave given as
(13). If condition (14) holds for all $\rho$ together with any of the
conditions

  * •

$\beta<0$ or

  * •

$\beta=0$, $\alpha>0$,

then the solution to problem (5), (6) is classically smooth at least until the
time $t_{*}=\min\limits_{\rho\in\mathbb{R}}T(\rho)$, $t_{*}\geq 2\pi$ (each
trajectory makes at least one revolution).

Let us denote

(16)
$T_{-}=\frac{1}{K_{-}-\nu^{2}/4}\left(\frac{\pi}{2}-\arctg\frac{\alpha/\beta+\nu/2}{\sqrt{K_{-}-\nu^{2}/4}}\right).$

######  Proposition 3.2.

Let $0<\nu<2\sqrt{K_{-}}$ (see (12)), the initial data (6) are bounded on
$\mathbb{R}$, belong to the class $C^{2}(\mathbb{R})$ and do not correspond to
a simple wave given as (13). Then if for some $\rho_{0}$ condition

(17) $\beta^{2}+\frac{\exp\\{\nu
T_{-}\\}}{K_{-}}(2\alpha-1)+\left(\frac{1}{K_{-}-\nu^{2}/4}-\frac{\exp\\{\nu
T_{-}\\}}{K_{-}}\right)\alpha^{2}>0,$

holds, then the solution to problem (5), (6) loses its smoothness before
$t_{*}=\min\limits_{\rho\in\mathbb{R}}T(\rho)$, $t_{*}\geq 2\pi$ (there is a
trajectory that does not make a revolution).

###### Proof.

We prove Propositions 3.1 and 3.2 simultaneously. First of all, note that
${Q}^{-1}=\lambda+u$ and this expression for $Q<0$ can be estimated from both
sides as

(18) $\psi_{-}(\theta)\leq\frac{1}{Q}\leq\psi_{+}(\theta),$

where

$\psi_{\pm}(\theta)=-\frac{\nu}{2}+{\tilde{K}}_{\pm}\tg\varphi_{\pm}+\exp\\{\nu\theta/2\\}\displaystyle{\lambda_{0}\sqrt{1+\tg^{2}\varphi_{\pm}}\over\sqrt{\displaystyle{(u_{0}+\nu/2)^{2}\over{\tilde{K}}_{\pm}^{2}}+1}},$
${\tilde{K}}_{\pm}=\sqrt{K_{\pm}-\nu^{2}/4},\quad\varphi_{\pm}={\tilde{K}}_{\pm}\,\theta+\arctg\frac{u_{0}+\nu/2}{{\tilde{K}}_{\pm}}.$

If in the domain of negative $Q$, that is, in quadrants III and IV, ${Q}^{-1}$
vanishes before it turns to infinity, then the solution loses its smoothness
in a finite time. Thus, we need to find the condition for $\psi_{-}(\theta)$
to vanish when $\psi_{-}(\theta)$ increases with respect to $\theta$.

On the other hand, if $\psi_{+}(\theta)\neq 0$ for $\theta\in(0,T_{+})$, where

$T_{+}=\displaystyle{1\over
1-\nu^{2}/4}\left(\frac{\pi}{2}-\arctg\displaystyle{u_{0}+\nu/2\over{\tilde{K}}_{+}}\right),$

then the trajectory returns to quadrant I and the solution remains smooth for
at least one revolution of the trajectory around the origin. Direct
calculation shows that

$\psi_{+}(\theta)=\frac{\exp\\{\nu\theta/2\\}}{\sqrt{\displaystyle{(u_{0}+\nu/2)^{2}\over{\tilde{K}}_{+}^{2}}+1}}\,F_{+}(\theta),$

where $F_{+}(\theta)$ characterizes the denominator $Q$ (see [15]) and is
defined as

$F_{+}(\theta)=1-\alpha+\left(\displaystyle{\nu\alpha+2\beta\over
2{\tilde{K}}_{+}}\,\sin{\tilde{K}}_{+}\theta+\alpha\,\cos{\tilde{K}}_{+}\theta\right)\,\exp(-\nu\theta/2).$

Thus, the conditions guaranteeing the smoothness of the solution during the
first oscillation correspond to the conditions of non-vanishing of
$F_{+}(\theta)$. It follows from the results of [15] that non-vanishing of
$F_{+}(\theta)$ for $\nu=0$ guarantees non-vanishing of this function for an
arbitrary $0<\nu<2$, therefore the inequality $\beta^{2}+2\,\alpha-1<0$ leads
to the desired result. Conditions $\beta<0$ or $\beta=0$, $\alpha>0$ ensures
that the trajectory falls into quadrant III. Note that if this condition is
not imposed, then in the absence of an upper bound for the trajectory for
$Q>0$ (in quadrants I and II), it cannot be guaranteed that the condition (14)
will be satisfied when the trajectory is in quadrant III. If the trajectory
starts, for example, in quadrant I, then it can go to infinity in quadrant IV
and not make a full revolution. Note that the results of the next section
allow one to obtain two-sided estimates of the position of the trajectory. In
conclusion, note that in [13] it is shown that a turn along a trajectory for
$K\geq 1$ takes time greater than or equal to $2\pi$. Proposition 3.1 is
proved. $\Box$

Let us proceed to the proof of Proposition 3.2. We will restrict ourselves
here to a rather rough but simple construction. More accurate construction of
a domain on the plane corresponding to the initial data leading to the loss of
smoothness is described below, in Remark 3.1. Here we note that since
$\lambda_{0}<0$, then

$\psi_{-}(\theta)>\bar{\psi}_{-}(\theta)=-\frac{\nu}{2}+{\tilde{K}}_{-}\tg\varphi_{-}+\exp\\{\nu
T_{-}/2\\}\displaystyle{\lambda_{0}\sqrt{1+\tg^{2}\varphi_{-}}\over\sqrt{\displaystyle{(u_{0}+\nu/2)^{2}\over{\tilde{K}}_{-}^{2}}+1}},$

where $T_{-}$ is defined by (16), and $\alpha/\beta=u_{0}$.

It is easy to calculate that $\bar{\psi}_{-}(\theta)$ has a root on
$(0,T_{-})$, if

$K_{-}-\exp\\{\nu
T_{-}\\}\,\frac{{\tilde{K}}_{-}^{2}}{(u_{0}+\nu/2)^{2}+{\tilde{K}}_{-}^{2}}\,\lambda^{2}_{0}>0,$

which can be rewritten as (17). Proposition 3.2 is proved. $\Box$

######  Remark 3.1.

One can see that

$\psi_{-}(\theta)=\frac{\exp\\{\nu\theta/2\\}}{\sqrt{\displaystyle{(u_{0}+\nu/2)^{2}\over{\tilde{K}}_{-}^{2}}+1}}\,F_{-}(\theta),$
$F_{-}(\theta)=1-\alpha+\left(\frac{\nu\alpha+2\beta}{{\tilde{K}}_{-}}\sin({\tilde{K}}_{-}\theta)+\alpha\cos({\tilde{K}}_{-}\theta)\right)\exp(-\nu\theta/2),$

therefore it is possible to construct on the plane $(\alpha,\beta)$ the curve
$\Phi_{\nu}(\alpha,\beta)=0$ (the procedure is described in detail in [15]).
The initial data corresponding to a solution that loses its smoothness already
during the first revolution of the phase trajectory are to the right of this
curve. Since we are dealing only with an estimate of $Q$, we cannot divide the
plane $(\alpha,\beta)$ into parts corresponding to a smooth and losing
smoothness solution exactly (to obtain a criterion for the formation of a
singularity).

For the case $\nu>2$, the estimate (18) is also valid, but here

$\psi_{\pm}(\theta)=-\frac{\nu}{2}+{\tilde{K}}_{\pm}\tanh\varphi_{\pm}+\exp\\{\nu\theta/2\\}\displaystyle{\lambda_{0}\sqrt{1+\tanh^{2}\varphi_{\pm}}\over\sqrt{\displaystyle{(u_{0}+\nu/2)^{2}\over{\tilde{K}}_{\pm}^{2}}+1}},$
${\tilde{K}}_{\pm}=\sqrt{\nu^{2}/4-K_{\pm}},\quad\varphi_{\pm}={\tilde{K}}_{\pm}\,\theta+{\rm
arcth}\frac{u_{0}+\nu/2}{{\tilde{K}}_{\pm}}.$

Based on these estimates, one can obtain statements similar to the previous
two, but here a statement of a different kind is also true.

######  Proposition 3.3.

For any initial data (6) (not corresponding to a simple wave given as (13))
belonging to the class $C^{2}(\mathbb{R})$ and bounded at $\mathbb{R}$ there
exists a value $\nu>0$ such that the solution to problem (5), (6) remains
classically smooth for all $t>0$ and stabilizes to a stationary state for
$t\to\infty$.

###### Proof.

Arguing as in the proof of Proposition 3.1, we obtain that for $\nu>2$

$\psi_{+}(\theta)=\frac{\exp\\{\nu\theta/2\\}}{\sqrt{\displaystyle{(u_{0}+\nu/2)^{2}\over{\tilde{K}}_{+}^{2}}+1}}\,F_{+}(\theta),$
$F_{+}(\theta)=1-\alpha+\left(\displaystyle{\nu\alpha+2\beta\over
2{\tilde{K}}_{+}}\,\sinh{\tilde{K}}_{+}\theta+\alpha\,\cosh{\tilde{K}}_{+}\theta\right)\,\exp(-\nu\theta/2)$

(see [15]), and the conditions guaranteeing the smoothness of solution
correspond to the conditions of non-vanishing of $F_{+}(\theta)$. It is easy
to see that for any fixed $\alpha$ and $\beta$ by increasing $\nu$ one can get
non-vanishing of $F_{+}(\theta)$ for all $t>0$. Moreover, $D$ and $Q$ tend to
zero as $t\to\infty$. Proposition 3.3 is proved. $\Box$

###  3.2. Calculation of the guaranteed number of revolutions

In this section, we describe a procedure that allows us to obtain an estimate
from below for the number of revolutions that makes projections onto the plane
$(D,Q)$ of the trajectory going from a particular point $\rho_{0}$ before the
formation of a singularity, and thus an estimate of the time of formation of a
singularity from below as infimum over $\rho_{0}\in\mathbb{R}$. The procedure
consists in constructing the constraints of the projection of the trajectory,
that is, the curves between which this projection is guaranteed to be located.
Especially, it is important to construct the upper limiter, since if the
condition (14) is not satisfied for it, then this condition cannot be
guaranteed for the projection of the trajectory itself. The equations of the
curves constituting the limiter are found explicitly. As mentioned above, this
procedure makes sense for small $\nu$.

Note that (11) implies the non-autonomous equation

(19) $\frac{1}{2}\displaystyle{dQ^{2}\over
dD}=-\nu\frac{Q}{(1-D)K(P)}-\frac{D}{(1-D)K(P)}-\frac{Q^{2}}{(1-D)}:=\Psi(D,Q,P).$

Since the density is positive, it follows that $1-D>0$. Let us denote

$\displaystyle\Sigma_{ij}$ $\displaystyle=$
$\displaystyle-\nu\frac{Q}{(1-D)K_{i}}-\frac{D}{(1-D)K_{j}}-\frac{Q^{2}}{(1-D)},$

where $i$, $j$ take the values $m$ or $p$, and $K_{m}=K_{-}$, $K_{p}=K_{+}$.
Thus,

$\displaystyle\Sigma_{pm}(D,Q)\leq\Psi(D,Q,P)\leq\Sigma_{mp}(D,Q),$
$\displaystyle\quad(Q>0,D<0),\,I\,\mbox{ },$
$\displaystyle\Sigma_{pp}(D,Q)\leq\Psi(D,Q,P)\leq\Sigma_{mm}(D,Q),$
$\displaystyle\quad(Q>0,D>0),\,II\,\mbox{ },$
$\displaystyle\Sigma_{mp}(D,Q)\leq\Psi(D,Q,P)\leq\Sigma_{pm}(D,Q),$
$\displaystyle\quad(Q<0,D>0),\,III\,\mbox{ },$
$\displaystyle\Sigma_{mm}(D,Q)\leq\Psi(D,Q,P)\leq\Sigma_{pp}(D,Q),$
$\displaystyle\quad(Q<0,D<0),\,IV\,\mbox{ }.$

In what follows, we denote the integral curves of the equations
$\displaystyle{dQ^{2}\over dD}=2\Sigma_{ij}(D,Q)$, (i.e. (19)) as
$\Lambda_{ij}$, where $ij$ takes indices $mm$, $mp$, $pm$, $pp$. Since in
quadrants I and II the function $D$ increases as it moves along the projection
of the phase curve, and in quadrants III and IV it decreases, then it follows
from Chaplygin’s theorem that the projection of the phase curve is bounded

  * •

in quadrant I from below by $\Lambda_{pm}$, from above by $\Lambda_{mp}$,

  * •

in quadrant II from below by $\Lambda_{pp}$, from above by $\Lambda_{mm}$,

  * •

in quadrant III from below by $\Lambda_{mp}$, from above by $\Lambda_{pm}$,

  * •

in quadrant IV from below by $\Lambda_{mm}$, from above by $\Lambda_{pp}$.

Thus, starting from a point, we switch the upper and lower boundaries of the
trajectory, moving from one quadrant to another (see Figure 1). As the top
limiter, we get a curve that unfolds around the origin. While this curve
intersects axis $Q=0$ at a point $D_{*}$, $0<D_{*}<\frac{1}{2}$, then, as
follows from Proposition 3.1, the projection of the phase trajectory makes
another revolution. This process can be continued until the point of
intersection with axis $Q=0$ in quadrant II is greater than $\frac{1}{2}$. If
$n$ is the number of complete revolutions, then the lifetime of the smooth
solution is estimated from below as $2\pi n$. Note that the fact that the
curve $\Lambda_{mm}$ intersects axis $Q=0$ at $D>\frac{1}{2}$ does not mean
that the projection of the phase trajectory itself also intersects $Q=0$ at
such a point. Numerical calculations show that the resulting estimate is
rather rough.

We can also construct the lower bounds of the projection of the phase
trajectory from the side of the origin using the curves $\Lambda_{pm}$ and
$\Lambda_{pp}$. However, it is easy to see that they form a spiral that turns
to zero, so if the conditions of Proposition 3.2 are not met, they will not be
met at the next turn. That is, no new information can be obtained.

Example. The curves $\Lambda_{mm}$ and $\Lambda_{mp}$, used as upper limiters
for the projection of the phase trajectory, can be found explicitly by the
formulas specified in [15]. They are rather cumbersome, so we consider the
case $\nu=0$ for illustration.

The curves $\Lambda_{mp}$ and $\Lambda_{mm}$ for $\nu=0$ have the form
$C(1-D^{2})=K_{+}Q^{2}+D^{2}$ and $C(1-D^{2})=K_{-}Q^{2}+D^{2}$, respectively.
The constant $C$ is first found from the initial data, and then in each next
quadrant from the coordinates of the point of intersection of the curve in the
previous quadrant with the coordinate axis.

If we choose the initial data (32) used in the calculations in Section 4
(Subsection 4.2), that is

$E_{0}(\rho)=0.4761\,\rho\,\exp\left\\{-0.0987654321\,{\rho^{2}}\right\\},\quad
P_{0}(\rho)=V_{0}(\rho)=0,$

we get $n=3$ and the guaranteed lifetime of a smooth solution is $6\pi$. In
this case, the smallest number of revolutions can be guaranteed for
trajectories starting from points $|\rho_{0}|\in(0.69,1.2)$.

From the results of numerics presented in Figure 2 for this case, we see that
in reality, the number of oscillations before blowup is 4.

Note that, as $\nu$ increases, the curves $\Lambda$ approach the origin, then
for each fixed starting point of the trajectory, the guaranteed number of
revolutions increases with $\nu$.

Figure 1. Projection of the phase curve of system (5), (11) onto the plane
$(D,Q)$ (solid). Bilateral limiters in each quadrant (dashed). 1 - curve
$\Lambda_{mp}$, 2 - curve $\Lambda_{mm}$, 3 - curve $\Lambda_{pm}$, 4 - curve
$\Lambda_{pp}$. The intersection of the upper limiter with the straight line
$D=\frac{1}{2}$ means that the continuation of the solution for the next
revolution is not guaranteed.

###  3.3. Small deviations from the trivial state

Previously we obtained the following result.

######  Theorem 3.1.

[13] For $\nu=0$, any solution to the Cauchy problem (5), (6) which is an
arbitrarily small deviation from the equilibrium $P=0$, $E=0$, such that
$2\sqrt{1+P_{0}^{2}(\rho)}+E_{0}^{2}(\rho)$ is not identically equal to a
constant, the derivatives turn to infinity in a finite time.

We want to show that for an arbitrarily small collision rate, this fact is no
longer true.

Note that the case $2\sqrt{1+P_{0}^{2}(\rho)}+E_{0}^{2}(\rho)=C$ for all
$\rho\in\mathbb{R}$ ($C_{\rho}=0$) corresponds to the simple wave for $\nu=0$,
see (13).

Let us show that the following theorem holds.

######  Theorem 3.2.

For any $\nu>0$ there exists $\epsilon>0$ such that the solution to the Cauchy
problem (5), (6) with initial data satisfying the conditions
$E_{0}^{2}(\rho)+P_{0}^{2}(\rho)+\alpha^{2}+\beta^{2}\leq\epsilon^{2}$
uniformly in $\rho$ remains smooth for all $\theta>0$.

We need the following lemma.

######  Lemma 3.1.

We denote

$\displaystyle{\alpha\over\beta}=u_{0},\quad\displaystyle{1-\alpha\over\beta}=\lambda_{0},\quad\beta\neq
0.$

Let $z(\theta)$ be a solution to the Cauchy problem

(20) $\displaystyle{d^{2}z\over d\theta^{2}}-\nu\displaystyle{dz\over
d\theta}+K(\theta)z=0,\quad z(0)=1,\,z^{\prime}(0)=-u_{0}.$

If for at least one point $\rho_{0}$ there is a moment $\theta_{*}>0$ such
that

(21) $z^{\prime}(\theta_{*})=\lambda_{0}e^{{\nu}\theta_{*}},$

then the derivatives of the solution of (5), (6) that do not correspond to a
simple wave given as (13) blow up in a finite time. Otherwise, the solution
remains smooth for all $\theta>0$.

###### Proof.

The change $u=-z^{\prime}/z$ reduces the first of the equations (15) to (20).
This equation is homogeneous, so we set $z(0)=1$, whence
$z^{\prime}(0)=-u_{0}$ follows. From the second equation (15) we have
$u=\lambda^{\prime}/\lambda-\nu$, and
$\lambda=\frac{\lambda_{0}}{z}\exp\left({\nu}\theta\right)$.

Thus,

$\displaystyle{u\over\lambda}=\displaystyle{D\over
1-D}=-\,\displaystyle{z^{\prime}\over\lambda_{0}}e^{-{\nu}\theta}.$

If $D$ goes to infinity at $\theta=\theta_{*}$, then
$-\,\displaystyle{z^{\prime}\over\lambda_{0}}e^{-{\nu}\theta_{*}}=-1$, that
is, the condition (21) holds. If there is no such moment, then $D$ and $Q$
remain bounded. The lemma is proved. $\Box$

Proof of Theorem 3.2. We give the proof for $\nu\in(0,2)$, although it can be
done similarly for $\nu>0$.

We start with case of general initial data and then prove the theorem for the
simple wave. In order to apply the lemma, we need an explicit form
$K(\theta)$. Let us make an assumption about the smallness of the initial
perturbation. If $P_{0}^{2}(\rho_{0})+E_{0}^{2}(\rho_{0})\leq\epsilon^{2}$,
$\epsilon\ll 1$, then (9) implies that $P$ and $E$ remain small as long as the
solution remains smooth, in particular, $P^{2}\leq\epsilon^{2}\ll 1$. The
function $P(\theta)$ satisfies the equation

${P}^{\prime\prime}+\nu P^{\prime}+\frac{P}{\sqrt{1+P^{2}}}=0,$

therefore, up to $o(\epsilon)$ we get
$P=\epsilon\,e^{-\frac{\nu}{2}\theta}\,\sin(\omega\theta+\theta_{0})$,
$\omega=\sqrt{1-\frac{\nu^{2}}{4}}$. Without loss of generality, we set
$\theta_{0}=0$, substitute $P$ in $K$ and expand the result in a series in
$\epsilon$. We obtain

(22)
$K(\theta)=1-\frac{3}{2}\,e^{-\nu\theta}\,\sin^{2}\omega\theta\,\epsilon^{2}+O(\epsilon^{4}).$

We neglect the terms of order higher than the second in (22), substitute the
result into (20) and seek a solution in the form
$z(\theta)=z_{0}(\theta)+\epsilon^{2}z_{1}(\theta)+o(\epsilon^{2})$. In the
standard way we get

$\displaystyle
z(\theta)=e^{-\frac{\nu}{2}\theta}\left(\frac{\nu/2-u_{0}}{\omega}\sin\omega\theta+\cos\omega\theta\right)+$
$\displaystyle\left[e^{\frac{\nu}{2}\theta}A_{1}\sin(\omega\theta+\theta_{1})+e^{-\frac{3\nu}{2}\theta}(A_{2}\sin(\omega\theta+\theta_{2})+A_{3}\sin(3\omega\theta+\theta_{3})\right]\,\epsilon^{2}+o(\epsilon^{2}),$

where $A_{i}=A_{i}(\nu,u_{0}),$ $\theta_{i}=\theta_{i}(\nu,u_{0}),$ $i=1,2,3$.

Let us fix $\nu>0$. For $\theta_{*}$ to exist (see Lemma 3.1), there must be a
positive root of the equation

(23)
$\lambda_{0}e^{\nu\theta}=z^{\prime}(\theta)=e^{-\frac{\nu}{2}\theta}\Theta_{1}+\epsilon^{2}[e^{\frac{\nu}{2}\theta}\Theta_{2}+e^{-\frac{3\nu}{2}\theta}\Theta_{3}]+o(\epsilon^{2}),$

where $\Theta_{i}=\Theta_{i}(\nu,u_{0},t)$, $i=1,2,3$, are bounded functions.
The absolute value of the left hand side of (23) grows as $\theta\to\infty$
faster than the absolute value of the right hand side, so a possible root can
be only for $\theta$, belonging to the first period of $\Theta_{1}$. The
amplitude of $\Theta_{1}$ is
$A_{\Theta_{1}}=\left(\frac{\left(\frac{\nu}{2}u_{0}-1\right)^{2}}{\omega^{2}}+u_{0}^{2}\right)^{1/2}$.
Thus, if condition $A_{\Theta_{1}}<|\lambda_{0}|$ holds, then as for
$\theta>0$ the amplitude of the $\epsilon$-independent expression on the right
hand side decreases, then by decreasing $\epsilon$, we can always ensure that
the graph of the growing exponential function on the left side does not
intersect with the graph of the function on the right side. This ensures that
there is no root of equation (23). The condition
$A_{\Theta_{1}}<|\lambda_{0}|$ implies

$\omega^{2}(2\alpha^{2}-1)+(\beta-\frac{\nu}{2}\alpha)^{2}<0,$

therefore for sufficiently small $\alpha^{2}+\beta^{2}$ this inequality holds.
Thus, the theorem is proved for the case of general initial data.

For the case of a simple wave (11) and (13) result in

(24) $\displaystyle{dQ\over
d\theta}=-Q\left(K(P)Q+\nu-\frac{P}{\sqrt{1+P^{2}}(E(P)+\nu P)}\right).$

We can see that $K(P)=1+O(P)$, $\frac{P}{\sqrt{1+P^{2}}(E(P)+\nu P)}=O(P)$,
therefore for any fixed $\nu>0$ there exist sufficiently small
$P^{2}_{0}(\rho_{0})+Q^{2}_{0}(\rho_{0})$ such that the solution to (24) tends
to zero as $\theta\to\infty$. Now Theorem 3.2 is proved for all cases. $\Box$

##  4\.  Numerical illustrations

The approach to the analysis of characteristics used to obtain analytical
results can be transformed into a high-precision approximate method for
calculating relativistic oscillations of cold plasma with allowance for
collisions. The specificity of the method lies in the identification of the
displacement function of particles (electrons) when using Lagrangian variables
with a function that characterizes the electric field when using Eulerian
variables.

In fact, the method presented below is a method of characteristics; however,
it is convenient to present its construction in terms of the Lagrangian
description of the medium, i.e. using the concepts of particles and their
trajectories.

Note that the simulated medium is represented by charged particles located in
the background field formed by stationary ions. Therefore, for each particle
there is an ”equilibrium” position when the background field is canceled and
which is convenient to take as values of its Lagrangian coordinate. In this
case, the displacement function relative to the equilibrium position is
responsible for the formation of the electric field:

(25) $\rho(\rho^{L},\theta)=\rho^{L}+R(\rho^{L},\theta),$

where $\rho^{L}$ is the ”equilibrium” position of the particle when it does
not contribute to the formation of the electric field, $R(\rho^{L},\theta)$ is
its displacement function generating an electric field at the point of the
trajectory $\rho(\rho^{L},\theta)$. According to [7], in the flat one-
dimensional case, there is a simple relationship between the field and the
displacement

(26) $R(\rho^{L},\theta)=E(\rho,\theta)\equiv
E(\rho^{L}+R(\rho^{L},\theta),\theta),$

which makes the functions of the electric field and displacement on the
particle trajectory indistinguishable. Note that this approach additionally
makes it possible not to solve simultaneously two identical differential
equations with different initial conditions (to determine the Eulerian
trajectory of a particle and the electric field along it). In this case, the
formal transition from the Euler coordinates $(\rho,\theta)$ to the Lagrangian
coordinates $(\rho^{L},\tau)$ is carried out by the usual transformation

$\tau\equiv\theta,\quad\rho^{L}=\rho-\int\limits_{\tau_{0}}^{\tau}d\tau^{\prime}V(\rho^{L},\tau^{\prime})\,,$

but instead of $\tau$ the old notation $\theta$ is used .

###  4.1.  Numerical algorithm

To find a numerical solution on the line $\rho\in(-\infty,\infty)$, we define
an arbitrary grid at the initial moment of time $\theta=0$:

$\rho_{1}(0)<\rho_{2}(0)<\dots<\rho_{M}(0),$

consisting of $M$ nodes. To each node $\rho_{k}(\theta=0),\;1\leq k\leq M,$
place a particle, marked with the Lagrangian coordinate $\rho^{L}_{k}$.
Moreover, for each particle from (5) and (26), there follow equations
describing the dynamics of particles in Lagrangian variables:

(27) $\begin{array}[]{c}\displaystyle{d\,P(\rho^{L}_{k},\theta)\over
d\,\theta}=-R(\rho^{L}_{k},\theta)-\nu
P(\rho^{L}_{k},\theta),\quad\displaystyle{d\,R(\rho^{L}_{k},\theta)\over
d\,\theta}=V(\rho^{L}_{k},\theta),\\\ \quad
V(\rho^{L}_{k},\theta)=\displaystyle{P(\rho^{L}_{k},\theta)\over\sqrt{1+P^{2}(\rho^{L}_{k},\theta)}}\,,\quad
k=1,2,\dots,M.\end{array}$

Let us use the equalities (25) and (26) to obtain the missing initial
conditions for systems of the form (27). At the node $\rho_{k}$ for
$\theta=0$, $E_{0}(\rho_{k})$ is defined, i.e.
$\rho_{k}=\rho^{L}_{k}+E_{0}(\rho_{k})$ holds, whence the initial conditions
follow for the particle with the number $k$:

(28) $P(\rho^{L}_{k},\theta=0)=P_{0}(\rho_{k}),\quad
R(\rho^{L}_{k},\theta=0)=E_{0}(\rho_{k}),\quad k=1,2,\dots,M,$

as well as the equilibrium value of the Lagrangian coordinate $\rho^{L}_{k}$
itself. The obtained relations allow, instead of problem (5), (6), written in
Eulerian variables, solve numerically the problem (27), (28) formulated in
Lagrangian variables.

The solution of problem (27), (28) determined in this way does not allow
finding the spatial derivatives of the required functions $P$ and $E$, which
makes it impossible to determine the electron density function $N$ in
accordance with (4). To avoid this drawback, it is enough to reformulate the
problem (7), (8) from Eulerian variables to Lagrangian variables, using the
reasoning presented above. Formal transformations give the equations

(29) $\begin{array}[]{l}\displaystyle{d\,Q(\rho^{L}_{k},\theta)\over
d\,\theta}=-D(\rho^{L}_{k},\theta)-\displaystyle{W^{2}(\rho^{L}_{k},\theta)\over\left(1+P(\rho^{L}_{k},\theta)^{2}\right)^{3/2}}-\nu
Q(\rho^{L}_{k},\theta),\vspace{1.0 ex}\\\
\displaystyle{d\,D(\rho^{L}_{k},\theta)\over
d\,\theta}=(1-D(\rho^{L}_{k},\theta))\,\displaystyle{Q(\rho^{L}_{k},\theta)\over\left(1+P(\rho^{L}_{k},\theta)^{2}\right)^{3/2}},\end{array}$

and initial conditions

(30) $Q(\rho^{L}_{k},\theta=0)=Q_{0}(\rho_{k}),\quad
D(\rho^{L}_{k},\theta=0)=D_{0}(\rho_{k}),$

corresponding to individual particles with numbers $k=1,2,\dots,M.$

Now, from the solution of (29), (30) and the relation (25) for each particle

(31) $\rho_{k}(\theta)=\rho^{L}_{k}+R(\rho^{L}_{k},\theta),\quad
k=1,2,\dots,M,$

it is possible at a Eulerian point of space $(\rho_{k}(\theta),\theta)$ to
determine the value of the electron density

$N(\rho_{k},\theta)=1-D(\rho^{L}_{k},\theta).$

Thus, the proposed numerical algorithm consists in finding solutions of
equations (27) with conditions (28) at the nodes of the Eulerian grid (31), as
well as equations (29) with conditions (30) for each of the particles numbered
$k=1,2,\dots,M.$ Note that the search for an approximate solution presupposes
the existence, uniqueness, and smoothness of the exact solution of the problem
due to Propositions 3.1 \- 3.3, and, as noted in Section 1, general theorems
on the existence and uniqueness of a locally smooth solution to the Cauchy
problem for systems of hyperbolic equations [5], pp.221-225. If the solution
has sufficient smoothness in the variable $\theta$, it seems very convenient
to use the classical Runge - Kutta method of the fourth order of accuracy
[12], otherwise (of less smoothness) one should use schemes of a lower order
of accuracy up to the Euler method. It is necessary to clarify that in this
situation the accuracy of the obtained approximation is determined exclusively
by the smoothness of the solution.

Note also that the stability of the time integration of the equations for the
momentum and the electric field is completely determined by the inequality
(10). There is no analogous property for equations describing the dynamics of
their derivatives; therefore, under initial conditions that do not lead to the
intersection of Lagrangian trajectories in the non-relativistic case, in the
relativistic case these trajectories can intersect. This situation leads to a
discontinuity in the function $E(\rho,\theta)$, and, as a consequence, the
singularity electron density functions in accordance with (4). Therefore, in
the process of calculations, it is necessary to constantly monitor the
preservation of the starting order of particles, i.e. conditions
$\rho_{k}(\theta^{n})<\rho_{k+1}(\theta^{n})$ for all $k=1,2,\dots,M-1.$
Violation of this condition for at least one value of $k$ means a blow-up
(impossibility of further application of the model), denoted by the term
”breaking”. The above means that observing or investigating a solution with
specific initial data only makes sense until the moment of blow-up.

For this reason, further we will keep in mind that the numerical method is
applied on a limited time interval $\theta\in[0,\theta_{br}]$, while the
desired solution to the problem (4) - (6) exists , is unique and bounded
together with the necessary derivatives. In other words, it is further assumed
that the order of Lagrangian particles does not change on the segment
$[0,\theta_{br}]$, that is, the inequality

$h_{k}(\theta)=\rho_{k+1}(\theta)-\rho_{k}(\theta)>0\quad\forall
k=1,2,\dots,M-1,$

holds together with the sufficient condition for a local in the time existence
(see Proposition 3.1)

$\left(\displaystyle{\partial
P(\rho,\theta)\over\partial\rho}\right)^{2}+2\,\displaystyle{\partial
E(\rho,\theta)\over\partial\rho}-1<0.$

Summarizing the above description of the numerical algorithm, we emphasize
that for solving the original problem (4) - (8) written in Eulerian variables,
it is convenient to go over to Lagrangian variables and use approximate
methods of integration over time to calculate the required functions at the
points $\rho_{k}(\theta^{n}),\quad k=1,2,\dots,M,$ belonging to the
trajectories of particles. This is quite enough for studying most of the
problems. However, it is quite possible that situations arise when it is
required to determine the solution at the given Eulerian points
$(\rho,\theta)$, which do not have to belong to the calculated trajectories of
the particles.

In this case, at the moment of time $\theta^{n}$ it is necessary to determine
first the interval $[\rho_{k}(\theta^{n}),\rho_{k+1}(\theta^{n})]$, to which
the given value $\rho$ belongs, and then use the interpolation procedure to
find the approximate value.

Taking into account that at the nodes of the Euler grid
$\rho_{k}(\theta^{n}),\quad k=1,2,\dots,M,$ not only the functions
$P(\rho,\theta)$ and $E(\rho,\theta)$ are used, but also their spatial
derivatives, it seems very convenient to use the Hermitian cubic interpolation
for this purpose. The derivation of the necessary formulas and error estimates
is given in [18] pp.24-39, the practical details of use, including the
necessary programs, are well described in [12], pp.110-114.

###  4.2.  Calculation results

As the initial conditions, we choose the functions

(32)
$E_{0}(\rho)=\left(\displaystyle{a_{*}\over\rho_{*}}\right)^{2}\rho\exp\left\\{-2\displaystyle{\rho^{2}\over\rho_{*}^{2}}\right\\},\quad
P_{0}(\rho)=V_{0}(\rho)=0.$

According to (32) and (4), initially the density minimum is at the origin, and
the density deviations from the equilibrium equal to one decay exponentially.
This choice is the most natural from the point of view of physics, since it
simulates the effect on a rarefied plasma of a short powerful laser pulse when
it is focused into a line (this can be achieved using a cylindrical lens), see
details in [17].

As the working area of localization of oscillations, we choose the segment
$[-d,d]$ at $d=4.5\,\rho_{*}$. In this case, the amplitudes of particle
oscillations at the boundaries will coincide in order of magnitude with the
machine accuracy, due to $\exp\\{-2d^{2}/\rho_{*}^{2}\\}\approx 2.58\cdot
10^{-18}$, that is, the near-boundary characteristics will coincide with the
straight lines $\rho\approx\pm d$ with a high degree of accuracy. It is
convenient to choose a uniform initial grid
$\rho_{k}(\theta=0)=kh-d,\;k=0,1,\dots,M,\,h=2d/M.$

As an illustration of the algorithm, we present the calculations of the blow-
up effect that can be observed after several periods of oscillations. We fix
the parameters $a_{*}=3.105,$ $\rho_{*}=4.5$ in order to satisfy the condition
of a local in time existence of a solution for system (5) from Proposition
3.1, and also to preserve the continuity with the results of numerical and
asymptotic analysis from [9]. The calculations for the extended system (5) -
(8) were carried out by the classical Runge - Kutta method of the 4th order of
accuracy at $\tau=h=0.01$ until the intersection of the trajectories of
Lagrangian particles, which in the Euler formulation corresponds to an
infinite electron density. In order to control the accuracy, calculations were
additionally carried out with grid parameters two times smaller than the main
(working) ones.

The results of calculations for the electron density as a function of time at
various parameters $\nu$, which characterizes the damping of oscillations due
to collisions of electrons and ions, are shown in Figures 2 - 4. Figure 2
shows the time dependence of the density maximum in a collisionless plasma at
$\nu=0$, i.e. in the absence of damping. It follows from Figure 1 that for the
given calculation parameters in the fourth oscillation period, a density
maximum is formed outside the oscillation axis, which already in the next
period at $\theta^{(0)}_{wb}\approx 29.5$ has an infinitely large value.

Figure 2. Time dependence of the electron density in nonlinear oscillations of
a collisionless plasma ($\nu\theta^{(0)}_{wb}=0$). The maximum density in the
entire computational domain (solid) and the density at the origin (dashed).
Figure 3. Time dependence of the electron density in nonlinear plasma
oscillations taking into account collisions ($\nu\theta^{(0)}_{wb}=0.2$). The
maximum density in the entire computational domain (solid), the density at the
origin (dashed).  Figure 4. Time dependence of the electron density in
nonlinear plasma oscillations with frequent collisions
($\nu\theta^{(0)}_{wb}=0.5$). The maximum density in the entire computational
domain (solid), the density at the origin (dashed).

When the collisions of electrons are taken into account, the time for the
appearance of the density singularity increases. For example, for
$\nu\theta^{(0)}_{wb}=0.2$, the off-axis maximum (see Figure 2) has an
infinite value only in the third period after its formation, and not in the
second period as it was in collisionless case. Calculations show that the
density singularity for given initial parameters arises only in the case of
relatively rare collisions, when the inequality $\nu\theta^{(0)}_{wb}\leq
0.422$. Note that the dependence of the breaking time on the parameter
$\nu\theta^{(0)}_{wb}$ with allowance for collisions was established by
asymptotic methods in [9]. When the equality $\nu\theta^{(0)}_{wb}=0.422$
holds, the time for the density singularity to appear is $\theta_{wb}\approx
75.22$, and it is about 2.5 times the breaking time in a collisionless plasma.
When the condition $\nu\theta^{(0)}_{wb}>0.422$ is satisfied, the density
singularity does not arise anymore. This scenario of the evolution of the
electron density maximum is shown in Figure 8 for $\nu\theta^{(0)}_{wb}=0.5$.
In this case, the off-axis maximum, after its formation, first increases, but
then it decreases due to the strong damping of oscillations. Figure 5 shows
the spatial distributions of the momentum $P$ and the electric field $E$ of a
collisionless plasma at the moment of blow-up $\theta\approx 29.5$, when the
absolute maximum of density outside the origin became unbounded, and the
trajectories of a pair of neighboring particles intersect. Note that, due to
the structure of equations (5), their solution will remain odd functions of
coordinates if the initial data have this property. The initial data (32) is
exactly that, therefore, graphs of functions are shown only on the positive
semiaxis.

It is easy to see that in the vicinity of the density maximum, the velocity
component has a discontinuity in the derivative (a weak discontinuity), but
not in the function itself, while a strong discontinuity is formed in the
electric field function. It is these qualitative characteristics of $P$ and
$E$ that ensure the breaking (blow-up) of oscillations. Let us emphasize that
blow-up has the character of a ”gradient catastrophe”, i.e. the functions $P$
and $E$ themselves remain bounded.

Figure 5. Spatial distribution of momentum and electric field at the moment of
overturning in a relativistic calculation of a collisionless plasma
($\nu\theta^{(0)}_{wb}=0$).

Note that the above results of calculations by the particle method were
completely reproduced by additional calculations according to the splitting
scheme in Euler variables and according to the time hopping scheme (”leap -
frog”) in Lagrange variables (see details in [4]).

Let us present estimates for the breaking time of nonlinear oscillations for
some typical plasma parameters [9]. If plasma oscillations are driven by
electric field (32) with parameters $a_{*}=3.105,$ $\rho_{*}=4.5$, then the
blowup time in a collisionless plasma is $\theta^{(0)}_{wb}\approx 29.5$, as
follows from the calculations. Hence, in accordance with the numerical result
for the threshold value of the dimensionless collision frequency
$\nu\theta^{(0)}_{wb}\leq 0.422$, it follows that the breaking effect takes
place when the following inequality is satisfied $\nu\leq 1.43\cdot 10^{-2}$.
In a fully ionized plasma, the dimensionless frequency of electron-ion
collisions is determined by the formula [2]

(33) $\nu=Z\displaystyle{\sqrt{8}\over 3}\eta^{3/2}\ln\Lambda,$

where $Z$ is the ion charge number, $\ln\Lambda$ is the Coulomb logarithm, and
the parameter $\eta$ is equal to the ratio of the interaction energy of
electrons $e^{2}N_{0e}^{1/3}$ to their kinetic energy $T_{e}$

(34) $\eta=\displaystyle{e^{2}N_{0e}^{1/3}\over T_{e}}.$

Let in a rarefied, fully ionized plasma with an ion charge $Z=5$, whose
electrons have a density $N_{0e}=10^{18}{\rm cm}^{-3}$ and a temperature
$T_{e}=50$ $eV$, a laser pulse propagates with a wavelength $\lambda=1.24$
$\rm\mu m$ (frequency $\omega_{0}\approx 1.5\cdot 10^{15}{\rm s}^{-1}$),
duration $\tau\approx 36$ fs and dimensionless electric field amplitude
$a_{0}\approx 2.5$. If a laser pulse is focused by a cylindrical lens into a
line with a transverse size $L_{x}\approx 24$ $\rm\mu m$, then in the wake
wave behind the pulse, mainly plane one-dimensional oscillations of electrons
with parameters $a_{*}=3.1,$ $\rho_{*}=4.5$ are excited, which are close to
ones used in the calculations above. Note that when a laser pulse with a
moderately relativistic intensity $a_{0}$ (basically between 1 and 3)
propagates in a rarefied plasma, the condition for optimal excitation of
plasma waves ($\tau_{*}=2$) is approximately conserved, and the amplitude of
plasma oscillations is related to the laser field by the same relation
$a_{*}^{2}\approx 1.52\,a_{0}^{2}$, as in the non-relativistic limit. For the
reduced plasma parameters, from formulas (33), (34) we find that the
dimensionless collision frequency $\nu\approx 0.5\cdot 10^{-2}$ is less than
the threshold value $1.43\cdot 10^{-2}$. Therefore, in this case, the breaking
effect takes place, and the breaking time is approximately equal to
$\theta_{wb}\approx 33$, since the parameter $\nu\theta^{(0)}_{wb}\approx
0.15$. If we consider the propagation of a laser pulse with the above
parameters in a plasma with the same density, but with a temperature of
$T_{e}=20$ eV, then calculations using formulas (33), (34) give the
dimensionless collision frequency $\nu\approx 1.8\cdot 10^{-2}$ (parameter
$\nu\theta^{(0)}_{wb}\approx 0.52$) that exceeds the threshold. Therefore, in
this case, strong damping suppresses the breaking of plasma oscillations in
the laser pulse wakefield.

##  5\.  Conclusion

In the present work, the influence of collisions of electrons on the breaking
(blow-up) of plane nonlinear plasma oscillations is investigated analytically
and numerically. If there are no collisions of electrons in the plasma, then
the breaking of plasma oscillations occurs due to the formation of a maximum
electron density outside the oscillation axis, which increases with time and
goes to infinity after several periods of oscillations. It is shown that, when
electron collisions are taken into account, the breaking is also associated
with an increase in this off-axis density maximum, but in a collisional plasma
its growth occurs more slowly, and the breaking time increases with an
increase in the collision frequency. It has been established analytically and
numerically that for any fixed initial data there is a certain threshold value
for the frequency of electron collisions, upon reaching which the density
singularity does not arise. It is shown that, if the collision frequency above
the threshold value, the density maximum outside the oscillation axis after
its formation increases only for a certain time, but then it decreases due to
strong damping of oscillations, i.e. the solution is stabilized to the trivial
state.

## Acknowledgment

Partially supported by the Moscow Center for Fundamental and Applied
Mathematics.

## References

  * [1] Alexandrov, A.F., Bogdankevich, L.S., Rukhadze, A.A.: Principles of plasma electrodynamics, Springer series in electronics and photonics (Springer: Berlin Heidelberg, 1984).
  * [2] Bellan, P. M.: Fundamentals of Plasma Physics (Cambridge University Press, Cambridge, 2006).
  * [3] Booker, H. G.: Cold Plasma Waves (Springer-Verlag, Berlin, 2004).
  * [4] Chizhonkov, E.V.: Mathematical Aspects of Modelling Oscillations and Wake Waves in Plasma (CRC Press, 2019).
  * [5] Dafermos, C.M.: Hyperbolic Conservation Laws in Continuum Physics. The 4th Edition, (Berlin-Heidelberg: Springer, 2016).
  * [6] Davidson, R.C.: Methods in nonlinear plasma theory (New York: Academic Press, 1972).
  * [7] Dawson, J.M.: Nonlinear electron oscillations in a cold plasma. Phys. Review 113(1959) No. 2, 383–387.
  * [8] Engelberg, S., Liu, H., Tadmor, E.: Critical thresholds in Euler-Poisson equations, Indiana University Mathematics Journal, 50 (2001) 109-157.
  * [9] Frolov, A.A., Chizhonkov, E.V.: Influence of electron collisions on the breaking of plasma oscillations. Plasma Physics Reports 44 (2018) 398-404
  * [10] Ginzburg, V. L.: Propagation of Electromagnetic Waves in Plasma ( Pergamon: New York, 1970).
  * [11] Infeld, E., Rowlands, G., Skorupski, A.A.: Analytically solvable model of nonlinear oscillations in a cold but viscous and resistive plasma. Phys. Rev. Lett. 102(2009) 145005(1–4)
  * [12] Kahaner, D., Moler, C., Nash, S.: Numerical Methods and Software (New York: Prentice-Hall International, Inc., 1989).
  * [13] Rozanova, O.S., Chizhonkov, E.V.: On the conditions for the breaking of oscillations in a cold plasma. Z. Angew. Math. Phys. 72(2021), 13, doi: /10.1007/s00033-020-01440-3
  * [14] Rozanova, O.S., Chizhonkov, E.V.: On the existence of a global solution of a hyperbolic problem, Doklady Mathematics, 101(2020) 254-256
  * [15] Rozanova, O.S, Chizhonkov, E.V., Delova, M.I.: Exact thresholds in the dynamics of cold plasma with electron-ion collisions. AIP Conference Proceedings (2020) 2302, No. 1, 060012. doi: /10.1063/5.0033619
  * [16] Verma, P.S., J.K.Soni, S. Segupta, P.K. Kaw Nonlinear oscillations in a cold dissipative plasma. Physics of Plasmas 17 (2010) 044503(1-4).
  * [17] Sheppard, C. J. R.: Cylindrical lenses — focusing and imaging: a review [Invited]. Applied Optics. 52(2013) 538-545
  * [18] Schultz, M.H.: Spline Analysis (New York: Prentice-Hall International, Inc., 1973).