A Ruzzon

Y Elskens
email: yves.elskens@univ-provence.fr

F Doveil

A symplectic, symmetric algorithm for spatial evolution of particles in a time-dependent field 6

Keywords: symplectic integrator, wave-particle interaction, hamiltonian chaos, resonance overlap, traveling wave tube PACS: 52.65.Cc, 05.45.Pq, 41.75.-i, 05.45.Ac, 52.65.Yy, 02.70.Ns, 52.20.Dq

A symplectic, symmetric, second-order scheme is constructed for particle evolution in a time-dependent field with a fixed spatial step. The scheme is implemented in one space dimension and tested, showing excellent adequacy to experiment analysis.

Introduction

Particle motion in a space-time dependent field is a classical fundamental problem of dynamics. It is generally well solved in most settings under most kinds of requirements. However, most solutions deal with formulations of the dynamics where the independent variable is time, viz. they aim at computing a function y such that the motion reads x = y(t). Yet in some settings one actually describes motion by the reciprocal function τ = y -1 , so that the motion reads t = τ (x).

One such instance is the propagation of electrons in a traveling wave tube, where it is natural to record particles when they pass at a fixed probe location, instead of getting a snapshot of their locations at a given time. Similar physical contexts are met in other particle beam devices, such as accelerators, klystrons, free electron lasers, electronic tubes for wave amplification,... [START_REF] Kartikeyan | Gyrotrons -High power microwave and millimeter wave technology[END_REF] To some extent, this description is somewhat analogous to eulerian descriptions of flows in hydrodynamics.

If one were interested in the evolution of a single particle, one could merely compute its motion as x = y(t) and deduce its "schedule" function t = τ (x) and related quantities as functions of spatial position. In this respect, many symplectic methods are available (see e.g. Refs [START_REF] Hairer | Geometric numerical integration : Structurepreserving algorithms for ordinary differential equations[END_REF][START_REF] Leimkuhler | Simulating hamiltonian dynamics[END_REF] for an overview), especially for separable hamiltonians of the form H(p, x, t) = K(p) + V (x, t). However, to describe a beam of many particles during their spatial progression, it is reasonable to follow them consistently in space, to generate numerical data sampled at the same (possibly many) space positions. It then becomes awkward to first evolve them in time and afterwards reconstruct their progression in spatial terms by interpolations.

For this purpose we reformulate in Section 2 the particle equations of motion, using the streaming variable x as independent variable (see figure 1). Since the original particle dynamics is hamiltonian, we ensure that the new description be symplectic by first expressing the action principle in terms of the timetable function τ . In the corresponding hamiltonian picture, the variable conjugate to τ is the energy ζ, and the generator of motion is momentum P.

In Section 3 we stress our requirements on the scheme and consider alternative strategies. Then we construct a first order symplectic scheme for the particle motion. The implicit part of the step can be performed either through algebraic solution of a cubic equation, or through a Newton iteration : we compare both procedures. Next we construct the adjoint, first order symplectic scheme, which also requires a Newton iteration, and we check its accuracy. Finally, we combine the direct and adjoint schemes to obtain a second order symmetric, symplectic, fixed ∆x scheme.

In Section 4 we benchmark our algorithm by analysing the particle motion in the field of a single harmonic wave, viz. we solve the pendulum motion in a galilean frame. Numerical simulations for realistic beam data generate beam deformations shown in Section 5.

Section 6 focuses on the evolution of a beam launched in presence of two harmonic waves. Simulations are confronted with experimental observations of the beam collected at the device outlet. Special attention is paid to the reproduction of a devil staircase structure, characteristic of the chaotic behaviour of the system, taking into account the finiteness of the experimental device.

In summary, experimental data often relate to limited interaction times, while numerical evidences and theoretical discussions of chaotic dynamics often deal with trajectories followed for long times in a compact domain of phase space. The agreement of our simulations with experimental evidence assesses the relevance of our algorithm to such experimental settings.

Evolution with respect to space

Rewriting the equations of motion in hamiltonian form with respect to space is straightforward in the symplectic formalism (see Section 2.2). However one may wish first a more pedestrian derivation, from the classical action principle.

Lagrangian viewpoint

The action for a non-relativistic particle with mass m moving along a one-dimensional axis Ox in a time dependent potential V (x, t) reads

S[y; t 0 , t 1 , x 0 , x 1 ] = t 1 t 0 L(y(t), ẏ(t), t)dt ( 1
)
where y is a continuously differentiable function of time t ∈ [t 0 , t 1 ] ⊂ R, subject to the constraints y(t 0 ) = x 0 and y(t 1 ) = x 1 , and the dot denotes derivative with respect to t. The lagrangian is

L(x, v, t) = mv 2 2 -V (x, t) . (2)
In the following we restrict all trajectories to the class of strictly monotone, increasing functions, viz. inf t 0 <t<t 1 ẏ(t) > 0. For these functions the reciprocal function τ : [x 0 , x 1 ] → [t 0 , t 1 ] : x → τ (x) exists ; τ is unique and also strictly monotone, increasing, continuously differentiable with

τ ′ (x) = dτ dx (x) = dy dt (τ (x)) -1 = 1 ẏ(τ (x)) (3)
where the prime denotes derivative with respect to x.

To rewrite (1) as a space-integral, we introduce the new lagrangian

L(t, u, x) := L(x, 1 u , t)u (4)
so that

S[τ ; x 0 , x 1 , t 0 , t 1 ] := S[τ -1 ; t 0 , t 1 , x 0 , x 1 ] = x 1 x 0 L(τ (x), τ ′ (x), x)dx . (5)
It is convenient to introduce the opposite to the (usual definition of) canonical momentum conjugate to τ ,

ζ = - ∂L ∂u (6)
and to perform the Legendre transform of -L, defining

P(ζ, t, x) := ζτ ′ + L(t, τ ′ , x) (7)
so that in the new variables the canonical Hamilton equations read

dτ dx = ∂P ∂ζ , (8)
dζ dx = - ∂P ∂τ . (9)
For the classical lagrangian (2) the new variables are the usual energy and linear momentum,

ζ = m 2τ ′ 2 + V (x, τ (x)) , (10)
P = 2m(ζ -V (x, τ (x))) . (11)

Hamiltonian viewpoint

The hamiltonian formulation of dynamics provides a direct path to the latter equations. Indeed it suffices to consider the symplectic 2-form dω := dpdx -dHdt [START_REF] Hairer | Geometric numerical integration : Structurepreserving algorithms for ordinary differential equations[END_REF] where p = m ẏ is conjugate to x and H = ζ is conjugate to t = τ , and to consider x as the independent variable along trajectories instead of t. The minus sign we introduced in ( 6)-( 7) ensures the usual signs in dω.

Of course, if the potential does not depend explicitly on time, ζ is a first integral. The above requirement, that ẏ nowhere vanishes, will be strengthened below by imposing P ≥ p min > 0 for some p min to ensure appropriate numerical accuracy.

Discrete time canonical transformations

