Christophe Cheverry

CAN ONE HEAR WHISTLER WAVES ?

Keywords: Rotating fluids, dynamical and hamiltonian systems, gyro-kinetic equations, geometric optics, WKB calculus, oscillatory integrals, wave-particle interaction, plasma physics, relativistic Vlasov-Maxwell equations, whistler waves

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction.

Whistlers [START_REF] Park | Whistler observations during a magnetospheric sudden impulse[END_REF][START_REF] Helliwell | Whistlers and Related Ionospheric Phenomena[END_REF] are very low frequency (VLF) electromagnetic (radio) waves generated by lightning [START_REF] Green | Lightning Effects on Space Plasmas and Applications[END_REF]. They were described long ago, by Storey [START_REF] Storey | An investigation of whistling atmospherics[END_REF] in 1953, as the sound of the dawn chorus. It was then suggested by Dungey [START_REF] Dungey | Loss of van allen electrons due to whistlers[END_REF] in 1963 that the natural or artificial emissions of whistler waves may be due to the oscillations in density of the charged particles trapped in the Van Allen belts [START_REF] Braun | Mathematical remarks on the Van Allen radiation belt: a survey of old and new results[END_REF]. This intuition has demonstrated its usefulness over the years [START_REF] Li | Modulation of whistler mode chorus waves: 2. role of density variations[END_REF][START_REF] Pandey | Pitch angle loss-cone anisotropic magneto plasma in presence of parallel electric a.c. field[END_REF][START_REF] Serra | Nonlinear shift of wave parameters[END_REF]. This field is currently the subject of intensive research [START_REF] Bortnik | The unexpected origin of plasmaspheric hiss from discrete chorus emissions[END_REF][START_REF] Ciurea-Borcia | Oblique whistler waves generated in cold plasma by relativistic electron beams[END_REF][START_REF] Santolik | New results of investigations of whistler-mode chorus emissions[END_REF][START_REF] Piel | Plasma physics: an introduction to laboratory, space, and fusion plasmas[END_REF][START_REF] Weia | Generation mechanism of the whistler-mode waves in the plasma sheet prior to magnetic reconnection[END_REF][START_REF] Yamaguchi | Ray tracing of whistler-mode chorus elements[END_REF] stemming from the physics community.

Many works have been put forth to explain what happens. Arising in connection with progresses in satellite observations (whistler waves are regularly detected by Cluster and Themis during substorms), there is now a number of contributions on this topic. The articles come from various scientific fields. They give rise to different interpretations of the phenomena. For a comprehensive overview, please refer to the recent PhD thesis in applied mathematics of S. Le Bourdiec [START_REF] Bourdiec | Méthodes déterministes de résolution des équations de Vlasov-Maxwell relativistes en vue du calcul de la dynamique des ceintures de Van Allen[END_REF], in astrophysics of A. Tenerani [START_REF] Tenerani | Dynamics of ion-scale coherent magnetic structures and coupling with whistler waves during substorms[END_REF] and in plasma physics of P. Pfannmöller [START_REF] Pfannmöller | Whistler Wave Propagation in Inhomogeneous Plasmas[END_REF]. These documents provide interesting investigations by means of observations, theoretical modelings and numerical studies. The references cited therein furnish up-to-date information.

The outer Van Allen belt is a layer of energetic electrons surrounding Earth. It is almost toroidal in shape. It is built with a cold magnetized plasma. It may be described by the relativistic Vlasov-Maxwell system. More precisely, the evolution in time τ ∈ R + (with τ ∼ 1 representing 10 3 s), in position x ∈ R 3 , and in velocity v ∈ R 3 of particles and plasma waves (in particular of whistler waves) is governed by: (1.1)

       ∂ τ f + ε -1 v • ∇ x f + ε -2 v × Be (x) • ∇ p f + (E + v × B) • ∇ p f = -ε -1 κ ∇ p M θ (p) • E , f |τ =0 ≡ f , ∂ τ E -ε -1 curl B = -ε -1 κ j , div E = κ ρ , E |τ =0 ≡ 0 , ∂ τ B + ε -1 curl E = 0 , div B = 0 , B |τ =0 ≡ 0 .
These equations are in dimensionless variables. The reader is referred to Paragraph 3.1.1 for a complete justification of the model (1.1). Keep in mind that the function M θ (•) is the Maxwell-Jüttner distribution, see (3.5)- (3.6); the number κ measures the strength of the coupling between the Vlasov and Maxwell parts, see (3.9); the small parameter ε := 10 -5 is the inverse of the electron gyrofrequency (relative to the reference frequency), see (3.9). The velocity v is limited by the speed of light, 1 in dimensionless variables. It is linked to the relativistic impulsion p through the relations:

(1.2) v(p) = p 1 + |p| 2 -1/2 , p(v) = v 1 -|v| 2 -1/2 .
The electric current density j and the electric charge density ρ are given by:

(1.3) j := R 3 f v dp , ρ := R 3
f dp , dp := (1 -|v| 2 ) -5/2 dv .

From now on, the symbol Be stands for the dipole model of the Earth's magnetic field. It bears noting that Be (•) is a non-constant function given by formulas (2.7) and (3.1). In practice, the contribution v × Be (x) reflects how, at position x, the external geomagnetic field Be (x) deflects the charged particles through the Lorentz force. The physical quantities to determine are the distribution function f (which indicates the particle density), together with the self-consistent electromagnetic field (E, B).

The large factor ε -1 = 10 5 is omnipresent in (1.1). Its occurrence in multiple places induces a full range of oscillations. These oscillations can interact through complex processes, whose mathematical study is interesting. Indeed, this should provide a solid foundation for a better understanding of wave-particle interactions. The corresponding program of research is still in its early stages. An important step is to exhibit the key mechanisms related to the production of plasma waves, like whistler waves.

To move in this direction, our strategy is to solve the Vlasov part independently, with a given field (E, B), typically with (E, B) = (0, 0). This can be done by the characteristic method. The solutions f ε thus obtained develop oscillations which take a very long time to appear. It starts with a transition from t := ε -1 τ ∼ ε to t ∼ 1 with fast rotations; and it goes on from t ∼ 1 to t ∼ ε -1 with again fast rotations (but of a different kind). During this maturation process, the features of the field Be (•) have a full impact. They result in special oscillating structures, based on kinetic phases, say ϕ 0 (t, x, v). Retain that ∇ v ϕ 0 ≡ 0, and that the oscillating wave front set induced by ϕ 0 implies specific geometries.

The above expressions f ε generate electric currents j ε . It suffices to apply (1.3). As this happens, the asymptotic behaviours of f ε are converted into time-space oscillations, say along T r(ϕ 0 )(t, x). These new oscillations are viewed as source terms at the level of (1.1). They can be absorbed or, on the contrary, they are propagated in the form of wave packets. Since ∇ v T r(ϕ 0 ) ≡ 0, after modelling, the selection process can be identified by looking at equations of Maxwell type, which are issued from (1.1). This is the Maxwell's system (3.20), where the relative permittivity ǫ r reflects the influence of Be . This leads to the study of a class of Fourier integral operators. Then, what matters is the crossing over time, in the cotangent space T * (R 3 ), of two geometrical objects: the Lagrangian submanifold generated by T r(ϕ 0 )(t, •) and the characteristic variety associated with (3.20).

This article is organized in accordance with the preceding guidelines. Part 2 is devoted to the study of the dynamical system (1.5) issued from the Vlasov equation. The main achievement is the derivation of long time gyro-kinetic equations. As a result, it becomes possible to determine effective conditions to impose on a strong inhomogeneous exterior magnetic field to obtain the long time dynamical confinement of plasmas. These conditions are met in the case of Be (•). However, they are not sufficient. Indeed, some self-consistent electromagnetic field (E, B) is produced at the level of (1.1). This takes the form of plasma waves, whose impact can destroy the dynamic stability. It is therefore also crucial to understand the basic mechanisms underlying the creation of such wave packets. In Part 3, this issue is investigated in the case of a weak coupling, and only concerning a specific kind of light. We fix j ε as above, and we consider (3.20). We focus on the whistler dispersion relation, and we recognize the essential features of whistler-mode chorus emissions.

1.1. Long time gyro-kinetic equations. The Vlasov part of (1.1) is the first equation of (1.1), the one on f . It controls the dynamics of particles. Assuming that κ = 0 and that the function (E, B)(ε, τ, x) is known, it can be solved by the characteristic flow:

(1.4) f ε (τ, x, v) = f x ε (-τ, x, v), v ε (-τ, x, v) , f ε (0, x, v) := f(x, v) ,

where (x ε , v ε )(τ, x, v) ∈ R 3 × R 3 is obtained through:

(1.5)

∂ τ x ε = ε -1 v ε , x ε (0, x, v) = x , ∂ τ v ε = ε -2 c(|v ε |) v ε × B + Ñ (v ε , E, B) , v ε (0, x, v) = v ,
with |v| < 1 and coefficients adjusted by taking into account (1.2):

(

1.6) c(v) := (1-|v| 2 ) 1/2 , Ñ (v, E, B) := (1-|v| 2 ) 1/2 E -(E •v) v +v ×B .
The Hamiltonian system (1.5) is proving to be a small perturbation of an integrable system. The asymptotic study of (1.5) falls under the scope of KAM or gyrokinetic theory [START_REF] Braun | Mathematical remarks on the Van Allen radiation belt: a survey of old and new results[END_REF][START_REF] Brizard | Foundations of nonlinear gyrokinetic theory[END_REF][START_REF]Topics in kinetic theory[END_REF]. There are three adiabatic invariants, providing a qualitative description of the motions. Practical experience reveals the following.

⋄ Physical observations. Recall that t := ε -1 τ . Moreover (Subsection 3.1.1), expressed in seconds, the observation time is t, and we have t := T t with T ≃ 10 -2 s. The trajectories associated with (1.5) are the combination of three basic motions:

(P1) Fast times (t ∼ ε or t ∼ 10 -7 s): The particles oscillate at the cyclotron frequency.

They move approximately in circles, which are contained in the orthogonal planes to the field lines, and which have a gyroradius (or Larmor radius) of size ∼ ε ≪ 1. The centers of these circles, sometimes called gyrocenters, are guided along the field lines;

(P2) Intermediate times (t ∼ 1 or t ∼ 10 -2 s): The centers follow the field lines. The particles bounce back and forth between two mirror points located just above the Earth's poles, crossing each time the equatorial plane;

(P3) Long times (t ∼ ε -1 or t ∼ 10 3 s): A longitudinal drift (initiation of a rotation around the Earth) becomes perceptible. It is associated with the creation of a ring current. Positive particles go towards west, negative towards east. ⋄

These motions have long been of interest to both mathematicians and physicists. Their description raises difficulties known under the name of the Störmer problem. This subject has been intensively studied; see the survey article [START_REF] Braun | Mathematical remarks on the Van Allen radiation belt: a survey of old and new results[END_REF] to have an overview. There are actually continuations in the direction of chaos [START_REF] Dilão | Differential equations, chaos and variational problems[END_REF]. The whole picture gives rise to toroidal regions surrounding the earth, the so-called Van-Allen belts, where the particles are trapped and oscillate with a lifetime from a few hours to 10 years. The hamiltonian approach is more concerned with the Lagrangian point of view. It does not offer much insight about the variations which can occur (at fixed times and at given spatial scales) when passing from a trajectory to one another. Such information is encoded elsewhere, in the Eulerian concept of what is the phase of an oscillation. Now, to capture the phases generated by (1.5) during long times τ ∼ 1, it is necessary to conveniently separate the scales. This means to revisit the Störmer problem by the way of modern averaging methods. It is this program that is achieved in Part 2. Section 2.1 starts with Assumption 2.1 allowing to straighten out in an orthogonal way the field lines associated with Be (•). Once that is accomplished, and as long as t ∼ 1, the asymptotic analysis of (1.5) relies on adaptations of the filtering method [START_REF] Chemin | An introduction to rotating fluids and the Navier-Stokes equations[END_REF][START_REF] Schochet | Fast singular limits of hyperbolic PDEs[END_REF], in the spirit of [START_REF] Gallagher | On pressureless gases driven by a strong inhomogeneous magnetic field[END_REF]. The situation during intermediate times t ∼ 1 is similar to the usual reports in rotating fluids [START_REF] Chemin | An introduction to rotating fluids and the Navier-Stokes equations[END_REF][START_REF] Frénod | Long time behavior of the two-dimensional Vlasov equation with a strong external magnetic field[END_REF][START_REF] Gallagher | On pressureless gases driven by a strong inhomogeneous magnetic field[END_REF][START_REF] Grenier | Pseudo-differential energy estimates of singular perturbations[END_REF][START_REF] Schochet | Fast singular limits of hyperbolic PDEs[END_REF]. The related phenomena are discussed in almost all plasma physics books [START_REF] Howard | Introduction to Plasma physics[END_REF][START_REF]Topics in kinetic theory[END_REF][START_REF] Piel | Plasma physics: an introduction to laboratory, space, and fusion plasmas[END_REF], in a separate chapter, about charged particle motion. Changing t into τ produces new problematical terms of size ε -1 . To partially absorb them, a rectification procedure is adopted in Section 2.2. This requires a consideration of the homological equation (2.31). This amounts to solve the first modulation equation governing the evolution during intermediate times. There is a corresponding flow map Ξ.

In a convenient choice of coordinates, the phase space (of dimension 6) is foliated by planes which are globally invariant under the action of Ξ. Each of these planes is equipped for Ξ with an energy function H, called the reduced hamiltonian (Definition 2.7).

Denote by C e (x) the mean curvature at x of the surface containing x and orthogonal to the field line passing over x. When the function C e (•) changes its sign along the field lines, the reduced hamiltonian H is a potential well. This provides a first confinement criterion. This can take place only if Be (•) undergoes variations in directions. This is especially so in the case of the geomagnetic field Be (•), when passing the equatorial plane.

The first modulation equation (2.31) is a differential equation. Its solutions follow the level curves of H, which are diffeomorphic to circles (at least when H is a potential well). The corresponding trajectories are therefore periodic. This is the origin of the notion of extended Van Allen belt EBēlt, which is clarified in Paragraph 2.2.4. Here is an explanation for the step (P2). Here lies also a key to derive uniform estimates in t (eliminating, for instance, problems of secular growth). We refer the reader to the introductory paragraph 2.2.1 for more comments, and to the whole Section 2.2 for a detailed analysis. After rectification, the system (1.5) can be written in the form (2.22). There remains terms of size ε -1 , but they act with a sort of decoupling. In Section 2.3, the discussion falls within Wentzel-Kramers-Brillouin (WKB) analysis [START_REF] Métivier | The Mathematics of Nonlinear Optics[END_REF][START_REF] Rauch | Hyperbolic Partial Differential Equations and Geometric Optics[END_REF]. More precisely, it combines techniques involving three scale expansions (inspired from diffractive geometric optics) and tools from [START_REF] Cheverry | Cascade of phases in turbulent flows[END_REF]. This includes formal computations, beginning in Paragraph 2.3.2 with the derivation of the second modulation equation (2.84), playing the part of long time gyro-kinetic equations. This is again a differential equation. The result is a second confinement criterion, requiring that all trajectories are bounded. This condition is met when dealing with Be (•). It is only then that a precise description of the step (P3) becomes accessible.

After that, there is need for more complete formal expansions in powers of ε.

E, B)(ε, τ, x) ∈ C ∞ ([0, 1] × R + × R 3 ; R 6 ) is a known given function.
Consider the centered dipole model (2.7) for the geomagnetic field Be . There is a time T ∈ R * + and a domain EBēlt of the phase space (x, v) ∈ R 3 × R 3 , which is relatively compact in (R 3 \ {0}) 2 , such that :

1) the projection of EBēlt onto the x coordinate is a toroidal region of R 3 which can be made conform to the shape of the Van Allen belt, designated by Bēlt; 2) for all (x, v) ∈ EBēlt, the solution

(x ε , v ε )(τ, x, v) to (1.5) remains for all τ ∈ [0, T ] in a neighbourhood of EBēlt, that is compact in (R 3 \ {0}) 2 ; 3) there are profiles X j (τ, x, v, ϑ, θ) and V j (τ, x, v, ϑ, θ) in C ∞ ([0, T ]×R 3 ×R 3 ×T 2 ; R 3 )
with j ∈ N such that, for all N ∈ N * , the asymptotic behaviour when ε goes to zero of the solution

(x ε , v ε )(τ, x, v) to (1.5) can be described with a precision O ε N -2 in L ∞ ([0, T ] × EBelt; R 6
) through the following multiscale and multiphase expansion:

(1.7) x ε v ε (τ, x, v) ∼ ε → 0 N j=0 ε j X j V j τ, x, v, ψ e s (τ, x, v) ε , ψ e pε (τ, x, v) ε 2 .
The secondary phase ψ e s (τ, x, v) does not depend on (E, B). It is specified in Definition 2.6, together with (2.115). The principal phase ψ e pε (τ, x, v) is an oscillation of small amplitude:

(1.8) ψ e pε (τ, x, v) := ψ e p0 (τ, x, v) + ε ψ e p1 τ, x, v, ψ e s (τ, x, v) ε .
Explicit formulas for the scalar functions ψ e p0 (τ, x, v) and ψ e p1 (τ, x, v, ϑ) are provided in the Definitions 2.6 and 2.7. For all k ∈ N, the C k -norms of the profiles (X j , V j ) with j ∈ N, of ψ e p0 and of ψ e p1 can be controlled by the C k -norm of (E, B). Beyond the case of the earth's magnetic field, the analysis of Part 2 should be the entry point into many new, interesting and challenging problems in physics, since it furnishes a systematic method to study very long-time effects induced by the cross product with a non constant vector field.

1.2. Wave-particle interaction in cold magnetized plasmas. A whistler wave is a very low frequency electromagnetic wave generated by lightning [START_REF] Green | Lightning Effects on Space Plasmas and Applications[END_REF][START_REF] Helliwell | Whistlers and Related Ionospheric Phenomena[END_REF][START_REF] Santolik | New results of investigations of whistler-mode chorus emissions[END_REF]. Whistler waves can be detected through spectrograms, giving a visual representation of the spectrum of frequencies in the signal as they vary with time. ⋄ Physical observations. It is worth pointing out that the whistler-mode signals recorded by satellite observations inside spectrograms share the following main characteristics: (C1) A broad (VLF) frequency range (from 100 Hz to 10 kHz) is implied. Within this range, chorus waves can occur in distinct frequency bands; (C2) The precise description of a chorus element can show a succession of wave packets with rising frequencies, see the sample produced in figure 1 p. 623 of [START_REF] Santolik | New results of investigations of whistler-mode chorus emissions[END_REF]. (C3) Whistler-mode chorus emissions consist of a succession of chorus elements, occurring with a repetition period of about or less than 1 second, marked by light-coloured vertical lines in the related spectrograms.

⋄ The present analysis is based on this concrete information. From (C1), we will deduce that some semiclassical regime is at stake, with an adequate dispersion relation. From (C2), we will infer the presence of harmonics, indicating nonlinear effects. But our foremost concern is to explain (C3). Though the generation mechanism of chorus alluded in (C3) has been extensively studied in the past, it is not yet well understood. Please refer to [START_REF] Santolik | New results of investigations of whistler-mode chorus emissions[END_REF] for recent results obtained during investigation in the Earth's magnetosphere. Part 3 gives a plausible mathematical interpretation of what happens concerning (C3). Subsection 3.1 is devoted to the modelling of whistler waves. It begins in Paragraph 3.1.1 with the nondimensionalization of the original equations, those considered in [START_REF] Bourdiec | Méthodes déterministes de résolution des équations de Vlasov-Maxwell relativistes en vue du calcul de la dynamique des ceintures de Van Allen[END_REF][START_REF] Pfannmöller | Whistler Wave Propagation in Inhomogeneous Plasmas[END_REF][START_REF] Tenerani | Dynamics of ion-scale coherent magnetic structures and coupling with whistler waves during substorms[END_REF]. Paragraph 3.1.2 outlines the way whistler waves are usually tackled. It is followed by the presentation of a quite different perspective. Unusual aspects come from the introduction of the parameter ε, from the variations of Be (•), and from the localization of the data (Paragraph 3.1.3). Most important of all, however, is the change in perspective subsequently introduced. Paragraphs 3.1.4, 3.1.5 and 3.1.6 propose to explain wave-particle interaction in terms of oscillations. They outline a mechanism of creation of light. This is the passage from ϕ 0 (t, x, v) to T r(ϕ 0 )(t, x), and then from T r(ϕ 0 )(t, x) to the production of wave packets.

On the one hand, the fast dynamics underlying the system of ODEs (1.5) serve to force the oscillating structures of f ε . They result in some macroscopic organization of particles, that is reflected in ϕ 0 (t, x, v). Computing the electric current density j ε , this yields T r(ϕ 0 )(t, x). On the other hand, as explained in Paragraph 3.1.6, the frequencies thus microlocalized are managed by a dispersion relation [START_REF] Helliwell | Whistlers and Related Ionospheric Phenomena[END_REF][START_REF] Piel | Plasma physics: an introduction to laboratory, space, and fusion plasmas[END_REF] which must take into account the presence of a non-homogeneous magnetized plasma. Combining these two aspects leads to a class of oscillatory integrals, introduced in Paragraph 3.1.7, involving phases Ψ(±, q; t, x; s, y, ς, ξ) whose principal characteristics can be viewed at the level of (3.39).

At this stage, the main question to address is to obtain stationary phase expansions. As noted in the introduction of Section 3.2, this relies on a study of the (interior and boundary) critical points associated with the scalar function Ψ(±, q; t, x; •). The behaviour of Ψ with respect to the time variable s has a remarkable property. Indeed, the preliminary estimates of Paragraph 3.2.1 reveals that Ψ(±, q; t, x; •, y, ς, ξ) is composed of a linear part in s (with a nonzero coefficient) plus some periodic part.

From that, it is possible to show the absence of boundary critical points (Paragraph 3.2.2). The number of interior critical points of Ψ(±, q; •) detected for different values of q near some (adequate) time t is finite (Paragraph 3.2.3). These critical points correspond to the emission of signals. In line with (C1), they are probably at the source of the wave packets composing a chorus element. In Paragraph 3.2.4, we focus on the spatial phase Ψ(±, q; t, x; s, y, π/2, ξ) by pointing out that the signals can only emanate from the mirror points. Then, a distinction must be drawn between two different situations. In the end, what matters are the respective positions of two quantities, which are relevant from a physics point of view. The first quantity is the maximum amplitudes b M of the geomagnetic field along trajectories (see Paragraph 3.2.1 for a precise definition); the second quantity is the threshold limit λ M issued from the dispersion relation (see Assumption 3.1).

The case c(w) b M < λ M is investigated in Paragraph 3.2.5. The other case c(w) b M > λ M is examined in Paragraph 3.2.6. In both situations, the aim is to identify the successive times s j (with j ∈ N) when the signals composing a chorus element are emitted. Depending on the situation, the episodes of whistler-mode chorus emissions may be shorter or longer. But, they all have one thing in common. As reflected in (C3), the signals appear again and again, over consecutive periods of almost equal length. This aspect is roughly highlighted below, through some Informal Statement 1. It is present in the background of Propositions 3.1 and 3.2, which contain the complete and rigorous Theorems. Further information about the interpretation of all features (C⋆) are available in Paragraph 3.2.7.

Informal Statement 1 (whistler-mode chorus emissions). The phase Ψ involves several critical points which can be classified by order of appearance s j , with s j ≤ s j+1 , the majority of which are separated by a uniform interval of time.

Part 3 must be viewed as a preliminary step towards more refined studies. There remain a number of important mathematical issues that need to be addressed. This includes, among other things, the propagation of the wave packets whose creation has been explained (Paragraph 3.2.8), and the problem of nonlinear interactions (Paragraph 3.2.9).

