Diagnostics of a gradient catastrophe for a class of quasilinear hyperbolic systems

1 Introduction

The base of many mathematical models is formed by quasilinear systems of hyperbolic equations leading to invariants, i.e., representable in the form

The typical formulation of problems for system (1.1) is the Cauchy problem

If in this case the functions f_k increase not very fast with the growth of r = (r₁, r₂,…,r_n)^T, for example, for any x, t, r we have

then the solutions r(x,t) to systems of such type remain bounded (in their absolute value) for any t. However, their derivatives can infinitely increase in absolute value when ξ_k(x,t,r) essentially depend on r (see [22], p. 82).

The effect of formation of infinite derivatives when the solution to a system of quasilinear hyperbolic equations is bounded is usually called the gradient catastrophe.

There exist classes of problems where the appearance of a gradient catastrophe means the application limit for the used mathematical model. In other words, when beginning with smooth initial data (1.2), for example, continuously differentiable ones, we get a discontinuous solution to system (1.1), this indicates the loss of physical sense of original functions. Let us clarify these arguments on the example of plasma oscillations.

Mathematical modelling of processes in collision-free cold plasma most often uses two approaches of Lagrange and Euler, i.e., the method of particles allowing one to trace their individual trajectories and a hydrodynamic description on the basis of partial differential equations (see, e.g., [1, 10, 15]). The effect of breaking of oscillations is well known [9] within the frameworks of both these approaches. In the first case the criterion of oscillations breaking is the intersection of electron trajectories, in the second case this is the fact that the function describing the electron density tends to infinity. Monograph [24] (see also Appendix in [13]) presents a rigorous justification of the appearance of a singularity under intersection of trajectories of particles.

The physical meaning of the breaking effect is quite simple. Electrons are moved from the equilibrium position due to one or other action, i.e., an electric field of charge separation is formed artificially. Further, under the action of Coulomb forces, particles tend to return to their equilibrium positions. However, moving with acceleration like a pendulum, electrons regularly pass past the specified position. Since the equations are nonlinear, the difference between the frequencies of oscillations of different particles becomes essential with time which leads to intersection of close trajectories. If two particles occupy the same position in space and time, then further tracing of their movement requires more complicated model than the classical electrodynamics can present because the infinite concentration of the electric charge requires special interpretation.

It is worth to note that the situation described above differs essentially from modelling in gas dynamics where discontinuous solutions are physically natural. The corollary from this nature is the fact that exact or approximate solution of the Riemann problem (Cauchy problem with a piecewise-discontinuous initial data) forms a basis of most modern numerical solution algorithms for the corresponding formulations of problems [18]. However, in problems related to plasma oscillations the formulation of the Riemann problem has no physical sense because a discontinuous initial function of the electric field already assumes an infinite charge concentrations at the points of discontinuity.

In the present paper we propose a two stage approach to detect the appearance of a gradient catastrophe.

At the first stage we propose to associate spatial gradients of initial functions (1.2) with the initial conditions for a special system of ordinary differential equations. This leads to consideration of a more simple Cauchy problem where we have to determine the condition of existence of a vertical asymptote (‘blow-up’ property). Since the auxiliary solution determines the local gradients of original problem (1.1)–(1.2), the fulfillment of this property causes a gradient catastrophe, i.e., the loss of continuity of the solution to the original problem. This stage essentially differs from other known approaches to the study of quasilinear hyperbolic systems such as the extended system ([22], p. 32), the majorant system ([22], p. 62), the method of nonlinear capacity [20], etc. It is worth to indicate a certain similarity with the ideas of [17], although the analysis proposed in the paper seems to be more simple. The first stage focused on the study of properties of the auxiliary Cauchy problem for a system of ordinary differential equations is in essence related to the study of the original formulation in Euler variables.

The second stage is aimed to more ‘fine’ analysis of the original problem in the equivalentLagrange formulation. It is applied under the assumption that, according to the results of the first stage, the initial data do not lead to a quick gradient catastrophe, i.e., to an breaking of oscillations appearing already at the first stage. The essence of the stage is in construction of a uniformly suitable asymptotic expansions of the original system in a weakly nonlinear approximation. Here we derive the dependence of the frequency of oscillations on the initial amplitude for equations describing trajectories of particular particles. Taking this into account, we make the conclusion concerning intersections of trajectories of particles and may estimate if necessary the asymptotics of this event depending on the time of oscillations.

2 Scalar Burgers equation and the system generalizing it

We begin the consideration of the proposed approach with the well-known problem

and analyze the conditions causing a gradient catastrophe for it.

Let us consider the so-called axial solution to the original equation, i.e., the sample solution of the form