To preserve the symplectic structure while integrating ( 10)-( 11) we use a sequence of canonical transformations. Note that [START_REF] Escande | Renormalization method for the onset of stochasticity in a hamiltonian system[END_REF] does not allow using a splitting method such as leap-frog [START_REF] Hairer | Geometric numerical integration : Structurepreserving algorithms for ordinary differential equations[END_REF][START_REF] Leimkuhler | Simulating hamiltonian dynamics[END_REF], as it is not the sum of integrable generators. In order to sample the spatial evolution regularly (as desired e.g. to follow many particles in parallel), we use a fixed spatial step ∆x, a strategy recommended e.g. by Hénon [START_REF] Hénon | On the numerical computation of Poincaré maps[END_REF] to obtain simple and accurate Poincaré maps (for a space periodic potential, the step ∆x is best taken as a fraction of the wavelength). A fixed spatial step also avoids generating noisy-like "spurious frequencies" in the simulated dynamics. Note that we could also use a time-integrator (leap-frog or other) with an adapting time step, ∆t = m∆x/P [START_REF] Hénon | On the numerical computation of Poincaré maps[END_REF], but it is not manifestly symplectic as P is not constant. Another approach, using a fixed time step and interpolating the trajectories for the spatial mesh, would raise the issue of constructing a symplectic interpolation scheme (of desired order). In this work we settle for the manifestly symplectic, fixed step approach, which can process many particles in parallel.

Our first integrator F ∆x : (τ, ζ, x) → (τ , ζ, x + ∆x) is chosen to provide an explicit, first-order approximation for the time increment. It is generated by

F (τ, ζ) = τ ζ + P(τ, ζ, x)∆x (13)
so that the system

ζ = ζ + ∂ τ P(τ, ζ, x)∆x , (14)
τ = τ + ∂ ζ P(τ, ζ, x)∆x , (15)
is a first-order, symplectic approximation to ( 8)-( 9). As the second equation is explicit with respect to time τ , the first equation is implicit with respect to the energy ζ. For the special case of momentum [START_REF] Escande | Renormalization method for the onset of stochasticity in a hamiltonian system[END_REF], it actually leads to a cubic equation,

ζ3 -(V + 2ζ) ζ2 + (2V ζ + ζ 2 ) ζ -V ζ 2 - m 2 (∂ τ V ) 2 ∆x 2 = 0 , (16)
which can be solved algebraically for ζ. Here V and ∂ τ V are computed at (x, τ ). Equation ( 16) usually has three real solutions, two of which being ∆x-close to ζ : the relevant root is such that

(ζ -ζ) ∂ τ V > 0.
It is advantageous to express [START_REF] Kartikeyan | Gyrotrons -High power microwave and millimeter wave technology[END_REF]

in the form ζ = ζ + (ζ -V )σ, so that σ = -a (1 + σ) -1/2 (17)
with the single parameter

a = m/2 (ζ -V ) -3/2 ∂ τ V ∆x ( 18
)
calculated at (x, τ ). The fixed point equation ( 17) is easily solved by the Newton method, which selects the "good" cubic root automatically (as aσ < 0) and converges very fast, especially if a is small. Analytically, this method stresses the small parameter a controlling the accuracy of our scheme : it involves a balance of the potential evolution ∂ τ V and the spatial step ∆x against (ζ -V ), viz. against the particle velocity P/m. For small velocity, the algorithm deteriorates -at worst it will miss turning points where P changes sign (which is forbidden by our assumptions on trajectories in the action principle). Let Z C ∆x : (τ, ζ, x) → (τ, ζ, x) and Z N ∆x : (τ, ζ, x) → (τ, ζ, x) denote respectively the cubic and Newton solvers. For perfectly accurate computations, they coincide and may be denoted identically Z ∆x . Note that Z ∆x is not symplectic, as

det DZ ∆x (τ, ζ, x) = [1 + ∂ ζ ∂ τ P(τ, ζ, x)∆x] -1 . ( 19
)
With the new energy ζ, [START_REF] Leimkuhler | Simulating hamiltonian dynamics[END_REF] immediately provides the new time, defining the map T ∆x : (τ, ζ, x) → (τ , ζ, x). This map is not symplectic either, as

det DT ∆x (τ, ζ, x) = 1 + ∂ ζ ∂ τ P(τ, ζ, x)∆x . (20)
Finally, we advance position, with I ∆x : (τ, ζ, x) → (τ, ζ, x + ∆x). The resulting integration scheme F ∆x = I ∆x • T ∆x • Z ∆x is symplectic by construction, within machine accuracy, as the planar map

T ∆x • Z ∆x is area-preserving : det D(T ∆x • Z ∆x ) = (det DT ∆x )(det DZ ∆x ) = 1.
A = 0.1, k = 0.2, v φ = 1, φ = π/4.
The variables advanced with ( 14) and ( 15) are in the form

τ ∼ = τ + τ ′ ∆x + τ ′′ 2 ∆x 2 + O(∆x 3 ) , (21)
ζ ∼ = ζ + ζ ′ ∆x + ζ ′′ 2 ∆x 2 + O(∆x 3 ) , (22)
where τ ′ and ζ ′ have been approximated with -∂ τ P(τ, ζ, x) and ∂ ζ P(τ, ζ, x). It follows that the most influential theoretical error, in every step of integration, is given by τ ′′ ∆x/2 for τ and ζ ′′ ∆x/2 for ζ. Table 1 compares the upper estimated theoretical errors, related to τ ′′ and ζ ′′ , and the maximum real simulation errors for five particle initial velocities. Slower particles are found affected by larger errors as expected.

To obtain a second order, symmetric scheme, we consider the adjoint map [START_REF] Hairer | Geometric numerical integration : Structurepreserving algorithms for ordinary differential equations[END_REF], which is also symplectic, generated by the function

F * (τ , ζ) = τ ζ -P(τ , ζ, x + ∆x)∆x (23) so that τ = τ -∂ ζ P(τ , ζ, x + ∆x)∆x , (24)
ζ = ζ -∂ τ P(τ , ζ, x + ∆x)∆x . (25)
For momentum [START_REF] Escande | Renormalization method for the onset of stochasticity in a hamiltonian system[END_REF], both equations involve the potential V (x + ∆x, τ ), which implies that one first solves (24) with respect to the new time τ by a Newton algorithm, and then computes the new energy ζ by (25). This defines the symplectic map

F * ∆x = Z * ∆x • T * ∆x • I ∆x : (τ, ζ, x) → (τ , ζ, x + ∆x).
One easily checks that I ∆x is self-adjoint, while Z * ∆x • Z -∆x and T * ∆x • T -∆x reduce to identity up to machine numerical tolerance (typically 10 -15 ).

Finally the composition

F (2) ∆x = F * ∆x/2 • F ∆x/2 = Z * ∆x/2 • T * ∆x/2 • I ∆x • T ∆x/2 • Z ∆x/2 (26)
is its own adjoint. It is thus symmetric and therefore second order [START_REF] Hairer | Geometric numerical integration : Structurepreserving algorithms for ordinary differential equations[END_REF].