Multiscale analysis of the Vlasov part.

Many works are devoted to the analysis of fast rotating fluids. In dimension three, the penalization matrix can always be represented by the cross product with a vector field B. The most studied case is when B is a constant vector. For a discussion and for a number of references, we can refer to the books [START_REF] Chemin | An introduction to rotating fluids and the Navier-Stokes equations[END_REF] and [START_REF]Topics in kinetic theory[END_REF] (for the gyrokinetic point of view). Look also at the contribution [START_REF] Frénod | Long time behavior of the two-dimensional Vlasov equation with a strong external magnetic field[END_REF] for specific development related to rarefied plasmas. A physically more relevant situation is when B undergoes variations in the spatial position. Decompose such a B(x) into B(x) = b(x) r1 (x), where b(x) and r1 (x) are respectively the amplitude and the unitary direction of B(x). As long as b ≡ 0 and t ∼ 1, without inducing special difficulties, a rotation brings back to the case of a vector field B having a constant direction. This means that, without limiting the generality of the foregoing, the discussion can be achieved with r1 ≡ e, where e ∈ R 3 is a fixed vector of norm 1. Then, the analysis is already more complicated than in the constant case. There are special phenomena at stake, see for instance the articles [START_REF] Cheverry | Semiclassical and spectral analysis of oceanic waves[END_REF][START_REF] Gallagher | On pressureless gases driven by a strong inhomogeneous magnetic field[END_REF]. When looking at the long time evolution, that is for t ∼ ε -1 instead of t ∼ 1, the changes in the unit directions generated by B, that is concerning r1 , also appear to play a crucial role. The reduction to r1 ≡ e is no more possible. Other tools must be developed, and other effects are expected. The aim of this Part 2 is precisely to fill this gap. The starting point is (1.5), where it is assumed that (E, B) is a known smooth function. The matter is to obtain a good description of the asymptotic behaviour of (x ε , v ε )(τ, x, v) when ε goes to 0. Keeping in mind the physical observations (Pj), the discussion will be systematically tested in the earth's context. This will appear throughout the text inside paragraphs that will be introduced by the word: • Application. The subscript "e" will be added each time the geomagnetic field is implied. The preliminary stage is to understand well what happens during intermediate times t ∼ 1. This is achieved in Section 2.1. This allows to explain the contents of Step (P1), and to catch basic features about (P2). Then, in Section 2.2, we consider longer times t ∼ ε -1 . There, the approach is completely new. The related results are in coherence with the observations reported in Step (P2) and beyond, they give a good insight into Step (P3).

2.1.

The averaging method during intermediate times t ∼ 1. Changing the time variable τ into t := ε -1 τ , the dynamical system (1.5) still contains fast rotating terms:

(2.1)

∂ t x ε = v ε , x ε (0, x, v) = x , ∂ t v ε = ε -1 c(|v ε |) v ε × B + ε Ñ (v ε , E, B) , v ε (0, x, v) = v .
Before applying at the level of (2.1) the classical filtering method (Paragraph 2.1.3), it is necessary to straighten out the field lines (Paragraph 2.1.2). To this end, Assumptions on the vector field B are needed (Paragraph 2.1.1). These hypotheses provide significant advantages with a view to achieving the above objectives. Most importantly, they enable to argue globally in space. Technically, their implementation allows in the next Paragraph 2.1.2 to straighten out the field lines in an orthogonal way. Also, they give a direct access to a geometrical interpretation of quantities that will be introduced later on.

2.1.1. Assumptions on the external magnetic field. We work under the following hypotheses, which are both intrinsic [START_REF] Weatherburn | On Lamé families of surfaces[END_REF] and quite general.

Assumption 2.1. [structure of the three-dimensional skew-symmetric matrix] Given some open set Õ ⊂ R 3 , the vector field B : Õ -→ R 3 satisfies the following conditions: (i) The function B is a smooth nonzero function: B ∈ C ∞ Õ; R 3 \ {0} . More precisely, there is a constant c ∈ R * + and a unit vector field r1 ∈ C ∞ ( Õ; S 2 ) such that:

(2.2) B(x) = -b(x) r1 (x) , |r 1 (x)| = 1 , b(x) := | B(x)| > c , ∀ x ∈ Õ .
(ii) There is a global smooth diffeomorphism y : Õ -→ O, defined through the curvilinear coordinates y j : Õ -→ R with j ∈ {1, 2, 3}, such that ∇ x y 1 is everywhere collinear to B, and giving rise to a triply orthogonal system.

Let δ j k be the Kronecker symbol (δ j k = 0 if j = k and δ j k = 1 if j = k). The two functions y 2 and y 3 are sometimes called Clebsch variables. Locally, the condition (ii) says that:

(2.3) ∇ x y j = dj t rj , t rj • rk = δ j k , dj := |∇ x y j | , ∀ (j, k) ∈ {1, 2, 3} 2 .
Exchanging the indices j = 2 and j = 3 if necessary, we can always make sure that the basis (r 1 , r2 , r3 ) is right-handed. The contravariant metric tensor (g ij ) associated with the change y is diagonal, with diagonal coefficients g jj = ( dj ) 2 . Recall that dj = ( dj ) -1 , where the dj are the Lamé coefficients. Adopt the conventions:

(2.4)

r j := rj • y -1 , d j := dj • y -1 , ∀ j ∈ {1, 2, 3} , b := b • y -1 .
As a consequence of (2.3)-(2.4), we have:

(2.5)

D x y • y -1 = D O , D =   d 1 0 0 0 d 2 0 0 0 d 3   , O = t O -1 =   t r 1 t r 2 t r 3   .
We will note y the current variable in the domain O. This variable y must not be confused with the diffeomorphism y : Õ -→ O. From now on, a symbol V representing a vector in R d with d ∈ N * will be denoted by

V = t (V 1 , • • • , V d ).
The geometric interpretation of (ii) is discussed in [START_REF] Weatherburn | On Lamé families of surfaces[END_REF]. The condition (ii) can be tested through differential equations on B.

It allows to incorporate most physical situations, for instance the geomagnetic case.

• Application 1 [the geomagnetic field]. The Van Allen belts are two layers of charged particles (plasmas) that are held in place [START_REF] Braun | Mathematical remarks on the Van Allen radiation belt: a survey of old and new results[END_REF] around the planet Earth by the geomagnetic field Be . Whistler mode waves are principally issued from the outer radiation belt, denoted by Belt in what follows, and simply referred to as the Van Allen belt. In the vicinity of Belt, the function Be can be adequately approximated as a dipole field, whose dipole axis D e is tilted at an angle of 10 degrees with respect to the earth's rotational axis. Let R ≃ 6, 371 × 10 6 m be the earth radius. The cartesian coordinate system x, with dimensionless spatial variable x := x/R, may be adjusted in such a way that D e and the (magnetic) equatorial plane Eq are given by:

D e := x ∈ R 3 ; x 1 = x 2 = 0 , Eq := x ∈ R 3 ; x 3 = 0 .
The domain Belt does not intersect D e . It can be viewed as a volume of revolution with axis D e , which is symmetric with respect to Eq. Moreover, the region Belt is contained in a spherical shell which extends from an altitude of about three to ten Earth radii. This last property leads to the following framework:

(2.6) Belt ⊂ Õe , Õe := x ∈ R 3 ; 3 < |x| < 9 , (x 1 ) 2 + (x 2 ) 2 = 0 .
Let b0 ≃ 3, 1 × 10 -5 T be the size of | Be |, as it is measured at the surface of the Earth. In rescaled units, the expression Be is replaced by Be := Be / b0 . In the conventional cylindrical coordinates ρ, φ and z issued from x, the map Be (•) does not depend on φ ∈ R. It can be described by [START_REF] Braun | Mathematical remarks on the Van Allen radiation belt: a survey of old and new results[END_REF]:

(2.7) Be (ρ, z) = ∇ × Ã = 1 (ρ 2 + z 2 ) 5/2   -3 ρ z 0 ρ 2 -2 z 2   , Ã = -ρ (ρ 2 + z 2 ) 3/2   0 1 0   .
The amplitude be of Be is the function: be (ρ, z)

:= (ρ 2 + z 2 ) -5/2 (ρ 2 + z 2 ) 2 + 3 z 2 (ρ 2 + z 2 ) 1/2 , be (ρ, 0) = ρ -3 .
It is clearly homogeneous of degree -3. This remark combined with (2.6) gives rise to:

(2.8) 0 < b m := 9 -3 ≤ be (x) ≤ b M := 2 × 3 -3 , ∀ x ∈ Õe .
In view of (2.8), the part (i) of Assumption 2.1 becomes obvious. On the other hand, we have div Be = 0, but also curl Be = 0. Thus, the Helmholtz decomposition of Be can be achieved in the form Be = ∇ x y 1 e . This special choice of y 1 e implies that d1 e ≡ be . To see why the remaining conditions in (ii) are verified, it suffices to produce:

(2.9)

y 1 e := z (ρ 2 +z 2 ) -3/2 ∈ R , y 2 e :=ρ 4 (ρ 2 +z 2 ) -3 ∈ R * -, y 3 e := φ ∈ R . With these conventions, the coordinates y 1 , y 2 and y 3 can be interpreted respectively as a sort of latitude, (negative) altitude and longitude. From (2.6), we can infer that:

O e = y e ( Õe ) ⊂ (y 1 , y 2 , y 3 ) ∈ R × R * -× R ; |y 1 | < 3 -1 , -9 -1 < y 2 < 0 . From (2.9), we can also compute ∇ x y j e = dj e t rj e with dj e := |∇ x y j e | for all j ∈ {1, 2, 3}, and check that the conditions t rj e • rk e = δ j k are satisfied. It appears here that a convenient choice for y e : Õe -→ O e can be obtained as the composition of two orthogonal change of variables. First, the one from x to (ρ, φ, z). Then, the one from (ρ, φ, z) to the toroidal coordinates y e = t (y 1 e , y 2 e , y 3 e ) defined in (2.9). In the coordinate system y, the field lines become a family of parallel straight lines:

Fl(ŷ) := (y 1 , ŷ) ∈ R × R * -× R ; y 1 ∈ R , ŷ := (y 2 , y 3
) . The inclusion (2.6) is only indicative. The minimal distance ρ m and the maximal distance ρ M , which can be measured from a particle to the earth's center knowing that this particle is located in Belt ∩ Eq, can be measured more accurately, with 3 < ρ m < ρ M < 9. Remark that y 2 = -ρ -2 when y 1 = 0. Accordingly, define y 2 m := -ρ -2 m and y 2 M := -ρ -2 M . The altitude y 2 being fixed, mark the maximal latitude h e (y 2 ) which can be reached by a particle trapped along the field line Fl(y 2 , y 3 ). The value h e (y 2 ) corresponds also to the extremal position of a mirror point having the altitude y 2 .

All ingredients y 2 m , y 2 M and h e (•) can be adjusted on the basis of concrete observations. They allow to determine precisely the shape of the Van Allen belt, leading to:

(2.10)

Belt = y ; |y 1 | ≤ h e (y 2 ) , h e ∈ C ∞ 0 (R * -; R + ) , supp h e = [y 2 m , y 2 M ] .
Lemma 2.1. The functions r j e := rj e • y -1 e , d j e := dj e • y -1 e and b e := be • y -1 e satisfy: ∂ y j r 3 e = 0 , ∀ j ∈ {1, 2} , (2.11a)

∂ y 3 d j e = 0 , ∀ j ∈ {1, 2, 3} , (∂ y 3 r 1 e ) • r 2 e = 0 , (2.11b) ∂ y 3 b e = 0 , ∂ y 1 b e (0, y 2 ) = 0 , ∂ 2 y 1 y 1 b e (0, y 2 ) > 0 , ∂ y 2 b e < 0 . (2.11c)
The function b e (•, y 2 ) is strictly increasing on R * + . On the other hand, for all y 2 ∈ R * -, the partial applications b e (•, y 2 ) and d j e (•, y 2 ) with j ∈ {1, 2, 3} are even in y 1 ∈ R. Proof. By construction r 3 e = t (-sin y 3 , cos y 3 , 0), explaining (2.11a). From the definitions in (2.9), we can easily deduce that d1 e and d2 e depend only on (ρ, z), whereas d3 e = ρ -1 . This gives rise to the first part of (2.11b). Knowing what is r 1 e for φ = 0, the values of r 1 e for φ = 0 are obtained by a rotation of angle φ. The vector ∂ y 3 r 1 e is collinear to r 3 e , and therefore (∂ y 3 r 1 e ) • r 2 e = 0. Obviously, b e does not depend on y 3 ≡ φ. In view of (2.9), we have the implicit relation:

(2.12) R 4 (y 1 ) 2 + R -y 2 = 1 , R := (ρ 2 + z 2 ) 1/2 .
With (2.7) and (2.12), we can compute:

(2.13) b 2 e = 1 + 3 (y 1 ) 2 R 4 R 6 , ∂ y 1 R = -2 y 1 R 5 1 + 3 (y 1 ) 2 R 4 , ∂ y 2 R = R 2 2 -y 2 1 + 3 (y 1 ) 2 R 4 .
From (2.13), it is easy to verify that ∂ y 2 b e (y) < 0, and that ∂ y 1 b e (y) > 0 when y 1 > 0.

The information about ∂ 2 y 1 y 1 b e (0, y 2 ) comes from ∂ 2 y 1 y 1 (b e ) 2 (0, y 2 ) = 18 R -2 . The symmetry with respect to Eq guarantees that b e (•, y 2 ) and d j e (•, y 2 ) are even in y 1 ∈ R.

The requirements on B made in Assumption 2.1 are fulfilled by other magnetic fields, such as Jupiter's (Voyager spacecraft detected whistler-like activity in the vicinity of Jupiter). Even in the case of the earth, they are compatible with more accurate models, which for instance incorporate solar effects. However, for the sake of simplicity, these aspects will not be taken into account.

2.1.2. Straightening out the field lines. Change (x, v) into (y, w) := y(x), O•y(x) v . Write the system (2.1) in terms of these new variables (y, w) to get:

(2.14)

∂ t y ε = D(y ε ) w ε , y ε (0, y, w) = y , ∂ t w ε = ε -1 B(ε, y ε , w ε ) × w ε + ε N (y ε , w ε , E, B) , w ε (0, y, w) = w .
By construction, we have B(ε,

•) = B 0 (•) + ε B 1 (•) with: (2.15) B 0 =   B 1 0 B 2 0 B 3 0   =   c(|w|) b 0 0   , B 1 =   B 1 1 B 2 1 B 3 1   =   [(D w • ∇ y ) t r 3 ] • r 2 [(D w • ∇ y ) t r 1 ] • r 3 [(D w • ∇ y ) t r 2 ] • r 1   .
The nonlinear perturbation N can be deduced from Ñ as follows:

(2.16) N (y, w, E, B) := O(y) Ñ t O(y) w, E • y -1 (y), B • y -1 (y) .

The dynamical system (2.14) is of the same type as (2.1) except that the singular part of the penalization matrix, which is represented above by B 0 , points now in the constant direction e 1 = t (1, 0, 0). Thus, the effect and the interest of the step 2.1.2 is to straighten out the direction of the vector field B(•). This is done only at leading order. Since the diffeomorphism y : Õ -→ O is an orthogonal change of variables, this operation creates at the level of the equation on w ε a skew-symmetric O(1) contribution that is materialized by the cross product B 1 ×. Observe that the action of B 1 is absent if and only if r 1 (and therefore r 2 and r 3 ) can be taken as a constant function.

2.1.3. Filtering out the main oscillations. The filtering method has been first introduced by S. Schochet [START_REF] Schochet | Fast singular limits of hyperbolic PDEs[END_REF], and then by E. Grenier [START_REF] Grenier | Pseudo-differential energy estimates of singular perturbations[END_REF]. It is today well-known. It has been applied in the article [START_REF] Gallagher | On pressureless gases driven by a strong inhomogeneous magnetic field[END_REF] to deal with a two-dimensional version of the present framework. Below, it is implemented to start the analysis. Introduce:

Λ := e 1 × =   0 0 0 0 0 -1 0 1 0   , e θΛ =   1 0 0 0 cos θ -sin θ 0 sin θ cos θ   , θ ∈ T := R/(2πZ).
We find that t Λ = -Λ. The matrix e θΛ is orthogonal. Consider a smooth (real-valued) scalar function φ ε (t, y, w), and replace w ε by u ε := e -φε Λ/ε w ε , so that

|w ε | = |u ε |. Define z ε := t (y ε , u ε ) ∈ R 3 × R 3
, which is a function of the variables (t, y, w). When expressing the preceding quantities in terms of z, oscillations are introduced. For instance:

(2.17) and: (2.18) B 0 (z) := c(|u|) b(y) , B 1 (θ, z) := B 1 1 (y, e θ Λ u) . Definition 2.1 (eikonal equation). The phase φ ε is adjusted in such a way that:

B 1 (ε, y ε , w ε ) = c(|w ε |) b(y ε ) + ε B 1 1 (y ε , w ε ) = B ε, ε -1 φ ε (t, y, w), z ε , with B(ε, θ, z) := B 0 (z) + ε B 1 (θ, z),
(2.19) ∂ t φ ε = B ε, ε -1 φ ε (t, y, w), z ε , φ ε (0) = 0 .
The right-hand side of (2. [START_REF] Joly | Nonlinear oscillations beyond caustics[END_REF]) is based on B 0 . Note that, for subsequent important reasons, the apparently small term ε B 1 is also incorporated inside B. The equation (2.19) is not self-contained. It must be coupled with a condition allowing to determine z ε .

Lemma 2.2 (the dynamical system after straightening and filtering). The time evolution of t (z ε , φ ε ) ∈ R 7 is governed by a differential equation which can be formalized according to:

(2.20) ∂ t z ε φ ε = A B ε, φ ε ε , z ε , z ε φ ε (0) = t (y, w) 0 ,
where the two applications

B(•, θ, •) and A(•, θ, •) are smooth in (ε, z) ∈ [0, 1] × R 6
, whereas the two functions B(ε, •, z) and A(ε, •, z) are smooth, and periodic in θ with period 2 π.

Proof. It suffices to compute B and A, and to check the listed properties. Just refer to the line (2.18) concerning B. On the other hand, we have:

∂ t y ε = D(y ε ) w ε = D(y ε ) e φε Λ/ε u ε , ∂ t u ε = ε -1 ∂ t φ ε + B 1 Λ u ε + e -φε Λ/ε   0 B 2 1 B 3 1   × w ε + ε e -φε Λ/ε N .
In view of (2.17), the eikonal equation is adjusted to make disappear the term with ε -1 in factor. Finally, we find that

A(ε, •) = A 0 + ε A 1 with:   A 1 0 A 2 0 A 3 0   (θ, z) = D(y) e θΛ u =   d 1 (y) u 1 d 2 (y) (cos θ u 2 -sin θ u 3 ) d 3 (y) (sin θ u 2 + cos θ u 3 )   , (2.21a)   A 4 0 A 5 0 A 6 0   (θ, z) = e -θΛ   0 B 2 1 B 3 1   × e θΛ u (2.21b) =   (B 2 1 sin θ -B 3 1 cos θ) u 2 + (B 2 1 cos θ + B 3 1 sin θ) u 3 -(B 2 1 sin θ -B 3 1 cos θ) u 1 -(B 2 1 cos θ + B 3 1 sin θ) u 1   ,   A 1 1 A 2 1 A 3 1   (θ, z) =   0 0 0   ,   A 4 1 A 5 1 A 6 1   (θ, z) = e -θΛ N (y, e θΛ u, E, B) . (2.21c)
In (2.21b), the functions B 2 1 and B 3 1 must be evaluated at the point (y, e θΛ u).

We are faced with (2.20). In comparison with (2.14), the source terms inside (2.20) are no more of size ε -1 . The applications B(•, θ, z) and A(•, θ, z) are smooth in ε ∈ [0, 1] and, due to the periodicity in θ, the source terms of (2.20) involve many oscillations at frequencies of the order ε -1 . Observe that the contribution issued from the part t (A 4 0 , A 5 0 , A 6 0 ) is nonlinear in z. Still, some skew-symmetric property has been preserved in so far as:

t (A 4 0 , A 5 0 , A 6 0 )(θ, z) = ∆(θ, z) u , t ∆(θ, z) = -∆(θ, z) , ∀ (θ, z) ∈ T × R 6 .
Define Ā0 as being the mean value of A 0 with respect to θ. Note that Ā0 is composed of the diagonal part t d 1 (y)u 1 , 0, 0 and of the skew-symmetric action ∆(z) u.

2.2.

The approach to reach long times t ∼ ε -1 . Change the variable t into τ := ε t. Then, the system (2.20) is simply translated into:

(2.22) ∂ τ z ε φ ε = 1 ε A B ε, φ ε ε , z ε , z ε φ ε (0) = t (y, w) 0 .
The long time evolution (t ∼ ε -1 or τ ∼ 1) is a delicate question for many reasons. There is the presence of singular terms, like ε -1 (A 1 0 , A 2 0 , A 3 0 ) and ε -1 B 0 , which at first sight inherit no special structure. And, even concerning the skew-symmetric part ε -1 ∆(θ, z) u, the filtering method cannot be repeated along the same lines as before. Objections come from the presence already of oscillations (∂ θ ∆ ≡ 0) and from the nonlinearity of the penalization matrix (∇ z ∆ ≡ 0). To go further, other arguments must be implemented.

Preliminaries and guidelines.

Looking at (2.22) gives the impression that the fast dynamics in long times are controlled by the unworkable term (A 0 , B 0 ). But it would be without paying attention to the presence of two scales, playing different roles. In fact, the oscillations inside (A 0 , B 0 ) must be handled separately. As we will see, what principally matters is the influence of Ā0 . ⋄ A first idea when dealing with (2.22) is to adequately separate the different scales. To this end, we have to introduce notations. Given T ∈ R * + , note T T := R/(T Z) and, as usual, note simply

T := T 2 π ≡ R/(2 π Z). Let Z(τ, t, θ) ∈ C 0 (R × T T × T). Introduce: (2.23) Z(τ, t) := 1 2 π 2 π 0 Z(τ, t, θ) dθ , Z (τ ) := 1 2 π T T 0 2 π 0 Z(τ, t, θ) dt dθ .
We can further decompose the Fourier series of Z according to:

(2.24) Z(τ, t, θ) = Z(τ,t) c (0,0) (τ ) Z (τ ) + p∈Z * c (p,0) (τ ) e ipt/T Z⋆ (τ,t) + (p,q)∈Z×Z * c (p,q) (τ ) e i(pt/T +qθ) Z * (τ,t,θ)
.

To better set the scene, we comment the above conventions. The symbols • and • mean that we extract the mean value of Z(τ, t, θ) in θ, and of Z(τ, t) in t. The marks • * and • ⋆ at the top right of Z indicate that we deal with the oscillating part of Z(τ, t, θ) in θ, and of Z(τ, t) in t. The actions ∂ -1 θ and ∂ -1

t will correspond to the inverse operators of ∂ θ and ∂ t , that is ∂ -1 θ ∂ θ = Id and ∂ -1 t ∂ t = Id,
with values in the set of functions with zero mean:

∂ -1 t Z⋆ (τ, t) := t 0 Z(τ, s) ds - 1 T T 0 t 0 Z(τ, s) ds dt , ∂ -1 t Z⋆ = 0 . Lemma 2.3. Let Z ∈ C 0 (T T ).
Retain that:

