Solutions to the cold plasma model at resonances

. — Little is known on the mathematical theory of hybrid and cyclotron solutions of Maxwell’s equations with the cold plasma dielectric tensor. These equations arise in magnetized plasmas to model an electromagnetic wave in a Tokamak. The solutions can behave extremely diﬀerently from those in vacuum. This work contributes to the local theory of the hybrid resonance by means of an original representation formula based on special functions, a certain eikonal equation and with a careful treatment of the singularity. Résumé (Solutions singulières résonantes pour un plasma). — La théorie mathématiques des équations de Maxwell avec le tenseur du plasma froid en présence d’une résonance cyclotron ou hybride est peu développée. Ces équations sont présentes pour modéliser la propagation d’une onde électromagnétique dans un plasma magnétique tel que celui d’un Tokamak et les solutions peuvent être très diﬀérentes de la propagation dans le vide. Ce travail contribue principalement à la théorie locale de la résonance hybride avec des formules originales de représentation à partir de fonctions spéciales. Ces formules sont obtenues au moyen d’une équation eikonale et d’un traitement spéciﬁque de la singularité.


Introduction
The goal of this work is to identify singular solutions of the cold plasma model at different resonances with an original mathematical method.Singular means here that a function u satisfies |(x − x * )u(x)| c > 0 uniformly with respect to x in a neighborhood of a point x * .In a more general sense, a singularity can also be a Dirac mass or a principal value.The problem comes from electromagnetic waves modeling in strongly magnetized plasmas and from the so-called resonant heating phenomenon.We more specifically consider the one species cold plasma model which represents a collection of zero-temperature electrons immersed in a uniform static magnetic field.Resonances correspond to limits of this model, and will be studied in this work thanks to a term representing friction of electrons on the background ions.With respect to the recent work [9], which was based on certain singular integral kernels and representations, the novelty of the present one lies in the use of a local representation of the solutions in terms of Bessel functions, and the introduction of a new stretching function satisfying an eikonal equation.Moreover, with this new simpler technique we are able to prove that the resonant heating is independent of the tensorial nature of the dissipation regularization term.However, a restriction of this approach is the local character of the results.In particular it is not possible to discuss boundary conditions at finite distance or at infinity with the techniques developed below.The reader interested by this topic may refer to [9], see also a recent alternative technique in [6].
In time harmonic regime and slab geometry in dimension two, for (x, y) ∈ R 2 , there is no resonance in the equation for the transverse magnetic polarization (also denoted as O-mode equations in plasma physics), while the equations for the transverse electric polarization (also denoted as X-mode) have two resonances.With the normalization ω 2 /c 2 = 1, the latter equations reduce to (1.1) The bulk magnetic field B 0 is oriented along the third axis, the electric field (E 1 , E 2 ) is transverse to the bulk magnetic field and the third component of the electric field is zero.The pulsation of the harmonic wave is ω > 0. As described in [18,19,22,9], a first order expansion in the cold plasma model yields the following dielectric tensor , where the cyclotron frequency is defined by ω c = e|B 0 |/m e while the plasma frequency is defined by ω p = e 2 N e /ε 0 m e .Here e is the electron charge, m e is the electron mass and ε 0 is the vacuum permittivity.The physical validity of this model is discussed in [4] are the parallel and perpendicular components-with respect to the bulk magnetic field B 0 -of the wave vector, v th is the electron thermal speed, and n is an integer.Note that the model is not valid if the frequency ω matches an integer multiple of the cyclotron frequency.Assume that the magnetic field and the electronic density depend on the horizontal variable x only: that is B 0 = B 0 (x) and N e = N e (x).
It corresponds to a wall facing a stratified plasma as depicted in Figure 1.1, as in a Tokamak where an antenna sends an electromagnetic wave toward the plasma, in order to probe or heat the plasma.The heating of a magnetized fusion plasma with X-mode in slab geometry: the domain.In a real physical device an antenna on the wall sends an incident electromagnetic wave toward the plasma.Here, for simplicity, the domain is the half plane {(x, y) ∈ R 2 , x > 0}.The medium is filled with a plasma with dielectric tensor given by (1.2).
such devices is an important issue for the ITER project [12].Idealized coefficients representing the main features of the physical coefficients at the wall of a Tokamak are described in Figure 1.2.More realistic plasma parameters may be much more oscillating as illustrated in the numerical studies [16,17].
Resonances correspond to a limit of the geometric optics, as the local wave number approaches infinity.They require a careful examination since some of them are related to singular behavior and absorption.The interested reader can refer to the physics textbooks [18,19] or to [22] for a historical paper on the topic.A recent reference is [10]; it also provides a list of additional references.Singular solutions of wave equations appear in other applications, such as cloaking [21] and metamaterials [7,2].As announced earlier, it is then natural to model the effect of collisions with a bath of static ions by an additional friction term, since the energy is absorbed by the ions.The validity of this modeling has been discussed in [3].Describing the behavior of the solutions as the friction parameter ν > 0 approaches zero will evidence the singularity of the solutions to the limit problem for ν = 0. Indeed one obtains solutions which The singular solutions will be built by means of new explicit representation formulas, allowing a careful extraction of the singularities.The singular behavior of the limit solutions is strongly related to the structure of the dielectric tensor (1.2).This crucial structure will be highlighted in Section 3.
The complex pulsation ω ν = ω+iν shifts the exact pulsation in the upper half-plane by the factor ν > 0. With this shift, Equations (1.1)-(1.2) become with the regularized dielectric tensor (1.4) where ω p and ω c are still respectively the plasma and cyclotron frequencies, and the regularized dielectric tensor is independent of the y-direction.The physical justification of ε ν can be found in plasma physics textbooks [4,18].For a recent mathematical treatment of these equations by means of singular integral equations, refer to [9].In fusion plasmas, the value of the friction parameter can be extremely small.For example, relative values of ν/ω ≈ 10 −7 are common.This extremely small value shows that the frictionless limit regime ν → 0 is relevant for some fusion applications.In such regimes the formal first order approximation of the dielectric tensor reads (1.5) The dissipation tensor D accounts for the underlying physical dissipation.Its coefficients are As previously mentioned, and stressed in Figure 1.2, the coefficients of these tensors are a priori non-constant functions of the horizontal variable x.
The structure of the limit tensor (1.2) shows two different resonances in the limit model ν = 0, defined thanks to the dispersion relation, see [18].The first one appears as a singularity in the coefficients of the equation, and is commonly referred to, in the plasma physics community, as a cyclotron resonance where ω c is a function of the x variable, a priori non-constant.We will show that mathematical solutions of Maxwell's equations (1.3) actually have no singular behavior at the cyclotron resonance: they are bounded in this regime, despite the fact that the dielectric coefficients are not.The other interesting regime corresponds to the vanishing of the diagonal part of ε 0 , that is to say It defines the hybrid resonance.In this regime the dielectric tensor is non-singular since the denominator is non-zero, i.e., the coefficients of the dielectric tensor are bounded.However, the solutions are known to be highly singular: they are not bounded near the resonance.A local expansion of the horizontal electric field yields representation formulas near isolated roots , where g ν is a function which is proved to be regular enough so that one can pass to the limit ν → 0 + .One gets [9] This is a well-known behavior in the plasma physics community, but has not received full attention by the mathematics community.One goal of this work is to develop a local comprehensive analysis of this phenomenon.As stressed in [9], a useful quantity to characterize the singular behavior of the solutions is the resonant heating of the plasma.Up to a non-essential factor ω, the resonant heating [4,9] is defined, for any where the Hermitian product of two complex vectors A and B is denoted as (A, B) = A • B. This formula has been justified in [9] where it is shown that Q > 0 in case of resonant heating using the simplification ε ν = ε + iνI, i.e., in the ν expansion (1.5) of ε ν the dissipation tensor is D = I.This simplification is mathematically justified since it corresponds to the limiting absorption principle.In this work we will show that the resonant heating is well-defined and takes the same value for a large class of dissipation tensors that includes physically based tensors.The analysis method uses 1 which is, after linearization around the bulk magnetic field, the parallel component of the magnetic field.System (1.3) is then equivalent to the following first order system Since the coefficients ε ν ij depend only on x, one can perform a Fourier transform in y.Denoting by iθ the Fourier coefficient, it yields (1.9) The Fourier coefficient θ can be considered as a frozen parameter from now on, and ∂ x will be replaced by d/dx denoting the derivative with respect to x.In this system, the derivative with respect to x appears only on W ν and E ν 2 .The energy integral (1.8) can be expressed for any J.É.P. -M., 2017, tome 4 It leads to which expresses a balance of energy.
It is convenient to consider a differential system with two unknowns E ν 2 and W ν , while E ν 1 can be computed thanks to the first equation of System (1.9) by
Eliminating E ν 1 , the reduced system then reads (1.12) where the matrix is defined as follows.
-The matrix M ν (x, θ) is constructed from the dielectric tensor by , where the coefficient d ν is the determinant of the dielectric tensor ε ν : (1.14) The determinant D ν (•, θ) of the matrix can be simplified as follows (1.15) It is already foreseeable that the roots of ε ν 11 will play a crucial part in some cases.The analysis of singular solutions around the hybrid resonance x = x * relies on the analyticity of the coefficients of the dielectric tensor.This assumption is very useful to study the problem with a convenient shift in the complex plane.The main results of this paper stem from a precise local study of the solution of the reduced system (1.12) when the entries of M ν have a simple pole at x * , and the determinant of the matrix has a simple pole as well.The mathematical method of singularity extraction is completely new to our knowledge.
Our main results can be summarized as follows: Lemma 2. -In our geometry, the cyclotron singularity of the dielectric tensor, defined by ω = ω c , is only an apparent singularity in the sense that the electromagnetic field is bounded and the resonant heating is zero Q cyc = 0.
Theorems 1 and 2. -The hybrid resonance, defined by ω 2 = ω 2 c +ω 2 p , yields a singular horizontal component of the electric field, i.e., |(x − x * )E ν=0 1 | c > 0 uniformly with respect to x in a neighborhood of x * , and a generic positive resonant heating Q hyb > 0. The limit of the horizontal component of the electric field can be expressed as the sum of a Dirac mass, of a Principal Value and of a bounded term, as in [9].
Corollary 1. -The hybrid resonant heating value is independent of the local dissipation tensor, and moreover the local dissipation tensor can be non-physical as well.The limit value for ν = 0 + of the electromagnetic unknowns is independent of the local dissipation tensor.The uniqueness of the limit solution with respect to the local dissipation tensor is established.This is a new result since in [9] the limit was established only for D = I.
Numerical section.-We provide numerical evidence of our claims with the help of a numerical method that can be used to treat other problems with the same structure.
Extension to multi-species.-We finally provide some formulas extending the method of singularity extraction to multi-species models, to model ions as a fluid as well.These formulas can be compared with similar formulas [4,23] for the cyclotron resonance and the hybrid resonance.
The work is organized as follows.The cyclotron singularity of the dielectric tensor is studied in Section 2, where the electromagnetic field is shown to stay bounded in the neighborhood of this singularity.The hybrid resonance is studied in detail in Section 3, where the singularity of the solution is carefully extracted thanks to convenient representation formulas with Bessel functions.We prove both the generic resonant heating formula and the fact that its value is independent of the local dissipation tensor.Numerical results are provided in Section 4. Some technical material is postponed to the appendix together with a particular analytical solution and the generalization to multi-species.
Acknowledgements.-The authors deeply thank the anonymous reviewers for their remarks, which helped a lot to improve the present work.