Validation : particle dynamics in a single wave

We test our schemes with the time dependent potential of a wave

V (x, t) = A cos(kx -kvt + φ) (27)
where A, k, v, φ are respectively the amplitude, wavevector, phase velocity and phase of the wave. Rescaling energy (and amplitude), space and time enables one to set m, k and v to unity, and the choice of the origin of time or space eliminates φ. Its integrability makes this dynamics a good benchmark. The accuracy of the determination of the adjoint map is checked by iterating first F ∆x for ∆x = 0.01 from x = 0 to x = L = 300, and then F * -∆x from x = L to x = 0, for five particles. Figures 2 display the discrepancies ∆ζ and ∆τ as functions of x for each particle and confirm that F * -∆x = F -1 ∆x to numerical accuracy. The order of the algorithms and their accuracy is further analysed in figure 3, comparing the first order and second order schemes for the motion of a particle with initial velocity v in = 1.5 in the field of a wave with A = 0.1, φ = π, k = 0.2, v φ = 1 over a length L = 100.

Beam dynamics in a single wave

A second accuracy check is provided by the Poincaré section of the beam by the positions x mod L = 0. The return map for (τ, ζ) variables is symplectic. As the particle motion is integrable (it reduces to the pendulum by a Galileo transformation to the reference frame comoving with the wave), each orbit must generate section points on lines satisfying the algebraic relation

ζ = mv 2 φ + H 2 ± v φ 2m H -2mA cos[k(x -v φ τ ) + φ] ( 28
)
where H is a constant. In particular, the motion on the wave separatrix corresponds to H = A. For H > A, this relation defines two branches for all times τ , which correspond to faster or slower circulating particles, while for H < A the relation defines the upper and lower part of trapped motion inside the wave's cat eye. Figure 4 shows that numerical trajectories perfectly reproduce these lines.

To assess the relevance of the algorithm to experiment we also follow the deformation of a beam of electrons injected in a single wave, e.g. in a traveling wave tube. As particles are accelerated or decelerated by the wave, the beam velocity profile is deformed. We thus inject N particles at x = 0, equally distributed over one time period of the wave, and plot the histogram of particle velocities as a function of abscissa x.

In figure 5 a cold beam is made of particles injected resonantly with the wave velocity, v b = v φ = 1, with m = k = 1. The wave amplitude, A = 0.002, determines the bounce frequency ω b = √ kA so that particle oscillations in the wave trough have a spatial period L b = 2πv φ /ω b = 140.5. The particles bounce indeed and, in agreement with the rotating bar approximation [START_REF] Mynick | Soluble theory of nonlinear beam-plasma interaction[END_REF], most of them reconvene every L b /2. Only the ones injected at times close to (φ + 2πn)/(kv φ ) (for integer n) enter the wave close to the X point and follow closely the inner side of the cat eye separatrix. The time these particles need to overcome half a wavelength can be arbitrarily long, so that they mark the boundary of the wave resonant domain at velocities v φ ± 2 A/m.

For particles injected with a velocity outside the wave cat's eye, the beam is modulated. A significant difference between (x, v) plots for propagating beams and the more familiar (x, v) Poincaré sections of particle-in-wave dynamics (see e.g. [START_REF] Escande | Renormalization method for the onset of stochasticity in a hamiltonian system[END_REF][START_REF] Escande | Stochasticity in classical hamiltonian systems : universal aspects[END_REF][START_REF] Elskens | Microscopic dynamics of plasmas and chaos[END_REF]) is the asymmetry between faster and slower particles, obvious in Figure 6.

In particular, for particles injected with a velocity above the wave cat's eye, the beam is moderately asymmetric. But for a particle injection velocity within the capture range [v φ -2 A/m, v φ + 2 A/m], the picture gets strongly deformed, as part of the beam is trapped as in figure 5 while part of it moves outside the wave cat's eye.

Particle dynamics in two waves : resonance overlap and chaos

The motion of a particle in the field of two waves is a paradigm of hamiltonian chaos. In our formulation, Poincaré sections are given by x mod L = 0, and the return map is symplectic, hence area-preserving in conjugate variables (τ, ζ). Figure 7 displays this Poincaré section, showing the growth of the chaotic domain for increasing wave amplitudes, and the destruction of KAM tori [START_REF] Escande | Renormalization method for the onset of stochasticity in a hamiltonian system[END_REF].

The corresponding transition to large scale chaos by increasing the resonance overlap parameter s = 2(

√ A 1 + √ A 2 )/(|v φ2 -v φ1 | √ m)
is also observed by recording the particle velocities at a fixed traveled distance L 0 , after being injected at a fixed velocity v in . As seen in figures 6, a cold beam injected in a wave cat's eye spreads over the velocity interval spanned by this cat eye, and if the beam is injected outside cat eyes it remains confined between the velocities of KAM tori on either side. Beam velocity spreading (also called heating) has been used to diagnose resonance overlap, and our numerical scheme is compared with experimental data [START_REF] Doveil | Experimental observation of resonance overlap responsible for Hamiltonian chaos[END_REF] in figure 8.

Moreover, the transition to large scale chaos in phase space is known to occur stepwise. For increasing wave amplitudes, successive KAM tori get destroyed, so that the beam invades velocity domains resulting from the merging of capture regions corresponding to "secondary" resonances [START_REF] Escande | Renormalization method for the onset of stochasticity in a hamiltonian system[END_REF]. The accessible velocity interval for the beam injected in one wave then grows like a devil's staircase, the higher steps corresponding to the merging with major secondary resonances. Figure 9 compares these domains obtained both numerically and experimentally [START_REF] Macor | Electron climbing a "devil's staircase" in wave-particle interaction[END_REF][START_REF] Doveil | Direct observation of a "devil's staircase" in wave-particle interaction[END_REF][START_REF] Buchanan | Richness in simplicity[END_REF]. While experimental data are blurred due to recording accuracy, numerical data have limited resolution due to the large number of particles (here only 25000) needed for the sharp observation of a threshold. Nevertheless, the agreement is quantitatively satisfactory.

Conclusion

The scheme has proved its relevance to describe accurately particle motion in a given field, along a single space dimension. It also provides new pictures of known behaviours, which complement more familiar, usually more symmetric plots in (x, v) space.

Extending our approach to three space dimensions is rather straightforward provided the particles stream along a given coordinate, say x. Higher-order symplectic schemes can also be constructed by composing several maps F γ i ∆x and F * γ i ∆x with appropriate substeps γ i [START_REF] Hairer | Geometric numerical integration : Structurepreserving algorithms for ordinary differential equations[END_REF]. More challenging is the issue raised by the particle feedback on the wave field, calling for a self-consistent space-based model of particles and waves evolution, in the spirit of models used for weak plasma turbulence [START_REF] Elskens | Microscopic dynamics of plasmas and chaos[END_REF].