The notion of axial solution was introduced in [3] for nonlinear problems describing laser-plasma interactions and possessing axial symmetry. Such solutions have local linear dependence in a neighbourhood of the initial point for the spatial coordinate. In this case we certainly do not assume axial symmetry, but for the sake of convenience we keep this name for real solutions having the indicated structure.

It is not difficult to see that the factor W(t) depends on time and satisfy in this case the following ordinary differential equation:

which should be supplemented with the real initial condition

Now we clarify the conditions for existence and uniqueness of the solution to Cauchy problem (2.3)–(2.4). Since the problem formulated here has the analytic solution

the condition we are interested in has the form ⩾ ≥ 0. Otherwise, β < 0 and the solution to problem (2.3)–(2.4) has a vertical asymptote for t = − 1/β (‘blow-up’ property), which means the appearance of a gradient catastrophe for original problem (2.1).

This conclusion is in good agreement with the known fact (see, e.g., [19]) that the nondecreasing behaviour of the smooth function U⁰(x) is necessary and sufficient for the existence of a smooth solution to (2.1) for all t > 0. Note that in this simple case the class of smoothness can be extended up to continuous functions, although we are primarily interested in classical (i.e., continuously differentiable) solutions.

Completing the example, we specify that just the behaviour of the gradient of the solution interests us because just the ‘blow-up’ property for the gradient determines the formation of a discontinuous solution from smooth initial data of the problem.

Now we consider the following system generalizing Burgers’ equation:

and endowed with the corresponding initial conditions

Assuming axial solutions

to system (2.5) often called quasi-orthogonal [23], we write down the system of ordinary differential equations

and supplement it with the real initial data

Problem (2.7), (2.8) can be solved explicitly. The case α = β is reduced to Burgers’ equation considered above for which the condition β < 0 forms a discontinuous solution (gradient catastrophe) from an arbitrarily smooth initial function.

Consider the case α ≠ β. The equality D(t) = W(t) + α − β holds for different initial conditions (2.8). Excluding the function D(t) from system (2.7) and applying the change y(t) = W(t) − (β − α)/2, we come to the problem

Its analytical solution is well known [16] and given by the formula

Since the absolute value of the hyperbolic tangent does not exceed one, for positive values of the argument btthe denominator in the formula y(t) vanishes if the inequalities β < 0 and α < 0 hold simultaneously. Negative values of the initial data are sufficient for the ‘blow-up’ property in problem (2.7)–(2.8). Since the parameters β and α characterize gradients of the initial functions in (2.6), the fulfillment of the condition

for at least one value −∞ < x < + ∞ leads to a gradient catastrophe in problem (2.5)–(2.6). A clear illustration of a discontinuous solution formation in the formulation considered here was presented in [23].

3 Plane non-relativistic electron oscillations

We consider the following more informative example of a Cauchy problem for equations in dimensionless form describing plane one-dimensional non-relativistic oscillations of electrons in cold plasma [4]:

Here V is the velocity of electrons, E is the functions specifying the electric field.

Let us consider the axial solutions

to system (3.1) and study their properties.

In this case, factors depending on time satisfy the following system of ordinary differential equations:

Supplement these equations with arbitrary real initial conditions

and find out the existence and uniqueness conditions for the solution to Cauchy problem (3.2)–(3.3).

Note that this Cauchy problem is not trivial because it admits both regular 2π-periodic solutions (for example, for small α and β), and solutions having singularities in a finite time interval (‘blow-up’ solutions). In the general case, known results cannot determine an exact boundary between sets of initial data generating solutions of different types even for polynomial right-hand sides [14]. Therefore, another view on the system of differential equations considered here seems to be useful. We have the following assertion [4].

4 Plane relativistic electron oscillations

We consider a more complicated example related to the simulation of plane relativistic oscillations in cold plasma. The corresponding system of equations have the following dimensionless form [11]:

Here V and P are the velocity and momentum of electrons, respectively, E is the function specifying the electric field. Formulating the Cauchy problem, this system should be supplemented with the initial conditions

Taking into account the approximate representation of the relativistic dependence of the momentum of the velocity

we get that axial solutions to system (4.1) coincide with axial solutions to equations (3.1) up to small summands of third order with respect to the parameter x. In other words, the Euler stage of diagnostics of the gradient catastrophe leads in this case to condition (3.5) again. Therefore, for small intensity of the initial electric field and small initial motion velocity of electrons we can easily eliminate the formation of a fast gradient catastrophe at the first period of oscillations.

Below we analyze the cause of the new ‘slow forming’ gradient catastrophe. Transform system (4.1) from Euler variables to Lagrange ones. We have