(2.25) (2.22). This requires a rectification procedure whose principle is exposed in Paragraph 2.2.2, and whose effective construction is achieved from Paragraph 2.2.3 to 2.2.6. This rectification is based on two operations. First, a conjugation of the system (2.22) with a flow induced by Ā0 . This pertains to a nonlinear conjugation since Ā0 is a nonlinear function, see (2.39) and (2.41). Secondly, while bearing in mind the presence of oscillations, an absorption of the oscillating terms through a sort of elliptic equation. ⋄ ⋄ A third idea is to exploit the properties of the geomagnetic field to show that the flow generated by Ā0 is globally defined in time, and in fact periodic in time (Proposition 2.1). It is at this level that the variations of the vector field r1 , that is of the unit directions generated by B, come into play. Depending on the situation, these variations can induce destabilizing effects or not. The geometrical criterion giving rise to a confinement is exposed in Assumption 2.2. It is satisfied in the case of magnetic dipoles, see Lemma 2.6. ⋄

t 0 Z(s) ds = Z t + t 0 Z ⋆ (s) ds , sup t∈R | t 0 Z ⋆ (s) ds| < +∞ . ⋄ ⋄ A second idea is to remove A 0 from the system
The main purpose of this Section 2.2 is to account for the observations which have been reported in Steps (P2) and (P3). We also wish to emphasize the two new Mathematical interpretations 1 and 2 (given in Paragraph 2.2.4) of physical phenomena: the first explains the bouncing back through some geometrical argument; the second highlights the notion of extended Van Allen belt, from which the Van Allen belt is derived.

2.2.2. Approximate rectification of the vector field A, homological equation. What is called here "approximate rectification", it is a (time dependent) oscillatory transformation of z ε (given by a lift) allowing to get rid of A 0 . Concretely, consider two maps:

Ξ : R + × R 6 -→ R 6 (t, z) -→ Ξ(t, z) , Ξ : R + × T × R 6 -→ R 6 (t, θ, z) -→ Ξ(t, θ, z) .
Then, look at the system:

(2.26) ∂ τ z ε ϕ ε (τ ) = 1 ε A b ε, τ ε , ϕ ε ε , z ε , z ε ϕ ε (0) = z iε 0 ,
where the initial data z iε is adjusted such that Ξ(0, z iε ) + ε Ξ(0, 0, z iε ) = (y, w), and where, assuming that the matrix D z Ξ(t, z) is invertible for all (t, z), the applications A and b are given (for ε small enough) by the chain rules:

A(ε, t, θ, z) := D z Ξ(t, z) + ε D z Ξ(t, θ, z) -1 A ε, θ, Ξ(t, z) + ε Ξ(t, θ, z) -∂ t Ξ(t, z) -B ε, θ, Ξ(t, z) + ε Ξ(t, θ, z) ∂ θ Ξ(t, θ, z) -ε ∂ t Ξ(t, θ, z) , (2.27) b(ε, t, θ, z) := B ε, θ, Ξ(t, z) + ε Ξ(t, θ, z) = B 0 Ξ(t, z) + O(ε) . (2.28)
There is a link between (2.22) and (2.26). Indeed, any solution (z ε , ϕ ε ) to (2.26) gives rise to a solution (z ε , φ ε ) to (2.22), just by posing:

(2.29) z ε (τ ) := Ξ ε -1 τ, z ε (τ ) + ε Ξ ε -1 τ, ε -1 ϕ ε (τ ), z ε (τ ) , φ ε (τ ) := ϕ ε (τ ) .
Since the functions A and b are smooth with respect to ε ∈ [0, 1], they can be expanded in powers of ε, say A(ε,

•) = A 0 + ε A 1 + ε 2 A 2 + • • • and b(ε, •) = b 0 + ε b 1 + ε 2 b 2 + • • • with: (2.30) A 0 (t, θ, z) := D z Ξ(t, z) -1 A 0 θ, Ξ(t, z) -∂ t Ξ(t, z) -B 0 Ξ(t, z) ∂ θ Ξ(t, θ, z) .
In view of (2.26), there is little possibility to get on z ε a uniform control in ε ∈ ]0, 1], which could be valid for long times τ ∼ 1, except if A 0 ≡ 0. This last condition is similar to some homological equation. Looking at the mean part of (2.30), it consists in imposing the following autonomous differential equation, called the first modulation equation:

(2.31) ∂ t Z0 -Ā0 ( Z0 ) = 0 , Z0 (t) ≡ Ξ(t, z) .
The positivity of B 0 = c b is guaranteed by (2.2) and the definition of c(•). This allows to divide by B 0 . Then, looking at the oscillating part of (2.30), we can give a sense to:

(2.32)

Ξ(t, θ, z) = (∂ -1 θ A * 0 ) θ, Ξ(t, z) / B 0 • Ξ(t, z) .
Assuming that Ξ(t, z) can be identified through (2.31), the relation (2.32) determines Ξ ≡ Ξ * without ambiguity. By the way, we recover the predominance of Ā0 over A * 0 . The problem is about (2.31). It is to show that the equation (2.31) has global bounded solutions.

For clarity and to exhibit important underlying features, it is better to detail the three intermediate changes Ξ 1 (z), Ξ 2 (z) and Ξ 3 (t, z) leading to:

(2.33) Ξ(t, z) := Ξ 1 • Ξ 2 • Ξ 3 (t, z) .
There will be a Paragraph devoted to the construction of each change Ξ j , for j ∈ {1, 2, 3} : Paragraph 2.2.3 for Ξ 1 , Paragraph 2.2.5 for Ξ 2 , and Paragraph 2.2.6 for Ξ 3 .

2.2.3.

The autonomous change of variables Ξ 1 (z). Use spherical coordinates with radial distance u, azimuth angle ς and polar angle a to represent the velocity vector u according to the decomposition

u = t (u , u ⊥ ) with u ≡ u 1 ∈ R, u ⊥ ≡ t (u 2 , u 3 ) ∈ R 2
, and:

(2.34)

u := u cos ς , u ⊥ = u sin ς t (cos a, sin a) , (ς, u, a) ∈ R × R * + × R . Similarly, decompose w = t (w , w ⊥ ) with w ≡ w 1 ∈ R and w ⊥ ≡ t (w 2 , w 3 ) ∈ R 2 into: (2.35) w := w cos ς w , w ⊥ = w sin ς w t (cos a w , sin a w ) , (ς w , w, a w ) ∈ R × R * + × R . Definition 2.2. [definition of Ξ 1 ] Take z = t (y, ς, u, a) and define: Ξ 1 (z) = t (Ξ 1 1 , Ξ 2 1 , Ξ 3 1 , Ξ 4 1 , Ξ 5 1 , Ξ 6
1 )(z) := ( y 1 , y 2 , y 3 , u cos ς , u sin ς cos a , u sin ς sin a ) . Under the action of Ξ 1 , the identity (2.30) reduces to:

(2.36) A 0 (θ, z) = D z Ξ 1 (z) -1 A 0 θ, Ξ 1 (z) , Ā0 (z) = D z Ξ 1 (z) -1 Ā0 Ξ 1 (z) ,
and the equation (2.31) is translated into:

(2.37)

∂ t Z0 -Ā0 ( Z0 ) = 0 , Z0 (0) = Zs1 i0 := (y, ς w , w, a w ) , Z0 = Ξ 1 ( Z0 ) .
We have also:

(2.38)

A(ε, θ, z) = D z Ξ 1 (z) -1 A ε, θ, Ξ 1 (z) = A 0 (θ, z) + ε A 1 (θ, z) ,
and we can obtain:

(2.39)

  A 1 0 A 2 0 A 3 0   (θ, z) = u   d 1 (y) cos ς d 2 (y) cos(θ + a) sin ς d 3 (y) sin(θ + a) sin ς   ,   Ā1 0 Ā2 0 Ā3 0   (z) = u   d 1 (y) cos ς 0 0   .
From above, we deduce that Z2 0 and Z3 0 do not depend on t. This is in coherence with [3] and gyrokinetic theory [START_REF] Brizard | Foundations of nonlinear gyrokinetic theory[END_REF]. We recover here that the main spatial effect (in position x ∈ Õ) during intermediate times (t ∼ 1) is a drift along the field lines.

Lemma 2.4. [conservation of the kinetic energy]

Assume that E ≡ 0. Then, the quantities

|v ε |, |w ε |, |u ε | or u ε are constant for all times τ ∈ R + , equal to |v| = |w| = w.
Proof. The Definition 2.34 of u ε says that u ε = |u ε |. Since the transformations allowing to pass from v ε to w ε , and then from w ε to u ε , are orthogonal, we find that:

|v ε | 2 = |O • y(x ε ) v ε | 2 = |w ε | 2 = |e φε Λ/ε u ε | 2 = |u ε | 2 = u 2 ε .
It suffices to argue at the level of v ε . From (1.5) and (1.6), we get:

∂ τ |v ε | 2 = 2 v ε • Ñ (v ε , 0, B) = 2 (1 -|v ε | 2 ) 1/2 v ε • (v ε × B) = 0 .
The interpretation at the level of (2.37) of Lemma 2.4 is that Z5

0 (t) = w for all t ∈ R * + , or equivalently that A 5 (ε, •) ≡ 0, that is A 5 0 (•) ≡ 0 and A 5 1 (•) ≡ 0. Now, when E ≡ 0, the proof indicates that ∂ t |v ε | 2 = O(ε)
. This means that the conservation of the kinetic energy remains true at the main order ε 0 . As a consequence, we still have A 5 0 ≡ 0, as confirmed by a direct computation. By construction:

(2.40)

A 4 0 A 6 0 (θ, z) = -B 2 1 sin(θ + a) + B 3 1 cos(θ + a) -cotan ς B 2 1 cos(θ + a) + B 3 1 sin(θ + a)
.

In (2.40), there is some hidden dependence in θ because, in view of (2.15), the functions B 2 1 (y, w) and B 3 1 (y, w) must be evaluated at the point w = e θΛ u. Taking this remark into account and averaging the expressions A 4 0 and A 6 0 with respect to θ gives rise to:

(2.41)

Ā4 0 Ā6 0 (z) = u sin ς C (y) cos ς E (y) , C (y) := [(d 2 ∂ y 2 t r 2 + d 3 ∂ y 3 t r 3 ) • r 1 ](y)/2 , E (y) := [(d 2 ∂ y 2 t r 3 -d 3 ∂ y 3 t r 2 ) • r 1 ](y)/2 .
Consider the surfaces

S k c := {x ∈ Õ ; y k (x) = c}, with k ∈ {1, 2, 3}.
Lemma 2.5. [geometrical interpretation of the coefficients C and E ] We have E (y) = 0. Moreover, C (y) coincides with the mean curvature of the surface S 1 y 1 , computed at the point x such that y(x) = y.

Proof. Let γ j be some integral curve of r j . It can be parameterized by y j or by its arc length s, with the identity d j ∂ y j ≡ ∂ s . By Dupin's theorem (the surfaces of a triply orthogonal system intersect each other along lines of curvature), γ j is necessarily a line of curvature for the surfaces S k c with k = j. It follows that the quantities (d j ∂ y j t r j • r 1 )(y) with j ∈ {2, 3} can be interpreted as the principal curvatures of the surface S 1 y 1 , computed at the point x such that y(x) = y. By definition, the average of these two principal curvatures, that is C (y), is the mean curvature of S 1 y 1 at y -1 (y). Another argument prevails concerning E . The quantities d 2 (∂ y 2 t r 3 ) • r 1 and (d 3 ∂ y 3 t r 2 ) • r 1 can be viewed as the two extra-diagonal terms of the Gauss map Γ : S 1 c -→ S 2 . Due to the symmetry of the Weingarten endomorphism, we are sure that E ≡ 0.

From (2.37), it follows that Z0 (t, y, w) satisfies:

(2.42) Z2 0 = y 2 , Z3 0 = y 3 , Z5 0 = w , Z6 0 = a w , ∀ t ∈ R ,
together with a planar motion directed by:

(2.43)

∂ t Z1 0 = w d 1 ( Z1 0 , y 2 , y 3 ) cos Z4 0 , Z1 0 (0) = y 1 , ∂ t Z4 0 = w C ( Z1 0 , y 2 , y 3 ) sin Z4 0 , Z4 0 (0) = ς w .
2.2.4. The reduced hamiltonian. By shrinking the size of O if necessary, we can always assume that the domain O is given in the following form:

(2.44)