It will also be interesting to apply this scheme to model the many-waves regime of weak plasma turbulence, where particle velocity undergoes a chaotic transport over a wide range [START_REF] Cary | Non quasilinear diffusion far from chaotic threshold[END_REF][START_REF] Bénisti | Origin of diffusion in Hamiltonian dynamics[END_REF][START_REF] Elskens | Microscopic dynamics of plasmas and chaos[END_REF][START_REF] Elskens | Nonquasilinear evolution of particle velocity in incoherent waves with random amplitudes[END_REF][START_REF] Besse | Validity of quasilinear theory : refutations and new numerical confirmation[END_REF], as properties of this transport are still controversial.

Figure 1 :

 1 Figure1: Sketch of the algorithm principle in (x, t) plot. The continued line is a trajectory (identical for both schemes, with small steps). Crosses sample positions (every 50∆t) as obtained by leap-frog with constant time step, circles sample times (every 50∆x) as obtained by our approach with constant space step. Straight lines are guides to the eye.

Figure 2 :

 2 Figure 2: Discrepancy, between direct evolution F ∆x and backward evolution F * -∆x , for energy (left panel) and arrival time (centre) as functions of position, for five particles in the field of a single wave with A = 0.1, k = 0.2, v φ = 1, φ = π. Particles injected at the origin at t = 0, with velocities 1.1, 1.5, 2, 2.5 and 3.5, move as shown on right panel.

Figure 3 :

 3 Figure3: Numerical accuracy for the first order algorithm (red, upper lines) and second order (blue, lower lines) versus calculation time per spatial unit length. Darker colour for the Newton method solution to the cubic equation, lighter colour for polynomial solver (slightly slower than Newton).

Figure 4 :

 4 Figure 4: (left) Poincaré section in (τ, ζ) variables, of the return map after one wavelength L. Trajectories for five trapped particles injected with v in = 1.8, 1.85, 2, 2.1, 2.25, for five untrapped particles with v in = 1.4, 1.5, 1.6, 1.7, 2.3, and for one particle injected on each separatrix branch. Wave parameters A = 0.02, v φ = 2, k = 0.2, φ = π ; particle mass m = 1. (right) Exact section lines (28).

Figure 5 :

 5 Figure 5: Velocity distribution function of particles along the x axis, injected at the wave phase velocity, v b = v φ = 1, with m = k = 1 and wave amplitude A = 0.002.

Figure 6 :Figure 7 :

 67 Figure 6: Velocity distribution function of particles along the x axis, when injected at v b = 1.6 (upper left), v b = 1.5 (upper right), v b = 1.4 (lower left) and v b = 1.3 (lower right), above the wave phase velocity, v φ = 1, with m = k = 1 and wave amplitude A = 0.04. The cat's eye boundaries lie at v φ + 2 √ A = 1.4 and v φ -2 √ A = 0.6.

Figure 8 :

 8 Figure 8: Velocities of particles after interaction with two waves, versus injection velocity. (top) Experimental data with dotted lines marking cat eyes boundaries ; (bottom) numerical results. (left) Non-overlapping cat eyes, s = 0.63 ; (right) overlapping trapping domains, s = 1.5.

Figure 9 :

 9 Figure 9: Velocities of particles after interaction with two waves, versus amplitude of waves A 1 = A 2 . Cold beam injected at phase velocity of one wave. (left) Numerical results for v b = v φ1 = 1.1487, v φ2 = 0.8609, k 1 = 0.6529, k 2 = 1.7424. (right) Experimental data.

Table 1 :

 1 It is first order only. v in δζ estimate δζ simulation δτ estimate δτ simulation The largest single step error for five particles with different initial velocity launched in the field of a single wave with

	2.5 0.16 1.5 0.11 0.9 0.10	6 • 10 -4 9 • 10 -4 0.008	0.02 0.10 8	0.005 0.03 0.5
	0.6 20	0.1	12	5
	0.4 3	0.3	70	16

c 2011 The authors. Reproduction of this article, in its entirety, for noncommercial purposes is permitted.

Acknowledgements

The work of A. Ruzzon was made possible by an Erasmus mobility from Università degli studi di Padova to université de Provence. This work benefited from fruitful discussions with D.F. Escande.

Appendix A. Explicit derivation of [START_REF] Mynick | Soluble theory of nonlinear beam-plasma interaction[END_REF] To reduce [START_REF] Macor | Electron climbing a "devil's staircase" in wave-particle interaction[END_REF] to [START_REF] Mynick | Soluble theory of nonlinear beam-plasma interaction[END_REF], note that the derivatives of momentum [START_REF] Escande | Renormalization method for the onset of stochasticity in a hamiltonian system[END_REF] read

and consider the dimensionless (both physically and numerically relevant) quantity σ = ( ζζ)/(ζ -V ), which characterizes the changes in particle energy per step and must be small for an accurate calculation. Equation ( 14) reads, with arguments (τ, ζ, x) in P and (τ,

using the dimensionless a defined by [START_REF] Ruzzon | Studio numerico della propagazione spaziale d'un fascio di test in un tubo a onde progressive e confronto con l'esperimento[END_REF]. The fixed point equation ( 17) then follows.

Appendix B. Explicit flowchart of the second order algorithm

The symmetric scheme (26) may be expanded as follows:


Boris Andreianov

Time Compactness Tools for Discretized Evolution Equations and Applications to Degenerate Parabolic PDEs

Keywords: time translates, Kruzhkov lemma, order-preservation, finite volumes

We discuss several techniques for proving compactness of sequences of approximate solutions to discretized evolution PDEs. While the well-known Aubin-Simon kind functional-analytic techniques were recently generalized to the discrete setting by Gallouët and Latché [15], here we discuss direct techniques for estimating the time translates of approximate solutions in the space L 1 . One important result is the Kruzhkov time compactness lemma. Further, we describe a specific technique that relies upon the order-preservation property. Motivation comes from studying convergence of finite volume discretizations for various classes of nonlinear degenerate parabolic equations. These and other applications are briefly described.

Introduction

Let us think of evolution equations set on a cylindrical domain Q := (0, T ) × Ω ⊂ IR + × IR N . Proving convergence of space-time discretizations of such equations often includes the three following steps : constructing discrete solutions and getting uniform (in appropriate discrete spaces) estimates; extracting a convergent subsequence; writing down a discrete weak formulation (e.g., with discretized test functions) and passing to the limit in the equation in order to infer convergence.

For the first step, obtention of estimates is greatly simplified by preservation, at the discrete level, of the key structure properties of the PDE (such as symmetry, coercivity, monotonicity of the diffusion operators involved; entropy dissipation, for the nonlinear convection operators in the degenerate parabolic case; etc.). For getting discrete a priori estimates test functions are often used, as in the continuous case. Therefore, some analogues of integration-by-parts formulas and chain rules are instrumental for the first step. For the examples we give in this paper, "discrete Boris Andreianov CNRS UMR 6623, 16 route de Gray, Besanc ¸on, France, e-mail: boris.andreianov@univ-fcomte.fr 1 duality" type schemes (mimetic, co-volume, DDFV; see, e.g., [START_REF] Andreianov | On 3D DDFV discretization of gradient and divergence operators. II. Discrete functional analysis tools and applications to degenerate parabolic problems[END_REF] and references therein) can be used to guarantee an exact integration-by-parts feature. In contrast, chain rules for derivation in time or in space must be replaced by approximate analogues, often taking the form of convexity inequalities (see, e.g., [START_REF] Andreianov | Discrete duality finite volume schemes for doubly nonlinear degenerate hyperbolic-parabolic equations[END_REF], [START_REF] Andreianov | On 3D DDFV discretization of gradient and divergence operators. II. Discrete functional analysis tools and applications to degenerate parabolic problems[END_REF]Sect. 4]).

