DISPERSION RELATIONS IN COLD MAGNETIZED PLASMAS

. Starting from kinetic models of cold magnetized collisionless plasmas, we provide a complete description of the characteristic variety sustaining electromagnetic wave propagation. As in [4, 13, 17], the analysis is based on some asymptotic calculus exploiting the presence at the level of dimensionless relativistic Vlasov-Maxwell equations of a large parameter: the electron gyrofrequency. The method is inspired from geometric optics [29, 33]. It allows to unify all the preceding results [8, 12, 38, 31, 37, 40], while incorporating new aspects. Speciﬁcally, the non trivial eﬀects [5, 9, 10, 24] of the spatial variations of the background density, temperature and magnetic ﬁeld are exhibited. In this way, a comprehensive overview of the dispersion relations becomes avail- able, with important possible applications in plasma physics [7, 28, 30].

1. Introduction. The Appleton-Hartree equation, sometimes referred to as the Appleton-Lassen equation, is a mathematical expression that describes the refractive index for electromagnetic wave propagation in a cold magnetized plasma. Historically [15], it was developed in the context of the magneto-ionic theory. The form that is still widely exploited in plasma physics [12,30,37] is valid in the case of a uniform magnetic field. It resulted from the works of E. Appleton [1] and D.
Despite the abundance of literature on the subject, more refined models are still the focus of ongoing research. It is because numerous factors may be considered: the framework may be relativistic [8,40]; it may look at the oblique propagation [31,38]; it can select a special frequency range (to point out for instance whistler waves [7,36]); it may incorporate the influence of boundaries [2,10]; and, last but not least, it can take into account the presence of spatial inhomogeneities. The aim of the present article is to discuss this issue through a new comprehensive approach that encompasses all the foregoing elements, while adding significant additional information concerning the last one.
A complete mathematical formulation of the problem is available by coming back to the original relativistic Vlasov-Maxwell (RVM) system. The general setting is as in [16], or see directly at (8)- (9). It allows to take a principle-based approach. One challenge is to fix adequately the physical framework. This must be done with a close link to concrete situations and in accordance with the basic concepts of plasma physics. As a prototype of a cold magnetized plasma, we will consider the earth's magnetosphere. This choice allows access to exhaustive, reliable and verifiable data. Simplifications come from three basic requirements. As stated in Assumption 1, the temperatures and the densities of the species are supposed to undergo only slight variations. As prescribed in Assumption 2, partially or fully ionized space plasmas have almost the same number of positive (ions) and negative (electrons) charges, and therefore behave quasineutral [18,19]. Moreover, as indicated in Assumption 3, most particles are in a state of local cold thermodynamic equilibrium [8,31]. There remains, however, a limited but very significant fraction of the plasma which is out of equilibrium. The related nonthermal facets play a central role in many applications. They are usually tackled through some stability analysis [26].
Another key aspect of astrophysical plasmas is the decisive intervention of some exterior magnetic fieldB e . From this point of view, the situation is very similar to that of magnetic confinement fusion [4,13,14]. The charged particles are deflected byB e through the Lorentz force. The strong influence ofB e leads to the introduction of a small dimensionless parameter ε, with 0 < ε 1, defined at the level of (49), and coming from the inverse of the electron cyclotron frequency. In addition, the geomagnetic fieldB e is inhomogeneous. More precisely, it can be approximated by a dipole model. With position x in some open set Ω ⊂ R 3 , it takes the form B e ≡ ε −1 B e (x), where the function B e (·) is subjected to Assumptions 4 and 5. As shown in (41), it follows that the Vlasov equation contains the penalization term ε −1 p × B e (x) · ∇ p , where p ∈ R 3 is a momentum.
A work of preparation allows to extract the dimensionless version of the RVM equations. Then, the parameter ε can serve as a unit of measurement to which all the dimensionless quantities can be compared. By this way, we are led to a singular version of the RVM system. The corresponding asymptotic analysis is part of a long tradition of mathematical works in gyrokinetics [4,13], fast rotating fluids [6] and geometric optics [29,33]. The uniform case, that is when the function B e (·) is constant, has been extensively studied. However, the realistic frameworks [5,7] imply non-constant functions B e (·). Now, the variations of B e (·) raise numerous difficulties. Their effects have not been investigated so far; they remain poorly understood; and they do not fall under the scope of classical results. The complications stem from the penalization term ε −1 p × B e (x) · ∇ p that involves a large skew-adjoint differential operator with variable coefficients.
The present article is starting to address this problem. It shows that the variations of B e (·), both in amplitude and directions, play a crucial role in electromagnetic wave propagation. In a logical order, they impact the dispersion relations, the eikonal equations, and finally the ray tracing methods [25,41]. They have therefore important practical implications. To make our main results readable and usable, some notation is needed.
In the magnetosphere, signals split up into two characteristic components. In the case of oblique propagation ( = 0), this translates into a partition V(x, ) = V o (x, ) V x (x, ). The part V o (x, ) is associated to ordinary waves (o -waves), whereas the subset V x (x, ) is related to extraordinary waves (x-waves). O -waves and x-waves are commonly observed. A precise mathematical definition of V o (x, ) and V x (x, ) is given in the statement below.
The description of the characteristic variety V through (1) is directly applicable in a large number of situations. That is why it is highlighted in this introduction. However, this representation implies a special choice of parameterization, which is not appropriate in the case = 0. Precisely, one important contribution of the present article is also to provide a global and intrinsic view of V (Definition 3.4) that allows to study V in its different facets (Section 3.4). In particular, in Paragraph 3.4.2, we investigate what happens if we fix τ ∈ R * + , that is if we consider sets given by ξ ∈ R 3 ; (t, x, τ, ξ) ∈ V . By this way, according to the values of τ , we can recover various pictures. This could include one or two (more or less) nested spheres (Figures 12 and 13). This can also result in one or two connected unbounded sets with conic shape (Figures 10 and 11). Another significant aspect of our study is the topological decomposition (valid only in the oblique case = 0) of V into connected components. From this perspective, the case of parallel propagation ( = 0), which is often presented as being indicative of the general situation, appears rather to be a very singular situation. As revealed in Paragraph 3.3.3, as can be seen by comparing Figure  7 ( = 0) and Figure 8 ( → 0), and as you can guess by looking at Figure  9, it is a composite of ordinary and extraordinary waves. The classical physical nomenclature (in terms of Alfvèn, whistler, ... and electrostatic waves) which is recalled in Paragraph 3.3.1 does not take into account this mixture. As a matter of fact, it is based on other considerations. Theorem 1.1 gives access to rigorously justified eikonal equations (Lemmas 3.30 and 3.32) governing the propagation of the phases φ. To this end, it suffices to replace in (101) the variables τ and ξ respectively by ∂ t φ(t, x) and ∇ x φ(t, x). The variations of B e (·) come into play through the dependence of and g ± (·) on B e (·), see (99) and Definition 3.7. As a by product, as shown in Section 3.4, purely parallel propagation cannot occur.
In the following text, special emphasis will be placed on frequencies which are in a range around or below, but comparable to the electron cyclotron frequency ε −1 . The reasons for this are the following. First, this range is where the exterior magnetic field B e (·) has the most influence by a clear separation between ordinary waves and extraordinary waves, and by the appearance of exactly two cutoff frequencies and two resonance frequencies. Secondly, as is well-known [22,28,40], this is where the propagation can be responsible for the energisation of confined plasmas and particle loss [38] through a mechanism of wave-particle interaction [7]. Now, when dealing with the plasmasphere, looking at such frequencies means to focus on Very Low Frequency waves (VLF radio frequencies in the range of 100 Hz to 10 kHz). Experimentally (Figure 17), there are a lower band and a upper band (corresponding to the two resonance frequencies) where VLF emissions arise. The monochromatic elements forming the fine structure of chorus have been mathematically interpreted in [7] as coming from a mesoscopic caustic effect. It was done on the basis of the classical toy model of [12,30,37], which is derived from parallel propagation ( = 0) in the presence of a uniform magnetic field. Beyond this preliminary approach, to fully understand the morphological properties of chorus emissions, Theorem 1.1 is required. As outlined in Paragraph 3.4.5, a lot of information about chorus emissions [27,28,35,41]  In Paragraph 2.4, we present such a framework issued from near-Earth space plasmas. This application involves a strong external magnetic field, satisfying special assumptions (Paragraph 2.5). In Paragraph 2.6, we derive the dimensionless form of RVM equations. By this way, we get a basic model (Paragraph 2.7) which can serve for the description of electromagnetic waves (Paragraph 2.8). This model corresponds to a linear perturbative regime. It is built with the system of equation (53)-(54)-(55) involving the small parameter ε ≡ |ε 1 | 1.
2.1. Relativistic Vlasov-Maxwell equations. The topic of RVM equations has been widely discussed [4,7,13,16,26]. The corresponding framework is recalled hereafter. The speed of light is c 0 2, 99×10 8 m s −1 . Let L ∈ R * + be a characteristic spatial length. The original spatial variable isx ∈Ω, whereΩ is some non-empty open set of R 3 . We fix the observation time T ∈ R * + as T := L/c 0 . The original time variable ist ∈ [0, T ]. There are corresponding rescaled versions: The original space and momentum variables are (x,p) with: We consider a plasma which is confined insideΩ, and which consists of N distinct species labelled by α ∈ {1, · · · , N }. The particles of the α th species have charge e α and rest mass m α . The number α = 1 will be associated with electrons. Thus, the elementary charge is e ≡ −e 1 1, 6 × 10 −19 C and the electron rest mass is m e ≡ m 1 9, 1 × 10 −31 kg. Recall that the proton-to-electron mass ratio β 1836 is a dimensionless quantity, so that: On the other hand, the charge e α is an integer multiple of e. More precisely: The velocityṽ α of a particle of type α is limited by |ṽ α | ≤ c 0 , and it is linked to the momentump ∈ R 3 through: The kinetic distribution function (KDF) of the α th species is denoted byf k α (t,x,p). It is composed of: x,p). The density ratio ν ∈ R * + between these two populations is assumed to be small and independant of α: The charge densityρ and the current density are given by: We impose a (stationary, divergence and curl-free) external magnetic fieldB e : Ω −→ R 3 . We also take into account some collective self-consistent electromagnetic field (Ẽ,B)(t,x), which is created by all plasma particles. Then, neglecting the collisional effects, the time evolution of the KDF can be modelled through the Vlasov equation: On the other hand, the self-consistent electromagnetic field (Ẽ,B)(t,x) is subjected to the Maxwell equations: In (9), the physical constant 0 8, 8 × 10 −12 F m −1 stands for the vacuum permitivity.