O = y = t (y 1 , y 2 , y 3 ) = t (y 1 , ŷ) ; y 1 ∈ ]y 1 -(ŷ), y 1 + (ŷ)[ , ŷ := (y 2 , y 3 ) ∈ Ô .
In (2.44), we assume that:

The set Ô ⊂ R 2 is open and bounded , (2.45a)

The functions

y 1 ± ∈ C ∞ ( Ô; R) are such that y 1 -≤ y 1 + . (2.45b) Assumption 2.2.
[mean curvature with a unique oriented change of sign] The set O can be described as in (2.44)-(2.45), and there is a smooth function y 1 0 : Ô -→ R such that the zero set of C is given by the graph of y 1 0 . In other words:

(2.46) C -1 (0) = y ; C (y) = 0 = y 1 0 (ŷ), ŷ ; ŷ ∈ Ô . Moreover, the sign of C (•) changes from -to + when crossing C -1 (0) : (2.47) ∂ y 1 C y 1 0 (ŷ), ŷ > 0 , ∀ ŷ ∈ Ô .
• Application 2 [change of sign along the equatorial plane]. For k ∈ {1, 2, 3}, introduce the surfaces S k ec := {x ∈ Õe ; y k e (x) = c}. Let C e (y) be the mean curvature of the surface S 1 ec , computed at the point x such that y e (x) = y.

Lemma 2.6. The function C e (•) does not depend on the variable y 3 . It is given by:

(2.48) C e (y 1 , y 2 ) = 3 2 1+3 (R 2 y 1 ) 2 -3/2 3+5 (R 2 y 1 ) 2 R y 1 , R := (ρ 2 +z 2 ) 1/2 ,
with R as in (2.12). Assumption 2.2 is verified with y 1 e0 ≡ 0. For all y 2 ∈ R * -, the partial function C e (•, y 2 ) is odd on R, and strictly increasing on R * + . Proof. Since the set S 1 ec is a surface of revolution with axis D e , the function C e (•) does not depend on y 3 . The identity (2.48) comes from a computation based on (2.9) and on the classical identity C e (y) = -2 -1 (div x r1 e ) • y -1 (y). Let a := R 2 y 1 . Exploit the relation (2.12) to extract:

∂ y 1 C e (y) = 3 2 2 a 2 G ′ (a 2 ) + 1 + a 2 1 -a 2 G(a 2 ) R 2 -y 2 1 + 3 (y 1 ) 2 R 4 , G(λ) := (3 + 5 λ) (1 + 3 λ) 3/2 .
The quantity within the brackets is positive, so that ∂ y 1 C e (y) > 0 with in particular:

(2.49) ∂ y 1 C e (0, y 2 ) = 9 2 (-y 2 ) 1/2 > 0 , ∂ 2 y 1 y 2 C e (0, y 2 ) = 9 4 (-y 2 ) 3/2 > 0 , ∀ y 2 ∈ R * -.
It is worth mentioning how to derive Assumption 2.2 from geometrical considerations. In view of (2.9), the surface S 1 e0 is the plane {z = 0}. This gives C -1 e (0) ⊃ {z = 0}. When c = 0, the surface S 1 ec has everywhere two non-zero principal curvatures of the same sign, and therefore everywhere a non-zero mean curvature. This furnishes (2.46) with y 1 0 ≡ 0. Moreover, the surface S 1 ec is the mirror image of S 1 e-c for the reflection with respect to the plane S 1 e0 , indicating that C e (•, y 2 ) is an odd function. The convention adopted about the orientation implies that y 1 → C e (y 1 , y 2 ) is increasing, explaining the sign in (2.47).

• Definition 2.3. Introduce the applications:

(2.50)

P 1 : O -→ R y -→ y 1 y 1 0 (ŷ) C (s,ŷ) d 1 (s,ŷ) ds , P 2 : ]0, π[ -→ R ς -→ -ln (sin ς) .
The sum H of P 1 and P 2 is called the reduced hamiltonian:

(2.51) H(z) := P 1 (y) + P 2 (ς) = y 1 y 1 0 (ŷ) C (s,ŷ) d 1 (s,ŷ) ds -ln (sin ς) , z = t (y, ς, u, a) .
Retain that H depends in fact only on (y 1 , y 2 , y 3 , ς). In contrast with the full hamiltonian associated to (1.5)

(2.52) ∂ t (H • Z0 )(t, y, w) = 0 , ∀ (t, y, w) ∈ R + × O .
Proof. Just compute:

∂ t (H • Z0 ) = (∂ y 1 P 1 )( Z1 0 , y 2 , y 3 ) ∂ t Z1 0 + (∂ ς P 2 )( Z4 0 ) ∂ t Z4
0 , and use (2.42), (2.43) together with the definitions of P 1 and P 2 to find 0.

The function P 2 : ]0, π[-→ R is obviously strictly convex. It has a global minimum at ς = π/2, with P 2 (π/2) = 0. This is a potential well. Define:

(2.53) g(ς) := sgn ς - π 2 | ln (sin ς)| 1/2 , g 2 ≡ P 2 , ς ∈ ]0, π[ . The function g : ]0, π[ -→ R
∂ y 1 f (y) > 0 , f y 1 0 (ŷ), ŷ = 0 , y = (y 1 , ŷ) ∈ ]y 1 -(ŷ), y 1 + (ŷ)[× Ô . Proof. By definition: P 1 y 1 0 (ŷ), ŷ = 0 , (∂ y 1 P 1 )(y) = d 1 (y) -1 C (y) , d 1 (y) > 0 , ∀ y ∈ O .
As a consequence of (2.46) and (2.47), we find that:

(∂ y 1 P 1 ) y 1 0 (ŷ), ŷ = 0 , (∂ 2 y 1 y 1 P 1 ) y 1 0 (ŷ), ŷ = d 1 y 1 0 (ŷ), ŷ -1 (∂ y 1 C ) y 1 0 (ŷ), ŷ > 0 . Just take f (y) := sgn y 1 -y 1 0 (ŷ) P 1 (y) 1/2 . The function H = f 2 + g 2 has a strict minimum equal to zero, at (y 1 , ς) = (y 1 0 , π 2 )
. The level curves of H are invariant sets for (2.43), diffeomorphic to a family of concentric circles. This indicates that the trajectories are periodic. This will be confirmed further. It should be remembered that the original system (1.5) is nearly integrable [START_REF] Braun | Mathematical remarks on the Van Allen radiation belt: a survey of old and new results[END_REF]. However, what emerged at the level of (2.43), that is at main order during intermediate times, is integrability.

• Application 3 [reason of the bouncing back] To our knowledge, the following explanation has not yet been produced in connection with Step (P2). This is why it is put forward.

Mathematical interpretation 1 (of the bouncing back). During intermediate times, the guiding centers move in the direction y 1 , with ŷ = (y 2 , y 3 ) fixed. They follow the projections on the base (the line y 1 ) of bicharacteristics drawn in the phase space (y 1 , ς), and given by the level curves of the reduced hamiltonian H. When the mean curvature C e (•, y 2 ) of the surfaces which are orthogonal to the field lines changes its sign along the field lines Fl(y 2 , y 3 ), these level curves are Jordan curves, and their projections are segments. As a result, the particles bounce back and forth between two "mirror points", which are the two extremities of these segments.

• From now on, we will work under Assumption 2.2, with energy levels

H(z) below h(ŷ) 2 where h ∈ C ∞ ( Ô; R + ) is such that: (2.55) 0 ≤ h(ŷ) < min lim sup y 1 → y 1 -(ŷ)
|f (y 1 , ŷ)| ; lim sup

y 1 → y 1 + (ŷ) f (y 1 , ŷ) .
In view of (2.43), the amplitude u of u influences only the speeds at which the particles move along trajectories. It does not modify the form of the trajectories. When u = 0 or when H e (z) = 0, the particle is stationary. This situation has no dynamical interest. It can be avoided through the two conditions 0 < u m ≤ u and 0 < H m ≤ H e (z), imposed for some u m ∈ ]0, 1[ and for some H m ∈ R * + . The case u = 1 is also not pertinent, because even the most energetic particles which are implied move below the speed of light. We will impose u ≤ u M < 1 for some u M with u m < u M < 1. All this amounts to focus on a bounded domain of the phase space of the following type:

(2.56) O ≡ O(h) := z ; z = t (y 1 , ŷ, ς, u, a) , ŷ = t (y 2 , y 3 ) ∈ Ô , 0 < H m ≤ H(z) < h(ŷ) 2 , 0 < u m < u < u M < 1 .
Coming back in the variables (y, u), and then in the original variables (x, v), we get:

(2.57)

Ω := Ξ 1 (z) ; z ∈ O , Ω := (x, v) ; y(x), O • y(x) v ∈ Ω .
The set Ω is the domain alluded in Theorem 1. Introduce also:

(2.58)

Ǒ := ž = t (h, ŷ, u) ; 0 ≤ h < h(ŷ) , ŷ ∈ Ô , u m < u < u M .
• Application 4 [drawing the shape of Van Allen belts]. Consider a particle starting at time t = 0 from the position (y, w). Define h := √ H e (y, w). During intermediate times, the trajectory of this particle is completely determined by ž = t (h, ŷ, w). The value of H is constant along the motion, equal to h 2 , but the values of P e1 and P e2 are not. The potential well P e1 (•, y 2 ) is maximal when |y 1 | reaches its maximum value along the trajectory. This occurs when y 1 = ±y 1 m (y 2 ), with y 1 m adjusted such that y 1 m (y 2 ) = f -1 ey 2 (h). This is when Z1 0 (t) = ±y 1 m . Then, we have Z4 0 (t) = π/2. At this time t, there is no drift along the field line Fl(y 2 , y 3 ). The particle is bouncing back. The potential well P e2 (•) is maximal when |ςπ/2| reaches its maximum value along the trajectory. This occurs when ς = ς m (modulo π), with ς m adjusted such that ς m := g -1 (h). This happens when Z1 0 (t) = 0, that is when the particle crosses the equatorial plane Eq. The value of ς m is called the equatorial pitch angle by physicists.

To ensure that the particles will not strike the atmosphere (and therefore be lost), the velocity vectors Ū0 (t) must point outside the loss cone [START_REF] Pandey | Pitch angle loss-cone anisotropic magneto plasma in presence of parallel electric a.c. field[END_REF], meaning that | Ū ⊥ 0 (t)| = 0, and most importantly, that the quotient Ū 0 (t)/| Ū ⊥ 0 (t)| ≡ cotan ς is adequately bounded. This can be achieved through a control on the pitch-angles ς m (away from 0 and π) or through a control on the mirror points y 1 m . In fact, both y 1 m and ς m are completely determined by the value of h m := f ey 2 (y 1 m ) = g(ς m ). This brings to light the notion of extended Van Allen belt EBelt, which is the region in position-velocity defined by: (2.59)

EBelt := (y, ς, u, a) ∈ O e × T × [u m , u M ] × T ; √ H e (y, ς) ≤ f ey 2 • h e (y 2 ) ,
where the function h e (•) was introduced in Paragraph 2.1.1.

Mathematical interpretation 2 (of the Van Allen belts). The Van Allen belt Belt may be seen as the spatial projection of an extended Van Allen belt EBelt. The shape of EBelt can be drawn by prescribing the values of u and adequate bounds on the energy levels of the reduced hamiltonian H e , see (2.59).

• It is EBelt which is pertinent from a dynamical point of view, for instance when dealing with localization properties of the initial data (y, ς w , w, a w ).

2.2.5.

The autonomous change of variables Ξ 2 (z), second piece of (2.33). To avoid multiple notations, we will consider during this second step that the function A obtained before, at the level of line (2.38), is the new starting point, and therefore rename it A. Thus, in this paragraph 2.2.5, we note z = t (y, ς, u, a). For such a point z chosen in Ω, we can write:

(2.60) (y 1 , ς) = f -1 ŷ (h cos υ), g -1 (h sin υ) , h := √ H(y 1 , ŷ, ς) ∈ [0, h(ŷ)] .
Lemma 2.9. [the first modulation equation in action-angle coordinates] The system (2.43) reduces to ∂ t h = 0, together with the following scalar equation on υ:

(2.61)

∂ t υ = w µ ž(υ) , µ ž(υ) := - C (f -1 ŷ (h cos υ),ŷ) 2 h cos υ cos g -1 (h sin υ) h sin υ .
Proof. The definition of υ in (2.60) leads to:

(2.62)

g ′ ( Z4 0 ) ∂ t Z4 0 = w C ( Z1 0 , ŷ) g ′ ( Z4 0 ) sin Z4 0 = h cos υ ∂ t υ .
With the definition of P 2 ≡ g 2 in (2.50), the polar decomposition (2.60) yields also:

(2.63) g ′ ( Z4 0 ) sin Z4 0 = - cos Z4 0 2 g( Z4 0 ) = - cos g -1 (h sin υ) 2 h sin υ .
Combining (2.62) and (2.63), we can easily deduce (2.61).

Of course, the function µ ž(•) is not defined when h = 0, υ = 0 (mod π) or υ = π 2 (mod π), but it can be extended there by continuity. Then, it becomes a smooth function on the whole domain R, which is periodic in υ with period 2 π. Moreover: Lemma 2.10. The speed of rotation µ is bounded and nowhere zero. More precisely:

(2.64) ∃ (µ -, µ + ) ∈ (R * + ) 2 ; µ -≤ µ ž(υ) ≤ µ + , ∀ (ž, υ) ∈ Ǒ × R .
Proof. By construction, C f -1 ŷ (•), ŷ andcos g -1 (•) are two strictly increasing functions, taking the value 0 at the point 0, with non zero derivatives at all points. It follows that µ ž(υ) ∈ R * + for all (ž, υ) ∈ Ǒ × R. The uniform control in (2.64) comes from an argument of compactness. Indeed, ž ∈ Ǒ and Ǒ is relatively compact in view of (2.45a) and (2.58). On the other hand, υ ∈ T and the torus T is compact.

The lower bound inside (2.64) is very important in what follows. In particular, it allows to consider the function 1/µ ž(υ). This is again a bounded, strictly positive function, which is periodic in υ with period 2 π. For all ž ∈ Ǒ, the application:

(2.65)

K ž : R -→ R , K ž(υ) := υ 0 dr µ ž(r) ,
is a diffeomorphism. With Lemma 2.3 in mind, introduce the decomposition:

(2.66)

K ž(υ) = k ž υ + K p ž (υ) , k ž := 1 2 π 2 π 0 dr µ ž(r) ∈ [µ -1 + , µ -1 -] ,
where

K p ž (•) is a 2 π-periodic function. The inverse of K ž, in the sense that K ž • K -1 ž (ν) = ν, can be put in the form: (2.67) K -1 ž (ν) = k -1 ž ν + Kp ž (ν) , Kp ž (ν + 2 π k ž) = Kp ž (ν) , ∀ ν ∈ R .
Using K ž, it is possible to normalize the speed of rotation. To this end, it suffices to replace the variable υ by ν := K ž(υ). The motion on the invariant tori expressed in terms of the canonical coordinates (h, ν) becomes linear in the angle variable ν:

(2.68)

∂ t h = 0 , ∂ t ν = w . Definition 2.4.
[definition of Ξ 2 ] Take z = t (h, y 2 , y 3 , u, a, ν) and define:

Ξ 1 2 (z) := f -1 ŷ h cos K -1 ž (ν) , Ξ 2 2 (z) := y 2 , Ξ 3 2 (z) := y 3 , Ξ 4 2 (z) := g -1 u h sin K -1 ž (ν) , Ξ 5 2 (z) := u , Ξ 6 2 (z) := a . Transforming (2.31) through Ξ 1 • Ξ 2 yields: (2.69) ∂ t Z0 -Ā0 ( Z0 ) = 0 , Z0 (0) = Zs2 i0 , Z0 = Ξ 1 •Ξ 2 ( Z0 ) , Ā0 = t (0, • • • , 0, u) .

The initial data Zs2

i0 can be determined through the relation

Ξ 1 • Ξ 2 ( Zs2 i0 ) = z = t (y, w). With h z := √ H(y, ς w
) and ž = t (h z , ŷ, w), this means that:

(2.70)

Zs2 i0 = (h z , ŷ, w, a w , ν z ) , ν z := K ž • arccos h -1 z f ŷ(y 1 ) .
At this stage, the vector field Ā0 is rectified. The components Zj 0 with j ∈ {1, 2, 3, 4, 5} are fixed parameters, whereas Z6 0 (t) is simply linear: Z6 0 (t) = Z6 i0 + w t.

2.2.6. The time dependent change of variables Ξ 3 (t, z), third piece of (2.33). This is just a modification of the last component.

Definition 2.5. [definition of Ξ 3 ] Take z = t (h, y 2 , y 3 , u, a, b) and define:

Ξ j 3 (t, z) := z j , ∀ j ∈ {1, • • • , 5} , Ξ 6 3 (t, z) := b + w t .
This operation changes Ā0 into Ā0 ≡ 0. But it makes the variable t appear at the level of the function A defined by (2.27), with Ξ = Ξ 1 • Ξ 2 • Ξ 3 and Ξ ≡ 0. Remarkably, since Ξ 2 is periodic in ν with period 2 π k z , the linear growth of Ξ 3 in the variable t is converted to a periodic behaviour in t of Ξ(•, z). (2.71)

T ≡ T ( Zs1 i0 ) ≡ T (ž) := 2 π k ž w -1 .
The solution Z0 to (2.37), or equivalently to (2.43), remains for all times t ∈ R + in the energy level

H( Zs1 i0 ). With h z = √ H( Zs1 i0
), it can be described by: (2.72)

Z0 (t) = f -1 ŷ h z cos K -1 ž (b z +w t) , y 2 , y 3 , g -1 h z sin K -1 h (b z +w t) , w, a w .
It is periodic in time t, with a period T depending on ž as in (2.71). The domain Ω, see (2.57), is invariant under the flow induced by (2.43). Given F ∈ C 0 (O; R), we have:

(2.73) F • Z0 (•, Zs1 i0 ) = π -π F f -1 ŷ (h z cos υ), y 2 , y 3 , g -1 (h z sin υ), w, a P ž(υ) dυ , where P ž(υ) := (2 π k ž) -1 µ ž(υ) -1 is a probability density function on the interval [-π, π].
Succinct proof of Proposition 2.1. A preliminary remark is that all components of the vector ž = t (h z , ŷ, w) are constants of the motion. The explicit formula (2.72) is an outcome of the preceding analysis. Concerning the periodicity, use (2.67) to obtain:

K -1 ž b z + w (t + T ) = k -1 ž (b z + w t + 2 π k ž) + Kp ž (b z + w t + 2 π k ž) = K -1 ž (b z + w t) + 2 π .
To get (2.73), just perform the change of variables υ = K -1 ž (b + w t), and use the definition of k ž in (2.66) to see that P ž is indeed a probability.

The identity (2.73) says that the time integral of the periodic function F • Z0 (t, Zs1 i0 ) over its period T can be converted into a spatial integral. The resulting expression F • Z0 (•, Zs1 i0 ) depends only on ž. This allows to define a map:

(2.74)

A : C 0 (O; R) -→ C 0 ( Ǒ; R) F (z) -→ A (F )(ž) := F • Z0 (•, z) .
The integrand in (2.73) is the pull-back of F in the action-angle variables (h, υ). When this integrand does not depend on υ, we simply recover A (F ) ≡ F . Otherwise, the computation of A (F ) can be more complicated. For instance, if F is a function of the variable y 1 only, the action of A on the observable F produces a function A (F )(ž) which can depend on all the coordinates in ž. Below, we remark that, under symmetry conditions, simplifications can occur when computing A (F )(ž).

Lemma 2.11. Let F (y 1 , ŷ, ς, u, a) ∈ C 0 (O; R). Assume that the functions C (•, ŷ) and d 1 (•, ŷ) are respectively odd and even. Then, the formula (2.73) can be simplified into:

(2.75)

A (F )(ž) = - 1 π k ž h 0 r E(F ) f -1 ŷ ( h 2 -r 2 ), ŷ, g -1 (r), w, a w C f -1 ŷ ( h 2 -r 2 ), ŷ cos g -1 (r) dr ,
where E(F ) is the even part of F (•, ŷ, •, w, a w ) centered at the point (0, π/2), that is:

(2.76) E(F )(y 1 , ς) := F (y 1 , ς) + F (-y 1 , ς) + F (y 1 , π -ς) + F (-y 1 , π -ς) .
Proof. In (2.73), decompose the domain of integration into:

[-π, π] = [-π, -π/2] ∪ [-π/2, 0] ∪ [0, π/2] ∪ [π/2, π] .
On each of these subdomains, change υ into r = h z sin υ, with adequate formulas to express h z cos υ in terms of r. This yields four integral contributions. Under the hypothesis of Lemma 2.11, the functions f (•, ŷ) and C f -1 ŷ (•), ŷ are odd. This allows to recombine the four integrals into a single integral, with E(F ) as in (2.76).

• Application 5 [description of Step (P2)]. From Lemma 2.1, we can easily deduce that the potential P e1 does not depend on y 3 , so that H e (z) = H e (y 1 , y 2 , ς). For similar reasons, in the earth's context, the formulas (2.70) and (2.71) become simply:

T e (z) = T e (y 1 , y 2 , ς w , w) = 2 π k (h,y 2 ) w -1 , h = √ H e (y 1 , y 2 , ς w ) , (2.77a) ν e (z) = ν e (y 1 , y 2 , ς w ) = K (h,y 2 ) • arccos h -1 f y 2 (y 1 ) . (2.77b) Moreover, the function y 1 → P e1 (y 1 , y 2 ) is even in y 1 ∈ R. Fix (y 2 , u) ∈ R * + × [u m , u M ]. Select h ∈ R *
+ , and consider the energy curve C e (h) := (y 1 , ς); H e (y 1 , y 2 , ς) = h 2 . As already observed, the minimal and maximal values of y 1 along C e (h) correspond to mirror points. They are obtained for ς = π/2, when g(ς) is minimal, equal to 0. They are opposite, located at the positions ±y 1 m (y 2 ) with y 1 m (y 2 ) determined by the two conditions:

y 1 m (y 2 ) ∈ R * + , H e y 1 m (y 2 ), y 2 , π/2 = h 2 .
The phase portrait associated to the dynamical system (2.43), drawn in the plane {(y 1 , ς)} contains C e (h) as a trajectory. Starting at time t = 0 from a position (y

(T /2) = -y 1 m .
It is strictly increasing on the interval [T /2, T ] with Z1 0 (3T /4) = 0 and Z1 0 (T ) = y 1 m , and so on.

Lemma 2.12 gives a mathematical justification to the step (P2), from which reference was made in the introduction. Indeed, the component Z1 (2.78)

0 (t) = Ξ 1 1 • Z 0 (t) = Z1 0 (t)
z ε ϕ ε (τ ) = z a ε ϕ a ε (τ ) + ε N r z ε r ϕ ε (τ ), z a ε ϕ a ε (τ ) = Z ε Φ ε τ, τ ε , ϕ ε (τ ) ε ,
where r z ε and ε r ϕ ε are intended to be remainders (satisfying uniform estimates in the sup norm with respect to ε ∈ [0, 1]), whereas the expressions Z ε and Φ ε are given by asymptotic expansions of the following type:

(2.79) Z ε Φ ε (τ, t, θ) = ε -1 0 Φp (τ ) + N -1 j=0 ε j Zj Φj (τ, t) + N j=2 ε j Z * j Φ * j (τ, t, θ) ,
involving profiles Z j and Φ j depending smoothly on (τ, t, θ) ∈ R + × T T × T. Taking into account (2.26) to replace ∂ τ ϕ ε , the derivative ∂ τ applied with the substitution performed in (2.78) leads to the differential action:

Op(z; ∂) ≡ Op(ε, t, θ, z; ∂ τ , ∂ t , ∂ θ ) := ∂ τ + ε -1 ∂ t + ε -2 b ε, t, θ, z) ∂ θ .
Plug the ansatz (2.78) into (2.26) to infer the constraint:

(2.80)

Op(Z ε ; ∂) Z ε Φ ε (τ, t, θ) - 1 ε A b (ε, t, θ, Z ε ) = ε N R z R ϕ (ε, τ, t, θ) = O ε N .
Note that ϕ ε is not at all present at the level of (2.80). Thus, when discussing (2.80), the knowledge of ϕ ε is not needed to move forward.

2.3.2.

First and second modulation equations. The first equation of (2.80), about Z ε , is decoupled from the second equation, about Φ ε . Let us start with the formal analysis concerning Z ε . To this end, substitute the expansion (2.79) into (2.80) and collect the terms with ε -1 in factor. By averaging in θ, we obtain the first modulation equation:

(2.81) ∂ t Z0 -Ā0 (t, Z0 ) = 0 , Ā0 ≡ 0 .
We recognize here the rectified version of (2.31), implying that Z0 (τ, t) = Z0 (τ ). On the other hand, the examination of the oscillating part in θ leads to:

Z * 1 = b 0 (t, Z0 ) -1 (∂ -1 θ A * 0 )(t, θ, Z0 ) , b 0 (t, z) := c(u) b Ξ 1 (t, z), y 2 , y 3 .
Recall that A * 0 ≡ 0 as a consequence of (2.32). Therefore, we find Z * 1 ≡ 0, as it has been already anticipated in (2.79). Then, examine what appears with ε 0 in factor to find:

(2.82) ∂ τ Z0 -A 1 (t, θ, Z0 ) + b 0 (t, Z0 ) ∂ θ Z 2 + ∂ t Z1 = 0 .
To compute A 1 (t, θ, z), come back to (2.27) in order to extract:

(2.83)

A 1 (t, θ, z) = D z Ξ(t, z) -1 A 1 θ, Ξ(t, z) + Ξ(t, θ, z) • ∇ z A 0 θ, Ξ(t, z) -B 1 (t, θ, z) ∂ θ Ξ(t, θ, z) -∂ t Ξ(t, θ, z) -Ξ(t, θ, z) • ∇ z b 0 (t, z) ∂ θ Ξ(t, θ, z) ,
where Ξ can be replaced as indicated in (2.32). Take the average of (2.82) in both t and θ to exhibit the second modulation equation: w) . The non linear differential equation (2.84) allows to determine Z 0 (τ, t, θ) ≡ Z0 (τ ) on some interval [0, T ] with T ∈ R * + . To understand the long time dynamic (t ∼ ε -1 or τ ∼ 1), as it appears at leading order in the original variables z = t (y, u), the procedure is to pass through (2.84), and then to exploit (2.29), giving rise to:

(2.84) ∂ τ Z0 -Ā1 ( Z0 ) = 0 , Z0 (0) = z i0 ≡ Zs2 i0 , Ξ(0, z i0 ) = t (y,
(2.85) Z0 (τ, t) = Ξ t, Z0 (τ ) , Z0 (τ ) = Ξ •, Z0 (τ ) .
Be careful, nothing guarantees that the map z -→ Ξ(•, z) is a diffeomorphism. Thus, it may be that a self-contained evolution equation on Z0 is not available. Also, there is no evident link between Z0 (0) = Ξ(•, z i0 ) and (y, w) = Ξ(0, z i0 ).

• Application 6 [description of step (P3)]. Come back to the Vlasov part of (1.1). The influence of (E, B) on the second modulation equation is displayed through the term A 1 inside A 1 , see (2.21c) and (2.83). Even if (E, B) ≡ (0, 0), the source Ā1 is likely to bring a non trivial contribution, by a mechanism issued from the fast dynamics.

Lemma 2.13. [self-induced longitudinal drift in the case of the geomagnetic field] Consider the equation (2.84) with data issued from the earth's context. Assume that (E, B) ≡ (0, 0). Then, the corresponding solutions Z0 (τ, z i0 ) are defined for all times τ ∈ R + . More precisely, they are given by a transport at constant speed: ) can be expressed in terms of (y, w) according to: (2.87)

(2.86) Z0 (τ, z i0 ) = z 1 i0 , z 2 i0 , z 3 i0 + Ā3 1 (z i0 ) τ , z 4 i0 , z 5 i0 + Ā5 1 (z i0 ) τ , z 6
z 1 i0 = √ H e (y 1 , y 2 , ς w ) , z 2 i0 = y 2 , z 4 i0 = w .
Besides, the long time mean motion, observed in the original variables t (y, u), satisfies :

(2.88) ( Z1 0 , Z2 0 , Z3 0 , Z4 0 )(τ, y, w) = t ( 0 , y 2 , y 3 + Ā3 1 (z i0 ) τ , 0 . Proof.
To determine Z0 (τ, z i0 ), we have to solve (2.84), that requires to compute Ā1 . To this end, note that the formula (2.83) applies in any of the stages of the transformations. For instance, we can use it after the first transformation Ξ 1 . Remark that (E, B) ≡ (0, 0) means that A 1 ≡ 0. With (2.32), these overviews provide:

(2.89)

A 1 = b -1 0 (D z Ξ 2 ) -1 (∂ -1 θ A * 0 • ∇ z )A 0 -b 1 A * 0 -b 0 ∂ t (b -1 0 ∂ -1 θ A * 0 ) -(∂ -1 θ A * 0 • ∇ z )(ln b 0 ) A * 0 θ, Ξ 2 • Ξ 3 (t, z) ,
where the quantities in the right hand side are the expressions which have been obtained just before Lemma 2.5. Thus, we work here with b j := B j • Ξ 1 for j ∈ {0, 1}, and with A * 0 given by (2.39) and (2.40). In the earth's context, there are many simplifications coming from (2.11) and from the special forms of B 0 and Ξ. The term b 1 is simply:

b 1 (t, θ, z) = Ξ(t, θ, z) • (∇ z B 0 ) Ξ(t, z) + B 1 θ, Ξ(t, z) = u sin ς sin (θ + a) d 2 ∂ y 2 (ln b) + d 3 (∂ y 3 t r 3 ) • r 2 .
On the other hand, (2.40) reduces to:

A 4 * 0 A 6 * 0 = u d 1 (∂ y 1 t r 2 ) • r 1 cos ς cos(θ + a) + D sin ς cos 2 (θ + a) * -d 1 (∂ y 1 t r 2 ) • r 1 cos 2 ς (sin ς) -1 sin(θ + a) -D cos ς cos(θ + a) sin(θ + a)
with

D := d 2 (∂ y 2 t r 2 ) • r 1 -d 3 (∂ y 3 t r 3 ) • r 1
. This provides useful information for computing the contribution A 1 , as indicated in line (2.89). Many simplifications occur when averaging A 1 with respect to θ. For instance, since b 0 does not depend on θ, we find that:

(D z Ξ 2 ) -1 ∂ t (b -1 0 ∂ -1 θ A * 0 ) ≡ ∂ -1 θ (D z Ξ 2 ) -1 ∂ t (b -1 0 A * 0 ) ≡ 0 .
More generally, all contributions which are linear or cubic in cos(θ+a) and sin(θ+a), as well as the products cos(θ + a) sin(θ + a), are canceled. There remains Ā1 = t (0, 0, Ā3 1 , 0, Ā5 1 , 0) with (there is a similar formula concerning Ā5 1 ):

Ā3 1 (t, z) = u 2 1 2 d 2 d 3 d 1 ∂ y 2 ln d 3 d 1 -(d 3 ) 2 d 1 (∂ y 3 t r 3 ) • r 2 (y 1 , y 2 ) sin 2 ς + d 3 (∂ y 1 t r 2 ) • r 1 (y 1 , y 2 ) cos 2 ς | (y 1 , ς, u) = (Ξ 1 2 , Ξ 4 2 , Ξ 5
2 ) • Ξ3(t, z) . The expressions Āj 1 with j ∈ {3, 5} do not depend on (z 3 , z 5 ) ≡ (y 3 , a). We can exploit (2.73) to compute the average of Ā1 with respect to t. This allows to remove the dependence in z 6 ≡ b. The coefficients Ā3

1 and Ā5 1 depend only on (z 1 i0 , z 2 i0 , z 4 i0 ). Then, consider the equation (2.85). Since Ξ j (z) = z j i0 = y j for j ∈ {2, 3}, the result (2.88) about Z2 0 and Z3 0 is a direct consequence of the formulas in (2.85). To get the information (2.88) concerning Z1 0 and Z4 0 , it suffices to apply (2.75) with respectively F (y 1 , ς) = y 1 and F (y 1 , ς) = cos ς.

It is easy to show that Ā3

1 (t, z) is a sum of positive terms as long as |z|/ρ ≤ arctan √ 2. In practice, the expression Ā3 1 (z) is positive in a much larger region. The long time mean motion is a signed rotation around the earth, localized in the equatorial plane Eq, with a sign depending on the charge of the selected particles (electrons or protons). This must be correlated with the creation of a ring current. The contents of Lemma 2.13 are consistent with the observations reported in step (P3) of the introduction. The long time conservation of the "energy levels" h and of the "heights" y 2 , represented by Z1 0 (τ, z i0 ) and Z2 0 (τ, z i0 ) at the level of (2.86), further supports the stability of the extended Van Allen belt EBelt (which is delimited by values of h and y 2 ) and, consequently, of the Van Allen belt. (2.90)

Z⋆ 1 (τ, t) = (∂ -1 t Ā⋆ 1 ) t, Z0 (τ ) , Z * 2 (τ, t, θ) = b 0 (t, Z0 ) -1 ∂ -1 θ A * 1 (t, θ, Z0 ) .
We proceed to see whether it is possible to find the Z j for j ≥ 2. We argue by induction, with the inductive hypothesis:

(H j ) :

The expressions Z k are known on the domain

[0, T ] × T T × T for 1 ≤ k ≤ j -2 ,
The expressions Z⋆ j-1 and Z * j are known on their domains of definition. We have already checked the validity of (H 2 ). Assume (H j ) for all j ∈ {2, • • • , N }. Then, consider what comes with ε j-1 in factor. This is a linearized version of (2.82) :

(2.91)

∂ τ Zj-1 -( Zj-1 • ∇ z )A 1 + b 0 ∂ θ Z * j+1 + ∂ t Z⋆ j = G j-1
, where G j-1 is a known function. From there, the procedure to deduce (H j+1 ) is a repetition of what has been done in Paragraph 2.3.2. There is just a point that must be stressed. By averaging (2.91) in t and θ, we recover the linear equation :

(2.92) ∂ τ Zj-1 -( Zj-1 •∇ z ) A 1 = Ḡj-1 ,
Zj-1 (0) = -Z⋆ j-1 (0, 0)-Z * j-1 (0, 0, 0) . It follows that the Cauchy problem (2.92) has a solution Zj-1 on the whole interval [0, T ]. It means that the life span of all the Z j is the same as the one of (2.84). Stop the induction when j = N . With the Z j , build Z ε as indicated in (2.79). To identify the profile Φ ε (τ, t, θ), consider the asymptotic expansion proposed in (2.79), and perform a formal analysis at the level of (2.

80). The study of what comes in factor of

ε -1 furnishes with Z0 (τ, t) = Ξ t, Z0 (τ ) = t ( Ȳ0 , Ū0 )(τ, t) ∈ R 3 × R 3 : (2.93) Φp (τ ) = c(w) τ 0 b( Ȳ0 ) (s) ds , Φ⋆ 0 = c(w) ∂ -1 t b( Ȳ0 ) ⋆ .
The contribution with ε j in factor, with j ∈ {0, • • • , N }, gives rise to:

(2.94)

∂ τ Φ j + ∂ t Φ j+1 + j k=0 b k ∂ θ Φ 2-k+j -b j+1 = 0 .
This must be completed by what comes the constraint ϕ ε (0) ≡ 0 in (2.26), that is:

(2.95) Φj (0) = -Φ⋆ j (0, 0) -Φ * j (0, 0, 0) . All coefficients b k have been previously identified. Exploiting (2.94) and (2.95), the Φ j can be determined successively, for j = 1 to N . By piecing the Φ j together, we get Φ ε . 2.3.4. Justification of the WKB approach, weak coupling. Fix any integer N ∈ N * . As long as the solution to (2.26) exists, it can be put in the form (2.78). As a by-product of the formal analysis of Paragraph 2.3, we can rely on (2.80). Since (R z , R φ ) is a smooth function on the compact set [0, 1] × [0, T ] × T T × T, it is uniformly bounded. Taking into account the initial data of (2.20) and the preceding adjustments, we are faced with:

(2.96) ∂ τ r z ε r ϕ ε = R z R φ ε, τ, τ ε , ϕ ε ε , r z ε , r z ε r ϕ ε (0) = - Z * N Φ * N (0, 0, 0) , where: (R z , R φ )(ε, τ, t, θ, r) := ε -N -2 b(ε, t, θ, Z ε ) -b(ε, t, θ, Z ε + ε N r) ∂ θ (Z ε , Φ ε ) + ε -N -1 (A, b) ε, t, θ, Z ε + ε N r) -(A, b) ε, t, θ, Z ε ) -(R z , R φ )(ε, τ, t, θ) . Recall that ∂ θ Z ε = O(ε 2 ) and A 0 ≡ 0. Therefore, there is a constant C ∈ R * + such that: sup |R z (ε, τ, t, θ, r)| ; (ε, τ, t, θ, r) ∈ [0, 1] × [0, T ] × T T × T × B(0, 1] ≤ C (|r| + 1) .
This control translates the fact that the dependence in r of the source term R z is nice, with uniform estimates in (ε, θ) ∈ [0, 1] × T. Note T ε ∈ R * + the life span associated to (2.26), or (it is the same thing) to (2.96). If T ε > T , rename T ε = T . By applying Gronwall's lemma to the first line of (2.96), we get easily:

(2.97) sup |r z ε (τ )| ; (ε, τ ) ∈ ]0, 1] × [0, T ε ] < +∞ . Now, the system (2.96) is also strongly nonlinear in r ϕ ε . However, the right hand terms depend on r ϕ ε only through ϕ ε , that is after substitution of ϕ ε for θ, whereas uniform estimates in θ are available. This is the reason for what is known as a weak coupling. This prevents to be faced with a finite time blow-up, with a time T ε shrinking to zero when ε goes to zero. In practice, recall that ∂ θ Φ ε = O(ε 2 ). Thus, the information (2.97) can be exploited to infer that:

(2.98) sup |ε r ϕ ε (τ )| ; (ε, τ ) ∈ ]0, 1] × [0, T ε ] < +∞ .
This last estimate is not as good as (2.97) but it prevents the solution of (2.96) to explode. We can assert that T ε ≡ T for all ε ∈ ]0, ε 0 ], and that r z ε is indeed a remainder. Concerning the component r φ ε , by implementing (2.98) at the level of (2.96), we observe a loss of precision by a fixed factor ε -1 . Briefly: Lemma 2.14. For all ε ∈ ]0, 1], the solution to the ode (2.26) is defined for all τ ∈ [0, T ]. Moreover, with (Z ε , Φ ε ) as in (2.79), and with a precision in L ∞ that is uniform with respect to τ ∈ [0, T ] and initial data z iε in a compact set, we can assert that:

(2.99) z ε ε ϕ ε (τ ) - Z ε ε Φ ε τ, τ ε , ϕ ε (τ ) ε = O ε N .
2.3.5. Description of the phases in long times. We can now assert that the asymptotic behaviour (when ε goes to zero) of (x ε , v ε ) during long times is given by (2.78) and (2.79). This reveals the presence of two types of oscillations.

⋄ Oscillations at frequencies of the order ε -1 coming from the periodic behaviour in t of the profiles. When incorporating the (smooth) dependence on the initial data (y, w), the period T becomes a function of (y, w). To study the effects of the fast variations in (y, w), it is necessary to put them at a uniform scale, that is to normalize T (say at the value 2 π). This operation makes a phase appear.

Definition 2.6. [notion of secondary phase] The function:

(2.100)

ψ s : (τ, y, w) -→ ψ s (τ, y, w) := 2 π τ /T = τ w/k ž ,
where ž = t (h z , ŷ, w) and h z = √ H(y, ς w ), is called the secondary phase. The expressions H, k h and ς w are given by (2.51), (2.65) and (2.35). In (2.100), the (smooth) dependence on (y, w) has been incorporated.

Introducing Zε (τ, ϑ, θ) := Z ε (τ, T ϑ/2 π, θ), we find that:

(2.101) z ε (τ ) = Zε τ, ψ s (τ ) ε , ϕ ε (τ ) ε 2 + O ε N , Zε ∈ C ∞ ([0, T ] × T 2 ; R 6 ) . ⋄
⋄ Oscillations at frequencies of the order ε -2 associated with ϕ ε . The access to ϕ ε is possible through the differential equation (2.26), which may seem little descriptive, or through the relation (2.78), which may sound nested. In fact, apart from the existence of ϕ ε on [0, T ], the analysis has so far left aside a precise description of what is ϕ ε . Information can be obtained by freezing a leading part of ϕ ε , and by absorbing the rest inside the profiles. This can be done by seeking ϕ ε in the form:

(2.102)

ϕ ε (τ ) = Φp (τ ) ε + Φ0 τ, τ ε + ε Φc ε τ, τ ε , Φp (τ ) ε 2 + Φ0 (τ, ε -1 τ ) ε .
Then, the relation (2.78) on ϕ ε is equivalent to:

(2.103) Φc ε -Φ1 - N -1 j=1 ε j Φ j+1 (θ + Φc ε ) | θ= Φp ε 2 + Φ0 ε = ε N -1 r ϕ ε .
Consider the condition:

(2.104)

I(ε, τ, t, θ, Φε ) := Φε -Φ1 (τ, t) - N -1 j=1 ε j Φ j+1 (τ, t, θ + Φε ) = 0 . Since ∂ ΦI = 1 + O(ε),
ε -Φε = O ε N -2 , so that: (2.105) ϕ ε (τ ) = Φp (τ ) ε + Φ0 τ, τ ε + ε Φε τ, τ ε , Φp (τ )
ε 2 + Φ0 (τ, ε -1 τ ) ε + O ε N -1 .
At this stage, reasoning modulo O ε N -1 , resorting to ϕ ε is no more necessary.

Definition 2.7. [notion of principal phase] Introduce:

ψ p0 (τ ) := Φp (τ ) ∈ C ∞ ([0, T ]; R) , ψ p1 (τ, ϑ) := Φ0 (τ, T ϑ/2 π) ∈ C ∞ ([0, T ] × T; R)
, where the functions Φp and Φ⋆ 0 are given by (2.93), whereas the part Φ0 is determined by solving ∂ τ Φ0b 1 = 0 with initial data Φ0 (0) = -Φ⋆ 0 (0, 0). The expressions ψ p0 and ψ p1 are smooth with respect to (y, w). The application:

ψ pε : [0, T ] × Ω -→ R
(τ, y, w) -→ ψ pε (τ, y, w) := ψ p0 (τ, y, w) + ε ψ p1 τ, y, w, ψs(τ,y,w) ε , is called the principal phase.

The behaviour of the family {z ε } ε when ε goes to zero can be described through a multiscale and multiphase expansion. More precisely: Lemma 2.15. Define Ẑε (τ, ϑ, θ) := Zε τ, ϑ, θ + Φε (τ, T ϑ/2 π, θ) . Then:

(2.106) z ε (τ ) = Ẑε τ, ψ s (τ ) ε , ψ pε (τ ) ε 2 + O ε N , Ẑε ∈ C ∞ ([0, T ] × T 2 ; R 6 ) .
Proof. Replacing Φc ε with Φε gives rise to a loss of ε N -2 . But this is compensated by the fact that ∂ θ Zε = O ε 2 . We still have (2.106) in L ∞ .

In comparison with (2.101), the progress in (2.106) is that ψ pε is a well identified function. The structure of the phase ψ pε reflects, in accordance with (P1) and (P2), some overlapping of oscillations. Intuitively, the phases ψ s and ψ p0 are related to the gyration around the field lines (for respectively intermediate and long time effects), whereas ψ p1 comes from the bouncing back and forth between the two mirror points. Note that the long time derivative of ψ pε is a large amplitude oscillation (compared to the amplitudes of ψ pε or b). Indeed, recalling (2.93), we have:

(2.107) ∂ τ ψ pε (τ, y, w) = c(w) b • Ȳ0 τ, y, w, τ ε + ε (∂ τ Φ0 ) τ, y, w, τ ε .
Using again (2.93), and (2.95) with j = 0 (knowing that Φ * 0 ≡ 0), we get: (2.108) ψ pε (0, y, w) = Φp (0) + ε Φ0 (0) + ε Φ⋆ 0 (0, 0) = 0 . ⋄ • Application 7 [about the principal phase induced by the geomagnetic field]. From (2.107) and (2.108), it is easy to deduce that:

(2.109) ψ pε (ε t, y, w) = ε c(w) t 0 b e • Ȳ0 (0, y, w, s) ds + O ε 2 .
The principal contribution of (2.109), the one with ε in factor, has a specific role to play in Part 3. This is why it is examined below. Lemma 2.16. Recall that T e (h, y 2 , w) and ν e (h, y 1 , y 2 , w) are the two smooth functions which can be determined through the explicit formulas (2.77a) and (2.77b). There exists a scalar function Γ e (•) such that, for all (τ, t) ∈ [0, T ] × R + :

(2.110) b e • Ȳ0 (τ, y, w, t) = Γ e (h, y 2 , ν e + w t) , h = √ H e (y 1 , y 2 , ς w ) .
Proof. Recall (2.11c) and that Ξ j 1 (z) = y j for j ∈ {1, 2}. By construction:

b e • Ȳ0 (τ, t) = b e J1 0 (τ, t), J2 0 (τ, t) , ( J1 0 , J2 0 )(τ, t) := (Ξ 1 2 • Ξ 3 , Ξ 2 2 • Ξ 3 ) t, Z0 (τ, z i0
) where Z0 (τ, z i0 ) is given by (2.86). In the case of Be , the applications Ξ 1 2 • Ξ 3 and Ξ 2 2 • Ξ 3 , given by Definitions 2.4 and 2.5, do not imply the variables y 3 and a:

Ξ 1 2 • Ξ 3 (t, z) = f -1 y 2 h cos K -1 (h,y 2 ) (b + u t) , Ξ 2 2 • Ξ 3 (t, z) = y 2 , z = t (h, y 2 , y 3 , u, a, b
) . In view of (2.86), the four remaining components h, y 2 , u and b are stationary, so that:

J1 0 (τ, t) = J1 0 (0, t) = f -1 z 2 i0 z 1 i0 cos K -1 (z 1 i0 ,z 2 i0 ,z 4 i0 ) (z 6 i0 + z 4 i0 t) , J2 0 (τ, t) = J2 0 (0, t) = z 2 i0 . Finally, noting that z i0 = Zs2
i0 with Zs2 i0 as in (2.70), we can directly recover (2.110), with T e and ν e given by (2.77a) and (2.77b), and:

(2.111) Γ e (h, y 2 , θ) := b e f -1 y 2 h cos K -1 (h,y 2 ) ( θ) , y 2 .
As a by-product of Proposition 2.1, Lemma 2.11 (applied here with F = b e • Ξ 1 ≡ b e ) and Lemma 2.16, with ž = (ž 1 , ž2 , ž4 ) = (h, y 2 , w), we find that: To reach the exact flow (x ε , v ε )(τ, •), see (2.78), the extra terms of the expansion (2.79) and the remainder ε N r z ε must be incorporated. On the one hand, the profiles Zj and Z * j with j ∈ N * are bounded functions (with at least ε in factor). On the one hand r z ε is controlled as indicated in (2.97). This means that (x ε , v ε ) remains at a distance ∼ ε of EBēlt.

(2.112) b e • Ȳ0 (τ, •) = Γ e (h, y 2 ) = -4 π k (h,y 2 ) h 0 r be f -1 y 2 ( √ h 2 -r 2 ),y 2 dr C f -1 y 2 ( √ h 2 -r 2 ),y 2 cos g -1 (r)
To justify the part 3) of Theorem 1, the starting point is Lemma 2.15. The matter is to follow how the asymptotic expansion (2.106) is transformed when coming back to the original variables (x, v). This can be done first by passing from z ε to z ε = t (y ε , u ε ) through (2.29), and then by using the inverse transformation:

(2.114) (x ε , v ε )(τ, x, v) = y -1 (y ε ) , t O(y ε ) e ϕε Λ/ε u ε τ, y(x), O • y(x) v .
The function Ξ(t, •) + ε Ξ(t, θ, •) is smooth, with uniform estimates in (t, θ) ∈ T T × T. On the other hand, since the matrix Λ is skew-symmetric, the action of e ϕε Λ/ε is unitary. It is not at these levels that difficulties appear. The problem comes when replacing ϕ ε at the level of (2.114) by the expression of ϕ ε obtained line (2.105), in order to make ψ pε appear. This introduces a loss of precision by a factor ε -2 . Still, by choosing N large enough, the approximation (1.7) becomes valid for arbitrary powers of ε. The nonlinear geometric optics is therefore completely justified in L ∞ on the long time interval [0, T ].

The last point to highlight is how to keep track of the initial data t (x, v) in (1.5). As already mentioned, there is a smooth dependence of z ε on z iε , when solving (2.26). The aim of (2.87) was to express z iε in function of t (y, w), resulting in Definitions 2.6 and 2.7. Now, to recover (1.7), it suffices to take:

(2.115)

ψ e pε (τ, x, v) := ψ pε τ, y(x), O•y(x) v , ψ e s (τ, x, v) := ψ s τ, y(x), O•y(x) v .
Note that the source term ∇ p M θ (p) • E can be taken into account, through some integral term. As long as E is fixed, this does not change the oscillating structure.

3. On the creation of light inside plasmas.

Our aim in this Part 3 is to get a qualitative mechanism to explain what the spectrograms contain. We focus our attention on the aspects underlined in line (C3) of Subsection 1.2, see the figure 1, p. 623 of [START_REF] Santolik | New results of investigations of whistler-mode chorus emissions[END_REF]. Concerning this issue, the subsequent analysis is convincing. But it is still at a preliminary stage. It must be seen as a stepping stone towards more refined studies. To begin with, it is important to further discuss the model (1.1).

3.1.

From Vlasov-Maxwell equations to oscillatory integrals. The first thing to do is to detail the origin of (1.1), starting from standard formulations [START_REF] Bostan | The Vlasov-Maxwell system with strong initial magnetic field: guiding-center approximation[END_REF][START_REF] Bourdiec | Méthodes déterministes de résolution des équations de Vlasov-Maxwell relativistes en vue du calcul de la dynamique des ceintures de Van Allen[END_REF][START_REF] Pfannmöller | Whistler Wave Propagation in Inhomogeneous Plasmas[END_REF][START_REF] Tenerani | Dynamics of ion-scale coherent magnetic structures and coupling with whistler waves during substorms[END_REF].

3.1.1. Some weakly non linear model issued from space plasma physics. Consider in the earth's magnetosphere, a population of electrons whose density in the phase space is denoted by f(t, x, v). The particles interact mainly by electro-magnetic fields (E, B) which are created collectively (collisions can be neglected). On the other hand, one must take into account the presence of a strong exterior magnetic field Be (x), undergoing non trivial variations in x. The function Be (•) is given by the dipole model of the Earth's magnetic field, see Application 1. With Be (•) as in (2.7), it can be approximated by:

(3.1) Be (x) = b0 Be (R -1 x) , b0 ≃ 3, 1 × 10 -5 T , R ≃ 6, 37 × 10 6 m .
In this context, relativistic Vlasov-Maxwell equations can be written:

(3.2)        ∂ t f + v • ∇ x f + q e (E + v × B + v × Be )• ∇ p f = 0 , ∂ t E -c 2 curl B = - q e ǫ 0 R 3 f v dp , div E = q e ǫ 0 R 3 f dp -n , ∂ t B + curl E = 0 , div B = 0 .
In (3.2), the physical constants are: q e ≃ 1, 6 × 10 -19 C is the (absolute value of the) charge of the electron; c ≃ 3 × 10 8 ms -1 is the speed of light; and ǫ 0 ≃ 8, 8 × 10 -12 F m -1 is the vacuum permittivity. The relativistic framework comes from the the energetic electrons contained in the outer radiation belt. The velocity v has a magnitude |v| which must be limited by c, and it is linked to the relativistic momentum p ∈ R 3 through the relations:

(3.3) v(p) = p m e 1 + |p| 2 m 2 e c 2 -1/2 , p(v) = m e v 1 - |v| 2 c 2 -1/2
, where m e ≃ 9, 1 × 10 -31 kg is the mass of the electron. It follows that:

∇ p = 1 m e 1 - |v| 2 c 2 1/2 ∇ v -∇ v v c ⊗ v c , |D v p(v)| = m 3 e 1 - |v| 2 c 2 -5/2
.

In (3.2), the letter n stands for the concentration of a background ion distribution. The outer radiation belt is globally neutral, meaning that:

(3.4) R 3 f(t, x, p) dp ≃ n .
A common approach is to primarily consider a steady state solution of (3.2), having the form (f, E, B) ≡ (f 0 , 0, 0). To this end, it suffices to select a function f 0 which depends only on |v| (or |p|) and which satisfies (3.4).

From a physics point of view, the expression f 0 can represent the plasma in thermal equilibrium. In a relativistic context, it is usually chosen as the Maxwell-Jüttner distribution:

(3.5) f 0 (t, x, p) ≡ f 0 (p) := f0 M θ p m e c , f0 := n 4 π m 3 e c 3 θ K 2 (1/θ)
.

In (3.5), the function K 2 (•) is the modified Bessel function of the second kind, whereas the temperature θ and the function M θ (•) are given by:

(3.6) θ := k T /(m e c 2 ) , M θ (p) := exp -1 + |p| 2 /θ , p ∈ R 3 .
The description of magnetospheric plasmas through expressions like f 0 is much too vague. Indeed, density measurements reveal fluctuations around f 0 , see [START_REF] Li | Modulation of whistler mode chorus waves: 2. role of density variations[END_REF]. This may involve:

1a. Anisotropic features: variations of f according to the amplitudes of the components v and v ⊥ of the velocity v (which are respectively parallel and orthogonal to the field lines) can be observed during intermediate times (t ∼ 1); 1b. Two-stream effects: the particles are going back and forth between the two mirror points. Typically, above a spatial position x, the function f can be localized near two opposite values of the parallel velocity (v = ±β for some β ∈ R * + ). The two aspects 1a and 1b are not taken into account when dealing with f 0 . Now, a way to capture more precise effects is to investigate what happens in the vicinity of (f 0 , 0, 0). To this end, we can seek the solutions of (3.2) in the form:

(3.7) (f, E, B)(t, x, p) = (f 0 , 0, 0)(p) + ( f0 a f f , c b0 a m E , b0 a m B )(t, x, p) ,
where we have used the nondimensional variables:

(3.8) t := t T , x := x R , v := v c , p := p m e c .
This way, all quantities are scaled relative to their characteristic unit of measure. Note that the multiplication of E by c allows to recover the classical symmetric form of Maxwell equations. The amplification factors a f and a m are introduced to make clear the comparison between the sizes of the perturbations and the amplitudes of f 0 and Be . Fix the observation time T such that T := c -1 R ≃ 10 -2 s. Introduce the two dimensionless quantities:

(3.9)

ω ce := q e b0 m e R c = q e b0 m e T , κ ≡ κ(θ) := q e R c n 4 π θ K 2 (1/θ) ǫ 0 m e .
The quantity q e b0 /m e is known as the electron gyrofrequency (or cyclotron frequency). The number ω ce is the electron gyrofrequency relative to the reference frequency T -1 . From now on, we take ε := 10 -5 ≃ ω -1 ce ≪ 1 as a small parameter governing the evolution. The value of κ is determined by the temperature θ. Recall that the average density in the outer Van Allen belt is of about n ≃ 50 protons/m 3 . It follows that κ(θ) ≃ θ -1/2 K 2 (1/θ) -1/2 . For θ = 1, we get κ ≃ 1. However, κ(•) is a very rapidly decreasing function. For (not so) large values of θ, the selection of κ ≪ 1, up to κ ≃ ω -1 ce , becomes consistent.

With p and v linked as indicated in (1.2), and j and ρ defined as in (1.3), the system (3.2) is transformed into:

(3.10)        ∂ t f + v • ∇ x f + ω ce a m (E + v × B) • ∇ p f + ω ce v × Be (x) • ∇ p f = -ω ce a m a -1 f ∇ p M θ (p) • E , ∂ t E -curl B = -(ω ce a m ) -1 a f κ 2 j , ∂ t B + curl E = 0 , together with: (3.11) div E = (ω ce a m ) -1 a f κ 2 ρ , div B = 0 .
In view of (3.10)- (3.11), the impact of the nonlinearity and the strength of the coupling are driven by a f , a m and κ. We can adjust a m in such a way that ω ce a m = ω -1 ce . We can fix a f to ensure a balance between the source terms in the Vlasov and Maxwell parts. Application of that criterion results in a f := ω ce a m /κ. Then, for θ = 1, we get κ ≃ 1 and a f ≃ ε. For larger values of θ, the selection of κ ≃ ε and a f ≃ 1 is more adapted. Note that the coefficients which are in factor of f , E and B at the level of (3.7) are respectively of size

a f f0 ≃ ε f0 (if κ ≃ 1) or a f f0 ≃ f0 (if κ ≃ ε), c a m b0 ≃ 0, 03 b0 , and a m b0 ≃ ε 2 b0
. The preceding values are in line with a perturbative regime. They allow to express (3.10)- (3.11) in terms of ε and κ, in order to obtain:

(3.12)        ∂ t f + v • ∇ x f + ε (E + v × B) • ∇ p f + ε -1 v × Be (x) • ∇ p f = -κ ∇ p M θ (p) • E , ∂ t E -curl B = -κ j , div E = κ ρ , ∂ t B + curl E = 0 , div B = 0 .
Now, it suffices to change t into τ := ε t in order to recover the scalings of the system (1.1).

According to the formulation (3.12), non linear effects can appear at leading order only during long times τ ≃ ε t ≃ 1. The same applies to the coupling if κ ≃ ε. And even if κ ≃ 1, the effective influence of the coupling should be smaller than it may at first seem. Indeed, as shall be seen later, the actual size of the electric current density j is obtained after applying a method of stationary phase. In addition, except for specific waves (like for instance whistler waves), the integration of (the oscillations contained in) E along the characteristics associated to the Vlasov equation induces cancellations.

The existence theory for systems like (3.2) has a very long history. We can cite for instance the result which has been established in 1988, in the pioneering work [START_REF] Glassey | Global existence for the relativistic Vlasov-Maxwell system with nearly neutral initial data[END_REF] of R. T. Glassey and J. W. Schaeffer. Since then, many improvements have been introduced. To have a view of the subject in a broader context, the reader can for instance refer to the book [START_REF]Topics in kinetic theory[END_REF]. The specificity of (3.12) comes from the presence of a small parameter ε and of a non constant exterior magnetic field Be (•). Beyond the existence theory, in connection with applications in physics, the special structures of the solutions to (3.12) are important aspects to develop.