where d /dt = ∂ /∂ t + V ∂ /∂ x is the total derivative in time.

Taking into account equality (3.7), in the case of plane one-dimensional oscillations, i.e., R(ρ^L, t) = E(ρ^L, t), basic system of equations (4.4) in Lagrange variables takes the form

If we suppose that the amplitude of oscillations is sufficiently small, equations (4.5) become weakly nonlinear and their approximate solutions can be constructed using the methods of the theory of perturbations [2]. We present here a brief derivation of the corresponding analytic formulas.

For definiteness sake, consider the following functions as initial conditions (4.2):

Such perturbation of the electric field initiating the oscillations is typical for a short laser pulse passing the plasma and having a Gaussian spatial distribution of intensity (see [6, 12]). Here a_* and ρ_* are parameters of the pulse.

Taking into account the approximate representation of relativistic dependence (4.3) of the momentum on the velocity and excluding the displacement R from system (4.5), we get the following equation for the velocity V:

Supply it with the initial conditions from (4.6)

We suppose the amplitude of oscillations A(ρ^L) is small, in this case we are interested in an approximate (asymptotic) solution to (4.7)–(4.8) uniformly bounded with respect to the variable t and differing from the exact solution by a value of third order of smallness, i.e., O(A³(ρ^L)). In this case we can assume

where, as above, ρ^L and x are related by the ratio ρ^L = x − E⁰(x).

The direct expansion of the solution to (4.7)–(4.8) with respect to powers of small parameter

contains, as usual, resonance summands of the form t sin t, i.e.,

Therefore, it is not applicable for long time intervals.

Using the standard methodology [2], one can easily obtain a correction term to the main frequency of oscillations dependent on the amplitude. As the result, the required bounded solution has the form (see [11])

Returning to formulation (4.7)–(4.8), we obtain the dependence of the frequency of oscillations on the spatial distribution of the initial amplitude A(ρ^L):

In this case, in contrast to the direct expansion, summands of third order of smallness do not increase in time. Apply the result obtained here to relations (4.5). In particular, we obtain that the trajectory of a particle has the following form up to terms of third order of smallness:

where ρ^L is the Lagrange coordinate of a particle in the equilibrium position not leading to generation of electric field. Since the initial amplitude of oscillations A(ρ^L) is not a constant value, trajectories of some neighbouring particles earlier or later have to intersect each other due to the difference in frequencies of oscillations.

This example illustrates the Lagrangian stage in prediction of the gradient catastrophe. In contrast to Euler stage, here we first analyze equations, but not initial conditions. For initial conditions we use the assumption that they do not generate the required effect on the first period of oscillations.

Let us describe the results of numerical experiments demonstrating the formation of a ‘slow’ gradient catastrophe. In order to satisfy condition (3.4), assume in (4.6) the following values of parameters: a_* = 2.07, ρ_* = 3.0.

Now we consider electron oscillation continuing for more than one period. It is interesting to track their formation with the use of numerical algorithms [11]. An important and convenient object of observation is the function of electron density

because just this function becomes singular at the gradient catastrophe.

In Fig. 1 red color indicates the spatial distribution of the electron density N at the initial time moment. In this case the density function is symmetric with respect to the origin, therefore, only the right-hand half of the graph is presented. The excess of positive charge at the origin leads to motion of electrons towards the center of the domain, which generates another distribution of the density function after a half of the oscillation period. This distribution is also presented in Fig. 1 (black color). Note that the electron concentration at the center of the domain might exceed the equilibrium (background) value equal to one many times. Fixed parameters lead to oscillations of small intensity when the oscillation amplitude is only about 10 times greater than the background value. If nonlinear plasma oscillations preserved their spatial form in time, then the electron densities presented in Fig. 1 were alternated regularly in each half of the period generating a strictly periodic sequence of extrema with the same amplitudes at the center of the domain.

Figure 1

Spatial distributions of electron densities varying each half of the period under regular oscillations, the maximum at the origin (black), the minimum at the origin (red).

However, we can see two trends in the process of oscillation. The first trend consists in the fact that off-axial oscillations are little ahead in phase than those on the symmetry axis (for x = 0) and this phase shift increases from period to period. The second trend is more visual. A gradual formation of the absolute maximum of density located out off the axis is observed and this maximum is comparable to the axial one. Figure 2 is a good illustration of these assertions. It shows variations of the electron density at the origin in time with red color and the dynamics of maximum values over the domain with black color. At first the oscillations are regular, i.e., global maxima and minima of the density over the domain alternate through a half of the period and are located at the origin. After the seventh regular (central) maximum at t ≈ 42.2 a new structure appears, this is an off-axial maximum of the electron density. The regular oscillations are still seen in this case in a neighbourhood of the origin. In its turn, the off-axial maximum increases in magnitude about two times at the time moment t ≈ 48.8 and in the next period, i.e., t ≈ 55.1, a singularity of the electron density appears at its position.