2.2.
A dominant stationary part in a state of local thermodynamic equilibrium. The unknowns in (8)-(9) are thef k α (·) and (Ẽ,B)(·). For ν = 0 and (Ẽ,B) = (0, 0), we simply find: This expression is assumed to satisfy (8) modulo some small term (to be specified later): It is also an equilibrium from the viewpoint of waves, meaning that: The aim of this subsection 2.2 is to give a detailed description of special functions f d α (·) satisfying (11) and (12), together with a number of other relevant physical constraints. As a matter of fact, the two paragraphs 2.2.1 and 2.2.2 will consider (12b). On the other hand, the paragraph 2.2.3 will give a precise description of f d α (·), complemented by the examination of (11) and the verification of (12a Denote simply by θ d α ∈ R * + and n d α ∈ R * + typical sizes ofθ d α (·) andñ d α (·). We require that the two quantitiesθ d α (x) andñ d α (x) do not deviate too far from θ d α and n d α . In other words:

Assumption 1. [possible but slight variations in temperatures and densities]
There is a constant c ∈ ]0, 1[ such that: Recall that k B 1, 38×10 −23 m 2 kg s −2 K −1 stands for the Boltzmann constant, and also retain the relationship 1 eV 1, 16 × 10 4 k B K. The electron temperature (T e ≡ T 1 ) and the ion temperatures (denoted by T α for α > 1) can be expressed either in kelvin (K) or in electronvolt (eV ). Because of the large difference in mass, the electrons will come to thermodynamic equilibrium amongst themselves much faster than they will come into equilibrium with the ions or neutral atoms. For this reason, the ion temperatures may be very different from (usually much lower than) the electron temperature: Based on the relative temperatures of the electrons, ions and neutrals, plasmas are classified as thermal or non-thermal. Introduce the thermal speed: It is linked to the dimensionless parameter θ d α through: These two quantities v th α and θ d α can also be viewed as measures of temperature. Combining (3) and (15), we get: As a rule of thumb, temperatures T α well below 100 eV are said cold ; those which are about 100 eV (θ d α 10 −2 ) are considered warm ; those with T α ranging from 100 eV to 10 keV (10 −2 θ d α 1) become progressively hot ; particles with higher energies (1 θ d α ) are termed energetic or relativistic.

Quasi-neutrality.
A plasma consists of approximately equal numbers of positively charged ions and negatively charged electrons. This property is expressed by (12b). In view of (7a) and (13), this can also take the following equivalent form.

Assumption 2.
The plasma is quasi-neutral in the sense that: The interpretation of (19) is the existence of a background neutralizing ion population. In view of (4), (14) and (19), we can infer that: In the scientific literature [23,40], the term cold plasma is sometimes associated with a Dirac mass in the momentump, in which case θ = 0 or equivalently M(r) ≡ (4 π) −1 r −2 δ |r=0 . However, the temperatures are not zero. The use of (21) is therefore more refined, while keeping track of concentration effects through the (possible) smallness of θ.
Electrons are lighter than ions and neutral atoms. This is why they can easier reach higher energies. This is especially true in the case of the outer Van Allen belt as it will be explained in Paragraph 2.4. It follows that the dominant populations can inherit more complex structures than (21). As long as θ d α 10 −2 , we can keep (21) to describef d α (·). But beyond, that is for hot plasmas (10 −2 θ d α 1) and all the more so for relativistic beams (1 θ d α ), it might be preferable to select a Maxwell-Jüttner distribution function. Such relativistic aspects will be examined in a forthcoming publication. Here, we will stay in the context of (21). Assuming (21), we find that: If the temperature and the density are constant, that is ifθ d α (·) ≡ θ d α andñ d α (·) ≡ n d α , the right-hand term of (24) disappears. Then, the expression (21) yields exact solutions to the equation (11). On the contrary, when ∇x(lnθ d α ) ≡ 0 or when ∇x(lnñ d α ) ≡ 0, this term is sure to contribute (more or less strongly). As revealed by the dimensionless version of the equations (see Paragraph 2.6), the induced effect depends on the variations ofθ d α (·) andñ d α (·) inx, but also on the size of the temperature θ d α . The expressionf d α (·) issued from (21) involves even functions. It makes no contribution to the current density in (7b). When dealing with (21), the constraint (12a) is sure to be satisfied.

Plasma phenomena out of equilibrium.
Collisionless plasmas are characterized by a Knudsen number K n larger than one. It follows that the possible discrepencies from the Maxwell-Boltzmann distributions are not immediately relaxed. Expressions like (10), where thef d α (·) are adjusted as in (21), are therefore not sufficiently exhaustive to fully describe realistic plasmas. In particular, in connection with observed phenomena, the following aspect may need to be incorporated: -Spatial variations in the positionx : As already taken into account, most concrete situations involve non constant functionsñ d α (·) − n d α andθ d α (·) − θ d α . A perturbative approach is required to absorb the remainders obtained when computing (8), see the right-hand term of (24).
-Anisotropic effects in the momentump : There is a plethora of factors that produce inhomogeneities according to the different directions |p| −1p ∈ S 2 ofp. Examples include lighting strikes, solar flares, ohming heating, neutral beam injections, high-frequency and radio-frequency waves. Obviously, such features cannot be included inside (21), because this expression depends only on |p|. -Dynamical aspects : The anisotropy inp can be enhanced by the external magnetic fieldB e (·) through the Lorentz force. This results in microscopic timedependent instabilities, known as anomalous transport. -Electromagnetic perturbations : In practice, fluctuations of the electric field and of the magnetic field (nearB e ) are recorded. In contrast with (10), the selfconsistent field (Ẽ,B)(·) is in general non-trivial. To grasp plasma instabilities or plasma processes which are not in thermal equilibrium, we can perform a stability analysis in the vicinity of (10). This amounts to take ν ∈ R * + with ν 1 in (6). From now on, the focus will be on (8)- (9) in the context of the perturbative regime (6). With a dominant stationary part adjusted as in (21), we look at the extra partf s α (·), which is governed by: together with: 2.4. Some concrete situation. The description through (21) of the underlying medium can be further specified with regard to the applications in physics. We will consider the case of Near-Earth space plasmas [24], inside the exosphere (where K n > 1). In the inner magnetosphere, we can distinguish two types of particle structures, which can involve special values of temperatures and densities, falling into the context of the description of cold plasmas from Paragraph 2.2.
• C1. The plasmasphere is a doughnut-shaped region which rotates with the Earth. It is a cold plasma ( 1 eV) made of H + (nominally about 80%), He + (10 − 20%) and O + (a few percent). As noted in [9], the average electronic temperature is T e 6 × 10 3 K, so that θ d 1 10 −3 . As a result, the selection of (21) for all α ∈ {1, · · · , N } is adequate. The plasmaspheric concentration is around a mean value which is close to n d 1 = 10 8 electrons/m 3 . The functionsñ d α (·) may vary slightly from place to place, especially when crossing the slot region. Global realistic profiles for n d α (·) are suggested in [9,24]. • • C2. The outer Van Allen belt partly overlaps with the plasmasphere. It is a concentric, tyre-shaped belt containing ionic and neutral species which may be globally cold, warm, or possibly hot (but certainly not beyond). This belt can also implicate highly energetic ( 1 MeV) electron fluxes which may be trapped [5,7] by the geomagnetic field. These fast charged parcticles travel along the field lines. They stay approximately at a fixed L-shell, where the density is almost constant. We can therefore takeñ d α (·) ≡ n d α for all α ∈ {1, · · · , N }. In particular, as in C1., we can fix n d 1 = 10 8 electrons/m 3 . Moreover, as long as the fraction of hot and energetic particles remains small, even for α = 1, we can still select (21) and incorporate the relativistic aspects insidef • In what follows, the two above contexts will be systematically tested. This will appear throughout the text inside paragraphs that will be preceded by the title Discussion.
2.5. Inhomogeneous magnetized plasmas. Most cold plasmas [4,7,13,18] are under the influence of a strong external magnetic field which can be prescribed through some adequate functionB e :Ω −→ R 3 , with amplitudeb e (x) := |B e (x)|. The functionb e (·) is assumed to be of the order b e ∈ R * + . More precisely, we can find c ∈ ]0, 1[ such that: In view of (2), we can consider the following rescaled version ofB e (·): Then, the condition (27) becomes: The function B e (·) generates a unit vector field: Complete e 3 (x) into some right-handed orthonormal basis (e 1 , e 2 , e 3 )(x), with: With the preceding ingredients, we can define some orthogonal matrix O(x) and some constant skew-symmetric matrix Λ according to:

Assumption 5. [divergence and curl-free external magnetic field] The function
. It is such that: Generally, magnetic field lines are curved. In most concrete situations, the function B e (·) is a well-known non constant function of x ∈ R 3 . As depicted below, cylindrical coordinates (ρ, ϕ, z) ∈ R * + × R 2 can be used to mark the position of x. Recall that, for any vector field A = A ρ e ρ + A ϕ e ϕ + A z e z , the divergence and the curl are given by the dipole model:

Discussion 1. [external magnetic field] In the case of Near-Earth plasmas, just
where R e 6, 3 × 10 6 m is the Earth radius. The plasmasphere and the two Van Allen belts occupy some area between altitudes of 2 R e and 8 R e . We can work with: The magnitude of the geomagnetic field, as it can be measured at the surface of the earth, is about b e 10 −5 T . The Earth's magnetic field B e ≡ B e e does not depend on the angle ϕ. As mentioned for instance in [5,7], it can be approximated by: Observe that the two functions B e e (·) and b e e (·) are homogeneous in (ρ, z) of degree −3. We can recover (29) with c = 8 −3 . Applying (34a) and (34b), we see easily that both conditions (33a) and (33b) are met. • 2.6. Dimensionless equations. The aim here is to write the system (25)- (26) in dimensionless form (Paragraph 2.6.1). By doing so, we have to take into account the effects induced by the spatial variations ofñ d . Then, it is useful to straighten the field lines (Paragraph 2.6.2 ). Moreover, in connection with the application to the earth's context, it is important to give (Paragraph 2.6.3) a precise description of (the relative sizes of) the various physical parameters.

Rescalings. Let us recall (28), and define rescaled versions ofñ
This says nothing about the comparison of the selfconsistent magnetic fieldB to the typical size b e ofb e (·). From the Ampère's law in (25), we can infer thatB ν θ d α b e . With this in mind, we can further define new unknowns by the relations: From now on, the time-spatial position is (t, x), with (t, x) ∈ M := [0, 1] × Ω. Let T * M be the cotangent bundle associated with M . With (37a), the vectors v α and p α are linked by the relations issued from (5), that is: Among the fundamental plasma parameters, we can mention (for α = 1) the electron gyrofrequency (or cyclotron frequency) ω ce ≡ ω c1 and the electron plasma frequency (or plasma oscillation) ω pe ≡ ω p1 . For α ∈ {2, · · · , N }, we could cite the ion gyrofrequencies ω cα and the ion plasma frequencies ω pα . For simplicity of presentation, we define below these quantities with an algebraic sign: There are corresponding dimensionless coefficients ε α and µ α , given by: A quick calculation indicates that the two relations of (42b) are propagated by the evolution equation (41)-(42a). Therefore, it suffices to check (42b) at the time t = 0, and then to focuss on (42a). (41) is not yet in a suitable form. Still, we need to straighten out the field lines. Recall (30)-(32) so that:

Straightening the field lines. Equation
In view of Discussion 1, the directions of the unit vector field e 3 (·), and therefore of B e (·), can vary with changes in x ∈ Ω. To remedy this situation, we replace simultaneously B e , B, E and p α according to: For the sake of simplicity, the subscript α that identifies the different momentum variables p α will be omitted. Concerning p ≡ p α ∈ R 3 , we can pass from cartesian to spherical coordinates, with orthonormal basis (e r , e , e ω ). This gives rise to: From now on, the spatial-velocity position is marked by y := (x, , ω, r) ∈ Ω × T 2 × R + . We modify f α (·) to fit with the preceding adjustements: sin ω sin , cos ) . As usual, the symbol S refers to the Schwartz space. We consider functions f (·) satisfying uniformly in (t, x, , ω) ∈ M × T 2 the conditions: The gradient ∇ p is converted into the spherical gradient ∇ p , with: The change of variables (x, p) → (x, p) on the right of (45) induces some extra term when transforming (v · ∇ x )f accordingly. Some application Q(·) does appear. This is a vector valued quadratic form in p, namely: Put aside the integral operators: The hierarchy between dimensionless parameters. For further analysis, it is crucial to produce values for the parameters ε α , θ d α and µ α which could be meaningful from a physical point of view. It is also important to compare these quantities to one another. To this end, the following dimensionless number (which comes from the inverse of the electron cyclotron frequency): will serve as a unit of measure.