As will be described in the next paragraphs, from this perspective, the introduction of ε and Be (•) can bring about a change in perspective.

3.1.2.

A brief foray into previous attempts. To incorporate the two characteristics 1a and 1b at the level of (3.2), a standard procedure in plasma physics [START_REF] Bourdiec | Méthodes déterministes de résolution des équations de Vlasov-Maxwell relativistes en vue du calcul de la dynamique des ceintures de Van Allen[END_REF][START_REF] Pfannmöller | Whistler Wave Propagation in Inhomogeneous Plasmas[END_REF][START_REF] Tenerani | Dynamics of ion-scale coherent magnetic structures and coupling with whistler waves during substorms[END_REF] is to modify at time t = 0 the initial data f 0 into f 0 + ν f 1 with ν ≪ 1, where f 1 is (in case 1a) a bi-Maxwellian distribution with thermal anisotropy:

(3.13) f 1 (v) := 1 β π 3/2 α 3 exp - v 2 α 2 exp - |v ⊥ | 2 β α 2 , β ∈ R * + ,
and (in case 0b) a bimodal Maxwellian distribution peaked around opposite speeds:

(3.14) f 1 (v) := 1 2 π 3/2 α 3 exp - (v -β) 2 α 2 + exp - (v + β) 2 α 2 exp - |v ⊥ | 2 α 2 .
There are practical reasons for preferring to change f |t=0 ≡ f 0 , without modifying E |t=0 ≡ 0 or B |t=0 ≡ 0. It is indeed expected that whistlers and related phenomena are initiated by changes of density, or strongly connected to density variations [START_REF] Li | Modulation of whistler mode chorus waves: 2. role of density variations[END_REF]:

2a. This can be caused by lightning effects [START_REF] Green | Lightning Effects on Space Plasmas and Applications[END_REF]; 2b. Looking upstream, it may also come from sudden changes [START_REF] Park | Whistler observations during a magnetospheric sudden impulse[END_REF][START_REF] Tenerani | Dynamics of ion-scale coherent magnetic structures and coupling with whistler waves during substorms[END_REF] in the solar wind pressure (as a consequence of solar flares).

The solutions f(ν, t, x, v) of (3.2), which are issued from initial data like

(f 0 + ν f 1 )(v) can always be put in the form f(ν, t, x, v) = f 0 (v) + ν f l 1 (t, x, v) + O ν 2 with f(0, t, x, v) = f 0 (v) where f l
1 is the solution to the equations (3.2) linearized around (f 0 , 0, 0). Numerical studies [START_REF] Bourdiec | Méthodes déterministes de résolution des équations de Vlasov-Maxwell relativistes en vue du calcul de la dynamique des ceintures de Van Allen[END_REF] indicate that the time behaviour of f l 1 might undergo weibel instabilities (in case 0a) or two-stream instabilities (in case 1b). These two types of instabilities are among the factors which are put forward to account for the sudden onset of whistler wave phenomena. However, there are some potential drawbacks to this approach: 2c. Contrary to the observations about the spatial location of the Van Allen belts, objects like f 1 are not well localized in space; 2d. Neither f 0 nor f 1 take into account the oscillations induced by the Vlasov part, and especially the occurrence of several different scales, revealed at the level of (1.1) by the nondimensionalization of the equations; 2e. In connection with the two preceding items 2c and 2d, the linearization procedure of (3.2) along f 0 emphasizes the instabilities induced by the term

( Ė + v × Ḃ) • ∇ p f 0 .
This aspect is found in (3.12) through the source terms. But the corresponding effects, although interesting, are not sure to prevail. 2f. This perturbation theory fails to really grasp the geometrical impact induced by the variations in x of the strong magnetic field Be (•). It just so happens that these variations are in fact at the core of the subject.

Other aspects deserve to be developed in addition to the foregoing. As already explained, it seems reasonable to associate (3.12) with some initial data having the form (f, E, B) | t=0 = (f 0 , 0, 0). Clarifications are needed as to the choice of f 0 (x, p). This must be done in connection with 2a and 2b.

• In the case of lightning discharges (2a), the support of f 0 (•, v) should be very localized.

That is why such events are so interesting from the practitioner's point of view. They give some experimental access to the properties of the characteristic flow (1.4).

• In the case of magnetic storms issued from solar flares (2b), the situation is less clear, since the perturbations come after complicated mechanisms. That said, less than half of the earth should be lit in that way.

In both cases 2a and 2b, the function f 0 (•, v) should be contained (uniformly in v) in a compact neighbourhood of the Van Allen belt. As already explained in Application 4, from a dynamical point of view, the propagation is organized by energy levels. At a constant energy level h, variations in y 1 are clearly linked with variations in ς. And finally, to get an adequate localization near the extended Van Allen belt, it is more convenient to argue with the variables h, y 2 , y 3 and u.

Another point is that the kinetic energy E c of the fast-moving particles (beams of electrons) is presumably close to a specific value. The values which are put forward concerning E c are typically of the order of ≃ 1 Mev. It should be recalled that E c depends on the amplitude of the velocity. To provide with some orders of magnitude, note that |v| ≃ 0, 95 c is equivalent to c(u) ≃ 0, 33, and this means that E c ≃ 1, 69 Mev. This furnishes a physical argument to work near the particular value |v| ≃ 0, 95 c. This can be done by assuming that f 0 (x, •) is of the form (3.14) with β = 0, 95 and α replaced by α/c (to take into account the change of scale v = c v). Since, it is expected that α/c ≪ 1, this corresponds to a highly concentrated distribution. In practice (3.17), we will take a Dirac mass δ w (•) at w = 0, 95. This choice avoids technicalities introduced by the mollified version of δ w (•). It is also consistent to some extent. Indeed, as long as E ≡ 0, the amplitude of the velocity is a conserved quantity, and we can examine separately what happens in the case of different values of |v|.

Recall that the control H e (y 1 , y 2 , ς, u) ≤ C gives rise to a bound on both y 1 and ς. Thus, given u ∈ [0, c], the above points can be taken into account by introducing adequate bump functions having the form:

(3.15) Ψ(y, ς, u, a) = Ψ √ H e (y 1 , y 2 , ς, u), y 2 , y 3 , Ψ ∈ C ∞ 0 (R + × R * -× R; R) .
To be more precise, it is convenient to select the values (h, y 2 ) coming from the identification of the extended Van Allen belt. To this end, it suffices to impose:

(3.16) supp Ψ(•) ⊂ K f := (h, y 2 , y 3 ) ; 0 ≤ h ≤ f ey 2 h e (y 2 ) , y 2 ∈ supp h e (•) , y 3 ∈ J .
Consider the system (1.1) written in the variables (y, ς, u, a) and finally, adjust the initial data f 0 of (3.12) in such a way that:

(3.17)

f 0 (y, ς, u, a) = Ψ(y, ς, u, a) f 0 (y, ς, u, a) δ w (u) ,
where δ w (•) is the Dirac mass at w (or, in the nonlinear case, a mollified version of it), whereas f 0 (•) is any smooth function. The system (3.12) is non linear. The parameter ε being fixed, it can be faced by an iterative method along the following lines:

• Preliminary step 0. Solve the Vlasov-Maxwell system (3.12) without the nonlinearities and without the source terms. There is no coupling. This amounts to take (E, B) ≡ (0, 0), and to consider the Vlasov equation, that is (1.5) in the case (E, B) ≡ (0, 0), together with (1.4) in the case f |t=0 = f 0 . This furnishes some distribution f 0 ≡ f 0 ε , and therefore some expression (f 0 , E 0 , B 0 ) = (f 0 , 0, 0).

In Part 2, we have achieved the step 0, yielding f 0 ε . Dealing with f 0 ε (t, x, p) in place of f 1 (v) is more significant. Indeed, in view of Section 2, we can assert that: 3a. The expression f 0 ε is well localized in space and velocity. The three quantities √ H e , y 2 and u are approximately or exactly preserved during the long time evolution. On the contrary, the coordinate y 3 undergoes some drift: the support of f 0 ε (τ, •) evolves to follow the ring current inside the Van Allen belt; 3b. The function f 0 ε inherits all the characteristics forced by the fast dynamics contained in the Vlasov part, including the various oscillations and, above each position x ∈ R 3 (with varying parameters along the field lines), some thermal anisotropy and also some two-stream effects; 3c. The geometric features of Be run through all the analysis of Part 2, leading finally to the refined notions of principal and secondary phases.

• Following steps n (with n ∈ N * ). Assume that (f n-1 , E n-1 , B n-1
) is a known expression.

Then, look at the Cauchy problem which is obtained by linearizing the equations of (3.12) along (f n-1 , E n-1 , B n-1 ). This is a system of linear PDEs which can be written as:

L(ε, t, x, p; ∂ t , ∇ x , ∇ p ) U = 0 , (t, x, p) ∈ R + × R 3 × R 3 , U = t (f, E, B) .
3.1.4. Connection to geometric optics. The fundamental object on which the propagation of singularities can be based is a subset of the cotangent bundle, denoted by Char ε L, and called the characteristic variety. Use the same letter L for the complete symbol associated with L, and note R n * := R n \ {0}. Define:

Char ε L := (t, x, p; s, ξ, η) ∈ R + × R 6 × R 7
* ; det L(ε, t, x, p; s, ξ, η) = 0 . It's not easy to make use of the information contained in Char ε L. It turns out that, due to the dependence on the parameter ε, the choice of the polarization, or according to the level of detail required, the discussion may take a wide variety of forms. We will limit ourselves here to pioneering efforts in the analysis of Char ε L. Before going forward, let us start with a brief overview of what is known.

The dynamical system (1.5) is also what allows to integrate pressureless gases driven by a strong magnetic field, that is:

(3.18) ∂ t u + (u • ∇)u + 1 ε Be (x) × u = 0 , u ∈ R 3 .
As long as t ∈ [0, T ] with T small enough, the asymptotic analysis of (3.18) can be completed by adapting to the three dimensional situation the method introduced in [START_REF] Gallagher | On pressureless gases driven by a strong inhomogeneous magnetic field[END_REF]. However, the approach of [START_REF] Gallagher | On pressureless gases driven by a strong inhomogeneous magnetic field[END_REF] cannot be extended to large values of T , and a fortiori to long times, due to a systematic crossing of the spatial characteristics. To go further, a more refined model is needed, based here on a kinetic description. All the steps of Part 2 can be translated in terms of manipulations involving the Vlasov equations. This can be done at the cost of many technical complications. Partial information on what happens is provided below.

Assume first that the coupling is inoperative, meaning that κ ≃ ε and that the terms of size ε can be removed. Then, the characteristic variety Char ε L contains in particular the set Char f ε L, which includes the polarizations along the particle density f . Recall the notations of (2.15) and of (2.34)-(2.35). After straightening out the field lines as in the step 2.1.2, the description of Char f ε L can be deduced from:

(3.19) ∂ t f + d j w j ∂ y j f + c(|w|) (B 1 × w) • ∇ w f + ε -1 λ ∂ a f = 0 , λ := c(|w|) b(y) |w ⊥ | .
Observe that ∂ a λ ≡ 0. Performing a Fourier analysis in the variable a ∈ T, the large term is transformed into large skew-symmetric terms of order zero, one for each Fourier mode. This appears to be in accordance with classical frameworks in geometric optics [START_REF] Métivier | The Mathematics of Nonlinear Optics[END_REF][START_REF] Rauch | Hyperbolic Partial Differential Equations and Geometric Optics[END_REF]. However, there is an obstacle that prevents to implement the above strategy: the other coefficients of (3.19) depend on a. To deal with this, the idea of Part 2 is to absorb the variations in a through the (partial: when w ⊥ = 0) ellipticity of λ ∂ a , and to extract from what remains some eikonal equation. This yields the phase ϕ 0 which is issued from (2.19) in the case ε = 0. Note that ∂ a ϕ 0 ≡ 0, but that ∇ w ϕ 0 ≡ 0. This means that the phase ϕ 0 which is forced by (3.19) is definitely a kinetic object. The same remarks apply to the principal and secondary phases.

It is worth noting that the approach of Section 2.1 to handle (3.19) during intermediate times is not standard. Now, the discussion in Sections 2.2 and 2.3 to reach long times τ ∈ [0, T ] with some fixed T ∈ R * + goes far beyond this. A key idea is the rectification procedure of Paragraph 2.2.2, through (2.29). The transformation (2.29) involves small amplitude oscillations. It must be differentiated in order to be interpreted in terms of wave front sets. This results in a large deformation (of size one, and not regular in ε ∈ ]0, 1]) of Char f ε L. It turns out that Part 2 is not based on a direct study of Char f ε L. There is a need for a work of preparation, without which it seems not possible to make progress.

The coupling with Maxwell equations adds yet another difficulty. The oscillating structure along ϕ 0 (t, x, p) is destroyed by the integration with respect to p, when computing the electric current j. Applying the stationary phase method, what remains are space-time oscillations implying basically the trace of ϕ 0 (t, x, •) at the critical points, which turn to be the mirror points (see Lemma 3.5). Call T r(ϕ 0 )(t, x) the phase thus generated. The general laws of geometric optics say that the corresponding oscillations: i) are absorbed if t, x, p; ∇ t,x T r(ϕ 0 )(t, x), 0 ∈ Char ε L for all p ∈ R 3 ;

ii) can be propagated if t, x, p; ∇ t,x T r(ϕ 0 )(t, x), 0 ∈ Char ε L for some p ∈ R 3 .

In the next Paragraph 3.1.5, the passage from ϕ 0 (t, x, p) to T r(ϕ 0 )(t, x) is interpreted from a physical point of view.

3.1.5. Wave-particle interactions seen from the perspective of oscillations. At this stage of the analysis, what's going on with oscillations seems to be of a two-prong nature. There are the dynamical features described in Part 2. The large term ε -1 v × Be (x) • ∇ p contained in the Vlasov part has the effect of forcing the presence of oscillations. The structures of these oscillations are encoded in the phases ψ s and ψ pε . The fact that ψ s does not depend on (E, B), and that the components ψ p0 and ψ p1 of ψ pε depend on (E, B) in a smooth way, is of particular interest. This is an indication of stability. This is a basic point because this justifies why the detected phases have a full meaning. Another issue is to discern how the forced oscillations are handled by the full system (3.12). We have seen that they coexist with space-time oscillations. This establishes a connection between features related to particles and waves. We emphasize that this is done through a linear effect: the coupling through the linear term j. As long as the nonlinearities can be neglected, the various oscillations can propagate independently. When dealing with linear hyperbolic equations, this stems from a principle of superposition of waves. All of the above draws on a mechanism of wave-particle interaction that relies first on linear aspects:

5a. The two phases ψ s (τ, x, v) and ψ pε (τ, x, v) are kinetic objects (they depend on v), emanating from a macroscopic vision of the motion of particles; 5b. The phase T r(ϕ 0 )(t, x) is a fluid object (from now on, "fluid" means something not depending on v), coming from the impact on the electric current j of the forced oscillations (involving ψ pε ). It leads potentially to the creation of plasma waves; before incorporating nonlinear facets:

5c. Nonlinear interactions of the kinetic and fluid oscillations, (in)stability analysis. The contrast between kinetic and fluid behaviors is reinforced by the fact that different polarizations are involved. This means that disjoint parts of Char ε L are used. Each of them requires a specific approach. We have already explained that a good strategy in the case 5a is to deform the part Char f ε L of Char ε L. But what about the case 5b ? 3.1.6. Dispersion relation for whistler waves. Subparagraphs i) and ii) of Section 3.1.4 suggest that the values of p are implied when discussing about T r(ϕ 0 ). This can produce a false impression. As a matter of fact, plasma physicists consider the phase T r(ϕ 0 ) as a fluid object, and they deal with it through the fluid model:

(3.20) ∂ τ (ǫ r E ε ) -ε -1 curl B ε = -ε -1 κ j ε , div (ǫ r E ε ) = κ ρ ε , ∂ τ B ε + ε -1 curl E ε = 0 , div B ε = 0 ,
where j ε and ρ ε are given by (1.3), whereas the symbol ǫ r ( = ε !) stands for a positive definite matrix, called the relative permittivity. However, how to compute ǫ r ? Most books [START_REF] Howard | Introduction to Plasma physics[END_REF][START_REF] Piel | Plasma physics: an introduction to laboratory, space, and fusion plasmas[END_REF][START_REF] Stix | Waves in Plasmas[END_REF] in plasma physics contain a chapter on waves in cold magnetized plasmas. There, details of ǫ r are provided. There are two main ways to determine ǫ r .

The first [START_REF] Howard | Introduction to Plasma physics[END_REF][START_REF] Piel | Plasma physics: an introduction to laboratory, space, and fusion plasmas[END_REF] relies on single particle models. The second way, in more specialized contributions as [START_REF] Ciurea-Borcia | Oblique whistler waves generated in cold plasma by relativistic electron beams[END_REF][START_REF] Stix | Waves in Plasmas[END_REF], is to linearize (3.2) along some distribution f 0 (p), like (3.5).

In both cases, on the Fourier side in (t, x), the Vlasov equation is used to express integrals of f with respect to p, like , in function of Ê. The remaining system involves only ( Ê, B).

It can be viewed as the Fourier side of (the left part of) the equation (3.20).

Having identified the matrix ǫ r , various dispersion relations can be derived from (3.20). This is usually achieved in the context of Appleton-Hartree equations, see [START_REF] Helliwell | Whistlers and Related Ionospheric Phenomena[END_REF]. For righthanded waves propagating parallel to the magnetic field, the relevant dispersion curve is drawn in [START_REF] Helliwell | Whistlers and Related Ionospheric Phenomena[END_REF][START_REF] Piel | Plasma physics: an introduction to laboratory, space, and fusion plasmas[END_REF]. This can be specified for wavenumbers k in a range corresponding to whistler modes, that is k ∼ 10 -5 m -1 when dealing with the typical frequency 1 -10 kHz.

Then, given some constant c ∈ ]0, 1/2], the following model law (see 5.3.3 or [START_REF] Helliwell | Whistlers and Related Ionospheric Phenomena[END_REF]) for ω(k) is proposed by physicists:

(3.21) c 2 |k| 2 ω(k) 2 = 1 + ε ω 2 pe ω(k) (1 -ε ω(k)) ∼ ε ω 2 pe ω(k) (1 -ε ω(k)) , c ε -1 ≤ ω(k) ≤ ε -1 ,
The formula (3.21) has the advantage of being robust and consensual [START_REF] Ciurea-Borcia | Oblique whistler waves generated in cold plasma by relativistic electron beams[END_REF][START_REF] Helliwell | Whistlers and Related Ionospheric Phenomena[END_REF][START_REF] Howard | Introduction to Plasma physics[END_REF][START_REF] Piel | Plasma physics: an introduction to laboratory, space, and fusion plasmas[END_REF] The symbol ω pe represents the electron plasma frequency. The ratio µ between the electron frequency and the plasma frequency is a small parameter, typically of the order 10 -1 . Please refer to [START_REF] Pfannmöller | Whistler Wave Propagation in Inhomogeneous Plasmas[END_REF]-p. 51 for a representative sample of experimental values. During geomagnetic disturbances, the number density of electrons n e and therefore ω pe (which depends linearly on √ n e ) increases inside whistler ducts. As a result, the overdense condition, which means that µ is small, can be reinforced up to the value µ ∼ 10 -2 or a little less, but still with ε ≪ µ. In view of this mechanism [START_REF] Helliwell | Whistlers and Related Ionospheric Phenomena[END_REF], the equality in (3.21) can be replaced by the equivalence ∼ on the right side of (3.21), leading to a simplified relation:

(3.23) ω(ξ) = ± 1 ε λ(ε µ |ξ|) , λ(z) := λ M z 2 1 + z 2 ≤ 1 .
The factor λ M will be called the threshold limit (of the dispersion law for whistler waves). It appears when the frequency ω(k) of (3.21), which is expressed in seconds, is changed to be in accordance with the scale of t, meaning that λ M ≃ T ∈ R * + . Recall the line (2.8), and that c(w) ≃ 3 -1 . The following quantities: In what follows, because this is more practical, we will work with (y, w), instead of the original variables (x, v). The change of variables performed in Paragraph 2.1.2 is orthogonal, whereas ω(•) depends only on the norm |ξ| of ξ. Thus, the formula (3.23) is not modified by such a transformation. Now, it is necessary to make a few comments with regard to the Assumption 3.1, in order to clarify its scope of validity but also its limitations:

(3.24) c(w) b m ≃ 3 -7 , λ M ≃ 10 -2 , c (
6a. The model (3.23) is sufficient to explain the creation of light inside a plasma, in the form of whistler-mode chorus emissions, as it is described in the introduction through the physical observation (C3) ; 6b. The law (3.23) is sure to be not specific enough to take proper account of the propagation of light, in the form of whistler waves. The affirmation 6a will be confirmed by the proofs of Propositions 3.1 and 3.2. We will see that, to obtain the Informal Statement 1, what counts most are the scalings of (3.23) and the majoration of λ(•), that is the existence of the threshold limit λ M .

To explain what is meant by 6b, the following remarks about (3.21) are in order. All cited references [START_REF] Ciurea-Borcia | Oblique whistler waves generated in cold plasma by relativistic electron beams[END_REF][START_REF] Helliwell | Whistlers and Related Ionospheric Phenomena[END_REF][START_REF] Howard | Introduction to Plasma physics[END_REF][START_REF] Piel | Plasma physics: an introduction to laboratory, space, and fusion plasmas[END_REF], as well as related works, deal with the case of a constant exterior magnetic field Be , and the law (3.21) only makes sense for a wave vector k which is parallel to the direction of this vector Be . Moreover, they start from formulations (kinetic equations or particle models) which are not in dimensionless form. The advantage is that this allows a high degree of generality. The downside is that this fails to grasp specific phenomena, which can be revealed by the presence of large or small dimensionless parameters, as in (3.12). There is still plenty of work to determine in the context of (3.12) what is the precise structure of the function λ(•), that is how it can depend on ε, x and ξ, with ξ not necessarily parallel to Be (x). Without a refined description of λ(•), it is not valid to capture the geometrical aspects underlying the propagation of whistler waves.

3.1.7. A class of oscillatory integrals. In this article, our aim is to give a straightforward access to the major aspects involved during the creation of whistler waves. With this in mind, we fix f ε as in (1.4), we compute j ε as in (1.3), and we examine the electromagnetic wave (E ε , B ε ) which is generated by the linear Cauchy problem made of (3.20) together with the initial data (E ε , B ε ) | t=0 = (0, 0). As usual when dealing with hyperbolic systems [START_REF] Joly | Nonlinear oscillations beyond caustics[END_REF], the solution (E ε , B ε ) can be represented as a finite sum, like:

(3.27) E ε B ε (t, x) = i τ 0 R 3
e i (τ -s) ω i (ξ)/ε e i (x-y)•ξ Π i (ξ) j ε (s, y) ds dy dξ .

In (3.27), the mapping Π i (•) stands for a spectral projector coming from (3.20), after some Fourier analysis [START_REF] Joly | Nonlinear oscillations beyond caustics[END_REF]. The scalar function λ i (ξ) is the associated eigenvalue.

Both Π i (•) and λ i (•) are C ∞ functions on R 3 \ {0}. They are positively homogeneous of degree zero and one respectively. The phenomena involving whistler waves can be recorded at the level of (3.27) by selecting the indice i giving rise to the dispersion relation ω i (ξ) satisfying (3.23). Then, the matter is to describe precisely, for data microlocalized at physically pertinent frequencies, the contribution brought by the integral corresponding to this choice of i. Particular care must be paid to the dependence on f ε via j ε .

At the level of (3.20), the constraint on the divergence is propagated. It can be forgotten. What matters is the source term j, which is now given by an integration with respect to w. with by construction:

Note (y ε , ς ε , u ε , a ε )(τ,
(3.29) fε (τ, y, ς, a) = w 2 sin ς (1 -w 2 ) 5/2 f 0 (y ε , ς ε , w, a ε )(-τ, y, w) × t O(y) w(ς, a) .
The oscillating structure of (y ε , ς ε , a ε ) which has been exhibited in Part 2 is preserved under nonlinear composition. Thereby:

(3.30) fε (τ, y, ς, a) = F ε, τ, y, ς, a, ψ s (τ, y, w) ε , ψ pε (τ, y, w) ε 2 + O ε N . The profile F ∈ C ∞ ([0, 1] × [0, T ] × R 3 × T 4
; R) can be expanded in a Fourier series to get:

(3.31) F = (p,q)∈Z 2 a pq (ε, τ, y, ς, a) e i (p ϑ+q θ) , a pq ∈ C ∞ ([0, 1]×[0, T ]×R 3 ×T 2 ; R) , with:

a pq (ε, τ, y, ς, a) = N j=0 ε j a pq j (τ, y, ς, a) + O ε N +1 , ∀ (p, q, N ) ∈ Z 2 × N .
Denote by "supp f " the support of the function f . Combining the definition (3.17) of f 0 and (3.29), we can obtain:

(3.32) supp fε (τ, •) ⊂ supp Ψ (y ε , ς ε , w, a ε )(-τ, •, w) .
Let us keep in mind the definition (3.15). Look at the (three) coordinates of Ψ after the substitution that is indicated in (3.32). As a byproduct of Section 2.3, they are smooth functions of τ , and they involve oscillations of small amplitude in ε. It follows that their derivatives with respect to both τ and t are uniformly bounded in ε ∈ ]0, 1]. Thereby, the function F (ε, τ, •, a) remains compactly supported in (y, ς) as long as τ is bounded. This property is transfered to the coefficients a pq j . Fix T ∈ R * + as in Paragraph 2.3.2. There is