The cyclotron resonance
The cyclotron singularity corresponds to a formal singularity in the denominator of (1.4).We consider the situation where there exists at least one x c in the domain such that ω c (x c ) = ω.For the simplicity of the analysis, we will assume that d dx ω c (x c ) = 0 so that x c is isolated.Clearly the singularity at x c in the dielectric tensor (1.4) in the limit case ν = 0 + disappears in the regularized case ν = 0.In order to focus on the cyclotron resonance, it has to be away from the hybrid resonance, that is ω p (x c ) = 0. Definition 2. -A point x c is referred to as an isolated cyclotron singularity if ω c (x c ) = ω and d dx ω c (x c ) = 0.
Note that such a singularity does not depend on ν, since ω c does not depend on ν.Let us now consider the reduced system (1.12),where the system's matrix M ν (•, θ) is defined with the physical dielectric tensor (1.4).The system is well-defined at a cyclotron resonance for all ν ∈ R, and one can pass to the limit in the following sense.
Lemma 1. -Assume that the magnetic field B 0 and the density N e are smooth.Assume an isolated cyclotron singularity occurs at x c , and assume there is no hybrid resonance on some neighborhood [x c − r, x c + r] for some r > 0. Assume (E ν 2 , W ν ) solves the reduced system (1.12) on [x c − r, x c + r] for all ν ∈ R, with a given Cauchy data which admits a finite limit when ν goes to 0 written as 2 , W ν ) tends to a bounded limit E 0 2 , W 0 as ν tends to zero, and (E 0 2 , W 0 ) solves the limit system (1.12) on [x c − r, x c + r] with ν = 0, with the same Cauchy data Proof.-Definition 1 with the physical tensor (1.4) yields after some careful calculations where The determinant (1.15) of this matrix can be recast as The limit value of α ν (•) as ν approaches 0 In this situation the matrix M ν (•, θ) is locally bounded and smooth.So the solutions of the reduced system (1.12) are smooth locally around x c and one can pass to the limit in this system.
Note that there is still a particular behavior in the neighborhood of this cyclotron resonance, known as "turning point behavior" [8,13].However, this does not induce any singular regime, except in the limit ω → +∞, which is not considered here.The consequences of the representation (2.1) on physical quantities are characterized in the following result.
Lemma 2. -Under the assumptions of Lemma 1, the component E ν 1 tends pointwise to a bounded limit E 0 1 as ν tends to zero and the value of the resonant heating is zero: Proof.-The horizontal component of the electric field can be evaluated thanks to (1.11).The coefficients in (1.11) can be expressed as Using (2.2), the limit ν = 0 value of these coefficients at So one can pass to the limit in (1.11) since the coefficients are bounded in the interval Therefore the horizontal part of the electric field E ν 1 admits a bounded and smooth limit E 0 1 .Consider now the heating term.On the one hand we notice that the electric field can be split into two parts the first vector being an eigenvector of the dielectric tensor:

The corresponding eigenvalue
p )E ν 2 /α ν , admits a finite limit for ν → 0, due to (1.11) and (2.4).So the second vector is O (ν + |x − x c |), which means that it counterbalances the singularity of the dielectric tensor.This shows that On the other hand the tensor ε ν tends almost everywhere (that is for x = x c ) to a Hermitian tensor, so that (ε ν E ν , E ν ) tends almost everywhere to zero.The Lebesgue dominated convergence theorem states that The proof is ended.
It evidences the fact that there is no resonant cyclotron heating in our model.This result is compatible with the literature [18,19,22] where a resonant cyclotron heating is possible, but only in another configuration with wave number parallel to the bulk magnetic field.This yields ∂ x = 0 and so is excluded from our discussion.

The hybrid resonance
This section gathers the main theoretical contributions of this work.The structural hypotheses on the general dielectric tensor for the following work to hold will be emphasized in Assumptions 1, 2 and 3.
The hybrid resonance corresponds to a singularity in the limit problem ν = 0.More precisely it corresponds to the vanishing of ε 0 11 .The singularity occurs when ε 0 11 appears in a denominator: in the expression of the x component of the electric field E 0 1 , see (1.11), as well as in the matrix of the (E 0 2 , W 0 ) system, see (1.13) with ν = 0.
Definition 3. -A point x h is referred to as a local hybrid resonance for the limit ε 0 12 (x h ) = 0. We will moreover assume that Going back to the physical dielectric tensor (1.4), one sees that ε 0 11 (x h ) = 0 corresponds to ).This case is referred to as hybrid regime in the literature, see [18,19,22].Unlike in the cyclotron resonance case, there exist singular solutions at the hybrid resonance.This is known since the seminal work of Budden [5], for an extremely particular case where ε 0 11 = 1 + x/x h .A recent mathematical analysis, performed in [9] with a singular integral equation, has evidenced the role of the limit resonant heating.Our main objective here is to provide additional understanding of the resonant heating by means of a local extraction of the singularity.
Hereafter follows the list of some assumptions on the general regularized tensor (3.1).They are sufficient to carry out the upcoming analysis.

Assumption 1 (Uniqueness and analyticity of the coefficients)
Assume x h is a local hybrid resonance.There exists Λ 1 > 0 such that Note that, for F satisfying this assumption of analyticity, One shall denote the derivative as usual by Assumption 2 (Structure of the dielectric tensor).-The dielectric tensor ε ν is such that As a result, the trace of the matrix M ν (x, θ) defined in (1.13) is identically zero on its domain.Moreover, ε 0 = (ε 0 ) * is a Hermitian matrix.
Also note that as a consequence of this assumption, the determinant of M ν is less singular than expected: since all the entries of the matrix are O((ε ν 11 ) −1 ), the determinant would be expected to be O((ε ν 11 ) −2 ); however (3.5) shows that it actually is only O((ε ν 11 ) −1 ).For the sake of clarity, consider the following notation: so that under Assumption 2, the matrix M ν introduced in Definition 1 reads These quantities are well-defined except possibly at the roots of ε ν 11 .
Remark 1. -In the model case ε ν = ε 0 + iνD with a diagonal dissipation tensor D = I, one can use the notations of [9] so that ε ν 11 = α(x) + iν and ε ν 12 = iδ(x).In this case the coefficients are So c ν does not depend on θ which can be deduced also in the general case from (3.6).
Assumption 3. -The regularized dielectric tensor ε ν (z) is locally analytic with respect to both variables z and ν.
As announced, it is straightforward to verify that these assumptions are satisfied by the physical tensor introduced in (1.4).
Stemming directly from Definition 3 and Assumption 2, a first property of the 2×2 differential system (1.12) is the following: Lemma 3. -Assume x h is a local hybrid resonance.Then for ν = 0 the determinant (1.14) of the dielectric tensor satisfies Proof.-Assumption 2 yields that ε 0 12 (x h ) = ε 0 21 (x h ) since the matrix is Hermitian.Using (3.2) and (3.5), one gets which is indeed a non-zero number in view of (3.3).
Note that moreover due to (3.5), one has The techniques developed in this work are direct consequences of these properties.The study that follows is performed in two steps: the analysis of the regularized problem, with ν = 0, then the study of the limit as ν approaches zero, leading back to the original problem.To this purpose it is crucial to develop uniform tools with respect to the regularization parameter ν, so that their properties still hold as ν goes to zero.