Discussion 2.
[about the size of ε] As indicated in (40), the number ε appears to be the ratio beween the reference frequency 1/T = c 0 /L and the gyrofrequency ω ce . This appears to be a small parameter. The plasmasphere begins above the upper ionosphere and extends outwards. It contains the inner Van Allen belt, which is located between 1 R e and 2 R e . Its outer boundary, known as the plasmapause, can vary with geomagnetic activity. During quiet period, it can expand outward to 7 R e or beyond. On the contrary, during magnetic storms, it moves at around 5 R e . The outer Van Allen belt covers altitudes of approximately 4 to 8 R e . In short, we can take the mean value L 5 R e , so that ε 10 −4 . • From now on, we take ε := 10 −4 1 as the small reference parameter to which all other quantities will be compared. For instance, with (3), keep in mind that:

Discussion 3.
[about the size of the coefficients θ d α ] The measurements recorded in the plasmasphere indicate that θ d 1 ε. The ratio θ d α /θ d 1 is given by (18). In view of (3) and (15), it is very small. The model (21) is therefore suitable for all α.
In the limit ε goes to zero, we can even say that the kinetic distribution functions are given by Dirac masses. This explains the cold plasma approximation, which is sometimes applied in geophysical research [23,40]. As regards the outer Van Allen belt, the situation is more problematic. It may still be considered that θ d α ε and that f d α (·) is as in (21) for all α ∈ {2, · · · , N }. However, the presence sometimes of a large amount of hot electrons might also point in favour of a Maxwell-Jüttner distribution. •