In this note, we give some insight into convergence proofs for different subclasses of degenerate elliptic-parabolic-hyperbolic PDEs under the general form1

u = b(v), w = ϕ(v), u t -div G(v) -a 0 ∇ ∇ ∇w + ψ(v) = f in Q = (0, T )×Ω (1)
with b(•), ϕ(•), ψ(•) continuous2 non-decreasing on IR, normalized by zero at zero, with a continuous convection flux G(•) and with a 0 : IR N → IR N of Leray-Lions type (see e.g. [START_REF] Alt | Quasilinear elliptic-parabolic differential equations[END_REF][START_REF] Andreianov | Discrete duality finite volume schemes for doubly nonlinear degenerate hyperbolic-parabolic equations[END_REF]; p-laplacian, with a 0 (ξ ) = |ξ | p-2 ξ is a typical example). For the sake of simplicity, homogeneous Dirichlet boundary condition on (0, T ) × ∂ Ω is taken.

But our main goal is to discuss the second step of the proofs3 , the one of getting compact4 sequences of discrete solutions. For linear problems, the two latter steps are somewhat trivial; indeed, mere functional-analytic bounds would lead to compactness in a weak topology, which is enough to pass to the limit from the discrete to the continuous weak formulation of the PDE. Thinking of nonlinear problems and passage to the limit in nonlinear terms, bounds in functional spaces can be sufficient when combined with basic compact embeddings; but this requires rather strong bounds involving e.g. some estimates of the derivatives. Regarding evolution PDEs of, say, porous medium type, L p bounds are available on the space derivatives but not on the time derivatives (those belong to some negative Sobolev spaces). In this situation, either compactness in an ad hoc strong topology is needed; or the weak compactness coming from uniform boundedness should be combined with some compactification arguments (compensated compactness, Young measures and their reduction, etc.) that exploit in a non-trivial way the particular structure of the PDE in hand (div-curl structure, pseudo-monotonicity, entropy inequalities, etc.).

In this note, we first recall in § 2 the fundamental techniques using only bounds in well-chosen functional spaces (see [START_REF] Aubin | Un théorème de compacité[END_REF][START_REF] Dubinskii | Weak convergence for elliptic and parabolic equations[END_REF][START_REF] Simon | Compact sets in the space L p (0, T ; B)[END_REF][START_REF] Amann | Compact embeddings of vector-valued Sobolev and Besov spaces[END_REF] for the continuous setting; see [START_REF] Gallouët | Compactness of discrete approximate solutions to parabolic PDEs -Application to a turbulence model[END_REF][START_REF] Emmrich | Doubly nonlinear evolution equations of second order: Existence and fully discrete approximation[END_REF] for the corresponding discrete results). In § 3, we present a collection of complementary techniques for estimation of time translates of families of functions that already possess some estimate of space translates. In § 4, we describe one indirect method for proving compactness and convergence of families of approximate solutions. The method heavily exploits the order-preservation property, required both for the PDE in hand and for the approximation scheme in use. Throughout the note, the exposition is motivated and illustrated by applications to approximation of several cases of problem (1) (different cases requiring different approaches).

2 Functional-analytic approach of Aubin-Lions-Dubinskii-Simon

In the continuous setting, one celebrated result is the Aubin-Lions or Dubinskii lemma ( [START_REF] Aubin | Un théorème de compacité[END_REF] and [START_REF] Dubinskii | Weak convergence for elliptic and parabolic equations[END_REF]) and its generalization by Simon [START_REF] Simon | Compact sets in the space L p (0, T ; B)[END_REF] (see also Amann [START_REF] Amann | Compact embeddings of vector-valued Sobolev and Besov spaces[END_REF]). To give an example relevant for the applications we have in mind, let us simply state here that a sequence (u h ) h bounded in L 1 (0, T ;W 1,1 (Ω )) and such that ((u h ) t ) h is bounded in L 1 (0, T ;W -1,1 (Ω )) is relatively compact in L 1 (Q), cf. [START_REF] Gallouët | Compactness of discrete approximate solutions to parabolic PDEs -Application to a turbulence model[END_REF]. More generally, compactness comes from an a priori bound on u h in some space L p (0, T, X) with X compactly embedded in L 1 (Ω ) (e.g., X = W 1,1 (Ω )) while the PDE brings information on boundedness of the time derivatives (u h ) t in some L q (0, T ;Y ) where Y can be a subspace of distributions on Ω equipped with a rather weak topology (e.g., Y = W -1,1 (Ω )). A discrete version of the Aubin-Simon lemma was recently proposed by Gallouët and Latché in [START_REF] Gallouët | Compactness of discrete approximate solutions to parabolic PDEs -Application to a turbulence model[END_REF]; it is based upon a careful reformulation of estimates in terms of "coherent" families (X h ) h , (Y h ) h of discrete spaces.

A related result taken from Simon [START_REF] Simon | Compact sets in the space L p (0, T ; B)[END_REF] and Amann [START_REF] Amann | Compact embeddings of vector-valued Sobolev and Besov spaces[END_REF] uses a bound on fractional time derivatives of u h . As it was demonstrated by Emmrich and Thalhammer in [START_REF] Emmrich | Doubly nonlinear evolution equations of second order: Existence and fully discrete approximation[END_REF], this version is quite appropriate in the time-discretized setting. Indeed, time fractional derivatives of order less than 1/2 exist even for piecewise constant functions. Technically, this method involves an indirect estimation of weighted time translates, under a form

T 0 T 0 |u h (t) -u h (s)| p
|t -s| 1+σ p ds dt with some p ≥ 1 and σ ∈ (0, 1/2).

These results only use bounds in functional spaces and very few of the underlying PDE properties. They offer a very wide spectrum of applications; yet they are difficult to apply on degenerate parabolic problems with non-Lipschitz nonlinearities. The difficulty comes from the fact that non-Lipschitz mappings make bad correspondence between linear functional spaces. Yet this difficulty is not a fundamental one; roughly speaking, it is settled by a careful use of translation arguments and of moduli of continuity. This is the object of the next section.

Direct estimation of time translates

In this section, the compactness question is studied using the one and only space5 L 1 (Q). By the Fréchet-Kolmogorov compactness criterion in L 1 (Q), uniform bounds on space and time translates of u h are needed. In the setting of the present note, the first ones are readily available. The difficulty lies in estimating the time translates as