3.1.
A convenient second order equation.-System (1.13) with matrix (3.7) can provide two different second order equations, by elimination of either of the two unknowns.In order to choose which of these two equations to study, consider the following properties of our problem.

1) is such that:
-There exists Λ 2 satisfying 0 < Λ 2 Λ 1 (with Λ 1 defined in Assumption 1) and C > 0 independent of the angle of incidence θ such that -The quantity The upper extra diagonal coefficient b ν does not dominate the diagonal part a ν for vanishing θ results in a singularity in this regime, and as a result the analysis is trickier for the equation on E ν 2 .Therefore the formulation considered is the one on W ν : -At the resonance, the limit equation (3.10) has a "regular singular point" [8,13].Indeed, thanks to Definition 3 and Lemma 4, the 1/c ν term is bounded while the a 2 ν /c ν and b ν terms behave as O((ε ν 11 ) −1 ) at the resonance.The term (a ν /c ν ) is locally bounded.
A consequence of Assumption 3 provides a unique complex root of the coefficient ε ν 11 in a neighborhood of the root of the limit coefficient ε 0 11 , thanks to the open mapping theorem: Lemma 5. -Under Assumptions 1 and 3, there exists A first order expansion is which shows that, choosing adequately Proof.-Define, in order to use the analytical implicit function theorem (Theorem 10.32 of [15]) the function There exists Λ 4 such that it is an analytic function in the ball B((x h , 0), Λ 4 ) ⊂ C 2 .One has ε(x h , 0) = 0, ∂ x ε(x h , 0) = ∂ x ε 0 11 (x h ) = 0, hence we can apply the implicit function theorem.There exists Λ 5 Λ 4 such that, for (x, ν) ∈ B((x h , 0), Λ 5 ), there exists a unique solution of ε(x, ν) = 0, which is denoted by X ν .In addition The assumptions of Definition 3 ensure that The case ε ν 11 = α(x) + iν mentioned in Remark 1 satisfies these hypotheses, with α(x h ) = 0, α (x h ) < 0. For this case one has, at first order, X ν = −iν/α (x h ) + O(ν 2 ), which means that the dominant part of the translated resonance is pure imaginary with non-zero imaginary part for ν = 0.An adequate scaling in (3.10) of the unknown W ν then provides an equation with no first order term.However, as the coefficient that is introduced for this scaling is the square root of c ν , we must use caution.We define the square root of x − X ν as the principal square root on complex numbers (which means that the branch cut is on x ∈ R − ).For the reader's convenience, the explicit expressions for the square root are: Check that −sc ν (x)(x − X ν ) has a strictly positive limit when ν → 0, hence locally in ν the square root of this quantity is well-defined.The rescaled unknown y is defined for x in a real neighborhood of x h by From (3.10) it satisfies on this neighborhood the equation In order to solve Equation (3.14), it is then crucial to isolate the singularity of the solution.We propose here to first understand the singular structure of the equation's coefficient.In the case of a local hybrid resonance, we notice that the most singular term of the coefficient is .
We then consider the following quantity.Definition 6. -Under Assumptions 1, 2 and 3, we define the coefficient-function R ν in B(x h , Λ 3 ) as: With this notation, Equation (3.14) reads Since the function R ν can also be defined by continuity at X ν , as stated in the following result, this form of the equation evidences the singularity of the coefficient.
Remark 3. -Note that the method described below is valid in a more general case than the hybrid resonance introduced in Definition 3, namely when ε 0 11 In this case one considers the limit k ν , for x → X ν , of which is finite.The value of k ν will generate the family of approximate solutions needed for the subsequent study.
Returning to the situation described in Definition 3, one has Lemma 6. -Under Assumptions 1, 2 and 3, there exists where X ν is the translated resonance.This makes sense by continuity for ν small since The proof of this Lemma is given in Section A. Since R ν (X ν ) is bounded, as expected, the singularity is explicit in Equation (3.15): there is a (x − X ν ) −2 singularity, and only if R ν (X ν ) = 0 then there is also a (x − X ν ) −1 singularity.
Illustration of the shift technique, with indication of the branch cuts for a correct definition of the square root.On the left: physical setting in the complex plane.On the right: mathematical setting for the study of the Bessel functions.
Since R ν is bounded and analytic, one can formally shift the equation and the unknown in the complex plane.This is illustrated in Figure 3.1.So we consider the function (3.16) which is a shit of (3.13) and write the equation for this new unknown.It can be justified that W ν , solution of (3.10), admits a convenient extension in a complex neighborhood of the hybrid resonance x h .It provides a rewriting of Equation (3.14) with these new variables.
-The shifted equation writes (3.17) Our strategy is first to construct an explicit quasi-solution of the shifted equation (7) which is more general than the original equation (3.14)-(3.15)because it is a complex equation, so as to provide an explicit quasi-solution of (3.14)- (3.15).The modifications to obtained exact solutions of the original equation (3.15) will be performed in a second step.