Discussion 4.
[about the size of the coefficients µ α ] In view of (40), the access to µ 1 requires to evaluate b e , m e and n d 1 . How to fix these quantities has already been explained. We find that µ := |µ 1 | 1. We remark that: We can also find that: In practice, the value of µ := |µ 1 | is of size 1. As indicated in Definition 3.21, it can be locally compared to b e (x) 1. The plasma is termed underdense when µ < b e (x) and overdense when b e (x) < µ. To track the influence of µ, this parameter will not be normalized in what follows. At all events, retain that the size of ε is always small, and far below µ. Discussion 5. [about the size of ν] We will adjust ν in such a way that ν ε. By this way, we can stay in a perturbative regime, even if θ d 1 1. Indeed: -Smallness ofB in comparison withB e : In view of (37c), the condition (6), the part νf s α (·) appears as a small modification off d α (·). This makes sense whatever the parameter θ d α is, small or large. Indeed, the amplitude off s α (·) as prescribed by (37b) with f α 1 is equivalent to the amplitude off d α (·) given by (21).
In concrete terms, we will impose ν ∼ ε.
2.7. The cold asymptotic regime. The discussion can be put in the context of some asymptotic analysis (when ε goes to 0). To this end, the coefficients ε α , θ d α and µ α must be viewed as functions of ε ∈ ]0, 1]. In view of the preceding study, when dealing with the plasmasphere, the three following assumptions can be retained: The perturbation is such that ν ∼ ε.
As already explained, the hypotheses (Cp1) and (Cp2) correspond to a cold thermal plasma. On the other hand, the hypothesis (Cp3) means that only a small fraction of the plasma is out of equilibrium. By combining information from Section 2.6 with (Cp1), (Cp2) and (Cp3), we can simplify (41)-(42a)-(42b) as indicated below. First, due to (Cp2), there is almost no distinction between the cases α = 1 and α ∈ {2, · · · , N }. The only difference is that |ε α | 1 for all α ∈ {2, · · · , N }, whereas |ε 1 | = ε 1. Thus, for all α ∈ {1, · · · , N }, we can impose: For α ∈ {2, · · · , N }, this equation (53) contains no singular term (of the order ε −1 ), and it is weakly nonlinear. On the contrary, for α = 1, we can note the presence of a fast rotating term (of the order ε −1 ), of a large source term (of size ε −1 ), and of some O(1) nonlinear contribution (of the form E · ∇ p f α ). The equation (42) does not need to be changed, except that the compatibility conditions contained in (42b) must be conveniently weighted: and except that the two expressions ρ(f α ) and (f α ) can be further specified by using (Cp1) and (52) in order to find: 2.8. Within the framework of geometrical optics. The plasmas can support a wide variety of wave phenomena. We refer to [10,12,28] for an overview. In most cases, these phenomena appear to be interconnected sets of particles and fields which can evolve in a periodically repeating fashion [7]. In order to capture the precise features underlying the propagation of such plasma waves, a good strategy is to start from the dimensionless version of the RVM system which has just been exhibited in Paragraph 2.7, made up of (53)-(54) together with (55). This allows to set the discussion within the coherent framework of geometric optics [29,33]. The corresponding asymptotic analysis is new for two main reasons. The first, which is well detailed in Sections 1 and 2.5, comes from the spatial variations of the magnetic field; the second is due to the mesoscopic precision of our model. It is important here to discuss this second aspect more thoroughly.
Being interested in the propagation of electromagnetic waves means to focus on oscillations of the self-consistent field (E, B)(·). Since the function (E, B)(·) depends only on (t, x), a key point is that only time-space oscillations can be involved at this level. With this in mind, we can introduce some smooth phase function φ ∈ C ∞ (M ; R), which depends on the macroscopic variable (t, x) ∈ M but certainly not on the kinetic variable p ∈ R 3 .

Assumption 6.
[non-stationary phase] The function φ is such that: The oscillating behaviour of (E, B)(·) may be viewed as a sum of monophase pieces which can be modelled by: Usually, the time evolution of (E, B)(·) is considered in the framework of MHD descriptions, through fluid models based on Maxwell's equations, involving the variables (t, x). This has the advantage of simplicity. But this also means various simplifying assumptions which are debatable when dealing with plasma phenomena out of equilibrium. To understand all subtleties induced by the underlying presence of p ∈ R 3 , it is necessary to come back to the original RVM system. To this end, we seek the complete solution u of (53)-(54) in the form of a basic monophase representation implying φ through: In (57), the profile U is assumed to be periodic in the fast variable θ ∈ T := R/(2 π Z). The coordinates inside (t, y) are considered as slow variables. When dealing with capital letters like U , the different font styles U , U and U will be used for expressions depending respectively on the variables (t, y, θ), (t, y) and (t, x, , ω). When studying (53)-(54), a first stage is to exhibit approximate solutions. Given N ∈ N * , we aim at a precision of the order O(ε N −1 ) through expansions like: In (58), the profiles U j (t, y, θ) are assumed to be smooth bounded real valued functions: with Fourier series: It is understood that the function F l j,α (·) and its derivatives at all orders satisfy (47). Plugg the real valued function u ε a (·) of (58) into (53)-(54) and into (55). Collect the contributions having the same power of ε in factor, sorted in increasing order. By this way, we get the condition: It turns out that the expressions G j (·) depend only on terms U i with i ≤ j + 1. In particular, for j = −1, we get the preliminary constraint: The equation (61) is interesting and difficult to solve. It contains polarization conditions on U 0 . Above all, it includes the so-called eikonal equation which allows to determine φ and which therefore governs the geometry of the propagation. An approximate solution u ε a (·) can be derived by solving successively the conditions G j ≡ 0 for j = −1, j = 0, and so on · · · up to j = N − 1. The corresponding WKB analysis would be meaningful, but it is postponed to other articles.
From now on, we focus on the initialization procedure, based on (61), which already requires a substantial amount of work.
3. Cold plasma dispersion relations. From now on, the matter is to solve the condition (61), which is inherited from the cold plasma model exposed in Paragraph 2.7, that is from (53)-(54)-(55). Given ∈ N, the Fourier coefficient G −1 depends only on U 0 . There remains: The expressions G −1 (·) are linear with respect to U l 0 with coefficients depending on the choice of φ. With (59) and (60), it follows that: The situation under study is very dispersive. After adjusting φ in order to obtain G l −1 ≡ 0 (and thereforeḠ −l −1 ≡ 0) for some l ∈ Z * , the other conditions G −1 ≡ 0 (with = |l|) are in general not verified (except for the trivial choice U 0 ≡ 0). This is why, at leading order, only one Fourier coefficient will be switched on.