a compact set K ⊂ R × R * -× R×]0, π[ such that: (3.33) supp a pq j (τ, •, a) ⊂ K , ∀ (τ, a, p, q, j) ∈ [0, T ] × T × Z 2 × N .
The frequencies implied in VLF signals are of the order of 1 kHz (in a range from 100 Hz to 10 kHz). To capture the related phenomena in terms of the present scalings, we have to introduce a localization of ξ in an annulus of size µ |ξ| ∼ ε -1 . This can be achieved by implementing a cutoff function:

(3.34) χ ∈ C ∞ 0 (R * + ; R + ) , χ ≡ 0 , supp χ ⊂ A := ξ ; 0 < r ≤ |ξ| ≤ R < +∞ .
To be sure that µ |ξ| ∼ ε -1 , we multiply by χ(ε µ |ξ|) on the Fourier side. The equation (3.20) on (E ε , B ε ) is linear. It follows that the contributions brought by the different projections Π i and the various terms in the expansion (3.31) can be studied separately. On this basis, the solution (E ε , B ε ) to (3.20) can be viewed as a linear combination:

(3.35) E ε B ε (ε t, x) = i ± p∈Z q∈Z j
I(i, ±, p, q, j; ε, t, x) , of oscillatory integrals having the following form:

(3.36) I(ε, t, x) := -κ ε 1-j ε t 0 K×T×R 3 e i Ψ(ε,ε t,x;s,y,ς,a,ξ)/ε 2 ã(s, y, ς, a, ξ) dV,
where dV is for the volume element dV := ds dy dς da dξ . For the indices i giving rise to (3.23), as explained in Paragraph 3.1.6, we find that (after specifying again the parameters): Ψ(±, p, q; ε, τ, x; s, y, ς, a, ξ) := ± (τs) λ(ε µ |ξ|) + (xy) • (ε 2 ξ) + ε p ψ s s, y, w(ς, a) + q ψ pε s, y, w(ς, a) , ã(p, q, j; s, y, ς, a, ξ) := χ(ε µ |ξ|) Π i (ξ) a pq j (s, y, ς, a) . At the level of (3.36), perform the change of variables s = ε s and ξ = ε -1 µ -1 ξ. Rename s and ξ as respectively s and ξ to find:

(3.37) I(ε, t, x) := κ ε 3-j t 0 K×T×R 3 e i Ψ(ε,ε t,x;εs,y,ς,a,(εµ) -1 ξ)/ε 2 ã εs, y, ς, a, (εµ) -1 ξ dV.

Look at the part ψ pε (εs, y, w) inside the function Ψ(•) of (3.37). Recall (2.109) and (2.110). Since w is fixed, and because y 3 and a do not appear at the level of b e • Ȳ0 (Lemma 2.16), we can skip the dependence on (y 3 , w, a) at the level of (2.109) which becomes:

(3.38) ψ pε (εs, y, w) = ε c(w) s 0 b e • Ȳ0 (0, y 1 , y 2 , ς, r) dr + O ε 2 .
It follows that the function Ψ(•) of (3.37) is of size O(ε). Moreover, since ψ s is linear in s, see (2.100), the contribution ε p ψ s (εs, •) is in fact of size O ε 2 . Isolating the main term, that is the one with ε in factor, makes a new phase Ψ(±, q; t, x; s, y, ς, ξ) appear, namely:

(3.39) Ψ ≡ Ψ(±, q; •) := ± (t -s) λ(|ξ|) + µ -1 (x -y) • ξ + q c(w) s 0 b e • Ȳ0 (0, y 1 , y 2 , ς, r) dr .
Since ∂ a Ψ ≡ 0, we can apply Fubini to put apart the integration with respect to a :

(3.40) I(ε, t, x) = t 0 K×R 3 e i Ψ(t,x;s,y,ς,ξ)/ε a(ε, s, y, ς, ξ) ds dy dς dξ , where a(p, q, j;

•) ∈ C ∞ (]0, 1] × [0, t] × R × R * -× R × T × R 3 ) is given by: a(ε, s, y, ς, ξ) := -κ ε 3-j χ(|ξ|) Π i (ξ) T e 2 i π p s
Te(y,w) e i Φ(ε,s,y,w) a pq j (ε s, y, ζ, a) da , with w ≡ w(ς, a) and, for the sake of completeness: Φ(ε, s, y, w) := q c(w)

s 0 ∂ τ Φ0 (εr, •) + ε -1 b • Ȳ0 (εr, •) -b • Ȳ0 (0, •)] | • ≡(y,w,r) dr .
From (3.33) and (3.34), it follows that the function a(•) is uniformly compactly supported. More precisely, with K := K × A compact, for all (p, q, j) ∈ Z 2 × N, we have:

(3.41) supp a(ε, s, •) ⊂ K , ∀ (ε, s) ∈ ]0, 1] × R .
We can recognize I(ε, •) as an oscillatory integral implying the phase Ψ given by (3.39) and the profile a. The considerations of Section 3.1 about (linear) wave-particle interactions, that is about the connection between the geometry (of forced oscillations) and the dispersion relation (of whistler waves) are now combined into the single object Ψ.

Stationary phase expansions.

The purpose of this Section 3.2 is to examine the asymptotic behaviour of I(ε, t, x) when ε ∈ ]0, 1] goes to zero. In this perspective, the influence of (p, q, j) inside a(•) plays no role. Moreover, we can normalize a by putting the multiplication by ε j-3 apart. In fact, what counts is only the smoothness of a(•) and to impose (3.41) on a(•). Nonetheless, it is important to take into account the sign ± and the choice of q ∈ Z, which are involved in the phase Ψ. For this reason, when needed, we will mark the dependence on (±, q). To describe I ≡ I(±, q; ε, t, x), we will adopt the perspective developed in [START_REF] Joly | Nonlinear oscillations beyond caustics[END_REF], with adaptations. Fix a position (t, x). A preliminary step is to locate the critical points C ≡ C(±, q; t, x) of the phase Ψ ≡ Ψ(±, q; t, x). When doing this, particular care must be given to the initial time s = 0 and to the final time s = t, which are boundary points. The set C can be separated into C = C i ∪ C b , where C i collects the interior critical points:

C i ≡ C i (±, q; t, x) := (s, y, ς, ξ) ∈ ]0, t[×K ; ∇ s,y,ς,ξ Ψ(t, x; s, y, ς, ξ) = 0 ,

whereas C b = C 0 b ∪ C t b
includes the boundary critical points of type r ∈ {0, t}, that is: C r b ≡ C r b (±, q; t, x) := (y, ς, ξ) ∈ K ; ∇ y,ς,ξ Ψ(t, x; r, y, ς, ξ) = 0 , r ∈ {0, t} . The non-stationary phase principle says that, when computing I, the contribution brought by integrating outside any fixed neighbourhood of C is of size O(ε ∞ ). This allows to restrict the discussion to points located in C. Now, the differential conditions involved in the definition of C can be decomposed into:

∂ s Ψ = 0 ⇐⇒ 0 = ∓ λ(|ξ|) + q c(w) b e • Ȳ0 (0, y 1 , y 2 , ς, s) , (3.42a) ∇ y Ψ = 0 ⇐⇒ 0 = -ξ + µ q c(w) s 0 ∇ y (b e • Ȳ0 )(0, y 1 , y 2 , ς, r) dr , (3.42b) ∂ ς Ψ = 0 ⇐⇒ 0 = s 0 ∂ ς (b e • Ȳ0 )(0, y 1 , y 2 , ς, r) dr , (3.42c) ∇ ξ Ψ = 0 ⇐⇒ 0 = ± µ (t -s) λ ′ (|ξ|) |ξ| -1 ξ + x -y . (3.42d)
To move forward, it is necessary to gather more information on the function b e • Ȳ0 , which represents how the amplitude of the geomagnetic field fluctuates along the field lines. This is done in Paragraph 3.2.1. The study of C b is based on the relations (3.42b) and (3.42c). In Paragraph 3.2.2, we show that C b = ∅ for t large enough. In Paragraph 3.2.3, we explain why all harmonics are not relevant (due to a dispersive effect), and why only mirror points can bring non trivial contributions. Then, it remains to examine C i . Note that the three conditions (3.42a), (3.42b) and (3.42c) do not depend on the variables (t, x), except through the condition s ≤ t. They select, independently of the choice of (t, x), the emission points (s, y, ς, ξ), that is the positions from which the signals are leaving. Their description is achieved in Paragraphs 3.2.5 and 3. with for example, in view of (2.110):

0 < b m (y 1 , y 2 , ς) ≡ b m (h, y 2 ) := inf Γ e (h, y 2 , w, •) , h = √ H e (y 1 , y 2 , ς) .
Again, due to the periodicity, there are optimal controls C j (y 1 , y 2 , ς) ∈ R * + with j ∈ {1, 2, 3} such that, for all s ∈ R + :

∇ y (b e • Ȳ0 )(0, y 1 , y 2 , ς, s) ≤ C 1 (y 1 , y 2 , ς) , (3.45a) s 0 ∇ y (b e • Ȳ0 ) ⋆ (0, y 1 , y 2 , ς, r) dr ≤ C 2 (y 1 , y 2 , ς) . (3.45b) s 0 ∂ h (b e • Ȳ0 ) ⋆ dr ≤ C 3 (y 1 , y 2 , ς) . (3.45c)
Then, compactness arguments allow to derive uniform bounds:

(3.46) C j (y 1 , y 2 , ς) ≤ C j < +∞ , ∀ (y 1 , y 2 , ς) ∈ K × T , ∀ j ∈ {1, 2, 3} .
This line, completed with (2.11c) and (2.49), gives rise to h -1 ∂ h Γ e (h, y 2 ) = c + O h 2 with c > 0. It suffices to select h ∈ R * + small enough to recover (3.49). It is even possible to quantify a maximal size of h allowing to get the lower bounds of (3.49). Since the remainders are O h 2 in place of O(h), the effect is enhanced. Lemma 3.2. [extra properties related to Γ e after substitution of h for √ H e (y 1 , y 2 , ς)] There are positive constants c 6 and C 6 such that, for all (y 1 , y 2 , y 3 , ς) ∈ K, we have:

(3.51) 0 < c 6 ≤ ∇ y Γ e √ H e (y 1 , y 2 , ς), y 2 ≤ C 6 < +∞ .
Proof. The domain K is compact. It suffices to show that the function under consideration never cancels. First, remark that:

∂ y 1 Γ e (h, y 2 ) = 1 2 ∂ h Γ e (h, y 2 ) H e (y 1 , y 2 , ς) -1/2 d 1 (y 1 , y 2 ) -1 C (y 1 , y 2
) . Taking into account Assumptions 2.2 and 3.2, this expression can be zero if and only if y 1 = 0. But then the derivative with respect to y 2 cannot also be zero. Indeed, the identity H e (0, y 2 , ς) = P 2 (ς) and (3.49b) give rise to:

∂ y 2 H e (0, y 2 , ς) = 0 , ∂ y 2 Γ e H e (0, y 2 , ς), y 2 = ∂ y 2 Γ e (h, y 2 ) = 0 .
3.2.2. Absence of boundary critical points for sufficiently large times. From now on, we suppose that Assumption 3.2 is verified. Knowing that ξ ∈ A, the relation (3.42b) written at time s = 0 or for q = 0 yields: ∇ y Ψ(±, q; t, x; 0, y, ς, ξ) = ∇ y Ψ(±, 0; t, x; s, y, ς, ξ) = ξ ≥ r > 0 .

Therefore C 0 b (±, q; t, x) = ∅ and C(±, 0; t, x) = ∅ for all (t, x). This preliminary result can be further elaborated. Lemma 3.3. [exclusion of boundary critical points] Given q ∈ Z * , define:

(3.52) t m (q, µ) := µ -1 |q| -1 C -1 1 r , t M (q, µ) := c -1 6 (µ |q| c(w) -1 R + C 2 .
For all t > t M (q, µ), we find that C b (±, q; t, x) = ∅ and that:

C i (±, q; t, x) ∩ (s, y, ς, ξ) ; s ∈ [0, t m (q, µ)[ = ∅ , (3.53a) C i (±, q; t, x) ∩ (s, y, ς, ξ) ; s ∈ ]t M (q, µ), t[ = ∅ . (3.53b) Proof. The part (3.53a) is a simple consequence of: ∇ y Ψ(±, q; t, x; s, y, ς, ξ) ≥ r -µ |q| C 1 s > 0 , ∀ s ∈ [0, t m (q, µ)[ .
To extract (3.53b), exploit first (3.47) to interpret (3.42b) according to:

∇ y Γ e (h, y 2 ) s = µ q c(w) -1 ξ - s 0 ∇ y (b e • Ȳ0 ) ⋆ (0, y 1 , y 2 , ς, r) dr , q = 0 .
Then, use a combination of (3.34), (3.45b) and (3.51).

For t > t M (q, µ), the study of C(±, q; t, x) is reduced to the one of C i (±, q; t, x). This information will be crucial in the next paragraph 3.2.3 since it allows to integrate by parts, separately with respect to s.

3.2.3.

Sorting of the harmonics. As usual when dealing with dispersive equations, not all harmonics q are relevant.

Lemma 3.4. [selection of the propagated harmonics] Fix any time t with t > t M (1, µ). The condition C(±, q; t, x) = ∅ can be realized only if:

(3.54) q ∈ H := q ∈ Z * ; |q| ≤ λ M b -1 m c(w) -1 , Card H < +∞ .
Proof. Fix any q ∈ Z. Remark that t M (1, µ) ≥ t M (q, µ) for all q ∈ Z * . From Lemma 3.3, for t > t M (1, µ), we have that:

C t b (±, q; t, x) = ∅ , C(±, q; t, x) ≡ C i (±, q; t, x) , ∀ q ∈ Z .
In this context, the function a(•) may as well be assumed to be compactly supported in the interval ]0, t[. Then, to make a preliminary selection of the critical points, we can look independently at (3.42a). Exploit (3.44) to get:

|∂ s Ψ(±, q; t, x; s, y 1 , y 2 , ς, ξ)| ≥ |q| c(w) b m (y 1 , y 2 , ς) + inf r≤|ξ|≤R λ(|ξ|) if ± q ∈ Z -, |q| c(w) b m (y 1 , y 2 , ς) -λ M if ± q ∈ N * .
The second lower bound invites to compare c(w) b m with λ M . Interesting situations are sure to appear because 1 < λ M b -1 m c(w) -1 ≃ 22, see (3.24). It should be emphasized here that the comparison of the two quantities c(w) b m and λ M makes sense from a physical point of view. On the one hand, the nondimensionalization of Paragraph 3.1.1 yields a function b e • Ȳ0 (0, y 1 , y 2 , ς, •) which is of size ≃ 1, with minimum and maximum values estimated in optimal manner as in (3.44). On the other hand, the threshold limit λ M , which was introduced in Assumption 3.1, is an experimental data that is also of size ≃ 1. The phase Ψ(±, q) for ± q ∈ Z -, or for ± q ∈ N * with |q| > λ M b -1 m c(w) -1 , has no stationary points, whereas the domain of integration is compact. As a result of nonstationary phase theorem, we find I(±, q) = O(ε ∞ ) in such cases. The quantity I(±, q) may be non-negligible modulo ε ∞ only if ± q ∈ N * and |q| ≤ λ M b -1 m c(w) -1 . Thus, the cardinal of H is finite. Moreover, H is a nonempty set if and only if -1 and +1 are in H . + such that, for all µ ∈ ]0, µ 0 [ and t > t M (1, µ), the condition (s, y, ς, ξ) ∈ C(±, q; t, x) implies that ς = π/2.

Proof. From Lemma 3.3, the discussion is reduced to times s such that s ≥ t m (q, µ). The condition ∂ ς Ψ = 0 is sure to be satisfied when ς = π/2. The value ς = π/2 corresponds to positions where w = 0, that is to mirror points. In this sense, when ς = π/2, we may say that the signals are generated from the two mirror points. It remains to show that the condition ς = π/2 is also necessary to obtain ∂ ς Ψ = 0. To this end, at the level of (3.55), consider the expression which is inside the brackets. It is the sum of two functions. Due to (3.49a), the first function is linear in s with a non zero slope. In view of (3.45c), the second function is bounded for all s ∈ R + . Concretely: λ -1 M , the right-hand term of (3.56), the minimizer, is sure to be not zero.

The level energy h ∈ R * + and the altitude y 2 ∈ R * -being fixed, the periodic function b e • Ȳ0 (0, y 1 , y 2 , ς, •) oscillates between a minimal value b m (h, y 2 ) and a maximal value b M (h, y 2 ). The following analysis depends on how c(w) b M (h, y 2 ) is positioned with respect to λ M . Given q with |q| c(w) b m (h, y 2 ) < λ M , two situations can be distinguished on the basis of the following chief characteristics:

3a. The case |q| c(w) b M (h, y 2 ) < λ M . The whistler episode is spread over a finite period of time t, still large for small values of µ. Namely, it is of size c µ -1 with a fixed constant c not depending on R ; 3b. The case |q| c(w) b M (h, y 2 ) > λ M . The whistler episode is of size c(R) µ -1 with a constant c(R) growing to infinity when the control R, on the norms |ξ| of frequencies, tends to +∞. From this perspective, it can be extended beyond any time t ∈ R * + . But, in practice, the size of |ξ| is limited. 3.2.5. Harmonics q implying variations below the threshold limit. This is the case 3a: the domain of definition of λ For s ≤ tM (q, µ), the relation (3.42b) and the upper bound (3.45a) give rise to: |ξ| ≤ R M (q, µ) := c -1 6 C 1 λ -1 (|q| c(w) b M ) + µ c -1 6 C 1 C 2 |q| c(w) . Since q ∈ H and µ ≤ µ 0 , a uniform bound on |ξ| is available:

|ξ| ≤ R M := c -1 6 C 1 λ -1 (λ M b -1 m b M ) + µ 0 C 2 λ M b -1 m .
From now on, we fix R with R > R M . This condition implies that there is no restriction coming from (3.42b) when working with (3.42a) on the time interval [0, tM (q, µ)]. Indeed, since the control on |ξ| is not too restrictive, the count of all possible values of ξ shall be exhaustive. The direction ξ can be freely adjusted in order to recover (3.42b). Then, we can substitute ξ in (3.42a) with the value extracted from (3.42b), and concentrate on the identity thus obtained.

Observe that the function Rsc(•) is bounded. Indeed, the triangle inequality states that: First of all, there is a need to explain why N y µ is finite. Since Osc(•) is periodic, and therefore bounded, from (3.63), we can easily deduce that S y µ is compact. Look at the square of (3.57). All expressions involved are real analytic: the functions b e and λ by definition, and the function Ȳ0 as a consequence of the Cauchy-Kovalevskaya Theorem (for ODEs). The zeros of a nonconstant analytic function are isolated. This is sufficient to conclude that N y µ < +∞. To help the understanding of what follows, please refer to the following picture 1. As a consequence of (3.79), two successive emission points are asymptotically separated by a uniform gap, in the sense that: Display of the solutions to F (s, µ) = 0.

The Informal Statement 1 expresses in summary form the contents of Propositions 3.1 and 3.2. It is put forward because it highlights a marker of whistler phenomena: the organization of spectrograms into light-coloured vertical lines. Now, consider two successive bands of light, recorded at some fixed position after an interval of about one second. The preceding analysis does not say that the wave packets on the right come from the propagation of those on the left. The plasma waves satisfy the Maxwell equations (3.20). As such, they are presumably observed only once (and it is over). Instead, the present study suggests that the wave packets on the right are issued from a (repeated) macroscopic kinetic reorganisation of particles (through the dynamics forced by Vlasov equations), and then from a mechanism of wave-particle interaction (see Subsection 3.1.5).

From this point of view, we can consider that all the wave packets found in spectrograms have been generated independently. Thus, as stated before, the matter of Section 3.1 is the creation of light inside plasmas.

The previous discussion leaves many questions unanswered. For instance, it says nothing about the spatial positions from which the plasma waves are emitted. As a matter of fact, the condition (3.42d) is not exploited. In the same vein, there are no comments about the propagation of waves, through (3.20), after their creation. The next Subsection 3.2.8 is inserted to partly remedy this.

3.2.8. About the propagation of whistler waves. This short paragraph makes a few comments on the propagation of singularities in the context of our toy model (Assumption 3.1). It is inspired from the text [START_REF] Joly | Nonlinear oscillations beyond caustics[END_REF], which goes much further. The Lagrangian manifold associated to the phase Ψ(±, q; •) is the set of points:

L(±, q; •) := t, x, ∇ t,x Ψ(±, q; t, x; s, y, ς, ξ) ; (s, y, ς, ξ) ∈ C(±, q; t, x) .

As a by-product of (3.74), (3.83) and (3.84), using the implicit theorem, it can be stated that the times s(±, q; y) giving rise to emission points by applying Propositions 3.1 and 3.2 are smooth functions of (y 1 , y 2 ), not depending on y 3 . From (3.42b), extract: ξ(±, q; y) := µ q s(±,q;y) 0 ∇ y (b e • Ȳ0 )(0, y 1 , y 2 , π 2 , r) dr . The knowledge of s(±, q; y) and ξ(±, q; y) furnishes global parametrizations of L(±, q; •), by implementing the following mapping ı(±, q) where, for simplicity, the reference to (±, q) is skipped on the right-hand side: The space-time projections of the bicharacteristic curves t → ı(±, q; t, y) give rise to the rays of geometric optics, which are here straight lines travelling at the group velocity. To apply the standard version of the stationnary phase theorem to the integral I(±, q; ε, t, x) near some interior critical point, say (s, y, ς, ξ) ∈ C i (±, q; t, x), degenerate situations must be avoided. In practice, the computation of the second derivative of Ψ(±, q; t, x, •) with respect to (s, y, ς, ξ) is needed, in order to avoid the caustic set.

ı(±, q) : (t, y) -→     t x ∂ t Ψ ∇ x Ψ     =    
Of course, straight lines do not reflect the real situation [START_REF] Bortnik | The unexpected origin of plasmaspheric hiss from discrete chorus emissions[END_REF][START_REF] Yamaguchi | Ray tracing of whistler-mode chorus elements[END_REF]. A more accurate model for the dispersion relation λ(•) is still lacking to implement ray tracing methods. This is also a prerequisite to explain the transformation of sferics into whistlers, and then hisses.

3.2.9. A wide variety of research perspectives. There are many other interesting issues which are more or less connected to the actual discussion. Just to mention two main directions: 8a. Devices (tokamaks, stellarators, • • • ) used to confine plasmas. Although there are differences [START_REF]Topics in kinetic theory[END_REF] (hot plasmas, collision processes, stronger influence of the electric field, other toroidal geometries, • • • ), the scales are similar [START_REF] Bostan | The Vlasov-Maxwell system with strong initial magnetic field: guiding-center approximation[END_REF][START_REF] Frénod | Long time behavior of the two-dimensional Vlasov equation with a strong external magnetic field[END_REF]. The dynamical aspects are again forced by a strong external inhomogeneous magnetic field, and the general Assumption 2.1 allows already to include most basic models; 8b. Nonlinear interactions of the kinetic and fluid oscillations. This is clearly a matter of nonlinear geometric optics, following on from [START_REF] Métivier | The Mathematics of Nonlinear Optics[END_REF][START_REF] Rauch | Hyperbolic Partial Differential Equations and Geometric Optics[END_REF]. A long-term objective might be a mathematical study of magnetic reconnection [START_REF] Weia | Generation mechanism of the whistler-mode waves in the plasma sheet prior to magnetic reconnection[END_REF]. In conjunction with 8a and 8b, there is an intriguing parallel between Van Allen Belts/magnetic islands, whistler waves/Mirnov oscillations, and magnetic reconnection/disruption.

Proposition 2 . 1 .

 21 [periodic properties of the solutions to the first modulation equation] Select Zs1 i0 = (y, ς w , w, a w ) ∈ O with O as in (2.56). Define Zs2 i0 = (h z , ŷ, w, a w , ν z ) as indicated in (2.70). Mark the position ž := t (h z , ŷ, w) ∈ Ǒ. Introduce b z := ν z as well as:

• 2 . 3 . 3 .

 233 The induction. Use the equation (2.82) to derive:

3. 1 . 3 .

 13 The localization and the oscillating structure of the data. Now, we move on to another interpretation of what occurs. We start with discussing the problems raised by 2c and 2d. The first stage is to describe what can be selected at the time t = 0.

  w) b M ≃ 2 × 3 -4 , are clearly comparable, with apparently c(w) b m < λ M < c(w) b M . At this stage, we can consider that whistler waves correspond to data involving frequencies of size µ |ξ| ∼ ε -1 , and the above can be summed up as follows. Assumption 3.1. [a toy model describing the dispersion relation for whistler waves] For right-handed waves propagating parallel to the magnetic field, the propagation is governed by a law of the form (3.23), implying a real analytic strictly increasing function λ(•) satisfying: λ : R * + -→ ]0, λ M ] , λ ′ : R * + -→ R * + , (3.25) λ(0) = 0 < c(w) b m < λ M := lim z→+∞ λ(z) . (3.26)

  •) the flow induced by (1.5) in the case (E, B) = (0, 0), interpreted after the transformation Ξ 1 of Paragraph 2.2.3. Since u ε (τ, •) ≡ w for all τ ∈ [0, T ], the localization induced by the Dirac mass δ w (u) is maintained at the initial position u = w. Now, the amplitude w of the velocity being fixed, the vector w becomes a function w(ς, a) of (ς, a), as indicated in (2.35), and sometimes simply denoted by w in what follows. The formula (1.3), with f replaced by f 0 as in (3.17), gives rise to: (3.28) j ε (τ, y) = T 2 fε (τ, y, ς, a) dς da ,

2 . 6 .

 26 Two situations can be distinguished according to the respective positions of the graph of c(w) b e • Ȳ0 (•) and λ M . The case where the supremum of c(w) b e • Ȳ0 (•) is below the threshold limit λ M is first examined, in Paragraph 3.2.5. The case of overlapping values is discussed next, in Paragraph 3.2.6. Finally, the propagation aspects are encoded inside (3.42d). They are briefly alluded in Paragraph 3.2.8. 3.2.1. Preliminary estimates. Recall the formulas (2.77b) and (2.110) which indicate that the functions ν e and Γ e , and therefore b e • Ȳ0 as well as derived quantities like (b e • Ȳ0 ) ⋆ , • • • , depend on ς only through h. The substitution of h is the one of (2.110). To avoid multiple notations, we will not change notations when such factorizations are implied. With this convention, we can for instance write: (3.43) b e • Ȳ0 (0, y 1 , y 2 , ς, r) = b e • Ȳ0 √ H e (y 1 , y 2 , ς), y 1 , y 2 , r . More specifically, the expression b e • Ȳ0 (0, •) ≡ Γ e can be viewed as a function of (h, y 2 ) only, as revealed by (2.112). Fix (y 1 , y 2 , ς) ∈ R × R * -× T. The function b e • Ȳ0 (0, y 1 , y 2 , ς, •) is periodic. It undergoes variations between a minimal value b m (y 1 , y 2 , ς) and a maximal value b M (y 1 , y 2 , ς). This means that (2.8) can be replaced by: (3.44) b m ≤ b m (y 1 , y 2 , ς) ≤ b e • Ȳ0 (0, y 1 , y 2 , ς, s) ≤ b M (y 1 , y 2 , ς) ≤ b M , ∀ s ∈ R + ,

3. 2 . 4 .

 24 Reduction to a spatial phase. From the perspective of the Maxwell equations (3.20), the variables x and v do not intervene in the same way. The velocity component v appears when integrating to identify j. On the basis of Paragraph 3.1.7, it is involved only through the component ς of w. Lemma 3.5. [the signals are necessarily emanating from the two mirror points] There is µ 0 ∈ R *

0 ∂

 0 Now, look at (3.42c) with ς in a compact set of ]0, π[. The formulas (2.77b) and (2.110) indicate that b e • Ȳ0 depends on ς only through h. Recall the convention (3.43). We have: (3.55) ∂ ς Ψ = -∂ h Γ e s + s h (b e • Ȳ0 ) ⋆ dr H -1/2 e c(w) 2 -1 cotan ς .

( 3 . 1 C - 1 3

 311 56)(∂ h Γ e s + s 0 ∂ h (b e • Ȳ0 ) ⋆ dr ≥ c 4 t m (q, µ) -C 3 ,∀ s > t m (q, µ) .For µ ∈ ]0, µ 0 [ with µ 0 = r c(w) b m C 4 C -1

6 µ

 6 -1 (•) contains [|q| c(w) b m , |q| c(w) b M ]. Select t > t M(1, µ). With Lemmas 3.3 and 3.5, the condition (s, y, ς, ξ) ∈ C(±, q) appears to be equivalent to ς = π/2, together with (3.42a) and(3.42b). On the one hand, from (3.42a), we can deduce that |ξ| < λ -1 (|q| c(w) b M ). On the other hand, from the equation (3.42b), together with (3.45b), (3.47) and (3.51), it is easy to infer that |ξ| ≥ µ |q| c(w) (c 6 s -C 2 ). It follows that:(s, y, ς, ξ) ∈ C(±, q) , ∀ s > tM (q, µ) := c -1 |q| c(w) -1 λ -1 (|q| c(w) b M ) + C 2 .

∇

  |Rsc(s)| ≤ s 0 ∇ y (b e • Ȳ0 ) ⋆ (0, y, π/2, r) dr . Now, the right hand term involves the integral of a (T /2)-periodic function, whose mean value is zero. It follows that: y (b e • Ȳ0 ) ⋆ (0, y, π/2, r) dr ≤ C 2 (h, y 2 ) , ∀ s ∈ R + , (3.63) lim sup s→±∞ |Rsc1(s)| = lim sup s→±∞ |η s + Rsc(s)| = +∞ .

Figure 1 .-1 2 .

 12 Figure 1. Below the threshold limit. Display of the solutions to F (s, µ) = 0. Now, we want to get a more accurate count of the elements contained in S y µ . First, come back to the relation (3.57). It can be realized only if s ≥ sm := O m C -1 1 µ -1 c(w) -1 . Define: µ 0 := min r c(w) b m C 4 C -1 1 C -1 3 λ -1 M ; c 6 O m c(w) -1 C -1 1 C -1 2 .For µ ≤ µ 0 and s as above, the part Rsc1(•) of Rsc(•) is sure to be non zero, with some positive lower bound that is uniform for µ near zero. It follows that the function Rsc1(•), and therefore Rsc(•), is analytic for s ≥ s m . The function Rsc1(•) has a particular form. It is the norm of a linear function plus a periodic function. Due to this specificity, the derivatives (at all orders) of Rsc1(•) are bounded.

(3.80) s 4 n -s 1 n

 41 = š4 -š1 + o(1) , s 1 n+1s 4 n = T /2 + š1 -š4 + o(1) .Proof. The function λ -1 (•) is not defined on the whole image of c(w) b e • Ȳ0 (h, y 1 , y 2 , •).The formulation (3.61) cannot be used here. Still, the system (3.42a)-(3.42b) amounts to:(3.81) F (s, µ) = 0 , F (s, µ) := c(w) b e • Ȳ0 (h, y 1 , y 2 , s)λ µ η s + µ Rsc(s) ,

Figure 2 .

 2 Figure 2. Case of overlapping values.

  t y ∓ µ ts(y) λ ′ (|ξ(y)|) |ξ(y)| -1 ξ(y) ±λ(|ξ(y)|) µ -1 ξ(y)

  But it is not until long times τ ∼ 1 that these objects operate as indicated in (2.100) and (1.8). The ingredients ψ e s and ψ e pε are the cornerstones of understanding whistler waves. The statement below, proven in Paragraph 2.3.6, highlights some key results of Part 2.

	Theorem 1. [global and long time modeling of radiation belt dynamics] Assume that the
	electromagnetic field (

This can be done by induction (Paragraph 2.3.3). The justification of the WKB approximations (with a loss of precision by a factor ε -2 ) comes into play after, in Paragraph 2.3.4. In the end of Section 2.3, see Paragraph 2.3.5, special emphasis is given to the important notions of secondary phase ψ s and of principal phase ψ pε . Both phases ψ e s and ψ e pε emerge already during intermediate times t ∼ 1.

Table of contents 1 .

 of1 Introduction. 1.1. Long time gyro-kinetic equations. 1.2. Wave-particle interaction in cold magnetized plasmas. 2. Multiscale analysis of the Vlasov part. 2.1. The averaging method during intermediate times t ∼ 1. 2.1.1. Assumptions on the external magnetic field. 2.1.2. Straightening out the field lines. 2.1.3. Filtering out the main oscillations. 2.2. The approach to reach long times t ∼ ε -1 . 2.2.1. Preliminaries and guidelines. 2.2.2. Approximate rectification of the vector field A, homological equation. 2.2.3. The autonomous change of variables Ξ 1 (z). 2.2.4. The reduced hamiltonian. 2.2.5. The autonomous change of variables Ξ 2 (z), second piece of (2.33). 2.2.6. The time dependent change of variables Ξ 3 (t, z), third piece of (2.33). 2.3. The long time WKB calculus. 2.3.1. Initialization. 2.3.2. First and second modulation equations. 2.3.3. The induction. 2.3.4. Justification of the WKB approach, weak coupling. 2.3.5. Description of the phases in long times. 2.3.6. Proof of Theorem 1. 3. On the creation of light inside plasmas. 3.1. From Vlasov-Maxwell equations to oscillatory integrals. 3.1.1. Some weakly non linear model issued from space plasma physics.

	3.1.2. A brief foray into previous attempts.
	3.1.3. The localization and the oscillating structure of the data.
	3.1.4. Connection to geometric optics.
	3.1.5. Wave-particle interactions seen from the perspective of oscillations.
	3.1.6. Dispersion relation for whistler waves.
	3.1.7. A class of oscillatory integrals.
	3.2. Stationary phase expansions.
	3.2.1. Preliminary estimates.
	3.2.2. Absence of boundary critical points for sufficiently large times.
	3.2.3. Sorting of the harmonics.
	3.2.4. Reduction to a spatial phase.
	3.2.5. Harmonics q implying variations below the threshold limit.
	3.2.6. Emission points in the case of overlapping values.
	3.2.7. Whistler-mode chorus emissions.
	3.2.8. About the propagation of whistler waves.
	3.2.9. A wide variety of research perspectives.
	References

  , the function H is related to the modulation equation(2.31). Lemma 2.7 below says that it governs what happens at leading order during intermediate times t ∼ 1. This is why H is called the reduced hamiltonian.Lemma 2.7. Fix ŷ ∈ Ô. The expression H(•, y 2 , y 3 , •, w, a) plays the part of an energy function for the planar motion (2.43), meaning that:

  1 , ς w ) on C e (h), the solution to(2.43) follows the curve C e (h) with a non zero derivative. It turns in a direct sense around the point (0, π/2). The function Z1 Lemma 2.12. [periodic properties and special features of the solutions to the geomagnetic first modulation equation] Let h ∈ R * + and y 1 m ∈ R * + such that H e (y 1 m , y 2 , π/2) = h 2 . Consider the solution ( Z1 0 , Z4 0 ) to (2.43), issued at time t = 0 from the position (y 1 m , π/2). The function Z1 0 (•) is periodic with the period T given in (2.71). It is strictly decreasing on the interval [0, T /2] with Z1 0 (T /4) = 0 and Z1 0

0 (•) is given by the first component in (2.72). It decreases up to -y 1 m if ς w ≥ π/2. It increases up to y 1 m if ς w ≤ π/2. Both positions -y 1 m and y 1 m are turning points. The informations given below explain what happens next, starting from (y 1 m , π/2).

  oscillates around the position y 1 0 (ŷ) = 0, which is in the equatorial plane Eq. • 2.3. The long time WKB calculus. From now on, we work with z = t (h, y 2 , y 3 , u, a, b). We use (2.29) with Ξ as in (2.33) and Ξ as in (2.32) to transform (2.22) into (2.26). As a consequence of Section 2.2, we can assert that A 0 ≡ 0. Knowing that A 0 ≡ 0, we can turn to the study of (2.26). What follows is inspired by multiscale methods, and more specifically by Wentzel-Kramers-Brillouin (WKB) techniques [6, 22, 27]. We recall the definitions (2.23) and (2.24) of the operators • and •. 2.3.1. Initialization. Given N ∈ N * , we seek the solution (z ε , ϕ ε ) to (2.26) in the form:

  Implicit function Theorem applied to (2.104) allows to uniquely determine the function Φε (τ, t, θ). It says that Φε is smooth in (ε, τ, t, θ) ∈ [0, 1] × [0, T ] × T T × T. Involved in the context of the relation (2.103), with the estimate (2.98 ) on r ϕ ε in mind, it also guarantees that Φc

  Proof of Theorem 1. Theorem 1 is a compilation of results obtained in Section 2. It makes reference to the domains Bēlt and EBēlt, which come from interpretations of (2.10) and (2.59) in terms of the variables (x, v).The Assertion 1) was discussed at the very end of Paragraph 2.2.4, see the Mathematical interpretation 2. How the set Bēlt can be adjusted to fit with concrete observations is explained in Paragraph 2.1.1, just before Lemma 2.1.The Assertion 2) brings information about the compactness and the stability of EBēlt. Since the map y(•) is a diffeomorphism and O • y(x) is an orthogonal matrix, it suffices to argue at the level of EBelt. As stated in Lemma 2.7, the set EBelt is preserved for all intermediate times t ∈ R by the flow associated to the first modulation equation (2.31), that is under the action of Ξ(t, •) for all t ∈ R. In view of (2.86), it is also invariant under the flow induced by the second modulation equation (2.84), that is under the action of Z0 (τ, •) for all τ ∈ [0, T ]. Considering (2.85), nothing is changing at this level.

	2.3.6. More precisely:
	(2.113)	Bēlt := x ; y(x) ∈ Belt ,	EBēlt := (x, v) ; y(x), ς, |v|, a ∈ EBelt .
			.	•

  . It is used below as a starting point for establishing a toy model describing in the context of (3.20) the dispersion relation for whistler waves. What is missing to fit with (3.20) is the dimensional analysis of (3.21). The speed of propagation is 1 at the level of (3.10). Thus, in dimensionless variables, the frequency c |k| ∼ 10 4 must be replaced by |ξ| ∼ ε -1 . This corresponds to a semiclassical regime. Introduce:

	(3.22)	ω pe :=	4 π n e q 2 e m e	1/2	,	µ :=	ω ce ω pe	.

Acknowledgments. This work received the support of the Agence Nationale de la Recherche: projet NOSEVOL (ANR 2011 BS01019 01) and ANR blanc DYFICOLTI.

Additional information may be obtained by applying to the expression b e • Ȳ0 (0, y 1 , y 2 , ς, •) the decomposition of Lemma 2.3. With h and Γ e as in (2.110) and (2.112), this leads to:

(3.47) s 0 (b e • Ȳ0 )(0, y 1 , y 2 , ς, r) dr = Γ e (h, y 2 ) s + s 0 (b e • Ȳ0 ) ⋆ (0, y 1 , y 2 , ς, r) dr .

In view of conditions (3.17) and (3.16), the part Γ e (h, y 2 ) of (3.47) only involves specific values of h, those which are related to the identification of the extended Van Allen belt. Since the function Γ e (•) does not depend on y 3 , it suffices to consider the projection of the compact K f , given by (3.16), onto the two first variables h and y 2 , that is:

(3.48) P r(K f ) := (h, y 2 ) ; 0 ≤ h ≤ f ey 2 h e (y 2 ) , y 2 ∈ supp h e (•) .

Assumption 3.2. [the mean value Γ e ≡ A (Γ e )(h, y 2 ) of Γ e is without critical points] There are positive constants c i ∈ R * + and C i ∈ R * + with j ∈ {4, 5} such that:

∀ (h, y 2 ) ∈ P r(K f ) , (3.49a)

∀ (h, y 2 ) ∈ P r(K f ) . (3.49b) Assumption 3.2 is realistic. Indeed: Lemma 3.1. [about the relevance of Assumption 3.2] For energy levels h ∈ R * + which are small enough, Assumption 3.2 is sure to be verified.

Proof. Since w is fixed, we can omit to mention w when it is not necessary. In the integral defining Γ e at the level of line (2.112), change r into h -1 r to get:

.

There are two ways of looking at (3.50). The first is to exploit the information (2.47) together with the definitions of f y 2 (•) and g(•) to derive:

This allows to justify the convergence of the integral (3.50), and to see that Γ e (0, y 2 ) = 0. The other view is to recall Proposition 2.1, where it was noted that the operation A (•) ≡ • of (2.74) amounts to integrate with respect to a probability density, so that:

The reason explaining the O h 4 instead of a O h 3 is that Γ e (•, y 2 ) is an even function. From this and (2.54), we can deduce that:

In view of the construction leading to f and g, see for instance (2.54), we have:

The situation 3a, if it happens, includes at least the value q = 1 with the sign +. Now, the discussion being similar for all other related values of q, it can be illustrated by investigating only the case of I ≡ I

We call S y µ the set of all solutions to (3.57). On the other hand, the cardinal of S y µ is denoted by N y µ := Card (S y µ ). Recall that the overdense condition corresponds physically to the selection of small parameters µ, with an order of magnitude of size ∼ 10 -2 . Proposition 3.1. [asymptotic time distribution of the emission points under the overdense condition; case where c(w) b M is below the threshold limit] The number N y µ of solutions to the equation (3.57) is finite. The elements s j µ of S y µ can always be ranked in order of increase, with 1 ≤ j ≤ N y µ . There is a constant

. This means that Ȳ0 (0, y, π/2, •) starts at initial time s = 0 from an extremal position ±y 1 m . Lemma 2.12 can therefore be applied. Since b e (•, y 2 ) is an even function, it follows that the function Osc(•) is periodic of period T /2. More precisely, it is strictly decreasing on the interval [0, T /4] from a maximum value O M ≤ λ -1 c(w) b M when s = 0, to a minimum positive value O m ≥ λ -1 c(w) b m that is reached when s = T /4. It is strictly increasing on the interval [T /4, T /2] from O m to O M , and so on. Briefly, we can assert that:

where, taking into account (3.51), we have to deal with the (strictly) positive number:

, and the auxiliary function Rsc(s) := Rsc1(s) + Rsc2(s) with:

The preceding property, together with (3.62), allows to introduce:

A preliminary remark is that the condition s ∈ S y µ ∩ [n T /2, (n + 1) T /2] can be realized only if the integer n is such that:

In view of (3.60), Osc ′ (s) = 0 if and only if Osc(s) = O m or Osc(s) = O M . After shifting the value of Osc(s) a little away from O m and O M , we can get a quantified version of this property. In fact, for C 9 ∈ R * + and µ 0 ∈ R * + respectively sufficiently large and small, we have for all µ ∈ ]0, µ 0 ] that:

Given n ∈ N, consider the condition:

The set S y µ can be decomposed into the disjoint union S y µ = Sy µ ∪ ( Sy µ ) c , where:

Recall that:

(∂ s F )(s, µ) = Osc ′ (s)µ ηµ Rsc ′ (s) , Osc ′ (s + T /2) = Osc ′ (s) . For µ = 0, the solutions of the equation (∂ s F )(s, 0) = Osc ′ (s) = 0 are of the form k T /4 with k ∈ Z, see (3.60b). In view of (3.60c), for s ≥ s m , we can apply the implicit function theorem to solve (∂ s F )(s, µ) = 0 for µ small. For all k ∈ Z, there is exactly one solution emerging locally from a position k T /4 ≥ s m . Using the periodicity of Osc ′ (•) and (3.64), we can check that, at least for µ ∈ ]0, µ 0 ] with µ 0 small enough (restrict µ 0 if necessary), there are no other solutions than the preceding ones. It follows that:

Combining (3.65), (3.67) and (3.68), it is easy to infer that:

Define the times t j µ with j ∈ {1, 2, 3, 4} in such a way that:

Consider also the translated versions of these times:

Select n as in (3.67). By construction, we have:

. For such n, observe that:

Secondly, deduce from (3.66) and (3.72) that:

On one side, we have (3.69) with a bound not depending on µ ∈ ]0, µ 0 ]. On the other side, from (3.73) and (3.74), it follows that, for all integer n satisfying (3.67), we have:

Let us now take a look at the second part of (3.59). Fix any constant

If necessary, diminish µ 0 to be sure that:

. Replace (3.67) with the more restrictive condition:

Define (see Figure 1) the times sj with j ∈ {1, 2, 3, 4} in such a way that:

The set of integers n satisfying (3.75) yields at least C µ -1 solutions in S y µ . Moreover, any of these solutions s j µ is such that:

Taking into account the variations of the function Osc(•), see (3.60d), this leads to:

, and therefore:

C) := min ( s3 -s2 ; s1 + (T /2) -s4 ) .

3.2.6. Emission points in the case of overlapping values. This is the case 3b, as is mentioned at the end of Subsection 3.2.4. Again, the related discussion is similar for all corresponding values of q. As before, it can be illustrated by investigating the case of I ≡ I 1 + , with

Choose µ 0 ∈ R * + to recover all the aforementioned conditions, including (3.78) as well as µ 0 ≤ (R 0 M ) -1 λ -1 (λ Mδ). Let (s, ξ) be a position that is adjusted as indicated in (3.42b). Define c 10 := 2 c -1 6 λ -1 (λ Mδ). We find that:

] . Select n as in Proposition 3.2. By construction:

Thus, there is at least one solution s