3.2.
Freezing and defreezing the coefficient-function. -This section focuses on Equation (3.17) for all ν = 0. We require that the function R(•) = R ν (x + X ν ) is locally analytic with coefficients uniform with respect to ν.We consider the more general from Equation (3.18) will be linked to Bessel's equation, replacing the coefficient function R by its value at the singularity.Indeed, consider the resulting new equation This equation is presented in the chapter on Bessel functions in the classical textbook [1], see Equation (9.1.50).According to this reference, solutions are then expressed in terms of the Bessel functions J 0 and Y 0 .Even if the equation is singular at the origin, its solutions are locally bounded as stressed below.
-If R(0) = 0, a pair of independent solutions of (3.19) can be expressed with Bessel functions as where λ = 2 −R(0).If R(0) = 0, a pair of independent solutions of (3.19) can be expressed with the simpler expression (3.20) and the function admits a infinite converging expansion in powers of z 2 where ψ(1) = −γ and ψ(k −1 for all k ∈ N such that k 2, and γ is the Euler constant. Therefore a more convenient pair of independent solutions of (3.19) also valid for λ = 0 reads Unlike the one proposed in Lemma 7, the pair of independent solutions (3.22) is uniform with respect to λ, for any function R continuous in the neighborhood of the origin.So when applied to the case it will be uniform with respect to ν. Indeed in this case, the parameter λ will depend on ν since then λ = 2 −R ν (X ν ).We will therefore denote by λ ν the value of λ = 2 −R ν (X ν ) to emphasize its dependence with respect to ν.So the second function from Lemma 7 has a logarithmic divergence when λ ν √ z ≈ 0, which means that it diverges if either λ ν ≈ 0 or √ z ≈ 0. On the other hand the logarithmic divergence of The functions J 0 and T 0 are odd, so do not need a branch cut.On the other hand a branch cut is needed for the logarithm function to be defined in a non ambiguous way.We take the branch cut The Bessel based functions, solutions of the equation with frozen coefficient (3.19), provide the formal dominant behavior of the solutions of Equation (3.18).Nevertheless the rigorous justification of this statement is not evident, mainly because the frozen J.É.P. -M., 2017, tome 4 coefficient R(0) is only a first order local approximation of the coefficient function R. Indeed the singularity of the 1/4z 2 term is of second order, so any naive iteration technique aiming at controlling the approximation error between both may diverge like O z −1 .
Overcoming this difficulty is the purpose of the next paragraph, which presents a technical process to control these errors.The idea is to explicit the link between Equation (3.18) and the frozen coefficient equation (3.19).Yet the singularity of the equation requires specific attention, and the process is two-fold: an intermediate equation is introduced to focus first on the singularity, before finally going back to the desired equation, namely (3.18).

Eikonal equation and stretching function. -Consider first another intermediate equation by reintroducing a varying coefficient in Equation (3.19):
Definition 8 (Stretching procedure).-Let ρ be an arbitrary smooth function defined in a complex neighborhood of the origin, such that ρ(0) = 0 and ρ (0) = 0. We will refer to it as the stretching function.This name comes, in particular, from McKelvey [14].
Let u be a given function.The stretching of u is the function u defined by No branch cut of the square root is needed for small |z| since we impose ρ (0) = 0.The equation for the stretched function is the following.Lemma 8. -Assume ρ is a stretching function and u solves the frozen coefficient Equation (3.19).Set s = (ρ ) −1/2 /(ρ ) −1/2 .Then the stretched unknown satisfies Proof.-The proof relies on the scaling of the stretched unknown: the coefficient of the first order derivative term in the computation of u is ((ρ ) −1/2 ) ρ + ((ρ ) 1/2 ) , which happens to be zero.
Note that the function s is a priori bounded around 0 since ρ (0) = 0, so that, like Equation (3.18), Equation (3.25) evidences the singularities of its coefficient.
At this point the stretching function is a tool that will be designed for the approximation of Equation (3.18).Indeed, although Equation (3.25) is different from the initial equation, namely (3.15), an adequate choice of stretching function leads the way back to it.We choose a stretching function as a solution of the following eikonal equation We call this equation the "eikonal" equation for the similar equation for φ if one seeks a solution of the wave equation (∆ − c −2 (x)∂ 2 t 2 )u = 0 under the form u(x, t) = a(x, k)e ik(φ(x)−t) .
Note that (3.26) implies for z, ρ(z) non-zero.As we choose ρ bounded in the neighborhood of 0, only the + sign is relevant.So one can write Use of the conjugate quantity of the numerator yields the equation where F is defined by Equation (3.27) is defined also at z = 0. Since F (0, a) = 2 [R(0)a − R(0)] is welldefined, then F (z, a) makes sense at least for small z and for bounded R.
We now reintroduce R = R ν (• + X ν ) and consider the more general equation in the complex plane where F ν is defined by Note that thanks to the fact that the denominator is 1 for (z, a) = (0, 0), F is analytic in a neighborhood of (0, 0) with coefficients uniform with respect to ν.
Lemma 9 (Solution of the eikonal equation).
The function σ ν = ρ ν /z is solution to the equation The value of F ν (0, a) = 4R ν (X ν ) (a − 1) is a well-defined (i.e., bounded, or nonsingular) real number.Equation (3.30) can be solved in the complex plane using a theorem from Coddington-Levinson [8, Th. 8.1, p. 34].For simplicity consider the initial data σ ν (0) = 1 which yields σ ν (0) = F ν (0, 1) = 0.This choice is arbitrary.Since σ ν (0) = 0, one has the Taylor expansion ρ ν (z) = z + O ν (z 3 ), where the constant in O ν (z 3 ) C|z| 3 is uniform with respect to ν.All other coefficients of the Taylor expansion of the stretching function ρ ν can be easily computed from the integral representation in the complex plane (3.29),where the path of integration is the straight line.The analyticity is provided by the Coddington-Levinson theorem.
Thanks to the stretching function described in Lemma 9, we can now express explicitly the solutions of Equation (3.25), where, in view of Lemma 7 and Equation (3.23), one needs to set Definition 9. -Under Assumptions 1, 2 and 3, consider Λ 7 defined in Lemma 9, and for all |ν| Λ 7 the stretching function ρ ν described in Lemma 9. A pair (U ν , V ν ) of independent solutions of Equation (3.25) is defined for all (ν, z) We adopt the convention (as one did before for the root of −sc ν ) that These functions are uniformly bounded with respect to ν.Then the functions U ν and V ν are bounded in a complex neighborhood of the resonance x h ∈ C, uniformly for small ν: there exists a constant Λ 8 ∈ (0, Λ 7 ) (with Λ 7 defined in Lemma 9) such that for all (ν, z) Moreover, Λ 8 can be determined such that there exists δ > 0 with in the same neighborhood.
Proof.-The property on ρ ν is immediate since σ ν (0) = 1.The uniform boundedness of U ν stems from the boundedness of J 0 (3.20).Concerning V ν , the log λ ν term has been carefully removed, see (3.21), so that the boundedness is achieved even for vanishing λ ν : the remaining log ρ ν (z) is controlled by the √ z contribution that comes from ρ ν (z).
It is then possible to apply the reverse shift.
are solutions of the equation on the real line where Proof.-By definition, all these functions are analytic in the ball B(x h , Λ 8 ) − {X ν } (that is except at X ν ) which is uniform with respect to ν.But, when ρ ν satisfies the eikonal relation (3.26) with R(z) = R ν (z + X ν ), the functions without the shift z → U ν (z) and z → V ν (z) are solutions of (3.25).Therefore the result holds for analytic functions on B(x h , Λ 8 )∩R for 0 < |ν| Λ 3 .So this relation is true everywhere by analytic continuation.The proof is ended.

3.2.2.
General solution with a Duhamel's principle and a limit process.-In order to go back to the equation with a variable coefficient, let us now consider the original problem (3.15).One can rewrite this equation as The last term is a perturbation with respect to Equation (3.32).As usual for such problems, this term is treated by means of the Duhamel's principle, where we will make major profit of the fact that the fundamental solutions of the singular equation (3.32) are bounded.
Lemma 12. -Suppose Assumptions 1, 2 and 3 hold.Consider for all ν a given reference point x * ∈ R close to x h in the ball of analyticity B(x h , Λ 8 ), where Λ 8 is defined in Lemma 10.Then the function is a solution of the integral equation where the matrix As usual for such Volterra integral equations, the initial point x * ∈ R can be chosen arbitrarily.
Proof.-This is the standard procedure of variation of parameters.Consider The solution y can be expressed as a combination of U ν (x − X ν ) and V ν (x − X ν ), with appropriate coefficients A ν (x) and B ν (x) To construct the coefficients A ν (x) and B ν (x), we first assume that So the first derivative of y is The second derivative reads Using (3.32), one gets Since y is solution of (3.38), one gets One can solve now the linear system made of (3.39) and (3.41).To compute the determinant of this linear system consider the Wronskian of U ν and V ν .Elementary calculations and the identity W (J0,Y0) (x) = 2/πx yield By analytic continuation the determinant of the linear system (3.39)-(3.41) is also equal to the same value.
Therefore the solution reads After integration, it yields the representation formula (3.36).The boundedness of the kernel comes from the properties of s ν , U ν and V ν in Lemma 10.It completes the proof.
Under the same conditions, one can define the integral operator K ν ), with a constant vector right hand side.Classically for this Volterra second type integral equation (see [20]), the solution is expressed with the resolvent kernel under the form where C ν (x * ) is the initial condition.Foreseeing the limit process, the following results state the limit properties of ρ ν and M ν , and of other useful quantities, as ν approaches zero.
Then the log and square root terms have limits which depend on the sign of ν: ) , and Proof.-The claim on ρ 0 is evident.The limit of the log term comes from the principal value of the logarithm: log z = a(z) + ib(z), where a(z) = log |z| and b(z) ∈ (−π, π].And the limit of the square root term comes from the principal value of the complex square root: , where θ(z) ∈ (−π, π] is the argument of z. Lemma 14. -Under Assumptions 1, 2 and 3, the matrix M ν defined in (3.37) has a limit in the following sense.For all Λ 9 ∈ (0, Λ 8 ) (where Λ 8 is defined in Lemma 10), one has that -The continuous kernels M + and M − are such that: Consider (for simplicity) a Cauchy data C ν (x * ) which admits a finite (given) limit in ν in the representation (3.44) of the solution of integral equation.Then the functions A ν and B ν also converge in C 0 [x h −Λ 9 , x h +Λ 9 ] to continuous limit functions A ± 0 and B ± 0 .
Proof.-The continuous limits of the kernel M ν are given by application of Lemma 13 since the kernel entries only depend on the functions U ν and V ν which are locally bounded, see (3.31), and continuous.The limit of functions A ν and B ν then directly stem from the limit of (3.44), where the kernel defined by (3.37)-(3.42)-(3.43)has a limit which is bounded.Since the limits are solutions of the limit integral equation with a bounded and continuous kernel, the limits A ± 0 and B ± 0 are continuous.It ends the proof.
It is away from the branch cut of the square root since ) and the value of s recalled in Definition 10.So W 1 is analytic in the neighborhood of x h .Mutatis mutandis the proof is the same for the function W ν 3 .
Considering the logarithmic part in Definition (3.21) of Y 0 , the situation is a little more involved concerning W ν 2 since the limits depend on the sign of ν.
Lemma 16. -Under Assumptions 1, 2 and 3, consider for all 0 < |ν| < Λ 9 the stretching function ρ ν described in Lemma 9.One has where the limit functions are and s defined in Definition 10.The limit functions W + 2 and W − 2 are equal on the interval (x h , x h + Λ 9 ), but are different on the interval (x h − Λ 9 , x h ), where they satisfy the jump relation Proof.-Formula (3.49) is an immediate consequence of the definition of W ν 3 .The second relation comes from (3.45): we note that a similar relation has been proved in [9].