Assumption 7.
[presence of a non-trivial monochromatic electromagnetic oscillation] There is some non-zero integer l ∈ Z * such that: Due to (63) and (64), the conditions inside (62) reduce to: This Section 3 is devoted to the analysis of (65).
3.1. The first step of the WKB calculus. With l ∈ Z * as in Assumption 7, introduce: Recall that, due to (56), we have (τ, ξ) = (0, 0). From (53), knowing that ε 1 = −ε and ε α 1 for α = 1, we can extract: together with: In view of (67), an electric oscillation (E l 0 ≡ 0) is correlated with an oscillation of the electron density distribution (F l 0,1 ≡ 0). On the contrary, the density distributions of ions contain no oscillations (at leading order). In order to satisfy (68), we must impose: With the f ε a,α as in (58), exploiting (69), the charge density ρ and the electric current  which are given by (55) can be expanded in powers of ε ∈ ]0, 1] according to: Coming back to (53)-(54), with the definitions of (48), we have to consider: The rest of the article is devoted to the study of the system (67)-(71) on (F l 0,1 , B l 0 , E l 0 ). The subscript α is not present at the level of (67)-(71). Thus, without any possibility of confusion, we can drop the subscript 1 at the level of θ d 1 , n d 1 and F l 0,1 . From now on, θ d , n d and F l 0 stand for θ d 1 , n d 1 and F l 0,1 .

Looking for reduced-form equations on the electric part.
The aim here is to extract from (67)-(71) necessary conditions involving only E l 0 . Let us go step-bystep.
• 3.1.1.a) Eliminating the presence of the variable r. The expression F l 0 ≡ F l 0,1 can always be factored into: . This allows to convert (67) into: On the other hand, we have: Coming back to (48b), it follows that: cos ω sin sin ω sin cos In the same way, we can compute: Eliminating the presence of the variable ω. This can be done by a Fourier analysis: Then, the condition (72) can be declined into: On the other hand, we are left with: together with: Now, from (75b), we can deduce that: Remark 1. For τ = 0, the first condition and the second condition inside (71c) are obvious consequences of respectively (78) and (71a). Thus, when τ = 0, we can forget the compatibility condition (71c).
• 3.1.1.d) Eliminating the presence of the density component F l 0 . The idea here is to exploit the equations of (75) in order to exhibit a relation between J 0 (F l 0 ) and E l 0 . To begin with, from (75), we can extract:

CHRISTOPHE CHEVERRY AND ADRIEN FONTAINE
where we have introduced the skew-Hermitian matrix: Multiply this identity by (sin ) 2 and integrate with respect to d . There remains: Then, by applying Σ to (79), we get: At this stage, we can state that a necessary condition to obtain non-trivial solutions U l 0 ≡ 0 of the system (65), satisfying E l 0 ≡ 0, is to impose det N(x, τ, ξ) = 0. The discussion about a sufficient condition requires to distinguish between the values of τ .
3.1.2. The stationary case (τ = 0). When τ = 0, we simply find det N(x, 0, ξ) = 0. In fact, the step 3.1.1.c) makes no sense since it removes the information contained in (71b). The same applies for the step 3.1.1.d) since the matrix Σ(x, 0) is not invertible. The situation τ = 0 must therefore be handled separately. Proof. The matter here is to deal with: The condition (83a) -or the condition (75) for τ = 0 -is equivalent to: together with F l,m 0 = 0 for all m ∈ Z \ {−1, 0, 1} and E l3 0 = 0. In the preceding lines, there is no condition on F l,0 0 . We can always adjust the component F l,0 0 so that the first condition of (83d) is satisfied, which amounts to impose: The part (83b) says that the direction E l 0 is parallel to ξ. In particular, we have ξ 3 = 0. Using (76) together with (84), we get ξ · J 0 (F l 0 ) = 0. This is all we need to adjust B l 0 through (83c) and the second relation in (83d). 3.1.3. The electron cyclotron resonance frequencies τ = ±b e (x) . The position x ∈ Ω being fixed, the electron cyclotron resonance frequency is given by the value |τ | = b e (x). It plays a crucial role in plasma physics (and also in condensed matter and accelerator physics). When |τ | = b e (x), the difficulty comes from the step 3.1.1.d). The matrix Σ x, b e (x) is not invertible. This is why this situation must be tackled separately. We will deal here with the case τ = +b e (the other case τ = −b e being very similar). Proof. From Paragraph 3.1.1, we know that the condition det N x, b e (x), ξ = 0 is necessary. To show that it is also sufficient, assume that det N x, b e (x), ξ = 0. Then, there is some E l 0 ≡ 0 such that: From (80), we can see that (1, i, 0) · Σ x, b e (x) = 0. Then, looking at (1, i, 0) · NE l 0 = 0, we can extract E l1 0 + i E l2 0 = 0. Therefore, the condition (75a) with τ = b e (x) is met. On the other hand, the relations (75b) and (75c) impose: The constraint (75d) implies F l,m In view of Remark 1, there remains to look at (71b), or equivalently at (79). First, with (86), we get: Recalling (76), we can define: By developing the content of the equation (79) with τ = b e (x), and by eliminating E l2 0 through the relation E l2 0 = i E l1 0 , we get: The subsystem (88a)-(88b) is overdetermined. It implies a compatibility condition, which can be exhibited by looking at the combination i (88a) + (88b) = 0, which is: Assuming (89), the equations (88a) and (88b) can be solved by adjusting J 1 0 (F l,−1 0 ) accordingly. It suffices to select adequately some F l,−1 0 (t, x) not depending on . This can be done because no condition on F l,−1 0 (·) has been imposed so far. To conclude, it suffices to observe that the system (85) and (88c)-(89) are equivalent. Indeed, the two first components of (85) correspond to the conditions (1, i, 0) · N E l 0 = 0 and (i, 1, 0) · N E l 0 = 0, and they turn to be E l1 0 + i E l2 0 = 0 together with (89) . On the other hand, the equation (88c) is the same as the third component of (85).
In the next Paragraph 3.2, we examine more carefully the content of this restriction.

The characteristic variety of cold magnetized plasmas.
Readers should refer to the next Paragraph 3.2.1 if there is doubt about notations or definitions. This information will be used repeatedly in what follows. It is also needed when applying Theorem 1.1.