∀h T -δ 0 Ω u h (t+δ ) -u h (t) ≤ ω(δ ) with lim δ →0 ω(δ ) = 0, (2)
ω(•) being a modulus of continuity, uniform in h. Here are two ways to obtain (2).

A discrete Kruzhkov lemma 6

Lemma 1 (Kruzhkov [16]). Assume that the families of functions

(u h ) h , (F α h ) h,α are bounded in L 1 (Q) and satisfy ∂ ∂t u h = ∑ |α|≤m D α F α h in D (Q).
Assume that u h can be extended outside Q, and one has 7

Q |u h (t, x+δ )-u h (t, x)| dxdt ≤ ω(δ ), with lim δ →0 ω(δ ) = 0, (3)
where ω(•) does not depend on h. Then (u h ) h is (relatively

) compact in L 1 (Q).
Clearly, this is an L 1 loc compactness result (one can apply the lemma locally in Q). For problem [START_REF] Alt | Quasilinear elliptic-parabolic differential equations[END_REF], the value m = 1 is relevant, because an L 1 bound is available for the flux G(v) -a 0 ∇ ∇ ∇ϕ(v) ; therefore we limit to this case the discussion of discrete analogues of Lemma 1. To give an idea of discrete versions of the Kruzhkov lemma 8 , assume we are given a family of meshes of Ω indexed by their size h and satisfying mild proportionality restrictions (e.g., for the case of two-point flux finite volume schemes as described in [START_REF] Eymard | Finite Volume Methods[END_REF], one needs for all neighbour volumes K,L, diam(K)+diam(L)≤const d K,L uniformly in h). Assume that on these meshes, spaces of discrete functions IR h and discrete fields (IR N ) h are defined (each element u h ∈ IR h or F h ∈ (IR N ) h is a piecewise constant on Ω function reconstructed from the degrees of freedom of the discretization method). Assume we are given discrete gradient and discrete divergence operators ∇ ∇ ∇ h and div h mapping between these spaces. Thus all discrete objects (functions, fields, gradient, divergence) are naturally lifted to L 1 (Q).

Let (δ h ) h be the associated time steps, let N h be the entire part of T /δ h . Assume that we are given an initial condition b 0 h and discrete evolution equations under the form