3.3.
Representation of the physical unknowns.-Being particularly careful about the various singularities encountered, one can now obtain meaningful representation formulas for the physical quantities (E 1 , E 2 , W ). Proposition 1. -Assumptions 1, 2 and 3 hold.Consider for all |ν| Λ 9 the stretching function ρ ν described in Lemma 9. Consider for all ν a given reference point x * ∈ R close to x h in the ball of analyticity B(x h , Λ 9 ), as well as a Cauchy data Bν (x * ) given through where the matrix G ν (x * ) is given by When ν goes to 0, this Cauchy data has a finite limit which can be any value in C 2 .Then the unique solution (E ν 1 , E ν 2 , W ν ) of the regularized system (1.9) on B(x h , Λ 9 )∩R is given by Proof.-The expression of is immediate from Lemma 12, since it gives the general solution (3.34) of the corresponding second order equation.It yields the first identity of (3.52).
Concerning the second identity of (3.52), one can start from by means of (1.12) and (3. To compute d dx W ν we directly differentiate all terms in the expression already obtained for W ν (first line of (3.52)).We observe that the derivatives of A ν and B ν vanish since It is then sufficient to plug this expression into (3.53), and to reorganize the sum to get the E ν 2 formula in the second line of (3.52).From (1.9), summing the first equation multiplied by ε ν 11 and the second equation by −ε ν 12 , and from Assumption 2, one gets Even if very simple, this algebra seems to be important since various cancellations of potential singular terms have been performed.It is then sufficient to plug the representation formulas for W ν and d dx W ν in (3.55) to obtain the last part of the claim.The Cauchy data C ν (x * ) is equal to Aν (x * ) Bν (x * ) , where A ν and B ν are known through (3.34)- (3.40).This gives ) .
This linear system can be inverted to obtain y(x * ), y (x * ).Then, use Definition (3.13) for y as well as the equation on W ν , E ν 2 given in Definition 1 to get y(x * ), y (x * ) in terms of W ν (x * ), E ν 2 (x * ), which leads to (3.51).Since A ν (x * ), B ν (x * ) is arbitrary, the Cauchy data is any vector in C 2 .
From Formulas (3.52) it is clear that the convergence with respect to ν is not the same for E ν 2 and W ν on the one hand, and for E ν 1 on the other hand.Indeed the most singular term in E ν 2 and W ν is the logarithm log ρ ν (• − X ν ): the last term in E ν 2 is non-singular since ε ν 11 /ρ ν (• − X ν ) is the ratio of two terms vanishing at order one at X ν -so is non-singular.
So E ν 2 and W ν are bounded in L p loc (x h − Λ 9 , x h + Λ 9 ) for p < ∞ and pass to the limit pointwise except at the singularity x h = lim ν→0 X ν .This result was already obtained for the slab geometry with a completely different technique in [9,Prop. 5.14] and is generalized here.So we state without detail the representation formulas for the limits.
The situation is different for the E ν 1 component of (3.52) since the division by ρ ν (• − X ν ) in the last term yields a singularity as ν approaches zero because This behavior is the same that already demonstrated in [9] using a singular integral equation technique.The difference is that we have now an explicit representation of this singular behavior for small ν.To express the limit of the singular term, the more efficient way is to use principal value and Dirac mass.Since the result is essentially the same as in [9] and all calculations are now evident starting from the representation formula (3.52) we state the result without details of the proof.The notations are that s = sign (i∂ ν ε ν 11 (x)/∂ x ε ν 11 (x)) |x=x h ,ν=0 , D is the Dirac mass and the principal value is defined by Theorem 1. -Suppose Assumptions 1, 2 and 3 hold.Consider for all |ν| Λ 9 the stretching function ρ ν described in Lemma 9, and C ν (x * ) (given by (3.51)) which goes to a finite limit when ν → 0. Consider the unique solution (E ν 1 , E ν 2 , W ν ) of the regularized system (1.9) on B(x h , Λ 9 ).Then the limit of E ν 1 , in D (B(x h , Λ 9 )), Of course the right hand side terms on the first line converge also in the space L p loc (x h − Λ 9 , x h + Λ 9 ) as in the previous proposition, and the pointwise limit holds away from the hybrid singularity x h .
For a, b ∈ B(x h , Λ 9 ) ∩ R\{x h }, one can pass to the limit with respect to ν: Proof.-There are many possibilities to compute this limit value from the previous representation formulas.
First case.-This is a direct consequence on the fact that a and b are on the same side of x h , that the solution is smooth away from x h and that ε 0 (x)X, X = 0 for any complex vector X since ε 0 (x) is Hermitian.
Second case.-Our proof of the result relies on three arguments which are a decomposition of the solution between a regular part and a singular part, an adapted decomposition of the resonant heating Q ν (a, b) where some terms naturally tend to zero, and a careful analysis of the remaining singular part.
Step 1: decomposition of the solution.-Let us decompose the electric field in two parts, a regular part R ν and a singular part S ν .The regular part R ν is regular enough, for example bounded in L 2 (−Λ 9 , Λ 9 ) uniformly with respect to ν and with a pointwise limit for x = x h .The singular part S ν is such that Step 2: decomposition of Q(a, b).-Let us introduce the Hermitian matrix This Hermitian matrix is uniformly bounded in view of all assumptions, that is for a constant uniform with respect to ν.The decomposition writes Step 3: estimate of all terms.-The two first terms of this expression tend to zero with ν in view of (3.60) and (3.62).Let us study the remaining term, which is These two terms are regular and have a natural limit as ν → 0 ± .One has Note that M ν (x) ∈ R. We thus deduce that As the limit exists in L 1 , we use the dominated convergence theorem to deduce that Therefore it is sufficient to pass to the limit for M