Notations and definitions.
With the rescaled version (36) of the electron density n d ≡ n d 1 and with µ ≡ |µ 1 | as in (51), introduce: With b e and τ as in (28) and (66), look at the skew-Hermitian matrix: With the preceding ingredients, consider the matrix: As will be seen later, all the values of (τ, ξ) ∈ R 4 \ {0}, including the singular cases 3.1.2 and 3.1.3, can be treated together. To this end, however, the variable τ must be factorized from the characteristic polynomial det N. As will be shown in Lemma 3.12, see the definitions (111) and (112), there exists an explicit function χ ∈ C ∞ Ω × (R 4 \ {0}) such that: The matrix O has been defined at the level of (32), while the matrix N is given by (94) with κ and Σ as in (92) and (93). For all x ∈ Ω, the function χ(x, ·) is polynomial with respect to the dual variables (τ, ξ) ∈ R 4 . Thus, given x ∈ Ω, the zeros of χ(x, ·) define some algebraic variety in T * (t,x) M .
The representation (96) of V offers a complete intrinsic vision of the characteristic variety. Now, we can decompose the direction ξ = t O(x) ξ ∈ R 3 of (66) in spherical coordinates, as we did for p ∈ R 3 in Paragraph 2.6.2. With the convention: S ph (r, ω, ) := t (r cos ω sin , r sin ω sin , r cos ) , we have: Note that the couple (ω, ) inside (98) differs from the the one of Figure 1. There will be no possible confusion since the decomposition (46) of p ∈ R 3 will no longer be used. Below, we fix (t, x) ∈ M and we introduce useful objects when investigating what happens in the cotangent space T * (t,x) M . According to the forthcoming presentation, the discussion will not directly involve the direction ξ but instead the variable ξ = t O(x) ξ through the spherical representation (r, ω, ) of ξ, see the definition (97)-(98) and Figure 2. In fact, due to some gyrotropic invariance, the angle ω will be almost absent. In particular, the function χ(·) of Lemma 3.12 does not depend on ω. It remains to deal with (r, ). Now, since the matrix O(x) is orthogonal, there is an obvious link between (r, ) and ξ, namely: With (r, ) as in (99) and χ x, (·) as in (111)-(112), retain that: As explained in the introduction, the characteristic variety V can be written: As indicated in (101), the subset V o (x, ) (for ordinary waves) and the subset V x (x, ) (for extraordinary waves) form some non-trivial partition: They can be recovered from Theorem 1.1 through ingredients τ ± 0 (·), τ ± ∞ (·) and g ± (·) that are clearly specified in the definitions below.

Definition 3.7. [dispersion relations for ordinary waves and extraordinary waves]
The links r 2 = g ± (x, , τ ) inside (1) involve the (generalized) Appleton-Hartree functions: where: Retain that all the expressions inside Definitions 3.6 and 3.7 can be interpreted as depending on ξ or ξ. It suffices to exploit (99).

The characteristic variety from the physical viewpoint.
In what follows, we will work in the non-singular case, with τ = 0 and with |τ | = b e (x). Introduce: In optics, the following basic notions often come into play.
Definition 3.8. The refractive index is the vector n ∈ R 3 given by (106).
Convections currents resulting from the plasma electrons and ions are usually accounted for the derivation of some dielectric tensor. Definition 3.10. The relative permittivity D(x, τ ), called sometimes the dielectric tensor, is the Hermitian matrix: In most plasma physics books [12,30,37], the presentation of V is achieved in a way that differs from (96). There, the construction of V is based on the Definitions 3.8, 3.9 and 3.10. Now, to recover (96) from n, σ and D, we can proceed as follows.
Proof. Let us start by computing the determinant: Thus, for τ ∈ R \ 0, ±b e (x) , the matrix Σ(x, τ ) is invertible. By construction, we have: The function M(·) of physicists is clearly not defined for the values τ = 0 and τ = ±b e (x). By contrast, the function χ(·) is defined and smooth everywhere, while containing as much information. That is why the use of χ(·) is mathematically more suitable. whereas χ x, (·) is a bivariate polynomial in (τ, |ξ|). More precisely, written as a polynomial in |ξ|, the function χ x, (τ, ·) is biquadratic: The coefficients A x, (·), B x, (·) and C x (·) are in R[τ ], given by the explicit formulas: Proof. First, note that: The computation of det M(x, τ, ξ) is completely similar to the calculus of R. Fitzpatrick, see the line (4.43) at page 94 of [12]: where: Then, with (113) and the definitions of n, γ and ς from (106), it suffices to write (114) as a rational fraction in powers of τ to get the result from Lemma 3.12.
The restriction χ x, (τ, |ξ|) = 0 inside (96) can be handled as a quadratic equation to be solved for |ξ| 2 , with discriminant ∆ := B 2 −4 A C. As a consequence of Lemma 3.13 below, the corresponding roots prove to be (positive, negative or zero) real numbers.
can be put in the form:

The function Q(·) is defined as in (105b). It is clearly non-negative. It vanishes if and only if
Proof. It follows from the second formulas listed in (112). When computing ∆, it suffices to identify the terms which are in factor of the different powers of sin 2 .

The physical motivations behind the study of the characteristic variety.
Looking at the set V ⊂ T * M is interesting because we can send or receive information only by using points (τ, ξ) in the phase space T * M which are microlocalized inside V .
Lemma 3.14. A necessary and sufficient condition to find solutions U l 0 to the system (65), with some non-trivial electromagnetic part (E l 0 , B l 0 ) ≡ 0 satisfying Assumption 7, is to impose the eikonal equation: Proof. When τ ∈ R * , in view of (95), the condition (118) is equivalent to: Now, from Paragraphs 3.1.4 and 3.1.3, we know that (119) is a necessary and sufficient condition. In the stationary case, when τ = 0, Paragraph 3.1.2 says that ξ 3 = 0 or equivalently that = π/2 is a NSC. Thus, it suffices to show that the restriction (96), that is: can be realized if and only if = π/2. With (112), just remark that: Since |ξ| = 0, we must have cos = 0, that is = π/2 as required.
The set V can be viewed as the union of disjoint subsets: Thereby, the set V *  (97), we can fix (x, ) ∈ Ω × [0, π], and define (we will not mark the star * at the level of V): With ξ as in (98), the coefficients A, B and C depend only on x, τ 2 and . Thus, starting from V(x, ), we can recover the whole section V * (x, ) by rotations: This rotational invariance can be clearly seen in Figures 10, 11, 12 and 13. By this way, the issue reduces to the study of V(x, ), which is generically a onedimensional subset of the quadrant R + × R + . Note also that the coefficients A and B satisfy: . This brief discussion can be summarized in a remark.