for n ∈ [1, N h + 1], b(v n h ) -b(v n-1 h ) δ h = div h [F n h ] + f n h in IR h , (4)
where families ( u n h n h , f n h n h (discrete functions) and F n h n h (discrete fields) are bounded in L 1 (Q). Assume that the discrete gradients ∇ ∇ ∇ h v n h n h are bounded in L 1 (Q) and that this bound implies a uniform translation bound in space of the family v h (this is true, e.g., when discrete Poincaré inequalities can be proved). Under these assumptions, reproducing at the discrete level the proof [START_REF] Kruzhkov | Results on the nature of the continuity of solutions of parabolic equations and some of their applications[END_REF] of Lemma 1 as it is done in [START_REF] Andreianov | Analysis of a finite volume method to solve a cross-diffusion population system[END_REF][START_REF] Andreianov | On 3D DDFV discretization of gradient and divergence operators. II. Discrete functional analysis tools and applications to degenerate parabolic problems[END_REF], one concludes that the family (b(v h )) h is relatively compact in L 1 (Q). Note that, the case m ≥ 2 would require more work. 6 There is a strong relation to the method of § 2. The Kruzhkov lemma allows for general moduli of continuity. E.g., for problem (1) with ϕ = Id, the Aubin-Lions-Dubinskii-Simon argument can be used if b(•) is Lipschitz continuous (with X = W 1,p 0 (Ω )) or Hölder continuous (with a fractional Sobolev space chosen for X), and the Kruzhkov lemma can be used for any continuous b(•). 7 In practice, space translation estimates of the kind (3) can be obtained via an estimate of some discrete gradients; notice that estimates of kind (3) are stable upon composing (u h ) h by a function b(•) which is uniformly continuous (as in (1), we mean that u h = b(v h )). 8 Here we give a rather heuristic presentation; see [START_REF] Andreianov | Analysis of a finite volume method to solve a cross-diffusion population system[END_REF] and [START_REF] Andreianov | On 3D DDFV discretization of gradient and divergence operators. II. Discrete functional analysis tools and applications to degenerate parabolic problems[END_REF] for two precise formulations covering, e.g., the two-point flux finite volume schemes ( [START_REF] Eymard | Finite Volume Methods[END_REF]) and DDFV schemes ( [START_REF] Andreianov | On 3D DDFV discretization of gradient and divergence operators. II. Discrete functional analysis tools and applications to degenerate parabolic problems[END_REF]

and ref. therein).

A classical technique for the "variational" setting

Following [START_REF] Alt | Quasilinear elliptic-parabolic differential equations[END_REF], by "variational" we mean a setting where the solution w is an admissible test function in the weak formulation of the PDE; e.g., (1) can be tested with w = ϕ(v). It typically comes along with a priori estimates that can be reproduced at the discrete level, provided the discretization is somewhat structure-preserving. 9The technique of [START_REF] Alt | Quasilinear elliptic-parabolic differential equations[END_REF] used, in its finite volume version, e.g., in [START_REF] Eymard | Finite Volume Methods[END_REF][START_REF] Eymard | Convergence of a finite volume scheme for nonlinear degenerate parabolic equations[END_REF][START_REF] Andreianov | Discrete duality finite volume schemes for doubly nonlinear degenerate hyperbolic-parabolic equations[END_REF], is to integrate 10 the equation in time from t to t +δ , take w h (t +δ )-w h (t) for test function, then integrate in (t, x). On problem [START_REF] Alt | Quasilinear elliptic-parabolic differential equations[END_REF], this leads to a uniform estimate

∀h > 0 T -δ 0 Ω b(v h )(t+δ ) -b(v h )(t) ϕ(v h )(t+δ ) -ϕ(v h )(t) ≤ ω(δ ). (5)
Then Lipschitz continuity of ϕ • b -1 (resp., of b • ϕ -1 ) can be used to infer uniform L 2 time translates of

w h = ϕ(v h ) (resp., of u h = b(v h )). Yet the L 1 time translates can be obtained in the case ϕ • b -1 (resp., b • ϕ -1
) is a merely continuous function.

• A technique for L 1 estimates involving non-Lipschitz nonlinearities (see [START_REF] Andreianov | Discrete duality finite volume schemes for doubly nonlinear degenerate hyperbolic-parabolic equations[END_REF]) Consider the case where ϕ := ϕ • b -1 is a uniformly continuous function (moreover, it is non-decreasing). Let π be a concave modulus of continuity for ϕ • b -1 , Π be its inverse, and set Π (r) = r Π (r). Let π be the inverse of Π . Note that π is concave, continuous, and π(0

) = 0. Set u δ = b(v h )(t+δ , x) and u = b(v h )(t, x). We have Q | ϕ(u δ )-ϕ(u)| = Q π • Π | ϕ(u δ )-ϕ(u)| ≤ |Q| π 1 |Q| Q Π | ϕ(u δ )-ϕ(u)| . Since | ϕ(u δ ) -ϕ(u)| ≤ π(|u δ -u|), we have Π (| ϕ(u δ ) -ϕ(u)|) ≤ |u δ -u| and Π (| ϕ(u δ ) -ϕ(u)|) = Π (| ϕ(u δ ) -ϕ(u)|)| ϕ(u δ ) -ϕ(u)| ≤ |u δ -u| | ϕ(u δ ) -ϕ(u)|.
Therefore, (5) implies an L 1 estimate of the kind (2) on w h = ϕ(v h ):

Q |w h (t+δ )-w h (t)| ≤ |Q| π 1 |Q| Q |u δ -u|| ϕ(u δ )-ϕ(u)| = |Q| π 1 |Q| ω(δ ) .
• Use of contraction arguments and absorption terms (see [START_REF] Andreianov | Convergence of approximate solutions to an elliptic-parabolic equation without the structure condition[END_REF]) Let us mention one more possibility for getting estimates of kind ( 2) for (1), which takes advantage of the monotonicity of ψ(•). Assume ϕ = Id in (1); to shorten the arguments, assume f = 0. Then L 1 translates in time of u h = b(v h ) can be estimated with every of the two preceding methods, the Kruzhkov lemma and a direct estimation of translates with variational techniques. This makes (b(v h )) h relatively compact; yet, when b -1 (•) is discontinuous, no information on compactness of (v h ) h is obtained this way. Now, let us use the translation (in time) invariance of the equation and the L 1 contraction property 11 natural for (1). This yields the estimate (see [START_REF] Andreianov | Convergence of approximate solutions to an elliptic-parabolic equation without the structure condition[END_REF])

Ω |b(v δ h ) -b(v h )|(T -δ ) + T-δ s Ω |ψ(v δ h ) -ψ(v h )| ≤ Ω |b(v δ h )(s) -b(v h )(s)| (6) for all s ∈ (0, T -δ ), where v δ h (t) = v h (t +δ ).
Integrating in s > α > 0, using the time translation bound for (b(v h )) h we get an L 1 ((α, T ) × Ω ) estimate of time translates of ψ(v h ). If ψ(•) is strictly increasing, this is enough for L 1 loc compactness of (v h ) h .

Applications to (1) and some other parabolic PDEs

• Application to a parabolic-hyperbolic PDE (see [START_REF] Andreianov | Discrete duality finite volume schemes for doubly nonlinear degenerate hyperbolic-parabolic equations[END_REF]) For problem (1) with b = Id, provided L p (Q) estimates of the discrete gradient of ϕ(v h ) are available, space translates of ϕ(v h ) (and the functions ϕ(v h ) themselves) can be estimated uniformly, and an estimate of the form ( 5) can be obtained. Then the above technique for exploiting (5) assesses the L 1 (Q) compactness of (ϕ(v h )) h , which is a first step of the convergence proof for this problem (see [START_REF] Andreianov | Discrete duality finite volume schemes for doubly nonlinear degenerate hyperbolic-parabolic equations[END_REF]) 12 .

• Application to an elliptic-parabolic PDE with the structure condition (see [START_REF] Andreianov | On 3D DDFV discretization of gradient and divergence operators. II. Discrete functional analysis tools and applications to degenerate parabolic problems[END_REF]) Assume ϕ = Id. Estimate (5) controls the L 1 time translates of b(v h ) similarly to what was described above 13 . If the structure condition

G(v) = F(b(v)) is satisfied, compactness of (b(v h ))
h is enough to pass to the limit, see [START_REF] Alt | Quasilinear elliptic-parabolic differential equations[END_REF] (cf. [START_REF] Ph | Sur un problème parabolique-elliptique[END_REF] and § 4).

• Application to a cross-diffusion system (see [START_REF] Andreianov | Analysis of a finite volume method to solve a cross-diffusion population system[END_REF])

The following kind of models comes from population dynamics:

u t -D 1 ∆ u -div (u + v)∇ ∇ ∇u + u∇ ∇ ∇v = u(a 1 -b 1 u -c 1 v), v t -D 2 ∆ v -div v∇ ∇ ∇u + (u + v)∇ ∇ ∇v = v(a 2 -b 2 u -c 2 v). (7)
Natural estimates for approximate solutions of ( 7) are L 2 bounds on

√ 1 + u + v ∇ ∇ ∇u, √ 1 + u + v ∇ ∇ ∇v;
this gives only an L 4/3 bound on the diffusion fluxes in [START_REF] Andreianov | Analysis of a finite volume method to solve a cross-diffusion population system[END_REF], thus we are not in a variational setting 14 . Therefore for a proof of convergence of finite volume approximations of the kind [START_REF] Eymard | Finite Volume Methods[END_REF], the Kruzhkov lemma was used in [START_REF] Andreianov | Analysis of a finite volume method to solve a cross-diffusion population system[END_REF] 15 .

• Application to convergence of some linearized implicit schemes (see [START_REF] Andreianov | Convergence of Discrete Duality Finite Volume schemes for the macroscopic bidomain model of the heart electric activity[END_REF]) In [START_REF] Andreianov | Convergence of Discrete Duality Finite Volume schemes for the macroscopic bidomain model of the heart electric activity[END_REF], discretization of the simplified version of cardioelectrical bidomain model:

v t -div M i (•)∇ ∇ ∇u i + H(v) = I ap (•), v t + div M e (•)∇ ∇ ∇u e + H(v) = I ap (•), v = u i -u e , (8)
was considered; here, the "ionic current" H(•) is a cubic polynomial. This nonlinear 11 Discrete version of (6) (see [START_REF] Andreianov | Convergence of approximate solutions to an elliptic-parabolic equation without the structure condition[END_REF]) assumes the L 1 contraction property (linked to order preservation via the Crandall-Tartar lemma) is preserved at the discrete level. Estimate ( 6) is exploited in § 4. 12 For Lipschitz ϕ(•), also the Aubin-Lions-Dubinskii-Simon and Kruzhkov lemmas could be used. For general ϕ(•), the author thinks that neither of these lemmas can replace the direct use of (5). 13 Alternatively, the Kruzhkov lemma can be used in a straightforward way, see [START_REF] Andreianov | On 3D DDFV discretization of gradient and divergence operators. II. Discrete functional analysis tools and applications to degenerate parabolic problems[END_REF]. 14 From the practical point of view, e.g. the first equation cannot be tested with u(t+δ ). 15 Alternatively, the discrete Aubin-Lions-Dubinskii-Simon lemma (see [START_REF] Gallouët | Compactness of discrete approximate solutions to parabolic PDEs -Application to a turbulence model[END_REF]) could be used here.

reaction term brings an estimate of vH(v) which bounds v in L 4 (Q). Time-implicit DDFV discretization of (8) preserves this structure; then the problem falls into the "variational" framework 16 and time translates can be estimated like in [START_REF] Alt | Quasilinear elliptic-parabolic differential equations[END_REF][START_REF] Eymard | Finite Volume Methods[END_REF][START_REF] Eymard | Convergence of a finite volume scheme for nonlinear degenerate parabolic equations[END_REF].

From the practical point of view, it is important to accelerate computations, and to consider a linearized method where the discretization of the reaction term H(v) is not fully implicit. Unfortunately, for theoretical analysis L 4 estimate for v h is not available any more; only a weaker estimate can be obtained with interpolation arguments. In [START_REF] Andreianov | Convergence of Discrete Duality Finite Volume schemes for the macroscopic bidomain model of the heart electric activity[END_REF], we applied the Kruzhkov lemma to exploit this weaker estimate 17 .

Advanced use of the underlying PDE features

Often mere functional-analytic bounds are not enough, but additional constraints coming from the particular structure of the approximated PDE may permit an indirect compactness/convergence proof. E.g., for the parabolic-hyperbolic PDE (1) (case b = Id) we proved the compactness of (ϕ(v h )) h in § 3. The two final steps (see [START_REF] Andreianov | Discrete duality finite volume schemes for doubly nonlinear degenerate hyperbolic-parabolic equations[END_REF]; see also [START_REF] Eymard | Convergence of a finite volume scheme for nonlinear degenerate parabolic equations[END_REF]) exploit fine PDE tools. First, the Minty argument (see, e.g., [START_REF] Alt | Quasilinear elliptic-parabolic differential equations[END_REF]) is used for (a 0 (∇ ∇ ∇ h ϕ(v h )) h ; second, the "nonlinear weak-* convergence" ( [START_REF] Eymard | Finite Volume Methods[END_REF][START_REF] Eymard | Convergence of a finite volume scheme for nonlinear degenerate parabolic equations[END_REF][START_REF] Andreianov | Discrete duality finite volume schemes for doubly nonlinear degenerate hyperbolic-parabolic equations[END_REF]) for (v h ) h is upgraded to strong convergence using entropy inequalities for (1).

Let us show how one very delicate case of (1), see [START_REF] Ph | Sur un problème parabolique-elliptique[END_REF], can be treated indirectly.

Compactness from monotone penalization and order-preservation

For getting (6), we already used the order-preservation structure for (1). Its further use, in conjunction with penalization, may lead to the following convergence proof.

• The structure needed for compactification Assume that one can prove uniqueness of a solution to a PDE (Eq 0 ) under study. Assume that (Eq 0 ) can be embedded "continuously" into a family (Eq ε ) of perturbed PDEs having the property that v ε 1 h ≤ v ε 2 h when ε 1 ≤ ε 2 , where v ε 1 h , v ε 2 h are the associated discrete solutions. Continuity in ε ∈ [-1, 1] means, we assume that limits as ε → 0 (if they exist) of exact solutions v ε of (Eq ε ) solve the limit equation (Eq 0 ).

Assume that for ε = 0, the corresponding sequence (v ε h ) h is well defined and it converges to an exact solution v ε of (Eq ε ). Then solutions (v 0 h ) h to the discretized equation (Eq 0 ) converge a.e., as h → 0, to the unique solution of (Eq 0 ). Indeed, write

v -1 h ≤ v -1/2 h ≤ ... ≤ v -1/m h ≤ ... ≤ v 0 h ≤ ... ≤ v 1/m h ≤ ... ≤ v 1/2 h ≤ v 1 h , (9)
and pass to the limit as h → 0 to define v ±1/m := lim h→0 v ±1/m h (up to extraction of a subsequence) solution to (Eq ±1/m ); then, (9) is inherited at the limit (except that (v 0 h ) h may not have a limit). By monotonicity, we can define v := lim m→∞ v -1/m and v := lim m→∞ v 1/m ; furthermore, we have v ≤ lim inf h→0 v 0 h ≤ lim sup h→0 v 0 h ≤ v. Both v, v solve (Eq 0 ). Thus, by uniqueness, (v 0 h ) h converges to v ≡ v the solution of (Eq 0 ).

• Application to an elliptic-parabolic PDE without the structure condition (see [START_REF] Andreianov | Convergence of approximate solutions to an elliptic-parabolic equation without the structure condition[END_REF])

We assume that ϕ = Id, ψ = 0 in (1). We have seen that compactness of (b(v h )) h can be established, e.g., with the Kruzhkov lemma. Under the structure condition G(v) = F(b(v)), this is enough to pass to the limit in the equation. But in general (see [START_REF] Ph | Sur un problème parabolique-elliptique[END_REF]) one lacks control of time oscillations of G(v h ), and the method of [START_REF] Alt | Quasilinear elliptic-parabolic differential equations[END_REF] fails. Yet it is enough to add penalization term of the form ψ ε (v) = ε(arctan v ∓ π 2 sign ε) to get into the setting where (6) can be exploited to control discrete solutions (v ε h ) h and to pass to the limit, as h → 0, for the ψ ε -penalized equation (1 ε ). The orderpreservation assumptions of the above method being fulfilled due to the choice of ψ ε , we get convergence of (v h ) h in the cases where uniqueness for (1) can be shown.

See[START_REF] Andreianov | Well-posedness results for triply nonlinear degenerate parabolic equations[END_REF] and references therein for well-posedness theory of such "triply nonlinear" equations. These are mathematical models for porous media, sedimentation, Stefan problem, etc..

Actually, we assume that either these functions are uniformly continuous, or v is bounded a priori.

When the compactification methods strongly utilize a particular structure of the underlying PDE, this step is in fact combined with the third step of passing to the limit.

Throughout the note, "compact" actually signifies "relavively compact".

Working in an h-independent space is an advantage for producing discrete versions of compactness arguments; yet the approach of[START_REF] Gallouët | Compactness of discrete approximate solutions to parabolic PDEs -Application to a turbulence model[END_REF] exhibits a simple and efficient use of h-dependent spaces.

Notice that for evolution PDEs governed by accretive in L 1 (Ω ) operators, of which (1) is an example, time-implicit discretizations are better suited for structure preservation. Use of numerical schemes in space that possess a kind of discrete duality (mimetic, co-volume, DDFV schemes, etc.) enables getting discrete estimates analogous to the continuous ones. For notions of solution involving some version of chain rule (e.g., entropy, renormalized solutions) orthogonality assumption on the meshes and isotropy assumption on the diffusion operator may be needed, see e.g.[START_REF] Andreianov | Discrete duality finite volume schemes for doubly nonlinear degenerate hyperbolic-parabolic equations[END_REF].

[START_REF] Ph | Sur un problème parabolique-elliptique[END_REF] Here, for the sake of simplicity, we stick to the terminology and notation of the continuous case.

Indeed, we have v h bounded in L 4 (Q) and H(v h ) is bounded in L 4/3 (Q) = (L 4 (Q)) * .

A discrete Aubin-Lions-Dubinskii-Simon argument could have been applied as well.

Acknowledgements The author thanks E. Emmrich for discussions the above techniques.