The change of variable
After integration where We claim that all terms have a natural limit in (3.66).One has When a < x h < b, the limit of (a − X ν )/ X ν is −∞ when X ν goes to 0 by positive values, and the limit of (a − X ν )/ X ν is +∞ when X ν goes to 0 by negative values.The sign of the imaginary part sign( X ν ) is given by (3.12) We thus deduce Then the value of the resonant heating is independent of D. The same for the pointwise limits of the electric field, the magnetic field and the numerical value of the resonant heating.
Proof.-It is sufficient to realize that the major assumption that was made, that is (3.2)-(3.4), is independent of the exact value of d 1 = 0, and that the other coefficients in D do not appear in the proof of the previous theorem.This is also clear in view of the value of the resonant heating (3.59).
An interesting consequence is that the initial dielectric tensor (1.4) can be replaced by the linear approximation ε 0 + iνD with D defined in (1.6).If one is interested only in the limit value, which is a reasonable assumption for fusion plasmas where ν ≈ 10 −7 may be encountered, this is a valid assumption.
One notices that the positivity condition d 2 1 > |d 2 | 2 that stems from the physics (1.7) is not mandatory for the mathematical results to hold, since this condition was never used in the proof.It is possible to modify even further the dissipation matrix which therefore becomes non-physical.One may consider for example an artificial tensor (3.68) Since only d 1 really matters in ε ν 11 which is the only quantity in the denominators visible in the matrix M ν defined in (1.13), the results of this work are still valid.We did not elaborate on such a matrix in this work since the condition d 3 = d 1 is non-physical and of less interest, nevertheless it must be stressed that the numerical results in next section with tensor D 2 validate this statement.

