Semi-discretization in time for Nonlinear Zakharov Waves Equations

In this paper we construct and study discretizations of a nonlinear Zakharov-wave system occurring in plasma physics. These systems are generalizations of the Zakharov system that can be recovered by taking a singular limit. We introduce two numerical schemes that take into account this singular limit process. One of the scheme is conservative but sensible to the polarization of the initial data while the other one is able to handle ill-prepared initial data. We prove some convergence results and we perform some numerical tests.


Setting of the problem and motivation
The strong Langmuir turbulence phenomenon is well-known problem in plasma physics [13].The key mechanism of this physical phenomenon is governed by the so-called Zakharov equations.These equations describe the slowly varying motions of the complex-valued envelope of the Langmuir electric field E 0 which nonlinearly interacts with large-scale fluctuation n 0 through the ponderomotrice force.They read in a dimensionless scalar form    (1.1) This system can be derived from the Klein-Gordon-wave system [4] or directly from the Euler-Maxwell system [15] using a small parameter ε = 1 ωpeT where ω pe is the electronic plasma pulsation and T a characteristic time of the pulse.The size of ε governs the regime of evolution of the solution.System (1.1) is only an approximation that can be sometimes not accurate enough to describe the full physics corresponding to the large variety of plasmas created by laser pulses.In [2], an intermediate system (between (1.1) and Euler-Maxwell in the hierarchy) has been introduced.It is less stiff than Euler-Maxwell but more precise than (1.1).In particular, it takes into account more deeply the oscillations of the electric field and allows "larger" values of ε.Its range of validity is expected to be wider than that of Zakharov.It reads, in a dimensionless form : endowed with the initial data System (1.2) is intermediate between Zakharov system (1.1) and Klein-Gordon-waves system, namely (1.4) Indeed, because of the presence of the parameter ε in (1.4), the electric field E is highly oscillatory whereas n is low frequency.Mathematically, these oscillations can be described through the change of function E(x, t) = E(x, t)e − it ε + c.c..For the low frequency ion fluctuation we set n(x, t) = N (x, t).By injecting (E, n) in system (1.4) and neglecting oscillatory terms in the nonlinear term of the second equation, one obtains system (1.2).Moreover in [2], it is shown that for any initial value ) with s large enough and ε|E 1 | → 0 as ε → 0, there exists an unique solution to (1.2) on a maximal existence time [0, T ε [ such that T ε > 0. Furthermore, if (E 0 , n 0 ) denotes the maximal solution of Zakharov equations (1.1) with initial value (E 0 , N 0 , N 1 ) defined on [0, T 0 [, then on one hand, lim ε→0 inf T ε ≥ T 0 , and on the other hand, for any T < T 0 , (E, N ) converges to (E 0 , n 0 ) on [0, T ] in C([0, T ]; Remark 1.1.The convergence result is easier than the one established in [4].Indeed the high frequency term in the non-linear force ∆|E| 2 have been filtered out.
In the non-compatible case (that is when oscillates, with a O(1) amplitude, at frequency 2 ε around the value ∂ t E 0 .In the compatible case (that is when E 1 = − i 2 ∆E 0 + E 1 + i ε N 0 E 0 ), the quantity E − E 0 oscillates at frequency 2 ε and these oscillations are of order of magnitude ε.In the last section, we will consider numerically larger initial data (cases where ∂ t E = O(1/ε) and show how to handle them.We can divide E into two components F and G involving the propagators e i ε (1∓ √ 1−ε∆)t respectively.Analogously, we decompose N into H + I.The functions (F, G, H, I) are solution of the split-system, namely (1.5) The initial data for F, G, H, I are recovered from E 0 , E 1 , N 0 , N 1 (see [2] for details).Note that this splitting is exact.Finally, system (1.2) owns the two following invariants where ψ is defined by ∂ t N = ∆ψ.
• Setting of the problem : It is clear that (E, N ) converges to (F, N ), as ε → 0, while G is the oscillatory part of E. The period of these oscillations is πε while their amplitude is of order ε.It implies that the more G oscillates, the more the contribution of G in E is negligible.From the numerical point of view, it means that we need a small time step in order to compute G precisely while it is not necessary since its contribution in E = F + G is small.However, for ε large enough it is essential to compute both the average and the oscillatory parts of the solution in view of the error estimate established in [4].
• Motivation : The aim of this paper is to propose a numerical approach for solving problem (1.2)-(1.3).
In fact, we investigate a numerical scheme which is precise and not expensive for any ε.In particular we study two schemes of Crank-Nicolson type for (1.2) using the splitting (1.5) for the associated discrete operators.We show that the non-centered in time scheme (called S 1 ) owns dissipative and polarization properties contrary to the centered in time scheme (called S 2 ).We will precise later what is polarization.Moreover we numerically check the error estimate established in [4] and confirm that system (1.2) compute precisely the propagation of short or ultra-short waves through a plasma.For a precise description of short or ultra-short waves, one can see [1].

Numerical schemes
We recall that the Crank-Nicolson scheme for (1.1) has been introduced by Glassey [6].It reads where E m j ≈ E 0 (t = mδt, x = x j ), n m j ≈ n 0 (t = mδt, x = x j ) and D + , D − are the forward and backward non-centered finite difference operators respectively.Several authors have studied numerical schemes for the Zakharov system (see for example [3], [14] and recently [7]).Analogously, we introduce the two following schemes for (1.2) : ⋆ the non-centered in time scheme The aim of this paper is to provide some mathematical results for semi-discretized version of these schemes, that is when D + D − is replaced by the Laplace operator ∆.The schemes then read : ⋆ the non-centered in time scheme (S 1 ) We prove existence, uniqueness and convergence of the discrete sequences (E m , N m ) for both schemes for fixed ε.We show that scheme (S 2 ) is conservative while scheme (S 1 ) is dissipative.Moreover scheme (S 1 ) polarizes the solution whatever the initial data we take.We will precise later on what we mean by polarize.The aim of this paper is to investigate the effect of the dissipation and the polarization phenomenon on the discrete solution given by (S 1 ).
The paper is organized as follows.In section 2, we state existence, uniqueness and convergence results of the discrete solution given by both schemes.We introduce the splitting of schemes (S 1 ) and (S 2 ).
Then we give some details of the proof of theorem 2.1, that is related to scheme (S 1 ).For scheme (S 2 ) we just check the key points of theorem 2.2.
In section 3, we discuss the behavior of the discrete solution given by both schemes.We study dissipative and polarization phenomenon of scheme (S 1 ).We compare the solution of (1.2) given by both schemes with the solution of Zakharov system (1.1).At the end we illustrate the transparency property [11] of system (1.4) established in [4].
Our mathematical result mainly depend on the study of the linear part of the schemes.Therefore a lot of technical tools are imported from [5].The main difference between this paper and [5] is the well known loss of derivative involve in the Zakharov system.This loss of derivative give rise to weaker convergence and stability results than in [5].

Theoretical results
Following [2] in the continuous case, we introduce an exact splitting for our schemes.Using these splitting, we prove existence, uniqueness and convergence results for the discrete solution for fixed ε.

Splitting of the schemes
We recall that in the continuous case, one divides E into two components F and G having the modes of propagation e i ε (1∓ √ 1−ε∆)t respectively.Analogously, we write N = H + I.The functions (F, G, H, I) are solution of the split system (1.5).One can see that G contains all the oscillatory part of the solution and the order of magnitude is of order |G(0)|.Moreover (F, N ) converges to the solution of (1.1) with the suitable initial data.The aim of this section is to present a similar splitting for (S 1 ) and (S 2 ) which will have the same behavior like the continuous case.In analogy with [5], we look for a splitting for (S 1 ) and (S 2 ) of the form where f m+1 , h m denote the right-hand side of (S 1 ) or (S 2 ).For the initial data of (2.1) we refer to [5].The operators α j , β j are given by α δt(1+Yj ) where X j , Y j are the roots of the characteristic equations of (S 1 ) and (S 2 ) for j = 1, 2. The nature of the operators α j and β j give us some information on the behavior of the schemes.We investigate the properties of each scheme separately.
• Splitting for (S 1 ) : One gets the following result (see [5] for details) : They are given by and .
Finally operators γ and ν are given by The proof is the same than the one of proposition 4.1.in [5].We omit it.Moreover it is clear that the Crank-Nicolson scheme (S 1 ) satisfies the two following approximate conservation laws : and where • Splitting for (S 2 ) : We have a similar result for (S 2 ).The function ).However, the operators (X 1 , X 2 ), which are solution of the following characteristic equation are different from the previous scheme.They read .
The discrete conservative quantities corresponding to (1.6) and (1.7) read and Remark 2.2.Operators X 1 and X 2 are unitary.It implies that iα 1 and iα 2 are skew-adjoint and generate unitary groups on any H s .Of course, it is not the case for (S 1 ).
We are now able to state the main theoretical results of this paper.

Existence and convergence
We give the two following results related to (S 1 ) and (S 2 ) respectively.For (S 1 ) one gets : ) and ε fixed in ]0, 1[.There exists δt 1 > 0 depending on initial data and ε such that for any δt ≤ δt 1 there exists T > 0 and a unique solution , where T = M δt.
• Let (F δt , G δt , H δt , I δt ) be the piecewise constant functions which value at time t ∈ [mδt, (m + 1)δt[ is (F m , G m , H m , I m )and (F, G, H, I) be the solution of (1.5) with the initial data F (0) = F 0 , G(0) = G 0 , H(0) = H 0 , I(0) = I 0 and T ∞ its existence time, then for ε fixed in ]0, 1[ and for all T < T ∞ (with For (S 2 ), one gets : . There exists δt 2 > 0 depending on the initial data such that for any δt ≤ δt 2 and ε > 0 there exists an unique solution G, H, I) be the solution of (1.5) with the initial conditions Remark 2.3.Note that these results are not "uniform" with respect to ε.This means that one can not perform a priori both limit ε → 0, δt → 0 contrary to the result of [5].This is certainly linked to the fact that, using the same kind of method in [2], the authors where not able to perform the limit ω → ∞ and c → ∞ in the Klein-Gordon-Wave system, namely ).Note that this problem has been solved in [12] using transparency type properties also used in [4].The transposition of such kind of properties in a semi-discrete framework is out of the scope of this paper.

Proof of theorems 2.1 and 2.2
We will give some details for the proof of theorem 2.1 related to (S 1 ).For theorem 2.2, we will just check the main steps.
Proof. of theorem 2.1.• We first prove existence and uniqueness of the discrete solution given by (S 1 ).We will use a fixed point procedure as in the continuous case.We write a Duhamel formula (2.8)-(2.9)with the four groups ∆ and defined below.The main point is to obtain uniform estimates with respect to δt.
We denote by (U δt 0 (t), U δt and (E δt (t), N δt (t)) the piecewise constant functions which value on [mδt, (m + 1)δt[ is (E m , N m ).One has where S δt 0 (mδt)f is the solution of and S δt 1 (mδt)g the one of The equivalent of (2.7) is where T δt 0 (mδt)p is the solution of and T δt 1 (mδt)q the one of In view of these notations, the discrete Duhamel formulation for (S 1 ) reads ) We introduce the spaces where c(K) is a constant to be determined later on.The time T is given by T = M δt, so that for fixed T , δt → 0 is equivalent to M → +∞.For (E m+1 ) 0≤m≤M−1 and (N m+1 ) 0≤m≤M−1 we construct the two following series (T (E m+1 , N m+1 )) 0≤m≤M−1 and (S(E m+1 , N m+1 )) 0≤m≤M−1 , namely Finally we define the application In view of the fixed point theorem, we will show that H is a contraction on B R .
In order to perform these estimate, we need to control the norm of the following operators : •X m+1 Proposition 2.5.Let K such that δt Remark 2.6.If S δt 0 (mδt)f is a solution of (S f 1 ), it implies that for all α and any regular function f , ∂ |α| S δt 0 (mδt)f is a solution of (S f 1 ) too.
For operators Y m+1 Proof.: We recall that operators Y 1 and Y 2 are given by ) corresponding to the nonlinear terms of (2.8)-(2.9),one gets : is given in the Fourier space by γ1 ( Xm+1 • Estimate of Q : One has By getting c(ε ii) We first give the relationship between Y m+1 and Y 1 , Y 2 , T δt 0 , T δt 1 where T δt 0 (mδt)p and T δt 1 (mδt)q are the groups corresponding to H + I introduced in section 1.As above, one gets In view of these previous formulas, operator ν1 ( Ŷ m+1 ) is given by So we get thanks to the expression of F (T δt 0 (mδt)) and F (T δt 1 (mδt)).This leads to result ii) of lemma 2.1.
In view of the fixed point procedure and the previous estimates of the operators that occur in the Duhamel formulation (2.8)-(2.9), the local existence of the solution of (S 1 ) in B R is straightforward.Note the well-known loss of derivative encountered in the Zakharov system is handled by i) of lemma 2.1.The consequence is a loss of uniformity with respect to ε.
For the convergence of the solution of (S 1 ) towards the continuous solution of (1.2)-(1.3),we write the Duhamel formula for (E(t), N (t)) solution of (1.2)-(1.3)between 0 and (m + 1)δt : and Substracting these identities from (2.8)-(2.9)respectively and taking the H s and H s−1 norms leads to and (2.13) In order to control the r.h.s. of (2.12), we decompose the integral into three terms, namely We decompose the integral term of (2.13) into three similar terms as In view of these previous notations, inequalities (2.12) and (2.13) are bounded by and (2.17) In the same way, as On the other hand, as the groups U 0,1 and V 0,1 are strongly continuous, For the first terms in the r.h.s. of (2.20) and (2.21), we write We recall that H s−1 (R d ) and H s (R d ) are algebras for s > d 2 + 1.Now according to the uniform estimates of (E, E δt ) and (N , N δt ) on [0, M δt], one has
with T = M δt.This finish the proof of theorem 2.1.
For scheme (S 2 ), in view of the previous results, the proof of the first point of theorem 2.2 is based on the following proposition : Proposition 2.8.Operators X 1 and X 2 corresponding to scheme (S 2 ) are unitary on H s (R d ).
Proof.Thanks to the formulas of X 1 and X 2 , it is clear that the groups X 3 Numerical results and comments

Diffusion and oscillations
The aim of this part is to investigate numerically the behavior for the Zakharov-wave system.We study the behavior of the solution of each scheme (S 1 ) and (S 2 ) when ε is small enough.Let us first recall the splitting for the continuous problem (1.5) (we do not rewrite the two last equations which do not depend on ε) : It is clear that It implies that (F, H + I) converge to the solution of Zakharov system (1.1) whereas G is an oscillatory part of E, solution of (1.2).Moreover the amplitude of these oscillations are of order |G(0)| H s while the frequency is 2 ε .The question is : how the discrete solution of each scheme behaves?To answer this question, we recall the discrete splitting of each scheme (2.1) (we rewrite only the two first equations) where δt(1+Xj ) for j = 1, 2. We do the asymptotic expansion of operators X 1 and X 2 .We investigate the discrete solution given by (S 1 ) for two values of ε.If δt is small enough with respect to ε, for example ε = δt 1/2 , one gets at the first order, when δt → 0, It follows that the leading part of G m can be written as It means that G m oscillates at frequency 2 ε and the amplitude of these oscillations damped exponentially in time for regular initial data G(0).On the other hand, if δt is "large" with respect to ε, for example when ε = δt, one has The leading part of G m becomes G M = √ 5 The amplitude of oscillations of G m are damped in exponentially in time even if their frequency is not the good one.In summary, when ε → 0, the discrete solution of (S 1 ) behaves like the continuous one.
For scheme (S2), one has estimates for the linear part that are independent of ε see [5].If we take as in the preceding case ε = δt, one gets 1− X1,2 , 1+ X1,2 as δt → 0 with X1,2 = (1∓i 1 2 ) −1 and so this means that the solution behaves like the continuous case.In other words, (F, H + I) converge to the solution of Zakharov system (1.1) whereas G oscillates at frequency 2 ε .Moreover its amplitude is of order |G(0)| H s which is taken small.We do some numerical experiments in order to illustrate the previous comments.We compute solution to (S 1 ) and (S 2 ) in 1D and we set in the periodic case.We use a finite difference discretization of −∆.We compute on the space-time interval [−5, 5] × [0, 2.10 −3 ] with 1024 points in space.We take f = 4e −20x 2 , N 0 = 0, N 1 = 0 and ε = 10 −5 .For g we choose two different values : • the first one, compared to the compatible case (that is • for the second one, we choose an arbitrary value independent of ε, for example g = e −20x 2 . In the two following figures, we represent the L 2 norm of E δt according to time given by • scheme (S 1 ) for two values of δt that are large with respect to ε (figure 1); • scheme (S 2 ), in the compatible case (left curve) and in the non-compatible case (right curve) for δt = 10 −6 (figure 2).In figure 1, we notice that the oscillations are damped when the time step increases : this is the dissipative phenomenon described above.In figure 2, the only difference between the two curves concerns the amplitude of the oscillations.Indeed, in the non compatible case, the amplitude is of order ε contrary to the compatible case where they are much smaller.

Polarization
In this part, we are going to show that scheme (S 1 ) is able to project the numerical solution on a polarized functions space.More precisely, we will compute the solution of the Cauchy problem (1.2) with "ill-prepared" initial data.Then we will show that these solutions converge, as ε tends to zero, to the solution of the Cauchy problem (1.2) with polarized initial data (which is asymptotically the same than the one of Zakharov system (1.1)).For a precise description of this type of initial data, see below.

The Klein-Gordon-Wave system
We use several notions introduced in [8], [9], [10].Let us first recall Klein-Gordon-wave system (1.4) in 1D : Introducing ϕ = ε∂ t E, ψ = √ ε∂ x E and U = (E, ϕ, ψ) such that (U, n) be a solution of the Cauchy problem for the following system : with Then we expand U (t, x) under the form ))e iθ + c.c. (3.4) with θ = − t ε and where c.c. denotes the complex conjugate.The quantity n is unchanged.Remark 3.1.The justification of these expansion has been established in [4].
In summary, any initial data which read with g bounded uniformly with respect to ε such that lim ε→0 εg = 0 are polarized.The associated initial condition of Zakharov system (1.1) reads • Case of ill-prepared initial data : we choose In opposition with the previous case, we can impose that lim ε→0 εg = 0 and in this case, the initial data is not bounded with respect to ε.It induces that |G(0)| H s (R) = O(1).In summary, any initial data which read ) are ill prepared.We take for example g = i f ε so U 0 (x) = f 2 (e + + e − ) + c.c..The associated initial condition of Zakharov system is It is recovered from U 0 (x) = f 2 (e + + e − ) + c.c. by getting |e − | = 0.
We investigate the convergence of the solution of each scheme when ε → 0 for polarized and ill prepared initial data.

Numerical experiments
We perform the computation on the space-time domain [−5, 5] × [0, 1] with 2048 points in space and for δt = 2, 5.10 −6 .We get some values of ε for which Zakharov model is available, that is for ε small enough.We will use the following errors involving the and e pola L 2 (t m ) := where F m is the mean value in time of E m .Therefore e L 2 (t m ) computes the error between the solution of each scheme (S 1 ) or (S 2 ) and the one of the Zakharov system whereas e pola L 2 (t m ) computes the error between the mean value in time of the solution of each scheme (S 1 ) or (S 2 ) and the one of the Zakharov system.For the four following figures, we represent e L 2 (t) with solid lines and e pola L 2 (t) with dotted lines.
In figures 5, 6, we represent e L 2 (t) (solid line) and e pola L 2 (t) (dotted line) obtained by scheme (S 1 ) (left curves) and scheme (S 2 ) (right curves).We get ε = 10 −4 for figure 5 and ε = 10 −6 for figure 6.We start with two different initial conditions E(t = 0) and E 0 (t = 0).However solid and dotted lines of figures 5 (a) and 6 (a) are superposed.It means that the solution of scheme (S 1 ) converges to the discrete solution of Zakharov system 1.1 whatever initial data we take.This is the polarization phenomenon.On the contrary only the part of F of the solution of scheme (S 2 )  converges to the solution of Zakharov system 1.1.We conclude that only scheme (S 1 ) exhibits the polarization property.

Large spectrum or short waves
The aim of this section is to show that system (1.2) models the propagation of short or ultra short waves [1] through a plasma contrary to Zakharov system (1.1).Moreover we illustrate the error estimate established in [4].In order to solve system (1.4) we use a finite difference method in 1D.The numerical scheme reads endowed with the initial conditions Operators D + and D − have been introduced in section 1.We introduce the errors in L .

1 and Y m+1 2 we have the following result : Proposition 2 . 7 .
Operators Y 1 et Y 2 are unitary in H s−1 (R d ).