Figure 2

Electron density dynamics, the maximum over the domain (black) and at the origin (red).

Figures 3 and 4 are presented for greater clarity of the off-axial maximum of the electron density. Figure 3 shows the spatial distribution of the density at t ≈ 48.8 when it has been completely formed and became comparable in absolute value to the regular (axial) maximum. Recall that in this case the density function is even relative to the origin and hence we present only the right-hand half of the graph. The graph of the density shown in Fig. 3 is a consequence of the distributions of the velocity V and the electric field E presented in Fig. 4. In their turn, the functions of velocity and electric field are odd relative to the origin in this case and we also present only right-hand halves of their graphs. Note that the velocity function tends to a jump of derivative in a neighbourhood of the density maximum and the function of the electric field takes a step form. Just these qualitative characteristics ensure the breaking of oscillations at the time t ≈ 55.1. It is worth noting that the breaking has a gradient catastrophe character, i.e., the functions V and E themselves remain bounded in this case.

Figure 3

Spatial electron density distribution at the moment of formation of the second off-axial maximum.

Figure 4

Spatial distributions of the velocity and electric field at the moment of formation of the second off-axial maximum.

The example of the diagnostics considered above has a pronounced combined structure, the Lagrange stage of analysis of equations naturally complements the Euler analysis of gradients of initial conditions.

5 Cylindrical non-relativistic electron oscillations

Below we consider the example when Lagrange’s analysis is more informative in the case of a gradient catastrophe than Euler’s one. The corresponding system of equations describing cylindrical one-dimensional non-relativistic oscillations of electrons in cold plasma has the following dimensionless form [12]:

Here V is the velocity of electrons, E is the function specifying the electric field.

Posing the Cauchy problem, we supplement this system with the initial conditions

As before, at the Euler stage of diagnostics we consider properties of axial solutions [21]. Recall that an axial solution to equations (5.1) means a real solution of the form

It is not difficult to see in this case that the factors dependent on time satisfy the following system of ordinary differential equations

Supplement the obtained equations with arbitrary real initial conditions

and clarify the conditions for existence and uniqueness of the solution to Cauchy problem (5.1)–(5.2). We have the following assertion [21]

Theorem 5.1

The necessary and sufficient condition for the existence and uniqueness of the smooth periodic solution to Cauchy problem(5.3)–(5.4)is the validity of the inequality

α<1/2.(5.5)

Theorem 5.1 implies that the sufficient condition for a gradient catastrophe for problem (5.1)–(5.2) is the validity of the inequality

∂E0(ρ)∂ρ⩾12(5.6)

for the initial function E⁰(ρ) at least for one value of ρ. This completes the Euler stage of diagnostics; further we assume that the inequality opposite to (5.5) holds true and we have to proceed to Lagrangian stage.

Transform equations (5.1) in Euler variables to equations relative to Lagrangian variables, i.e.,

dV(ρL,t)dt=−E(ρL,t),dE(ρL,t)dt+V(ρL,t)E(ρL,t)ρ=V(ρL,t)(5.7)

where d /dt = ∂ /∂ t + V ∂ /∂ ρ is the total time derivative.

In this case the function R(ρ^L, t) determining the displacement of the particle having the Lagrange coordinate ρ^L and trajectory ρ(ρ^L, t),

ρ=ρL+R(ρL,t)(5.8)

satisfies the equation

dR(ρL,t)dt=V(ρL,t).(5.9)

In this case of cylindrical oscillations there is no simple relation between the displacement and electric field (2.7), therefore, we need additional calculations. Expressing the velocity V through the displacement R by formula (5.9), write the second equation of (5.7) in the form

(ρL+R)dEdt+EdRdt=(ρL+R)dRdt.(5.10)

Equation (5.10) possesses the first integral

(ρL+R)E=12(ρL+R)2+C

where the constant C is determined from the condition that the electric field vanishes under the absence of displacement of particles. This gives the following expression for the electric field:

E(ρL,t)=12(ρL+R(ρL,t))2−(ρL)2ρL+R(ρL,t).

Using the obtained formula, we can rewrite equations (5.7), (5.9) in the more convenient form

dR(ρL,t)dt=V(ρL,t),dV(ρL,t)dt=−12(ρL+R(ρL,t))2−(ρL)2ρL+R(ρL,t).(5.11)