Remark 3.
All the geometrical information is provided by the intersection of V * x with the quadrant of an hyperplane. More precisely, we can focus on: and then recover the set V * x by changing τ ∈ R + into −τ , by replacing ξ 3 ∈ R + by −ξ 3 , and by applying rotations with axis R t (0, 0, 1) and angles ω ∈ [0, 2 π[. This means also that, in practice, we can restrict our attention to the case where (τ, ) ∈ R + × [0, π/2]. Lemma 3.13 indicates that the situation = 0 (existence of a double root when τ = κ) and the other cases ∈ ]0, π 2 ] must be distinguished. In plasma physics books, the two values = 0 and = π 2 are usually studied separately, and they are exploited to obtain a classification of waves. However, this viewpoint does not give access to a complete vision of the characteristic variety V . In what follows, we will rather consider that the generic situation is when ∈ ]0, π 2 [, while the two cases = 0 and = π 2 can be recovered as (very special) endpoints. 3.3.1. Parallel = 0(π) non-singular τ ∈ R * + \ {b e (x)} propagation. The matter here is to study the equation χ x,0 (τ, |ξ|) = 0, which is simply: Either we have τ 2 − κ 2 = 0 or Z = |ξ| 2 must be subjected to: , this quadratic equation has two roots Z 1 (x, τ ) and Z 2 (x, τ ) given by: The set V(x, 0) is usually (see for instance [12]-chapter 4) separated into connected parts, corresponding to the different possibilities marked above. a) Parallel left-handed circularly polarized wave |ξ| 2 = Z 1 (x, τ ) . Knowing that τ ∈ R * + , this can be achieved for some |ξ| ∈ R + if and only if: We find only one connected component: Parallel right-handed circularly polarized wave |ξ| 2 = Z 2 (x, τ ) . Knowing that τ ∈ R * + , this can be achieved for some |ξ| ∈ R + if and only if: Accordingly, the set V r (x, 0) can be split into two connected parts V − r (x, 0) and V + r (x, 0) which are distinguished below.
b.1) The first part V − r (x, 0) contains both the dispersion relation for Alfvèn waves (values of τ such that 0 τ ) and the dispersion relation for whistler (or electron cyclotron) waves (values of τ such that τ b e ): 2) The second is the standard right-handed circularly polarized wave: c) Parallel longitudinal wave or electrostatic plasma wave (τ = κ). With κ as in (92), this means a link between τ and x, whereas no condition is imposed on |ξ|.
For τ = τ ± ∞ , as a consequence of Lemma 3.13, the equation χ x, (τ, |ξ|) = 0 has two distinct real solutions which can be described by: When g ± (x, , τ ) < 0, there is no mathematical solution. The associated waves are usually considered to be evanescent. On the contrary, when 0 ≤ g ± (x, , τ ), the condition (131) is achieved for some |ξ| ∈ R + . This means that a wave can be propagated. To distinguish between these two situations, we have to study the functions g ± (x, , ·). To illustrate the following discussion, we can directly plot the graphs of the functions g ± (x, , ·). The content of Figures 4 and 5 is justified below.   For all x ∈ Ω, the function τ + ∞ (x, ·) is strictly increasing on the interval [0, π/2], and it satisfies: For all x ∈ Ω, the function τ − ∞ (x, ·) is strictly decreasing on [0, π/2], with: Proof. By direct calculation.

Lemma 3.17. [comparison of resonance and cutoff frequencies] If
Proof. Assume that √ 2 b e (x) < κ(x). In view of (133), to get (134), it suffices to note that: . The other comparisons follow from straightforward computations.
. The functions g + (x, , ·) and g − (x, , ·) of (104) are of class C ∞ respectively on the domains With the expressions given in Definition 3.7, the functions g ± (x, , ·) are clearly of class C ∞ on R + \ {τ + ∞ (x, ), τ − ∞ (x, )}. Moreover, using the right part of (131), they can be written in the form: The difficulty is to show that g ± (x, , ·) is well defined and smooth near τ ± ∞ (x, ). To this end, we have to examine more precisely the denominator in (136). First, combining (112b), (112c) and (130), the term ∆ = b 2 e κ 4 Q of (117) can be written in the form: In view of the second formulas in (105a) and (112a), we have: From (137) and Lemma 3.16, we can easily deduce that: And therefore: This means that the denominator of (136) does not vanish. This gives the conclusion.
We are now able to distinguish between the different regions of wave propagation. This can be done by selecting the connected parts of V(x, ). This means to isolate where the functions g − (x, , ·) and g + (x, , ·) are positive. We will adopt the following classification. a) Oblique ordinary wave |ξ| 2 = g + (x, , τ ) ≥ 0 . With Lemmas 3.18 and 3.19, together with the intermediate value theorem, we get: The set V o (x, ) consists of two connnected parts V − o (x, ) and V + o (x, ) which are defined by: The b) Oblique extraordinary wave |ξ| 2 = g − (x, , τ ) ≥ 0 . The set V x (x, ) is also composed of two connected parts V − x (x, ) and V + x (x, ) which are defined by: Retain that: The case of perpendicular propagation ( = π/2) can be incorporated within the general framework of oblique propagation. However, using (103), we get τ − ∞ (x, π/2) x (x, π/2) and V + x (x, π/2) which correspond to extraordinary waves. Some sets V o (x, ) and V x (x, ) are represented on Figure 6 below.    (150) and (151). This allows to provide a fresh perspective on what happens when = 0. When = 0, the formula (104) yields: .
Recall (143) which guarantees that κ < τ + 0 . As a consequence of (153), the definitions (128) of V + r (x, 0) and (151b) of V + x (x, 0) coincide. The function g ± (x, 0, ·) is not suitable for the description of the other components of V(x, 0). This is because is not continuous at the point τ = κ. For = 0, the identity (103) gives rise to: Under the overdense condition, we find: 0) . Under the underdense condition, we can further decompose V − x (x, 0) into the disjoint union of V −− x (x, 0) and V −+ x (x, 0), with: Then, we can recover: Finally, due to the disontinuity of the functions g ± (x, 0, ·), the case τ = κ has to be treated separetely. In fact, the corresponding study is exactly the same as in Paragraph 3.3.1c). We recover the case of parallel longitudinal waves. These

3.4.
Other parametrizations of the characteristic variety. Until now, to describe the set V , we have looked at the direction ξ ∈ R 3 in fact at |ξ| ∈ R * + for a given ∈ [0, π] as a function of the frequency τ ∈ R * + , with τ selected in convenient subdomains. The aim of this Section 3.4 is to investigate the other ways to proceed.
3.4.1. The characteristic variety at a fixed value of ω. The description of V * x through the functions g ± (·) does not involve ω and the sign of τ . As observed in Paragraph 3.2.4, this means that V * x is invariant under rotations with axis R t (0, 0, 1) and angles ω ∈ [0, 2 π[. Since we can recover the case ω = 0 from the case ω = 0 by such rotations, it suffices to consider the case ω = 0. To this end, as pointed in Remark 3, we can cut V * x with a well-chosen quadrant of an hyperplane. Moreover, we can consider ∈ [0, π] as a variable rather than as a parameter (like we did before in Paragraph 3.3.2). By doing so, we get a three dimensional representation of V * x , which could be materialized through: In accordance with the classification of Paragraph 3.3.2 into the categories a) and b), the set V * x can be separated into connected parts: V + x (x) := ( , τ, r) ; τ + 0 (x) ≤ τ , r 2 = g − (x, , τ ) . The sets V ± o (x) and V ± x (x) are drawn on Figure 9. We also include (in the plane r = 0) the graphs of the functions τ ± ∞ (x, ·) Figure 9. The characteristic variety in the variables (τ, , r).