Numerical illustrations
4.1.Numerics.-We show numerical results which illustrate some of the theoretical results.These results have been computed with the Matlab solver developed by Lise-Marie Imbert-Gérard during her PhD [11].The code solves the system under the form (1.12) with the ode23s subroutine of Matlab.This routine is adapted to stiff problems.
We use three different dielectric tensors that all have the same limit ε 0 for ν = 0.The three dielectric tensors are defined by three different dissipation tensors.The first one denoted as D 0 corresponds to the exact formula (1.4).The second one D 1 = 1 0 0 1 corresponds to the linear approximation ε 0 + iνD 1 .One can notice that D 1 = I corresponds to the limiting absorption principle that was studied in [9].The last is a Hermitian and positive matrix as is required on the physical basis (see (1.6) for example).These three methods satisfy the main hypothesis of Corollary 1 that is (D i ) 11 > 0 for i = 0, 1, 2. In each case, the matrix M ν (x, θ) is constructed and the problem solved numerically in the interval [−5, 5] with arbitrary Cauchy data E ν 2 (5) = 1 and W ν (5) = 0. Notice that the Cauchy data is prescribed on the right boundary inside the non-propagative zone, and is propagated from the right to the left by the Matlab solver.With this, the numerical method naturally captures an exponentially growing term (from the right to the left) which is therefore exponentially decreasing from the left to the right in the propagative region and so models correctly the condition at infinity.Table 4.1.Values of the heating calculated with Formula (3.58) with different ν and different dissipation tensors.As predicted by the theory, the limit resonant heating in the last column is independent of the dissipation tensor.
On the contrary if the numerical Cauchy data is prescribed on the left boundary in the propagative zone, it inevitably captures a non-physical exponentially increasing function at infinity.Finally, note that the Cauchy data E ν 2 (5) = 1 is extremely small with respect to the order of magnitude of the numerical solutions inside the domain.All calculations are performed with different but positive values of the regularization parameter ν > 0. The other parameters are chosen so that x h = 0 is a local hybrid resonance, x < 0 is the propagative region and x > 0 is the non-propagative region.
The convergence as ν goes to the 0 of the heating term appears clearly in Table 4.1.In Figure 4.1 we display the real part of E ν 1 , E ν 2 , W ν for four different values of ν, and for the linear models with D 1 and D 2 .We observe the numerical convergence: even if the results for ν = 0.5 are very different, the functions converge numerically to the same limit for ν = 0.5 × 10 −3 .
The traditional approach is to calculate v in terms of E (which introduces a singularity at x such that ω i c (x) = ω), then to deduce j (or v e ) in terms of E (which introduces a singularity at x such that ω e c (x) = ω) and finally to replace j in Maxwell equations.This structure is not the natural one in the sense that the dielectric tensor that one obtains has singularities when ν goes to 0. Instead of using this approach, we seek the TE field (E 2 (x), B 3 (x) := W (x))e iθy .Let (j 1 , j 2 ) be the components of the electric current transverse to b.The TE field being solution of we have thus to obtain (j 2 , E 1 ) in terms of E 2 , W .It can be checked that j 1 , j 2 , v 1 , v 2 , E 1 are solutions of (ω e c ) 2 = 0.This polynomial in X has two strictly positive roots, one root X − between p 2 and 1, the other root X + greater than 1.For any ω > 0, this gives rise to two hybrid resonances X h − and X h + such that ω e c (X h − ) = ωX + .They give rise to the behavior described in this paper.Other formulas can be sought for.

Figure 1 . 1 .
Figure1.1.X-mode in slab geometry: the domain.In a real physical device an antenna on the wall sends an incident electromagnetic wave toward the plasma.Here, for simplicity, the domain is the half plane {(x, y) ∈ R 2 , x > 0}.The medium is filled with a plasma with dielectric tensor given by (1.2).

Figure 1 . 2 .
Figure 1.2.Idealized parameters for the X-mode equations in slab geometry.The electronic density x → N e (x) is low at the boundary, and increases towards a plateau.The background magnetic field B 0 is here constant for simplicity.

Lemma 4 .
-Assume that x h is a local hybrid resonance.Then, under Assumption 2, the dielectric tensor ε ν (3.

Remark 4 .
-The function J 0 (z) admits a infinite converging expansion in powers of z 2 the second function of(3.22)  is by construction only for √ z ≈ 0. As a result Functions (3.22) are well-defined for z = 0 as ν goes to zero, as well as for ν = 0.The Bessel function (3.21) recasts as Y 0 (z) = T 0 (z) + 2 π log (z/2) J 0 (z) and the second function in (3.22) then reads
7).From (3.6) one has 1/c ν = −ε ν 11 /d ν and a ν /c ν = iθε ν 12 /d ν .So one can write (3.53) Consider on the one hand (3.3)-(3.4)and (3.12) which yield an estimate of the shift of X ν in the complex plane, and on the other hand of the representation formulas for E ν 2 and E ν 1 .Then one has the estimates in L 2(3.60) R ν L 2 (a,b) C and S ν L 2 (a,b) C |ν| ,for some constant C > 0 independent of ν.The regular part is estimated using the integrability in space of the logarithm function.The singular part is estimated thanks to the basic bound b a |1/x + iν| 2 dx π/|ν|.
The two equations are (using a = C 0 c, where C 0 is a given positive constant)W = cE 2 − C 0 cW , E 2 = C 0 cE 2 + bW .We replace E 2 by W /c + C 0 W in the second equation, to obtain (W /c) = (C 2 0 c + b)W .Note that W = J 0 (ρ) is a solution of this equation The external magnetic field is B 0 = B 0 (x)b, b is a unit vector.The so-called fluid equations are given through Newton's lawm e (−iω + ν)v e −e(E + B 0 (x)v e b), m i (−iω + ν)v = Ze(E + B 0 (x)v ∧ b)from which one deduces the equation for the electric current thanks to v e = v− 1 eNe(x) j:(−iω + ν)j = ε 0 ω 2 p (1 + p)(E + B 0 (x)v ∧ b) − eB 0 (x)