Relation (5.11) implies the following useful equation relative to R(ρ^L, t):

d2Rdt2=−12(ρL+R)2−(ρL)2ρL+R(5.12)

where the Lagrange coordinate ρ^L characterizing a particle participates as a parameter. This equation is not new, it was derived in [9] from qualitative reasons where it was also used for analysis of breaking of electron oscillations.

Using standard technique [2], we can easily obtain a correction to the principal frequency of oscillations dependent on the amplitude, i.e.,

ω=1+R02(ρL)12(ρL)2,R0(ρL)=R(ρL,0).

As the result, the required bounded solution to (5.12) has the form (see [12]):

R(ρL,t)=R0(ρL)cos⁡ωt+O(R03(ρL)).(5.13)

In this case all summands of the third order of smallness do not increase in time as this was in the case of plane relativistic oscillations.

The dependence of the frequency on the initial amplitude in formula (5.13) has a more complicated character in comparison with formula (4.10). In particular, it implies that under linear dependence of the initial amplitude on values of the Lagrange coordinate ρ^L we can get no intersection of trajectories of initial particles. On the other hand, if for initial conditions (5.2) for system (5.1) we take the following analogue of (4.6):

V0(ρ)=0,E0(ρ)=a∗ρ∗2ρexp−2ρ2ρ∗2

then numerical experiments from [21] convincingly demonstrate the occurrence of a ‘slow’ gradient catastrophe (see Figs. 5 and 6).

Figure 5

Spatial distribution of electron density, the regular maximum on the axis of symmetry (black), the first off-axial maximum (red).

Figure 6

Spatial distributions of the velocity and electric field at the moment of gradient catastrophe.

If, similar to the case of plane oscillations, we take the parameters of initial conditions not leading to a gradient catastrophe in the first period of oscillations, then each period will be characterized by the maximum of electron density appearing strictly on the axis of symmetry and not changing its magnitude. This distribution of density is shown in Fig. 5 with black color and corresponds to the values a_* = 0.365, ρ_* = 0.6.

After several periods of oscillations (in our case, six) an additional off-axial maximum of the density function is formed at some distance from the axis and it is a predecessor of intersection of trajectories. The density distribution having the first off-axial maximum is marked in Fig. 5 with red color.

After some time this distribution disappears, but after about a period it appears again with a greater maximal value exceeding the regular (axial) maximum almost twice. In contrast to the maximum of density positioned on the axis, the off-axial density maximum quickly increases after its appearance from period to period. Thus, with regular oscillations with the minimum and maximum of density located strictly on the axis, a gradient catastrophe of the solution occurs in a period at the off-axial maximum place. It is hard to construct such graph (with infinite values of the density), therefore, Figure 6 presents spatial distributions of the velocity and electric field corresponding to the breaking moment. The step-wise function of electric field and the kink of velocity derivative visually illustrate the formation of singularity in the electron density. Recall that in this case of axially symmetric oscillations we have the formula

N=1−1ρ∂∂ρρE.

6 Conclusions

A two-stage (combined) approach is proposed for diagnostics of a gradient catastrophe in equations of hyperbolic type written with respect to invariants.

The first stage which is convenient to call Euler is aimed, first of all, to analysis of initial data and consists in the study of the ‘blow-up’ property of the auxiliary system of ordinary differential equations. Such approach is traditional in the diagnostics of a gradient catastrophe for hyperbolic systems though the technical aspects may be essentially different. The present research demonstrates the application of the first stage both for model formulations (Burgers equation, quasi-orthogonal system), and for problems having meaningful physical sense.

The second stage essentially supplements the first one, but is less traditional. It is convenient to call it Lagrangian. At this stage, first of all, we analyze equations written relative to Lagrange variables. Here, using the methodology of perturbations for construction of solutions being uniformly suitable in long time intervals, we can reveal the dependence of the frequency of oscillations on the amplitude generated by the initial conditions. This predicts possible intersection of trajectories of particular particles, which also means the occurrence of a gradient catastrophe. Two stages of diagnostics are applied to hyperbolic systems having more complex structure of solutions, namely, plane relativistic and non-relativistic electron oscillations and also axially symmetric non-relativistic electron oscillations in cold plasma.

In essence, the proposed approach is based on the equivalent interpretation of oscillation breaking, i.e., in Lagrange variables this is an intersection of trajectories of particles and in Euler variables this is a singularity of the electron density function. The results of the study can be generalized to more complex formulation including the influence of ion dynamics of oscillating processes and to cases of two and more spatial dimensions [5, 7].