3.4.2.
The characteristic variety at a fixed value of τ . Our purpose here is to give a detailed description of the following set: For a given ξ ∈ R 3 , we use the notation ξ = S ph (|ξ|, ω ξ , ξ ) which comes from the convention of Figure 2 to represent ξ. In connection with what has been done in Paragraph 3.3.2, the set ν(x, τ ) can be decomposed into ν( The access to ν is achieved through the functions g ± (·). Now, given (x, τ ) ∈ Ω×R + , the domain of definition, the sign, the regularity and the behaviour of g ± (x, ·, τ ) : [0, π/2] −→ R depend strongly on the choice of τ . Recall that we have to deal with the additional condition |ξ| 2 = g ± (x, , τ ) ≥ 0. Thus, of particular interest are the intervals where g ± (x, ·, τ ) is non-negative. We can first deal with the case of g + (·). Proof. From the proof of Lemma 3.18, we know that g + (x, ·, τ ) is defined and continuous except if τ = τ − ∞ (x, ). Given τ ∈ R + , from Lemma 3.16, this can happen for some angle − ∞ (x, τ ) ∈ [0, π/2] if and only if we are in case i). Now, assuming i), we have: Then, coming back to the expression (104) and using the information (157), we can identify as in i) the interval where g + (x, ·, τ ) is non-negative, as well as the asymptotic behaviour (156). In fact, to get both i) and ii), it suffices to look at Figure 5.
The wave equation describes the propagation of electromagnetic waves in a vacuum. In its dimensionless version, it takes the form: The characteristic varieties associated with these two equations (on E and B) are the same. They are simply given by the light cone: Then, given τ ∈ R * + , we recover the sphere of radius τ : Now, in the presence of an exterior magnetic field, the picture (161) is strongly affected. The first effect, already encoded within (155), is the well-known separation of ν into the two parts ν o and ν x corresponding to ordinary and extraordinary waves. The second impact is a disappearance of ν o (x, τ ) and ν x (x, τ ) when τ falls into the gap between resonance and cutoff frequencies, and their resurgence in the form of cones below resonance frequencies. With this in mind, we can highlight the following distinctions. o (x, τ ), which are respectively above and below the plane {ξ; ξ 3 = 0} and given by: Proof. This statement is a direct consequence of Lemma 3.22, case i). In particular, the asymptotic behaviour (156) implies that neither ν a o (x, τ ) nor ν b o (x, τ ) are bounded. On the contrary, for large values of r, these two sets become close to the two cones: The distortion from ν l is here obvious. As long as τ remains sufficiently small, which means in practice that we deal with Very Low Frequency waves (at approximately 1 kHz), there is no similarity between ν o (x, τ ) and ν l (τ ). The same applies to ν x (x, τ ).
R t (0, 0, 1) + ∞ Figure 11. Extraordinary resonance cone. νx(x,τ) Proof. This statement is a direct consequence of Lemma 3.23, case ii).  The ordinary resonance cone can appear with an ordinary sphere if κ(x) < √ 2b e (x) (see Lemma 3.17). An extraordinary resonance cone can occur with an ordinary sphere which, in this case, is located inside. For the sake of clarity such configurations are not made apparent in Figure 10 and 11). For τ > τ + 0 (x), the extraordinary sphere and the ordinary sphere coexist (see Figures 12 and 13), the first being included in the second. At higher frequencies, that is concretely above 1 MHz, the distinction between ν o (x, τ ) and ν x (x, τ ) becomes almost invisible. For  (121), recalling the convention (98), we can assert that: As noted in the article [31], interpreted as a polynomial in the variable τ instead of |n|, the characteristic equation χ x, (τ, |ξ|) = 0 turns to be of fourth order in τ 2 .
More precisely, we have χ x, (τ, r) = P x, ,r (τ 2 ) with: Explicit algebraic expressions τ (x, |ξ|, ) are of course available, but they are not so easy to implement. We will proceed differently. Another way to extract information is to study the properties of the function g ± (x, , ·). We start by looking at the two extreme situations ( = 0 and = π/2). Then, we consider the oblique case ∈ ]0, π/2[. , , Proof. For L, remark that: From (124), we know that Then, looking at the righthand term of (169), we get the first result. Otherwise, just compute ∂ τ R(x, τ ).

When
= 0, there is a mix between ordinary and extraordinary waves. This effect is enhanced under the underdense assumption, see Figure 14. It follows that the independent use of L(x, ·) or R(x, ·) is not suitable to get inversion formulas that are globally defined in . The reason is that the expressions L(x, ·) or R(x, ·) (taken separately) do not allow a continuous connection to g ± (x, ·) when goes to 0. To remedy this, we have to blend the functions L(x, ·) and R(x, ·) accordingly. This induces the technicalities reproduced below. We start with the overdense case (Definition 3.21) which is simpler. For ordinary waves below the resonance frequency, just take: to recover: r) . For extraordinary waves located below the resonance frequency, we have to follow the red curve in Figure 14-(a). This means introducing: κ(x)) , in order to obtain: r) . For ordinary waves above the cutoff frequency, we have to reconstruct the blue curve in Figure 14-(a). Just define: that gives access to: V + o (x, 0) = (τ, r) ; 0 ≤ r , τ = F 2 +,ov (x, r) . We now move to the underdense case. As already noted in Paragraph 3.3.3, the situation is more intricated when b e (x) > κ. However, the logic is the same. We have to follow the coloured lines (blue and red) at the level of Figure 14-(b). Briefly define: if r ≥ R(x, κ(x)) , The situation of V + x (x, 0) is simpler because V + x (x, 0) just coincides with V + r (x, 0). It suffices to define: V + x (x, 0) = (τ, r) ; 0 ≤ r , τ = F 2 − (x, r) .
3.4.4. The eikonal equation. In Paragraph 3.1, we have approximated the RVM equations by using the WKB theory. The first outcome of this procedure is an eikonal equation which provides a simplified but pertinent representation of electromagnetic wave propagation in terms of geometric optics. relation, as found in [12,30,37]. Below, by combining the refined model exhibited in Theorem 1.1 with the information from [7], we can get a better understanding of what happens.
According to [7], the wave packets are emitted inside the characteristic variety V , all along the resonance cones (drawn at the level of Figures 10 and 11). Moreover, a signal which is generated at the time t from the position (t, x, τ, ξ) ∈ V starts to travel at the velocity |∇ ξ λ i −, (x, ξ)|. In view of the dispersion relations, the size of |∇ ξ λ i −, (x, ξ)| may be significant when ξ is located away from the resonance frequencies τ ± ∞ (x, ) with given by (99). Then, as described in Theorem 2 of [7], the electromagnetic signals should appear again and again over consecutive periods of almost equal length. They are marked by light-coloured vertical lines in the spectrograms, see Figure 17. On the contrary, the size of |∇ ξ λ i −, (x, ξ)| becomes small when ξ approaches τ ± ∞ (x, ). Then, there is almost no propagation. This could account for the creation of quasi-electrostatic waves. This difference in behaviour could explain why chorus emissions are usually composed of discrete narrowband wave elements accompanied by banded incoherent waves [27].  Moreover, chorus emissions usually occur in two frequency ranges, a lower band and some upper band, which are separated by a gap [27,32,35]. The high part of the lower band might correspond to the resonance frequency τ − ∞ . On the other hand, the upper band could be delimited from below by the cutoff frequency τ − 0 and at the top by the resonance frequency τ + ∞ . To well illustrate these correspondences, the comparison between Figures 16 and 17 is facilitated above. Note also that the lower and upper bands would be associated with respectively the ordinary and extraordinary modes of propagation.
The fact [32] that lower-band waves tend to be field-aligned ( 0) whereas upper-band waves seem to be highly oblique ( π/2) is another aspect that can be interpreted as a consequence of our analysis. Indeed, away from the emission points, the (non propagating) electrostatic waves should not be detected. Such waves are generated from the vertical lines in Figures 6 and 7. This yields the line τ ≡ τ + ∞ = κ in the perpendicular case = 0, and τ ≡ τ − ∞ = 0 in the parallel case = π/2. As a result, away from the emitting regions, one would expect to observe 420 CHRISTOPHE CHEVERRY AND ADRIEN FONTAINE principally a lower-band (related to τ − ∞ ) when 0, and a upper band (related to τ + ∞ ) when π/2.
In conclusion, the refined structures of the characteristic variety V help to improve our understanding of the morphological properties of chorus emissions.