\documentclass[JEP,XML,SOM,Unicode]{cedram}
\datereceived{2015-01-07}
\dateaccepted{2016-12-27 }
\dateepreuves{2017-01-10}

\usepackage{mathrsfs}
\let\mathcal\mathscr
\multlinegap0pt
\newenvironment{smallpmatrix}{\big(\begin{smallmatrix}}{\end{smallmatrix}\big)}
\newenvironment{enumeratei}
{\bgroup\def\theenumi{\roman{enumi}}\def\theenumii{\arabic{enumii}}\begin{enumerate}}
{\end{enumerate}\egroup}
\newenvironment{enumeratea}
{\bgroup\def\theenumi{\alph{enumi}}\begin{enumerate}}
{\end{enumerate}\egroup}
\newcommand{\RedefinitSymbole}[1]{%
\expandafter\let\csname old\string#1\endcsname=#1
\let#1=\relax
\newcommand{#1}{\csname old\string#1\endcsname\,}%
}
\RedefinitSymbole{\forall} \RedefinitSymbole{\exists}
\let\epsilon\varepsilon
\def\mto{\mathchoice{\longmapsto}{\mapsto}{\mapsto}{\mapsto}}
\let\dpl\displaystyle
\def\To#1{\mathchoice{\xrightarrow{\textstyle\kern4pt#1\kern3pt}}{\stackrel{#1}{\longrightarrow}}{}{}}

\newcommand{\psfrac}[2]{\sfrac{(#1)}{#2}}
\newcommand{\spfrac}[2]{\sfrac{#1}{(#2)}}

\newcommand{\bA}{\boldsymbol{A}}
\newcommand{\bB}{\boldsymbol{B}}
\newcommand{\bC}{\boldsymbol{C}}
\newcommand{\bE}{\boldsymbol{E}}
\newcommand{\bI}{\boldsymbol{I}}
\newcommand{\bM}{\boldsymbol{M}}
\newcommand{\bR}{\boldsymbol{R}}
\newcommand{\bS}{\boldsymbol{S}}
\newcommand{\bX}{\boldsymbol{X}}
\newcommand{\bk}{\boldsymbol{k}}

\newcommand{\PV}{\mathrm{P.V.}}

\DeclareMathOperator{\sign}{sign}

\newcommand{\be}{\begin{equation}}
\newcommand{\ee}{\end{equation}}

\theoremstyle{plain}
\newtheorem{theor}{Theorem}
\newtheorem{Lemma}{Lemma}
\newtheorem{cor}{Corollary}
\newtheorem{prop}{Proposition}
\newtheorem{assumption}{Assumption}

\theoremstyle{definition}
\newtheorem{df}{Definition}

\theoremstyle{remark}
\newtheorem{rem}{Remark}

\begin{document}
\frontmatter

\title{Solutions to the cold plasma model at~resonances}

\author[\initial{B.} \lastname{Després}]{\firstname{Bruno} \lastname{Després}}
\address{Sorbonne Universités, UPMC Université Paris 06, UMR 7598, Laboratoire Jacques-Louis Lions\\
75252 Paris Cedex 05, France}
\email{despres@ann.jussieu.fr}
\urladdr{https://www.ljll.math.upmc.fr/despres/}

\author[\initial{L.-M.} \lastname{Imbert-Gérard}]{\firstname{Lise-Marie} \lastname{Imbert-Gérard}}
\address{Courant Institute of Mathematical Sciences, New York University\\
251 Mercer Street, New York, NY 10012-1185, USA}
\email{imbertgerard@cims.nyu.edu}
\urladdr{http://www.cims.nyu.edu/~imbertgerard/}

\author[\initial{O.} \lastname{Lafitte}]{\firstname{Olivier}\nobreakauthor\lastname{Lafitte}}
\address{ Université Paris 13, Sorbonne Paris Cité, LAGA (UMR 7539)\\
93430 Villetaneuse, France\\
and\\
DM2S/DIR, CEA Saclay\\
91191 Gif-sur-Yvette cedex, France}
\email{lafitte@math.univ-paris13.fr}

\thanks{Lise-Marie Imbert-Gérard and Bruno Després acknowledge the support of ANR under contract ANR-12-BS01-0006-01. Moreover, this work was carried out within the framework of the European Fusion Development Agreement and the French Research Federation for Fusion Studies. It is supported by the European Communities under the contract of Association between Euratom and CEA. The views and opinions expressed herein do not necessarily reflect those of the European Commission. Lise-Marie Imbert-Gérard also acknowledges the support of the Applied Mathematical Sciences Program of the U.S. Department of Energy under contract DEFG0288ER25053.}

\begin{abstract}
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 differently 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.
\end{abstract}

\subjclass{34M03, 78A25, 34C11}

\keywords{Cold plasma model, resonant Maxwell's equations, hybrid resonance, eikonal equation, singular solutions of ODEs}

\alttitle{Solutions singulières résonantes pour un plasma}

\begin{altabstract}
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 diffé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écifique de la singularité.
\end{altabstract}

\altkeywords{Modèle du plasma froid, équations de Maxwell résonantes, résonance hybride, équation eikonale, solutions singulières d'EDO}

\maketitle
\tableofcontents
\mainmatter


\section{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)|\geq 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 \cite{jmpa}, 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 \cite{jmpa}, see also a recent alternative technique in \cite{campi}.

In time harmonic regime and slab geometry in dimension
two, for $(x,y)\in \mathbb 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 $\omega^2/c^2=1$, the latter equations reduce
to
\begin{equation}\label{eq:000}
\begin{cases}
{}-\partial^2_{y}E_{1} +\partial^2_{xy}E_{2} -\epsilon_{11}
E_{1} -\epsilon_{12} E_{2} =0, \cr
\phantom{{}-{}}\partial^2_{xy}E_{1} -\partial^2_{x}E_{2} -\epsilon_{21}
E_{1} -\epsilon_{22} E_{2} =0.
\end{cases}
\end{equation}
The bulk magnetic field $ \bB_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 $\omega>0$.
As described in \cite{stix,swanson,CW,jmpa}, a first order expansion
in the cold plasma model yields the following dielectric tensor
\be \label{eq:hyp01}\arraycolsep3.5pt
{\epsilon}=
\begin{pmatrix}
1-\dfrac{\omega_p^2}{
\left(\omega^2-\omega_c^2\right)}&i\, \dfrac{\omega_c \omega_p^2}{\omega \left(\omega^2-\omega_c^2\right)} \\[10pt]
-i\, \dfrac{\omega_c \omega_p^2}{\omega \left(\omega^2-\omega_c^2\right)}
&1-\dfrac{\omega_p^2}{
\left(\omega^2-\omega_c^2\right) }
\end{pmatrix},
\ee
where the
cyclotron frequency is defined by
$
\omega_c= \sfrac{e |\bB_0|}{m_e}
$ while the plasma frequency is defined by
$
\omega_p=\sqrt{\sfrac{e^2 N_e}{\epsilon_0m_e} }
$. Here $e$ is the electron charge, $m_e$ is the electron mass and $\epsilon_0$ is the vacuum permittivity.
The physical validity
of this model is discussed
in \cite{branbilla}[p. 166]
for waves such that
$
|\bk_\parallel | v_{th}\ll |\omega - n\omega_c|$ and $
|\bk_\perp |v_{th}\ll\omega_c,
$
where
$\bk_\parallel $ and~$\bk_\perp$ are the parallel and perpendicular components---with respect to the bulk magnetic field $\bB_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 $\omega$ 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
$ \bB_0= \bB_0(x)$ and $N_e=N_e(x)$.
It~corresponds to a wall facing a stratified plasma as depicted in Figure
\ref{figure:0},
as in a Tokamak where an antenna
sends an electromagnetic wave toward the plasma, in order to probe
or heat the plasma.
\begin{figure}[ht]
\begin{center}
\begin{tabular}{ccc}
\includegraphics[width=6.5cm,angle=0]{./despres-et-al_figs/dispo}
\end{tabular}
\end{center}
\caption{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)\in\mathbb R^2,\,x>0\}$.
The medium is filled with a plasma with dielectric
tensor given by \eqref{eq:hyp01}.}
\label{figure:0}
\end{figure}
The heating of a magnetized fusion plasma with such devices
is an important issue for the ITER project \cite{iter}.
Idealized coefficients representing the main features of
the physical coefficients at the wall of a Tokamak are
described in Figure~\ref{figure:1}.
More realistic plasma parameters may be much more oscillating
as illustrated in the numerical studies \cite{dscdh,gusakov}.

\begin{figure}[ht]
\begin{center}
\begin{tabular}{ccc}
\includegraphics[width=6.5cm,angle=0]{./despres-et-al_figs/physcoeff}
\end{tabular}
\end{center}
\caption{Idealized
parameters for the X-mode equations in slab geometry.
The electronic density $x\mto N_e(x)$ is low at the boundary,
and increases towards a plateau.
The background magnetic field $B_0$ is here constant for simplicity.}
\label{figure:1}
\end{figure}

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 \cite{stix,swanson} or to \cite{CW} for a historical
paper on the topic. A recent reference is~\cite{dumont}; it also provides a list of additional references.
Singular solutions of wave equations appear in other applications, such as cloaking \cite{weder2} and
metamaterials \cite{chen:lipton,bonnet2}.
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 \cite{brambilla2}.
Describing the behavior of the solutions as the friction parameter $\nu>0$ approaches zero will evidence the singularity of the solutions to the limit problem for $\nu=0$. Indeed one obtains solutions which may behave like
$$
E_1^\nu\approx \frac{c}{rx+i\nu} \quad\mbox{with } r,c\neq 0.
$$
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 \eqref{eq:hyp01}. This crucial structure will be highlighted in Section \ref{sec:hyb}.

The complex pulsation
$
{\omega}_\nu=\omega +i \nu$ shifts the exact pulsation in the upper half-plane by the factor
$\nu>0$.
With this shift, Equations \eqref{eq:000}--\eqref{eq:hyp01}
become
\begin{equation} \label{eq:00}
\begin{cases}
-\partial^2_{y}E_{1}^\nu +\partial^2_{xy}E_{2}^\nu -\epsilon_{11}^\nu
E_{1}^\nu -\epsilon_{12}^\nu E_{2}^\nu =0, \\[3pt]
\phantom{{}-}\partial^2_{xy}E_{1}^\nu -\partial^2_{x}E_{2}^\nu -\epsilon_{21}
^\nu
E_{1}^\nu -\epsilon_{22}^\nu E_{2}^\nu =0.
\end{cases}
\end{equation}
with
the regularized dielectric tensor
\be \label{eq:hyp1}\arraycolsep3.5pt
{\epsilon}^\nu=
\begin{pmatrix}
1-\dfrac{{\omega}_\nu\omega_p^2}{
\omega
\left(\omega_\nu^2-\omega_c^2\right)}&i\, \dfrac{\omega_c \omega_p^2}{\omega \left(\omega_\nu^2
-\omega_c^2\right)} \\[10pt]
-i\, \dfrac{\omega_c \omega_p^2}{\omega \left(\omega_\nu^2-\omega_c^2\right)}
&1-\dfrac{{\omega}_\nu \omega_p^2}{
\omega\left(\omega_\nu^2-\omega_c^2\right) }
\end{pmatrix},
\ee
where $\omega_p$ and $\omega_c$ are still respectively the plasma and cyclotron frequencies, and
the regularized dielectric tensor is independent of the $ y$-direction.
The physical justification of ${\epsilon}^\nu$ can be found in plasma physics textbooks \cite{branbilla,stix}.

For a recent mathematical treatment of these equations by means of
singular integral equations, refer to \cite{jmpa}.
In fusion plasmas, the value of the friction parameter can be extremely
small.
For example, relative values of $\sfrac{\nu}{\omega}\approx 10^{-7}$ are common.
This extremely small value shows that the frictionless limit regime $\nu\to 0$ is relevant
for some fusion applications.
In such regimes the formal first order approximation of the dielectric tensor reads
\be \label{eq:princip}
{\epsilon}^\nu=
{\epsilon}^0+
i\nu
\mathbb D + O(\nu^2),
\ee
with
\be \label{eq:dede}\arraycolsep3.5pt
\epsilon^0=
\begin{pmatrix}
1-\dfrac{\omega_p^2}{\omega^2-\omega_c^2}&i\, \dfrac{\omega_c \omega_p^2}{\omega \left(\omega^2-\omega_c^2\right)} \\[10pt]
-i\, \dfrac{\omega_c \omega_p^2}{\omega \left(\omega^2-\omega_c^2\right)}&1-\dfrac{\omega_p^2}{\omega^2-\omega_c^2} \end{pmatrix}
\qquad\text{and}\quad\arraycolsep1.5pt
\mathbb D= \left(
\begin{array}{ccc}
\lambda_1 &-i \lambda_2 \\
i \lambda_2 & \lambda_1
\end{array}
\right).
\ee
The dissipation tensor $\mathbb D$
accounts for the underlying physical dissipation. Its coefficients are
\be \label{eq:dede22}
\lambda_1=\frac{\omega_p^2\left(\omega^2 + \omega_c^2 \right) }{
\omega \left(\omega^2-\omega_c^2 \right)^2 }>0
\quad\text{and}\quad\lambda_2=\frac{2 \omega_c \omega_p^2 }{\left(\omega^2-\omega_c^2
\right)^2 }\in \mathbb R.
\ee
Since $\lambda_1-\lambda_2=
\sfrac{\omega_p^2 }{
\omega \left(\omega+\omega_c \right)^2 }>0$, the matrix
$\mathbb D=\mathbb D^*>0$ is indeed positive.
As~previously mentioned, and stressed
in Figure \ref{figure:1}, the coefficients
of these tensors
are a priori non-constant functions of the horizontal variable $x$.

The structure of the limit tensor \eqref{eq:hyp01} shows two different resonances in the limit model $\nu=0$, defined thanks to the dispersion relation, see \cite{stix}.
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
$$
\omega_c(x)=\omega,
$$
where $\omega_c$ is a function of the $x$ variable, a priori non-constant.
We will show that mathematical solutions of Maxwell's equations \eqref{eq:00} 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 $\epsilon^0$, that is to say
$$
\omega_p^2(x)+ \omega_c^2(x)=\omega^2.
$$
It defines the \emph{hybrid resonance}. In this regime the dielectric tensor is non-singular since the denominator
is non-zero, \ie 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 $
\omega_p^2(x_h)+ \omega_c^2(x_h)=\omega^2 $: typically
\[
E_1^\nu(x)=\frac{g_\nu(x)}{\omega_\nu-\sqrt{ \omega_p^2(x)+ \omega_c^2(x) }},
\]
where $g_\nu$ is a function which is proved to be regular enough so that one can pass to the limit
$\nu\to 0^+$. One gets \cite{jmpa}
\begin{multline*}
\label{eq:eess}
E_1^+=\lim_{\nu \to 0^+}E_1^\nu=
\alpha D(\cdot -x_h)+ \PV\frac{g_0(\cdot)}{\omega -\sqrt{ \omega_p^2(\cdot)
+ \omega_c^2(\cdot) }}, \\ D\text{ is the Dirac mass and }\alpha \in \mathbb C.
\end{multline*}
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 \cite{jmpa}, a useful quantity
to characterize the singular behavior of the solutions is
the \emph{resonant heating} of the plasma. Up to a non-essential factor $\omega$, the resonant heating
\cite{branbilla,jmpa}
is defined, for any $(a,b)$ such that $-L\leq a<x_h<b\leq \infty$,
by
\be \label{eq:Q}
\mathcal Q=\lim_{\nu \to 0^+}
\biggl(\Im
\int_a^b \left(\epsilon^\nu \bE^\nu, \bE^\nu \right) dx \biggr),\qquad \bE^\nu=\left(
\begin{array}{c}
E_1^\nu \\
E_2^\nu
\end{array}
\right),
\ee
where the Hermitian product
of two complex vectors $\bA$ and $\bB$ is denoted as $\left(\bA, \bB\right)= \bA\cdot \overline{\bB}$.
This formula has been justified in \cite{jmpa} where
it is shown that $\mathcal Q>0$ in case of resonant heating
using the simplification $\epsilon^\nu=\epsilon+i\nu \bI$, \ie in the $\nu$ expansion \eqref{eq:princip} of $\epsilon^{\nu}$ the dissipation tensor is $\mathbb D=\mathbb I$. This simplification
is mathematically justified since it corresponds to the \emph{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
$W^\nu=\partial_xE^\nu_{2}-\partial_yE^\nu_{1}$ which is, after linearization around the bulk
magnetic field, the parallel
component of the magnetic field. System \eqref{eq:00} is then equivalent to the following first order system
$$
\begin{cases}
\begin{array}{rcl}
\phantom{{}-{}}\partial_yW^\nu -\epsilon_{11}^\nu E_{1}^\nu -\epsilon_{12}^\nu E_{2}^\nu &=0, \cr
{}-\partial_xW^\nu -\epsilon_{21}^\nu E_{1}^\nu -\epsilon_{22}^\nu E_{2}^\nu &=0, \cr
W^\nu \,+\,\partial_yE_{1}^\nu -\partial_xE_{2}^\nu &=0.
\end{array}
\end{cases}
$$
Since the coefficients $\epsilon_{ij}^\nu$ depend only on $x$, one can perform a Fourier transform in $y$. Denoting by $i\theta$ the Fourier coefficient, it yields
\be\label{eq:sysinit}
\begin{cases}
\begin{array}{rcl}
\phantom{{}-{}d}i\theta W^\nu -\epsilon_{11}^\nu E_{1}^\nu -\epsilon_{12}^\nu E_{2}^\nu &=0, \cr
{}-\dfrac d {dx}W^\nu -\epsilon_{21}^\nu E_{1}^\nu -\epsilon_{22}^\nu E_{2}^\nu &=0, \cr
W^\nu\,+\,i\theta E_{1}^\nu - \dfrac d {dx}E_{2}^\nu &=0.
\end{array}
\end{cases}
\ee
The Fourier coefficient $\theta$ can be considered as a frozen parameter
from now on, and~$\partial_x$ will be replaced by $\sfrac d {dx}$ denoting the derivative with respect
to $x$.
In this system, the derivative with respect to $x$
appears only on $W^\nu$ and $E_{2}^\nu$.
The energy integral \eqref{eq:Q} can be expressed
for any $(a,b)$ such that $ -L\leq a<b<\infty$ using the identity
\begin{align*}
\int_a^b \left(\epsilon^\nu \bE^\nu, \bE^\nu \right) dx
&=\int_a^b \Bigl({ i\theta W^\nu } \overline{E_{1}^\nu} - \frac d {dx}W^\nu \overline{E_{2}^\nu} \Bigr) dx\\
&= \int_a^b \Bigl({ i\theta W^\nu } \overline{E_{1}^\nu} + W^\nu \frac d {dx} \overline{E_{2}^\nu} \Bigr) dx-
W^\nu (b) \overline{E_{2}^\nu} (b) + W^\nu (a) \overline{E_{2}^\nu} (a)\\
&= \int_a^b \left| { W^\nu } \right|^2 dx-
W^\nu (b) \overline{E_{2}^\nu} (b) + W^\nu (a) \overline{E_{2}^\nu} (a).
\end{align*}
It leads to
\be \label{eq:enbal}
\Im
\int_a^b \left(\epsilon^\nu \bE^\nu, \bE^\nu \right) dx
= - \Im\left(
W^\nu (b) \overline{E_{2}^\nu} (b) \right) + \Im\left(W^\nu (a) \overline{E_{2}^\nu} (a) \right),
\ee
which expresses a balance of energy.

It is convenient to consider a differential system with two unknowns $E_2^\nu$ and $W^\nu$, while $E_1^\nu$ can be computed thanks to the first equation of System \eqref{eq:sysinit} by
\be \label{eq:extract}
E_{1}^\nu=\frac{i\theta W^\nu-\epsilon_{12}^\nu E_{2}^\nu}{\epsilon_{11} ^\nu}.
\ee
Eliminating $E_1^{\nu}$, the reduced system then reads
\be \label{eq:fond1}
\frac d {dx}\begin{pmatrix}
E_{2}^\nu(x) \cr
W^\nu (x)
\end{pmatrix}=
\bM_\nu(x,\theta)
\begin{pmatrix}
E_{2}^\nu(x) \cr
W^\nu(x)
\end{pmatrix},
\qquad -L<x,
\ee
where the matrix is defined as follows.

\begin{df}\label{def:M}
The matrix $\bM_\nu(x,\theta) $
is constructed
from the dielectric tensor
by
\be \label{eq:mm}
\bM_\nu(\cdot,\theta)=
\begin{pmatrix}
-\dfrac{i\theta \epsilon_{12}^\nu}{\epsilon_{11}^\nu}&
\displaystyle
1-\dfrac{\theta^2}{\epsilon_{11}^\nu}
\vspace{2mm}
\cr
-\dfrac{d^\nu}{\epsilon_{11}^\nu}
&
-\dfrac{i\theta \epsilon_{21}^\nu}{\epsilon_{11}^\nu}
\end{pmatrix},
\ee
where the coefficient $d^\nu$ is the determinant of the dielectric
tensor $\epsilon^\nu$:
\begin{equation} \label{eq:dnunu}
d^\nu=\epsilon_{11}^\nu \epsilon_{22}^\nu-\epsilon_{12}^\nu\epsilon_{21}^\nu.
\end{equation}
\end{df}

The determinant $\mathfrak D_\nu(\cdot, \theta)$ of the matrix
can be simplified as follows
\be \label{eq:detdet}
\begin{split}
\mathfrak D_\nu(\cdot, \theta) =\det \bM_\nu(\cdot,\theta) &=-\Bigl(\frac{i\theta\epsilon_{21}^\nu }{\epsilon_{11}^\nu }
\Bigr)^2+\frac{d^\nu}{\epsilon_{11}^\nu}
\Bigl(1-\frac{\theta^2}{\epsilon_{11} ^\nu}\Bigr)\\
&=\frac{d^\nu}{\epsilon_{11}^\nu}-\theta^2 \frac{ \epsilon^\nu_{22} }{ \epsilon^\nu_{11} }= \frac{d^\nu}{\epsilon_{11}^\nu}-\theta^2.
\end{split}
\ee
It is already foreseeable that the roots of $\epsilon_{11}^\nu$ 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~\eqref{eq:fond1}
when the entries of $\bM_\nu$ 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:

\subsubsection*{Lemma \ref{heat:cyc}}
In our geometry, the cyclotron singularity of the dielectric tensor, \hbox{defined} by
$\omega=\omega_c$, is only an apparent singularity in the sense that the electromagnetic field
is bounded and the resonant heating is zero
$\mathcal Q_\mathrm{cyc}=0$.

\subsubsection*{Theorems \ref{theore:fin0} and \ref{theore:fin}}
The hybrid resonance,
defined by $\omega^2=\omega_c^2+\omega_p^2$, yields a singular horizontal
component of the electric field, \ie $|(x-x_*)E_1^{\nu=0}|\geq c>0$ uniformly
with respect to $x$ in a neighborhood of $x_*$, and a generic positive resonant heating
$\mathcal Q_\mathrm{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 \cite{jmpa}.

\subsubsection*{Corollary \ref{cor:last}}
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 $\nu=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 \cite{jmpa} the limit was established only for $\mathbb D =\mathbb I$.

\subsubsection*{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.

\subsubsection*{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 \cite{branbilla,zizi} for the cyclotron resonance and the hybrid resonance.

\medskip
The work is organized as follows.
The cyclotron singularity of the dielectric tensor is studied
in Section \ref{sec:cyc}, where the electromagnetic field is shown to stay bounded in the neighborhood of this singularity. The hybrid resonance is studied in detail in Section~\ref{sec:hyb},
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 \ref{sec:num}.
Some technical material is postponed to the appendix together with
a particular analytical solution and the generalization to multi-species.

We will use the notations $\mathbb R^*=\mathbb R-\{ 0\}$ and $\mathbb C^*=\mathbb C-\{ 0\}$.

\subsubsection*{Acknowledgements}
The authors deeply thank the anonymous reviewers for their remarks, which helped a lot to improve the present work.


\section{The cyclotron resonance}
\label{sec:cyc}

The cyclotron singularity corresponds to a formal singularity in the denominator
of \eqref{eq:hyp1}. We consider the situation where there
exists at least one $x_c$ in the domain such that
$
\omega_c(x_c)=\omega$.
For the simplicity of the analysis, we will assume that $\frac{d}{dx}\omega_c(x_c)\neq 0$ so that
$x_c$ is isolated.
Clearly the singularity at $x_c$ in the dielectric tensor \eqref{eq:hyp1}
in the limit case $\nu=0^+$ disappears in the regularized case $\nu\neq 0$. In order to focus on the cyclotron resonance, it has to be away from the hybrid resonance, that is $\omega_p(x_c)\neq 0$.

\begin{df}
A point $x_c$ is referred to as an isolated cyclotron singularity if
$\omega_c(x_c)= \omega$ and $\frac{d}{dx}\omega_c(x_c)\neq 0$.
\end{df}

Note that such a singularity does not depend on $\nu$, since $\omega_c$ does not depend on $\nu$.
Let us now consider
the reduced system \eqref{eq:fond1}, where the system's matrix $\bM_\nu(\cdot,\theta)$ is defined with the physical dielectric tensor \eqref{eq:hyp1}. The system is well-defined at a cyclotron resonance for all $\nu\in\mathbb R$, and one can pass to the limit in the following sense.

\begin{Lemma}\label{lem:E2W}
Assume that the magnetic field $\bB_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^\nu,W^\nu)$ solves the reduced system \eqref{eq:fond1} on $[x_c-r,x_c+r]$ for all $\nu\in\mathbb R$, with a given Cauchy data which admits a finite limit when $\nu$ goes to $0$ written as
$(E_2^\nu,W^\nu)(x_c-r)\to (\alpha,\beta)$.

Then $\left(E_{2}^\nu,W^\nu \right) $ tends to a bounded limit $\left(E_{2}^0,W ^0 \right) $
as $\nu$ tends to zero, and
$(E_2^0,W^0)$ solves the limit system \eqref{eq:fond1} on $[x_c-r,x_c+r]$
with $\nu=0$, with the same Cauchy data
$(E_2^0,W^0)(x_c-r)=(\alpha,\beta)$.

\end{Lemma}
\begin{proof}
Definition \ref{def:M} with the physical tensor \eqref{eq:hyp1} yields
after some careful calculations
\begin{multline}
\label{eq:matm0}
\bM_\nu(x,\theta)\\
= \dfrac{1}{
\alpha_\nu(x)
}
\begin{pmatrix}
\theta \omega_c\omega_p^2&
(1-\theta^2) \omega ({\omega_\nu}^2-\omega_c^2)-{\omega_\nu}\omega_p^2
\vspace{2mm}
\\
-\omega ({\omega_\nu}^2-\omega_c^2)+2{\omega_\nu}\omega_p^2
-\dfrac{\omega_p^4}{\omega}&
-\theta \omega_c\omega_p^2
\end{pmatrix}
(x,\theta),
\end{multline}
where $\alpha_\nu(\cdot)=\omega ({\omega_\nu}^2-\omega_c^2)-{\omega_\nu}\omega_p^2$.
The determinant \eqref{eq:detdet} of this matrix can be recast as
\[
D_\nu(\cdot, \theta)
= \frac{d^\nu}{\epsilon_{11}^\nu}-\theta^2= \frac{\omega({\omega_\nu}^2-\omega_c^2)-2{\omega_\nu}\omega_p^2 +\sfrac{\omega_p^4}{\omega} }
{
\alpha_\nu}(\cdot)-\theta^2.
\]
The limit value of $\alpha_\nu(\cdot)$ as $\nu$ approaches $0^+$
is
$
\alpha_0(\cdot)=\omega \left({ \omega}^2-\omega_c^2- \omega_p^2\right)$.
Since the cyclotron resonance is not a hybrid resonance, that is
$|\omega_p|\geq C>0$ uniformly on $(x_c-r,x_c+r)$, one has that
\be \label{eq:limalpha}
\alpha_0(x_c)=-\omega \omega_p^2 \neq 0,
\ee
since $\omega=\omega_c(x_c)\neq 0$ as consequence of the normalization $\omega^2/c^2=1$. In this situation the matrix $\bM_\nu(\cdot,\theta)$
is locally bounded and smooth. So the solutions
of the reduced system \eqref{eq:fond1}
are smooth locally around $x_c$ and one can pass to the limit
in this system.
\end{proof}

Note that there is still a particular behavior in the neighborhood of this cyclotron resonance, known as ``turning point behavior'' \cite{co:ed,LWZ}. However, this does not induce any singular regime, except in the limit $\omega\to +\infty$, which is not considered here.
The consequences of the representation \eqref{eq:matm0} on physical quantities are characterized in the following result.

\begin{Lemma} \label{heat:cyc}
Under the assumptions of Lemma \ref{lem:E2W},
the component $E_{1}^\nu $ tends pointwise to a bounded limit $E_{1} ^0$ as $\nu$ tends to zero and
the value of the resonant heating is zero:
\be \label{eq:heatcyc}
\lim_{\nu\to 0^+}
\left(
{\Im}
\int_{x_c-r}^{x_c+r}
\left(\epsilon^\nu \bE^\nu, \bE^\nu \right) dx \right)
=0.
\ee
\end{Lemma}

\begin{proof}
The horizontal component of the electric field can be evaluated thanks to
\eqref{eq:extract}. The coefficients in \eqref{eq:extract} can be expressed as
$$
\frac{i\theta }{\epsilon_{11} ^\nu}
=
\frac{i\theta \omega \left(\omega_\nu^2-\omega_c(\cdot)^2 \right) }{\alpha_\nu(\cdot, \theta)}
\quad\text{and}\quad
\frac{- \epsilon_{12}^\nu }{\epsilon_{11} ^\nu}
= \frac{- i \omega_c(\cdot)
\omega^2_p }{\alpha_\nu(\cdot, \theta)}.
$$
Using \eqref{eq:limalpha}, the limit $\nu=0$ value of these coefficients at $x_c$ is
\be \label{eq:limcoeff}
\left(\frac{i\theta }{\epsilon_{11} ^0}\right) \left(x_c \right)
=0
\quad\text{and}\quad
\left(\frac{- \epsilon_{12}^0 }{\epsilon_{11} ^0}\right) \left(x_c \right)
= i.
\ee
So one can pass to the limit in \eqref{eq:extract} since the coefficients are bounded
in the interval $[x_c-r, x_c+r]$. Therefore the horizontal
part of the electric field $E_1^\nu$ admits a bounded and smooth limit $E_1^0$.

Consider now the heating term.
On the one hand we notice that
the electric field can be split into two parts
$$
\bE^\nu=
E_2^\nu
\left(
\begin{array}{c}
i \\
1
\end{array}
\right)
+
\left(
\begin{array}{c}
E_1^\nu -i E_2^\nu
\\
0
\end{array}
\right),
$$
the first vector
being an eigenvector of the dielectric tensor:
$$
\epsilon^\nu \left(
\begin{array}{c}
i \\
1
\end{array}
\right)=
\mu
\left(
\begin{array}{c}
i \\
1
\end{array}
\right).
$$
The corresponding eigenvalue $\mu = 1-\sfrac{\omega_p^2}{\omega\left(\omega_\nu + \omega_c \right)}$ is bounded so non-singular at~$x_c$. The second vector is $E_1^\nu-i E_2^\nu= (\omega_\nu -\omega_c) F^\nu$,
with
\[
F^\nu= \left[i\theta(\omega_\nu+\omega_c) \omega W^\nu
-
i((\omega_\nu+\omega_c) \omega -\omega_p^2)E_2^\nu
\right ]/\alpha^\nu,
\]
admits a finite limit for $\nu\to 0$,
due to
\eqref{eq:extract}
and \eqref{eq:limcoeff}. So the second vector
is
$O\left(\nu +|x-x_c|\right)$, which means that it counterbalances the singularity of the dielectric tensor. This shows that $\epsilon^\nu \bE^\nu$ is bounded uniformly on $[x_c-r, x_c+r]$ with respect to~$\nu$. Since $ \bE^\nu$ is also uniformly bounded, it yields that $\left(\epsilon^\nu \bE^\nu, \bE^\nu \right)$ is bounded uniformly on $[x_c-r, x_c+r]$ with respect to $\nu$.

On the other hand the tensor
$ \epsilon^\nu$ tends almost everywhere (that is for $x\neq x_c$) to a Hermitian
tensor, so that
$\Im
\left(\epsilon^\nu \bE^\nu, \bE^\nu \right)$
tends almost everywhere to zero.
The Lebesgue dominated convergence theorem states that
$$
\lim_{\nu\to 0^+}
\left(
\Im
\int_{x_c-r}^{x_c+r}
\left(\epsilon^\nu \bE^\nu\cdot \bE^\nu \right) dx \right)
=0.
$$
The proof is ended.
\end{proof}

It evidences the fact that there is no resonant cyclotron heating in our model.
This result is compatible with the literature \cite{stix,swanson,CW}
where a resonant cyclotron heating
is possible, but only in another configuration
with wave number parallel to the bulk magnetic field.
This yields $\partial_x =0$ and so is excluded from our discussion.

\section{The hybrid resonance} \label{sec:hyb}

This section gathers the main theoretical contributions of this work.
The structural hypotheses on the general dielectric tensor
\be \label{eq:epsilonmu}
\epsilon^\nu(x) =
\begin{pmatrix}
\epsilon_{11} ^\nu(x)& \epsilon_{12} ^\nu(x)\\[3pt]
\epsilon_{21}^\nu(x) & \epsilon_{22}^\nu(x)
\end{pmatrix}
\in \mathbb C^{2\times 2}, \quad \forall \nu \in\mathbb R,
\ee
for the following work to hold will be emphasized in Assumptions \ref{ass:UnA}, \ref{ass:SDT} and \ref{ass:eps}.

The hybrid resonance corresponds to a singularity in the limit problem $\nu=0$. More precisely it corresponds to the vanishing of $\epsilon_{11}^0$. The singularity occurs when~$\epsilon_{11}^0$ appears in a denominator: in the expression of the $x$ component of the electric field~$E_1^0$, see \eqref{eq:extract}, as well as in the matrix of the $(E_2^0,W^0)$ system, see \eqref{eq:mm} with $\nu=0$.

\begin{df}\label{def:LR}
A point $x_h$ is referred to as a local hybrid resonance for the limit $\nu=0$ of
the system \eqref{eq:00}--\eqref{eq:epsilonmu} if the limit diagonal coefficient of
$\epsilon_{11}^\nu(x)$
is such that
\begin{equation} \label{eq:dnu2}
\epsilon_{11}^0(x_h)=0\quad\text{and}\quad\partial_x\epsilon_{11}^{0}(x_h)\in \mathbb R^{*},
\end{equation}
and if the extra-diagonal part is locally non-zero
\begin{equation} \label{eq:dnu2:rab}
\epsilon_{12}^0(x_h)\neq 0.
\end{equation}
We will moreover assume that
\begin{equation} \label{eq:dnu2.2}
i \partial_\nu\epsilon_{11}^{0}(x_h)\in \mathbb R^{*}.
\end{equation}
We denote by $s\in \{-1,1\}$ the sign of $-\partial_x\epsilon^0_{11}(x_h)$.
\end{df}
Going back to
the physical dielectric tensor \eqref{eq:hyp1}, one sees that $ \epsilon_{11}^0(x_h)=0$ corresponds to
$$
\omega^2=\omega_c^2(x_h)+ \omega_p^2(x_h).
$$
This case is referred to as hybrid regime in the literature, see \cite{stix,swanson,CW}. Unlike in the cyclotron resonance case, there exist singular solutions at the hybrid resonance.
This is known since the seminal work of Budden
\cite{budden}, for an extremely particular case where $\epsilon_{11}^0=1+\sfrac{x}{x_h}$. A recent mathematical analysis, performed in \cite{jmpa}
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 \eqref{eq:epsilonmu}. They are
sufficient to carry out the upcoming analysis.

\begin{assumption}[Uniqueness and analyticity of the coefficients]
\label{ass:UnA}
Assume $x_h$ is a local hybrid resonance. There exists $\Lambda_1>0$ such that
\begin{itemize}
\item the dielectric tensor $\epsilon^{\nu} $ is analytic in the ball
$x\in B(x_h,\Lambda_1)$ for all $\nu \in [-\Lambda_1,\Lambda_1]$,
\item $x_h$ is the unique root of $\epsilon_{11}^0$ in $B(x_h,\Lambda_1)$.
\end{itemize}
\end{assumption}

Note that, for $F$ satisfying this assumption of analyticity, $F(z)=\sum_{n\geq 0} a_n(z-x_h)^n$, with
$a_n\in \mathbb C^{2\times 2}$ a matrix with the summability condition
$\sum \vert a_n\vert R^n<+\infty$ for $R<\Lambda_1$. One shall denote the derivative as usual
by $F'(z)= \sum_{n\geq 0} na_n(z-x_h)^{n-1}$.
The chain rule reads
$\frac{d}{dx}(F(\phi(x)))=F'(\phi(x))\frac{d\phi}{dx}(x)$
for $\phi$ a function from an interval $I \subset \mathbb R$ to~$\mathbb C$.

\begin{assumption}[Structure of the dielectric tensor]
\label{ass:SDT}
The dielectric tensor $\epsilon^{\nu} $ is such that
\begin{equation}\label{sym:coeff:eps}
\epsilon_{11}^\nu= \epsilon_{22}^\nu \quad \text{and} \quad
\epsilon_{12}^\nu = -\epsilon_{21}^\nu, \qquad
\forall \nu \in [-\Lambda_1,\Lambda_1].
\end{equation}
As a result, the trace of the matrix $\bM_\nu (x,\theta)$ defined in \eqref{eq:mm} is identically zero on its domain.
Moreover, $\epsilon^0=(\epsilon^0)^*$ is a Hermitian matrix.
\end{assumption}

Also note that as a consequence of this assumption, the determinant of $\bM_\nu$ is less singular than expected: since all the entries of the matrix are $O((\epsilon_{11}^\nu)^{-1})$, the determinant would be expected to be $O((\epsilon_{11}^\nu)^{-2})$; however \eqref{sym:coeff:eps} shows that it actually is only $O((\epsilon_{11}^\nu)^{-1})$. For the sake of clarity, consider the following notation:
\begin{df}
Define for all $\nu\in [-\Lambda_1,\Lambda_1]$
\be\label{eq:cM}
a_\nu=
-\frac{i\theta \epsilon_{12}^\nu}{\epsilon_{11}^\nu}, \quad
b_\nu=
1-\frac{\theta^2}{\epsilon_{11}^\nu}, \quad
c_\nu=
-\frac{d^\nu}{\epsilon_{11}^\nu},
\ee
so that under Assumption \ref{ass:SDT}, the matrix $\bM_\nu$ introduced in Definition \ref{def:M} reads
\be \label{eq:matm}
\bM_\nu (x,\theta)=
\begin{pmatrix}
a_\nu(x,\theta)&b_\nu(x,\theta) \cr
c_\nu (x)&-a_\nu(x,\theta)
\end{pmatrix}.
\ee
These quantities are well-defined
except possibly at the roots of $\epsilon_{11}^\nu$.
\end{df}

\begin{rem} \label{rem:rem1}
In the model case $\epsilon^\nu=\epsilon^0+i\nu\mathbb D$ with a diagonal dissipation tensor $\mathbb D=\mathbb I$,
one can use the notations of \cite{jmpa} so that
$\epsilon_{11}^\nu=\alpha(x)+i\nu$ and $\epsilon_{12}^{\nu}=i\delta(x)$. In this case the coefficients are
\begin{align*}
a_\nu(x,\theta)&=\frac{\theta\delta(x)}{\alpha(x)+i\nu},\\
b_\nu (x,\theta)&= 1-\frac{\theta^2}{\alpha(x)+i\nu},\\
c_\nu (x)&=\frac{\delta^2(x)}{\alpha(x)+i\nu}-(\alpha(x)+i\nu).
\end{align*}
So
$c_\nu$ does not depend on $\theta $ which can be deduced also in the general case
from \eqref{eq:cM}.
\end{rem}

\begin{assumption}\label{ass:eps}
The regularized dielectric tensor $\epsilon^\nu(z)$ is locally analytic
with respect to both variables $z$ and $\nu$.
\end{assumption}

As announced, it is straightforward to verify that these
assumptions are satisfied by the physical tensor introduced in \eqref{eq:hyp1}.

Stemming directly from Definition \ref{def:LR} and Assumption \ref{ass:SDT}, a first property of the $2\times 2$ differential system \eqref{eq:fond1} is the following:
\begin{Lemma}
Assume $x_h$ is a local hybrid resonance.
Then for $\nu=0$ the determinant \eqref{eq:dnunu} of the dielectric tensor satisfies
\begin{equation} \label{eq:ddd}
d^0(x_h)=
-\left| \epsilon^0_{12}(x_h)\right| ^2 <0.
\end{equation}
\end{Lemma}

\begin{proof}
Assumption \ref{ass:SDT} yields that $ \epsilon^0_{12}(x_h)=\overline{ \epsilon^0_{21}(x_h)}$ since the matrix is Hermitian.
Using \eqref{eq:dnu2} and \eqref{sym:coeff:eps}, one gets
$$
d^0(x_h)= \epsilon^0_{11}(x_h)^2 -\left| \epsilon^0_{12}(x_h)\right| ^2 = -\left| \epsilon^0_{12}(x_h)\right| ^2<0,
$$
which is indeed a non-zero number in view of \eqref{eq:dnu2:rab}.
\end{proof}

Note that moreover due to
\eqref{sym:coeff:eps}, one has
$$
\epsilon^0_{12}(x_h)=- \overline{ \epsilon^0_{12}(x_h)} \in i \mathbb R^*.
$$
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 $\nu\neq 0$, then the study of the limit as $\nu$ approaches zero, leading back to the original problem. To this purpose it is crucial to develop uniform tools with respect to the regularization parameter $\nu$, so that their properties still hold as $\nu$ goes to zero.

\subsection{A convenient second order equation}

System \eqref{eq:mm} with matrix \eqref{eq:matm} 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.

\begin{Lemma}\label{ass:ac}
Assume that $x_h$ is a local hybrid resonance. Then,
under Assumption \ref{ass:SDT},
the dielectric tensor $\epsilon^\nu$
\eqref{eq:epsilonmu}
is such that:
\begin{itemize}
\item There exists
$\Lambda_2 $ satisfying $0<\Lambda_2\leq\Lambda_1$ (with $\Lambda_1$ defined
in Assumption \ref{ass:UnA})
and $C>0$ independent of the angle of incidence $\theta$ such that
\be \label{eq:ac}
\Bigl| \frac{a_\nu (x,\theta)}{c_\nu (x)} \Bigr| \leq C |\theta |
\quad\text{for $x\in B(x_h,\Lambda_2)$ and for $|\nu |\leq \Lambda_2$.}
\ee
\item The quantity
$$
\Bigl| \frac{a_\nu (x_h,\theta)}{b_\nu (x_h)} \Bigr| = \Bigl| \frac{\epsilon_{12}^\nu (x_h)}{\theta} \Bigr|
$$
is unbounded for $\theta\to 0$.
\end{itemize}
\end{Lemma}

The upper extra diagonal coefficient $b_\nu$ does not
dominate the diagonal part $a_\nu$ for vanishing $\theta$ results in a
singularity in this regime, and as a result the analysis is trickier for the equation on $E_2^\nu$. Therefore the
formulation considered is the one on~$W^\nu$:
\be \label{eq:w}
\frac d {dx} \Bigl(\frac1 {c_\nu} {\frac d {dx} W^{\nu}}\Bigr)=
\Bigl(\frac{a_\nu^2}{c_\nu}+b_\nu-
\Bigl(\frac{a_\nu}{c_\nu}\Bigr)'\Bigr)W^{\nu}.
\ee

\begin{rem}
At the resonance,
the limit equation \eqref{eq:w} has a ``regular singular point'' \cite{co:ed,LWZ}. Indeed, thanks to Definition \ref{def:LR} and Lemma \ref{ass:ac}, the $\sfrac1 {c_\nu}$ term is bounded while the $\sfrac{a_\nu^2}{c_\nu}$ and $b_\nu$ terms behave as $O((\epsilon_{11}^\nu)^{-1})$ at the resonance.
The term $\left(\sfrac{a_\nu}{c_\nu}\right)'$ is locally bounded.
\end{rem}

A consequence of Assumption \ref{ass:eps} provides a unique complex root of the coefficient~$\epsilon_{11}^\nu$ in a neighborhood of the root of the limit coefficient $\epsilon_{11}^0$, thanks to the open mapping theorem:

\begin{Lemma}\label{lemma:3lam}
Under Assumptions \ref{ass:UnA} and \ref{ass:eps},
there exists
$\Lambda_3$ satisfying $0<\Lambda_3\leq\Lambda_2$ (with $\Lambda_2$ defined in Lemma \ref{ass:ac}) such that
if $|\nu|\leq \Lambda_3$ there exists
a unique point \hbox{$X_\nu\in B(x_h,\Lambda_3)\subset \mathbb C$}
(in a complex neighborhood of $x_h$)
solution of
\be \label{eq:xmu}
\epsilon_{11}^{\nu}(X_\nu)=0.
\ee
A first order expansion is
\be \label{eq:xnuxnu}
X_\nu=x_h - \frac{\partial_\nu \epsilon_{11}^0(x_h) }{\partial_x \epsilon_{11}^0(x_h)}\,\nu +O(\nu^2),
\ee
which shows that, choosing adequately $\Lambda_3$, for $0<|\nu| \leq \Lambda_3$:
$X_\nu \not \in \mathbb R$ and $X_\nu\to x_h$ for $\nu \to 0$.
\end{Lemma}

\begin{proof}
Define, in order to use the analytical implicit function theorem (Theorem 10.32 of \cite{rudin}) the function
$$\epsilon(x,\nu)=\epsilon_{11}^{\nu}(x).$$
There exists $\Lambda_4$ such that it is an analytic function in the ball $B((x_h,0),\Lambda_4)\subset \mathbb C^2$. One has
$\epsilon(x_h,0)=0, \partial_x\epsilon(x_h,0)=\partial_x\epsilon_{11}^0(x_h)\not=0$, hence we can apply the implicit function theorem. There exists $\Lambda_5\leq \Lambda_4$ such that, for $(x,\nu)\in B((x_h,0), \Lambda_5)$, there exists a unique solution of $\epsilon(x,\nu)=0$, which is denoted by $X_\nu$. In addition $$X_\nu=x_h-\frac{\partial_\nu\epsilon(x_h,0)}{\partial_x\epsilon(x_h,0)}\,\nu+O(\nu^2).$$
The assumptions of Definition \ref{def:LR} ensure that $i\sfrac{\partial_\nu\epsilon(x_h,0)}{\partial_x\epsilon(x_h,0)}\in \mathbb R^*$.
\end{proof}

\begin{df}\label{def:Xnu}
The unique point $X_\nu\in\mathbb C$ solution of \eqref{eq:xmu} is called
the translated hybrid resonance.
\end{df}

The case $\epsilon_{11}^{\nu}=\alpha(x)+i\nu$ mentioned
in Remark \ref{rem:rem1} satisfies these hypotheses, with $\alpha(x_h)=0$, $\alpha'(x_h) < 0$. For this case one has, at first order, $X_\nu=-i \sfrac \nu{\alpha'(x_h)}+O(\nu^2)$, which means that the
dominant part of the translated resonance is pure imaginary with non-zero imaginary part
for $\nu\neq 0$.
An adequate scaling in \eqref{eq:w} of the unknown~$W^\nu$ 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_{\nu}$, we must use caution.
We define the square root of $x-X_{\nu}$ as the principal square root on complex numbers (which means that the branch cut is on $x\in \mathbb R_-$). For the reader's convenience, the explicit expressions for the square root are:
\begin{enumeratea}
\item
for $ x-\Re X_{\nu}>0$,
\begin{multline*}
\sqrt{x-X_{\nu}}=\sqrt{
\frac12\bigl(
((x-\Re X_{\nu})^2+(\Im X_{\nu})^2)^{\sfrac12}+x-\Re X_{\nu} \bigr)}\\
- i\frac{\Im X_{\nu}}{ \sqrt{ 2 \bigl(((x-\Re X_{\nu})^2+(\Im X_{\nu})^2)^{\sfrac12}+x-\Re X_{\nu} \bigr) } }\,;\hspace*{1.6cm}\mbox{}
\end{multline*}
\item
for $ x-\Re X_{\nu}\leq 0$,
\begin{multline*}
\sqrt{x-X_{\nu}}=\frac{\vert \Im X_{\nu}\vert}{
\sqrt{ 2 \bigl(((x-\Re X_{\nu})^2+(\Im X_{\nu})^2)^{\sfrac12}+\Re X_{\nu}-x\bigr)}}\\
- i \sign
\left(\Im X_{\nu} \right)
\sqrt{ \frac12 \bigl(((x-\Re X_{\nu})^2+(\Im X_{\nu})^2)^{\sfrac12}+\Re X_{\nu}-x
\bigr)}\,.
\end{multline*}
\end{enumeratea}
Check that $-sc_{\nu}(x)(x-X_{\nu})$ has a strictly positive limit when $\nu\to 0$, hence locally in $\nu$ 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
\be\label{rescaled}
y(x)=\frac{\sqrt{x-X_{\nu}}}{\sqrt{-sc_{\nu}(x)(x-X_{\nu})}}\,W^{\nu}(x)\in \mathbb C.
\ee
From \eqref{eq:w} it satisfies on this neighborhood the equation
\begin{multline}
\label{fundamental-equation0}
\frac{d^2 y(x)}{dx^2}
=\biggl(
a_\nu(x)^2 +b_\nu(x) c_\nu (x)- c_\nu(x)
\Bigl(
\frac{a_\nu(x)} {c_\nu(x)}
\Bigr) '\\
+\frac{ \sqrt{-sc_{\nu}(x)(x-X_{\nu})}}{\sqrt{x-X_{\nu}}}
\Bigl(\frac{\sqrt{x-X_{\nu}}}{ \sqrt{-sc_{\nu}(x)(x-X_{\nu})}} \Bigr)'' \biggr)
y(x).
\end{multline}
The complicated coefficient $\sfrac{\sqrt{x-X_{\nu}}}{\sqrt{-sc_{\nu}(x)(x-X_{\nu})}}$ is the principal root of $(-sc_{\nu})^{-\sfrac12}$
In order to solve Equation \eqref{fundamental-equation0}, 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
\begin{equation*}\frac{ \sqrt{-sc_{\nu}(x)(x-X_{\nu})}}{\sqrt{x-X_{\nu}}}
\biggl(\frac{\sqrt{x-X_{\nu}}}{ \sqrt{-sc_{\nu}(x)(x-X_{\nu})}} \biggr)''.
\end{equation*}
We then consider the following quantity.

\begin{df} \label{df:rmu}
Under Assumptions \ref{ass:UnA}, \ref{ass:SDT} and \ref{ass:eps}, we define
the coefficient-function~$R_\nu$ in $B(x_h,\Lambda_3)$ as:
\begin{multline*}
R_\nu(x)=(x-X_\nu)
\biggl(a_\nu^2(x) + b_\nu(x) c_\nu(x) -
c_\nu (x)
\Bigl(\frac{a_\nu}{ c_\nu} \Bigr)'(x)\\
+\frac{ \sqrt{-sc_{\nu}(x)(x-X_{\nu})}}{\sqrt{x-X_{\nu}}}
\biggl(\frac{\sqrt{x-X_{\nu}}}{ \sqrt{-sc_{\nu}(x)(x-X_{\nu})}} \biggr)'' (x) + \frac1{ 4(x-X_\nu)^2 } \biggr).
\end{multline*}
\end{df}

With this notation, Equation \eqref{fundamental-equation0} reads
\begin{equation}\label{fundamental-equation}
\frac{d^2 y(x)}{dx^2}=
\Bigl(-
\frac{1}{4(x-X_\nu)^2}+\frac{R_\nu(x)}{x-X_\nu}
\Bigr)y(x), \qquad x\in \mathbb R.
\end{equation}
Since
the function $R_\nu$ can also be defined by continuity at $X_\nu$, as stated in the following result,
this form of the equation evidences the singularity of the coefficient.

\begin{rem}
Note that the method described below is valid in a more general case than the hybrid resonance introduced in Definition \ref{def:LR}, namely when $\epsilon_{11}^0(x_h)=\partial_x\epsilon_{11}^0(x_h)=0$ and $\partial^2_{x^2}\epsilon_{11}^0(x_h)\in \mathbb R^*$. We denote by ${\tilde s}$ the sign of $\partial^2_{x^2}\epsilon_{11}^0(x_h)$. In this case one considers the limit $k_\nu$, for $x\to X_\nu$, of
$$(x-X_\nu)^2\biggl(a_\nu^2 + b_\nu c_\nu -
c_\nu
\Bigl(\frac{a_\nu}{ c_\nu} \Bigr)'+
\sqrt{ {\tilde s}c_\nu }
\left(\sfrac1{ \sqrt{ {\tilde s}c_\nu} } \right)'' \biggr), $$
which is finite. The value of $k_\nu$ will generate the family of approximate solutions needed for the subsequent study.
\end{rem}

Returning to the situation described in Definition \ref{def:LR}, one has
\begin{Lemma} \label{lemme:rana}
Under Assumptions \ref{ass:UnA}, \ref{ass:SDT} and \ref{ass:eps},
there exists $\Lambda_6 $ satisfying $0<\Lambda_6\leq \Lambda_3$
(with $\Lambda_3$ defined in Lemma \ref{lemma:3lam})
such that the coefficient function $z\mto R_\nu(z)$ is analytic on a ball $z\in B(x_h,\Lambda_6)$ uniformly for $|\nu | < \Lambda_6$.
Moreover,
$$
R_\nu(X_\nu)= \left[
-\frac{i\theta(\epsilon_{12}^\nu)'+(\epsilon_{12}^{\nu})^2}{(\epsilon_{11}^{\nu})'}-\frac{(\epsilon_{12}^{\nu})'}{\epsilon_{12}^{\nu}}+\frac14 \frac{(\epsilon_{11}^{\nu})''}{(\epsilon_{11}^{\nu})'}
\right](X_\nu),
$$
where $X_\nu$ is the translated resonance. This makes sense by continuity for $\nu $ small
since
$(\epsilon_{11}^{\nu})'(X_\nu)=
\partial_x \epsilon_{11}^{\nu}(X_\nu)\neq 0$ from
\eqref{eq:dnu2}.
\end{Lemma}
The proof of this Lemma is given in Section \ref{annex}.
Since $R_\nu(X_\nu)$ is bounded, as expected, the singularity is explicit in Equation \eqref{fundamental-equation}: there is a $(x-X_\nu)^{-2}$ singularity, and only if $R_\nu(X_\nu)\neq 0$ then there is also a $(x-X_\nu)^{-1}$ singularity.

\begin{figure}[h]
\begin{center}
\begin{tabular}[h]{cccccc}
\scalebox{0.45}{\input{./despres-et-al_figs/compi.pdf_t}} & & &&&
\scalebox{0.45}{\input{./despres-et-al_figs/compi4.pdf_t}}
\end{tabular}
\end{center}
\caption{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.}
\label{figpremier}
\end{figure}

Since $R_\nu$ is bounded and analytic, one can formally shift
the equation and the unknown in the complex
plane. This is illustrated in Figure \ref{figpremier}. So we consider the function
\be \label{be:be}
\widetilde{y}(z)=
\frac{\sqrt z } { \sqrt{ -sc_\nu(z+X_\nu) z } }\,W^{\nu}(z+X_\nu), \quad z\in \mathbb C^*.
\ee
which is a shit of \eqref{rescaled}
and write the equation for this new unknown.
It can be justified that $W^\nu$, solution
of \eqref{eq:w},
admits a convenient
extension in a complex neighborhood of the hybrid resonance $x_h$.
It provides a rewriting of Equation \eqref{fundamental-equation0} with these new variables.

\begin{df}
[Shifted equation]
\label{df:hnu}
The shifted equation writes
\be \label{eq:gg}
\frac{d^2 \widetilde y(z)}{dz^2}=
\Bigl(-
\frac{1}{4 z^2}+\frac{R_\nu(z+X_\nu)}{z}
\Bigr)\widetilde y(z), \qquad z\in \mathbb C^*.
\ee
\end{df}

Our strategy is first to
construct an explicit quasi-solution of the shifted equation
\eqref{df:hnu} which is more general than
the original equation \eqref{fundamental-equation0}--\eqref{fundamental-equation}
because it is a complex equation, so as to provide an explicit quasi-solution
of \eqref{fundamental-equation0}--\eqref{fundamental-equation}.
The modifications to obtained exact solutions of the original
equation \eqref{fundamental-equation} will be performed
in a second step.

\subsection{Freezing and defreezing the coefficient-function} \label{subsec:cbf}

This section focuses on Equation \eqref{eq:gg} for
all $\nu \neq 0$. We require that the function $\mathcal R(\cdot) = R_\nu(x+X_\nu)$ is locally analytic with coefficients uniform with respect
to $\nu$.
We consider the more general from
\be \label{eq:ggcal}
\frac{d^2 \widetilde y(z)}{dz^2}=
\Bigl(-
\frac{1}{4 z^2}+\frac{\mathcal R(z)}{z}
\Bigr)\widetilde y(z), \qquad z\in \mathbb C^*.
\ee

Equation \eqref{eq:ggcal} will be linked to Bessel's equation, replacing the coefficient function~$\mathcal R$ by its value
at the singularity. Indeed, consider the resulting new equation
\begin{equation} \label{eq:froz}
\frac{d^2 Y (z)}{dz^2}=
\Bigl(-
\frac{1}{4z^2}+\frac{\mathcal R (0)}{z}
\Bigr)
Y(z), \qquad z\in \mathbb C^*.
\end{equation}
This equation is presented in the chapter on Bessel functions in the classical textbook~\cite{AS}, 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.

\begin{Lemma} \label{lemme1}
If $\mathcal R(0)\neq 0$,
a pair of independent solutions of \eqref{eq:froz} can be expressed
with Bessel functions as
\[
\begin{cases}
\dpl
z \mto \sqrt{ z }
J_0
\left(\lambda \sqrt{ z } \right), & z\in \mathbb C,
\\
\dpl
z \mto \sqrt{ z }Y_0
\left(\lambda \sqrt{ z }
\right), & z\in \mathbb C,
\end{cases}
\]
where $\lambda=2
\sqrt{ -\mathcal R(0) }$.

If $\mathcal R(0)=0$,
a pair of independent solutions of \eqref{eq:froz} can be expressed with the simpler expression
\[
\begin{cases}
\dpl
z\mto \sqrt{ z}, & z\in \mathbb C,
\\
\dpl
z\mto \sqrt{ z} \log \left(z\right), & z\in \mathbb C.
\end{cases}
\]
\end{Lemma}

\begin{rem} \label{remark:1}
The function $J_0(z)$ admits a infinite converging expansion in powers of~$z^{2 }$
\be \label{eq:j0}
J_0(z)=\sum_{k=0}^\infty \frac{ \left(-\frac14 z^2\right)^k }{ (k!)^2 },
\ee
and the function
\[
Y_0(z)-\frac2\pi \log \left(\sfrac z 2 \right) J_0(z)
\]
admits a infinite converging expansion in powers
of $z^2$
\be \label{eq:y0}
T_0(z):=Y_0(z)-\frac2\pi \log \left(\sfrac z 2\right) J_0(z)=
-\frac2\pi \sum_{k=0}^{\infty }
\frac{ \psi(k+1)} { (k!)^2 } \Bigl(-\frac{ 1}4z^2 \Bigr)^k,
\ee
where $\psi(1) = -\gamma$ and $\psi (k) =-\gamma+ \sum_{\ell=1}^{k-1} \ell^{-1}$ for all $k\in\mathbb N$ such that $k\geq 2$, and $\gamma$ is the Euler constant.

Therefore a more convenient pair of independent solutions of
\eqref{eq:froz}
also valid for $\lambda=0$ reads
\begin{equation}\label{eq:rep0}
\begin{cases}
\dpl
z\mto \sqrt{ z}\,
J_0
\left(\lambda \sqrt{ z } \right), & z\in \mathbb C,
\\[3pt]
\dpl
z\mto \sqrt{ z}\,
\Bigl[Y_0
\left(\lambda \sqrt{ z }
\right)
-\frac2\pi \log \left(\sfrac \lambda 2 \right)J_0
\left(\lambda \sqrt{ z } \right)
\Bigl], & z\in \mathbb C.
\end{cases}
\end{equation}
\end{rem}
Unlike the one proposed in Lemma \ref{lemme1},
the pair
of independent solutions
\eqref{eq:rep0} is uniform with respect to $\lambda$, for any function $\mathcal R$ continuous in the neighborhood of the origin. So when applied to the case
\begin{equation} \label{eq:rnu0}
\mathcal R = R_\nu(\cdot+X_\nu),
\end{equation}
it will be uniform with respect to $\nu$.
Indeed in this case, the parameter $\lambda$ will depend on $\nu$ since then $\lambda=2\sqrt{ - R_\nu(X_\nu) }$. We will therefore denote by $\lambda_\nu$ the value of $\lambda=2\sqrt{ - R_\nu(X_\nu) }$ to emphasize its dependence with respect to $\nu$. So the second function from Lemma \ref{lemme1} has a logarithmic
divergence when $\lambda_\nu \sqrt{ z} \approx 0$, which means that it diverges if either
$\lambda_\nu\approx 0$ or $\sqrt{ z}\approx 0$. On the other hand the logarithmic
divergence of the second function of \eqref{eq:rep0} is by construction only for $\sqrt{z}\approx 0$. As a result Functions \eqref{eq:rep0} are well-defined for $z\neq 0$ as $\nu$ goes to zero, as well as for $\nu=0$.

The Bessel function \eqref{eq:y0} recasts as $Y_0(z)=T_0(z)+\frac2\pi \log \left(\sfrac{z}2\right) J_0(z)$
and the second function in \eqref{eq:rep0}
then reads
\begin{equation} \label{eq:rep0_bis}
z\mto \sqrt{ z}
\Bigl[T_0
\left(\lambda \sqrt{ z }
\right)+
\frac1\pi \log \left(z \right)J_0
\left(\lambda \sqrt{ z } \right)
\Bigl].
\end{equation}
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
$$
\mbox{Arg}(z)\in (-\pi, \pi)\Longleftrightarrow
z \in \mathbb C - (-\infty, 0].
$$

The Bessel based functions, solutions of the equation with frozen
coefficient \eqref{eq:froz},
provide the formal
dominant behavior of the solutions
of Equation \eqref{eq:ggcal}. Nevertheless
the rigorous justification of this statement is not evident, mainly
because the frozen coefficient $\mathcal R(0)$ is only a first order local approximation
of the coefficient function~$\mathcal R$. Indeed the singularity of the $\sfrac{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\left(z^{-1} \right)$.

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 \eqref{eq:ggcal} and the frozen coefficient equation \eqref{eq:froz}. 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 \eqref{eq:ggcal}.

\subsubsection{Eikonal equation and stretching function}

Consider first another intermediate equation by reintroducing a varying coefficient in Equation \eqref{eq:froz}:

\begin{df}[Stretching procedure]
Let $\rho $ be an arbitrary smooth function defined in a complex neighborhood of the origin, such that $\rho(0)=0$ and $\rho'(0)\neq 0$.
We will refer to it as the stretching function. This name comes, in particular, from McKelvey~\cite{McK}.

Let $u$ be a given function.
The stretching of $u$ is the function $\tilde u$ defined by
$$
\widetilde u (z)= (\rho'(z))^{-\sfrac12} u (\rho(z)).
$$
\end{df}
No branch cut of the square root is needed for small $|z|$ since we impose $\rho'(0)\neq 0$.
The equation for the stretched function is the following.

\begin{Lemma}
Assume $\rho$ is a stretching function and $u$ solves the frozen coefficient Equation \eqref{eq:froz}.
Set $s=\sfrac{
\left((\rho')^{-\sfrac12}
\right)''
}
{
(\rho')^{-\sfrac12}
}$.
Then the stretched unknown satisfies
\be \label{eq:g}
\frac{ d^2 \widetilde u(z)}{dz^2}=\biggl(
-\frac{ \rho'(z) ^2}{ 4 \rho(z)^2 }
+ \frac{ \mathcal R(0) \rho'(z)^2}{ \rho(z) }
+s(z)
\biggr) \widetilde u(z), \qquad z\in \mathbb C^*.
\ee
\end{Lemma}

\begin{proof}
The proof relies on the scaling of the stretched unknown: the coefficient of the first order derivative term in the computation
of $\widetilde u ''$ is $((\rho')^{-1/2}) ' \rho' + ((\rho')^{1/2})'$, which happens to be zero.
\end{proof}
Note that the function $s$ is a priori bounded around $0$ since $\rho'(0)\neq 0$, so that, like Equation \eqref{eq:ggcal}, Equation \eqref{eq:g} evidences the singularities of its coefficient.

At this point the stretching function is a tool that will be designed for the approximation of Equation \eqref{eq:ggcal}.
Indeed, although Equation \eqref{eq:g} is different from the initial equation, namely \eqref{fundamental-equation}, 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
\be\label{eikonal}
\rho' (z)^2
\biggl(-
\frac{1}{4\rho (z)^2}+\frac{\mathcal R(0)}{\rho(z)}
\biggr)=
-
\frac{1}{4z^2}+\frac{\mathcal R(z)}z,
\qquad
z\mbox{ in a neighborhood of }0.
\ee
We call this equation the ``eikonal'' equation for the similar equation for $\phi$ if one seeks a solution of the wave equation $(\Delta -c^{-2}(x)\partial^2_{t^2})u=0$ under the form
\[
u(x,t)=a(x,k)e^{ik(\phi(x)-t)}.
\]
Note that \eqref{eikonal} implies
\[
\frac{\rho'(z)}{\rho(z)}=\pm\frac{1}{z}\sqrt{\frac{1-4zR(z)}{1-4\rho(z)R(0)}}
\]
for $z, \rho(z)$ non-zero. As we choose $\rho$ bounded in the neighborhood of 0, only the $+$ sign is relevant.
So one can write
\begin{align*}
\Bigl(\frac{\rho(z)}{z} \Bigr)'(z)= \frac{\rho'(z)}{z}- \frac{\rho(z)}{z^2}&=
\frac{\rho(z)}{z^2} \biggl(\sqrt{\frac{1-4zR(z)}{1-4\rho(z)R(0)}} -1 \biggr)\\[3pt]
&=
\frac{\rho(z)}{z^2} \frac{\sqrt{1-4zR(z)} - \sqrt{1-4\rho(z)R(0)}}{\sqrt{1-4\rho(z)R(0)}}.
\end{align*}
Use of the conjugate quantity of the numerator yields
the equation
\begin{equation}\label{eikonal2}
\Bigl(\frac{\rho(z)}{z} \Bigr)'(z)=
\frac{\rho(z)}{z} F \Bigl(z,\frac{\rho(z)}{z} \Bigr),
\quad z\mbox{ in a neighborhood of }0,
\end{equation}
where $F$ is defined by
$$
F(z,a)
=\frac{
4\left[ \mathcal R(0)a- \mathcal R(z)\right]
}{
\sqrt{1-4z \mathcal R(z) }
\sqrt{ 1-4z \mathcal R(0) a }
+
\left(1-4z \mathcal R(0) a \right)
}.
$$
Equation \eqref{eikonal2} is defined
also at $z=0$. Since $F(0,a)=2\left[ \mathcal R(0)a- \mathcal R(0)\right]$ is well-defined, then $F(z,a)$ makes sense at least for small $z$ and
for bounded $\mathcal R$.

We now reintroduce $\mathcal R= R_\nu(\cdot+X_\nu)$ and consider
the more general equation in the complex plane
\begin{equation}\label{eikonal2:2}
\sigma_\nu'(z) =
\sigma_\nu(z) F_\nu\left(z,\sigma_\nu(z) \right),
\quad \rho_\nu(z)=z \sigma_\nu(z),
\qquad z\mbox{ in a neighborhood of }0,
\end{equation}
where $F_\nu$ is defined by
$$
F_\nu(z,a)
=\frac{
4\left[ R_\nu(X_\nu)a- R_\nu(z+X_\nu)\right]
}{\sqrt{
1-4 R_\nu(z+X_\nu) }\sqrt{
1-4 R_\nu(X_\nu) a }
+\left(1-4 R_\nu(X_\nu) a \right)
}.
$$
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 $\nu$.

\begin{Lemma}
[Solution of the eikonal equation]
\label{lm:ronu}
Under Assumptions \ref{ass:UnA}, \ref{ass:SDT} and \ref{ass:eps}, consider Equation \eqref{eikonal2:2}.

Then there exists a constant $\Lambda_7\in (0, \Lambda_6]$
(with $\Lambda_6$ defined in Lemma \ref{lemme:rana})
such that:
\begin{itemize}
\item For all $\vert\nu\vert\leq \Lambda_7$ there exists a solution to Equation \eqref{eikonal2:2}
on a maximal interval $I_\nu=]a_\nu, b_\nu[\subset ]-\Lambda_3, \Lambda_3[$
(with $\Lambda_3$ defined in Lemma \ref{lemma:3lam}). 

\item A solution of the eikonal equation \eqref{eikonal2:2} is
$\rho_\nu(z)/z$ for $z\in B(0, \Lambda_7)$ with
\begin{equation} \label{eq:stc2}
\rho_\nu(z)=z \exp \biggl(
z \int_0^1 F_\nu(zt,\sigma_\nu(zt))dt
\biggr).
\end{equation}
It is the unique solution such that $\rho_\nu'(0)=1$.
Note that $F$ is evaluated on the straight line from $0$ to $z$.

\item This function is analytic in the ball $z\in B(0, \Lambda_7)$.
\end{itemize}
\end{Lemma}

\begin{proof}
Note that $]-\Lambda_3, \Lambda_3[$ is not necessarily the maximal existence interval of the ordinary differential equation \eqref{eikonal}.

The function $\sigma_\nu= \sfrac{\rho_\nu}{z} $ is solution to the equation
\begin{equation} \label{eq:sigma}
\sigma_\nu'(z)=
\sigma_\nu(z) F_\nu(z,\sigma_\nu(z)), \qquad z\mbox{ in a neighborhood of }0.
\end{equation}
The value of $F_\nu(0,a)= 4 R_\nu(X_\nu)
\left(a -1\right)
$
is
a well-defined (\ie bounded, or non-singular) real number.
Equation
\eqref{eq:sigma} can be solved
in the complex plane using
a theorem from Coddington-Levinson
\cite[Th.\,8.1, p.\,34]{co:ed}.
For simplicity consider the initial data $\sigma_\nu(0)=1$ which yields
$\sigma_\nu'(0)=F_\nu(0,1)=0$. This choice is arbitrary.

Since $\sigma_\nu'(0)=0$,
one has the Taylor expansion $\rho_\nu(z)=z+O_\nu(z^3)$,
where the constant in $\left| O_\nu(z^3)\right| \leq C |z|^3$
is uniform with respect to $\nu$.
All other coefficients of the Taylor expansion of the stretching function $\rho_\nu$
can be easily computed from the \emph{integral representation in the complex plane}
\eqref{eq:stc2},
where the path of integration is the straight line.
The analyticity is provided by the Coddington-Levinson theorem.
\end{proof}

Thanks to the stretching function described in Lemma \ref{lm:ronu}, we can now express
explicitly the solutions of Equation \eqref{eq:g}, where, in view of Lemma \ref{lemme1} and Equation \eqref{eq:rnu0}, one needs
to set
$$
\lambda_\nu=
2\sqrt{- R_\nu(X_\nu)}
=2\sqrt{-\mathcal R(0)}.
$$

\begin{df}\label{df:UnuVnu}
Under Assumptions \ref{ass:UnA}, \ref{ass:SDT} and \ref{ass:eps}, consider $\Lambda_7$ defined in Lemma \ref{lm:ronu}, and for all $\vert\nu\vert\leq \Lambda_7$
the stretching function $\rho_\nu$ described in Lemma \ref{lm:ronu}. A pair $(U_\nu,V_\nu)$ of independent solutions of
Equation \eqref{eq:g} is defined for all $(\nu,z)\in (-\Lambda_7,\Lambda_7)\times B(0,\Lambda_7)$ by
\begin{equation*}
\left\{\begin{array}{rl}
U_\nu(z)&\displaystyle = \sqrt{\frac{\rho_\nu(z) }{\rho_\nu'(z)}}\,
J_0
\left(\lambda_\nu \sqrt{ \rho_\nu(z)} \right),\\[10pt]
V_\nu(z)&\displaystyle = \sqrt{\frac{\rho_\nu(z) }{\rho_\nu'(z) }}
\left[Y_0
\left(\lambda_\nu \sqrt{ \rho_\nu(z)} \right) -\dfrac2\pi \log\left(\sfrac{\lambda_\nu}2 \right)
J_0
\left(\lambda_\nu \sqrt{ \rho_\nu(z)} \right) \right], \\[10pt]
& \displaystyle = \sqrt{\frac{\rho_\nu(z) }{\rho_\nu'(z) }}
\left[T_0
\left(\lambda_\nu \sqrt{ \rho_\nu(z)} \right)
+\dfrac1\pi \log\left(\rho_\nu(z) \right)
J_0
\left(\lambda_\nu \sqrt{ \rho_\nu(z)} \right) \right].
\end{array}\right.\end{equation*}
We adopt the convention (as one did before for the root of $-sc_{\nu}$) that
\[
\sqrt{\frac{\rho_\nu(z) }{\rho_\nu'(z) }} =
\biggl({\frac{ \int_0^1 \rho_\nu'(tz)dt }{\rho_\nu'(z) }} \biggr)^{\sfrac12}\sqrt z
\]
is defined in $ \mathbb C- (-\infty,0]$, that is apart from the branch cut $(-\infty,0]$.
\end{df}

These functions are uniformly bounded with respect to $\nu$.

\begin{Lemma} \label{lemma:bound}
Under Assumptions \ref{ass:UnA}, \ref{ass:SDT} and \ref{ass:eps}, consider
the stretching function $\rho_\nu$ \eqref{eq:stc2}
for $\vert\nu\vert\leq \Lambda_7$ defined in Lemma \ref{lm:ronu}.

Then the functions $U_\nu$ and $V_\nu$ are bounded in a complex
neighborhood
of the resonance $x_h\in \mathbb C$, uniformly for small $\nu$:
there exists a constant $\Lambda_8\in (0,\Lambda_7)$
(with $\Lambda_7$ defined in Lemma \ref{lm:ronu})
such that for all $(\nu,z)\in [-\Lambda_8,\Lambda_8]\times \overline{ B(0,\Lambda_8) }$
\be \label{eq:bound_uv}
\left|U_\nu(z) \right|+ \left|V_\nu(z) \right| \leq C.
\ee
Moreover, $\Lambda_8$ can be determined such that there exists $\delta>0$ with
$$
\left| \rho_\nu'(z) \right|\geq \delta>0, \quad \forall \nu\in \overline{B(0, \Lambda_8)}
$$
in the same neighborhood.
\end{Lemma}

\begin{proof}
The property on $\rho_\nu$ is immediate since $\sigma_\nu(0)=1$.
The uniform boundedness of $U_\nu$ stems from the boundedness of $J_0$ \eqref{eq:j0}.
Concerning $V_\nu$, the $\log \lambda_\nu$
term has been carefully removed, see \eqref{eq:y0}, so that the boundedness is achieved
even for vanishing $\lambda_\nu$: the remaining $\log \rho_\nu(z)$ is controlled by the $\sqrt z$
contribution that comes from $\sqrt { \rho_\nu(z) }$.
\end{proof}

It is then possible to apply the reverse shift.

\begin{Lemma}
The functions $g(x)= U_\nu(x-X_\nu)$ or
$g(x)= V_\nu(x-X_\nu)$ are solutions of the equation
on the real line
\begin{equation} \label{eq:realshift}
\frac{d^2}{dx^2}g (x) =\Bigl[
-
\frac{1}{4(x-X_\nu)^2}+\frac{R_\nu(x)}{x-X_\nu}+
s_\nu(x-X_\nu)
\Bigr]
g(x), \ \ \forall x\in B(x_h,\Lambda_8)\cap \mathbb R,
\end{equation}
where
$s_\nu=\sfrac{
\left((\rho_\nu')^{-\sfrac12}
\right)''
}
{
(\rho_\nu')^{-\sfrac12}
}$.
\end{Lemma}

\begin{proof}
By definition, all these functions are analytic in the ball $B(x_h,\Lambda_8)-
\left\{X_\nu \right\}$ (that is except at $X_\nu$) which is uniform
with respect to $\nu$.
But, when $\rho_\nu $ satisfies the eikonal relation \eqref{eikonal} with $\mathcal R(z)=R_\nu(z+X_\nu)$, the functions without the shift $z\mto U_\nu(z)$ and $z\mto V_\nu(z)$
are solutions of \eqref{eq:g}.
Therefore the result holds for analytic functions
on $B(x_h,\Lambda_8)\cap \mathbb R$ for $0<| \nu |\leq \Lambda_3$. So this relation is true everywhere by analytic continuation.
The proof is ended.
\end{proof}

\subsubsection{General solution with a Duhamel's principle and a limit process}\label{ssec:DP}
In order to go back to the equation with a variable coefficient, let us now consider the original problem
\eqref{fundamental-equation}.
One can rewrite this equation as
\begin{multline}
\label{eq:encoreg}
\frac{d^2}{dx^2}y (x) =\Bigl[
-
\frac{1}{4(x-X_\nu)^2}+\frac{R_\nu(x)}{x-X_\nu}+
s_\nu(x-X_\nu)
\Bigr]
y(x)\\
- s_\nu(x-X_\nu) y(x), \quad x\in \mathbb R.
\end{multline}
The last term is a perturbation with respect to Equation \eqref{eq:realshift}. 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
\eqref{eq:realshift} are bounded.

\begin{Lemma}\label{lem:Cnu}
Suppose Assumptions \ref{ass:UnA}, \ref{ass:SDT} and \ref{ass:eps} hold.
Consider for all $\nu$ a given reference point $x^* \in \mathbb R$ close to
$x_h$ in the ball of analyticity $B(x_h,\Lambda_8)$, where $\Lambda_8$ is defined in Lemma \ref{lemma:bound}.
Then the function
\begin{equation} \label{eq:yw1}
y(x)= A_\nu(x) U_\nu(x-X_\nu)
+
B_\nu(x) V_\nu(x-X_\nu)
\end{equation}
is a solution of \eqref{eq:encoreg} for all $(\nu,x)\in (-\Lambda_8,\Lambda_8)\times B(x_h,\Lambda_8)\cap\mathbb R$ if and only if
$\bC_\nu$, defined by
\begin{equation} \label{eq:c40}
\bC_\nu(x)=
\begin{pmatrix}
A_\nu(x) \\
B_\nu(x)
\end{pmatrix}, \quad x\in \mathbb B(x_h,\Lambda_8)\cap\mathbb R,
\end{equation}
is a solution of the integral equation
\be \label{eq:c}
\bC_\nu(x)=
\bC_\nu(x^*)
+
\pi \int_{x^*}^x
\mathcal M_\nu(t)
\bC_\nu(t)dt,
\ee
where the matrix
\be\label{eq:defM}
\mathcal M_\nu(t)
=s_\nu (t -X_\nu)
\begin{pmatrix}
-U_\nu(t -X_\nu)V_\nu(t -X_\nu) & -V_\nu(t -X_\nu)^2 \\[3pt]
U_\nu(t -X_\nu)^2 & U_\nu(t -X_\nu)V_\nu(t -X_\nu)
\end{pmatrix}
\ee
is uniformly bounded in $L^\infty$ for $(\nu,t)\in (-\Lambda_8,\Lambda_8)\times B(x_h,\Lambda_8)\cap\mathbb R$.
\end{Lemma}

As usual for such Volterra integral equations, the initial point $x^*
\in \mathbb R$ can be chosen arbitrarily.

\begin{proof}
This is the standard procedure of variation of parameters. Consider 
\begin{equation}\label{eq:uv}
y'' - \Bigl[
-
\frac{1}{4(x-X_\nu)^2}+\frac{R_\nu(x)}{x-X_\nu}+
s_\nu(x-X_\nu)
\Bigl]
y = f,
\end{equation}
with
$f(x) = -s_\nu(x-X_\nu) y(x)$.
The solution $y$ can be expressed as a combination of $U_\nu
(x-X_\nu)$ and $V_\nu(x-X_\nu)$, with appropriate coefficients $A_\nu(x)$ and $B_\nu(x)$
$$
y(x)=A_\nu(x) U_\nu(x-X_\nu)
+
B_\nu(x) V_\nu(x-X_\nu).
$$
To construct the coefficients $A_\nu(x)$ and $B_\nu(x)$,
we first assume
that
\be \label{eq:ab1}
A_\nu'(x) U_\nu(x-X_\nu)
+
B_\nu'(x) V_\nu(x-X_\nu)=0.
\ee
So the first derivative of $y$ is
\begin{equation} \label{eq:ywprime}
y'(x)= A_\nu(x) U_\nu'(x-X_\nu)
+
B_\nu(x) V_\nu'(x-X_\nu).
\end{equation}

The second derivative reads
$$
y''(x)=
A_\nu(x) U_\nu''(x-X_\nu)
+
B_\nu(x) V_\nu''(x-X_\nu)
+
A_\nu'(x) U_\nu'(x-X_\nu)
+
B_\nu'(x) V_\nu'(x-X_\nu).
$$
Using \eqref{eq:realshift}, one gets
\begin{multline*}
y''(x)= \Bigl[
-
\frac{1}{4(x-X_\nu)^2}+\frac{R_\nu(x)}{x-X_\nu}+
s_\nu(x-X_\nu)
\Bigr]
y(x)\\
+
A_\nu'(x) U_\nu'(x-X_\nu)
+
B_\nu'(x) V_\nu'(x-X_\nu).
\end{multline*}
Since $y$ is solution of \eqref{eq:uv},
one gets
\begin{equation} \label{eq:ab2}
A_\nu'(x) U_\nu'(x-X_\nu)
+
B_\nu'(x) V_\nu'(x-X_\nu)
=f,
\end{equation}
with
\[
f(x)=
- s_\nu(x-X_\nu) \left(A_\nu(x) U_\nu(x-X_\nu)
+
B_\nu(x) V_\nu(x-X_\nu) \right).
\]
One can solve now the linear system made of
\eqref{eq:ab1} and \eqref{eq:ab2}.
To compute the determinant of this linear system
consider the Wronskian of $U_\nu$ and $V_\nu$.
Elementary calculations and the identity
$\mathcal W_{\left(J_0,Y_0 \right)}(x)=\sfrac 2{\pi x}$
yield
\begin{align*}
\mathcal W_{\left(U_\nu, V_\nu \right)}(x)&=
\frac{\rho_\nu(x-X_\nu) }{\rho_\nu'(x-X_\nu) }\;
\lambda_\nu\; \frac{\rho_\nu'(x-X_\nu)}{2\sqrt{\rho_\nu(x-X_\nu)}}\;
\mathcal W_{\left(J_0,Y_0 \right)}\left(\lambda_\nu\sqrt{\rho_\nu(x-X_\nu)}\right)\\[5pt]
&=
\frac{ \lambda_\nu \sqrt{\rho_\nu(x-X_\nu)}}2\; \mathcal W_{\left(J_0,Y_0 \right)}\left(\lambda_\nu\sqrt{\rho_\nu(x-X_\nu)}\right)=\frac{1} \pi
\end{align*}
where the convention is that
\[
\sqrt{ \rho_\nu(x-X_\nu) }=
\sqrt{ \int_0^1 \rho'_\nu(t(x-X_\nu)) dt }\;
\sqrt{x-X_\nu }
\]
is defined apart from the branch cut $x-X_\nu\in (-\infty,0] \subset \mathbb C$ thanks to $\rho'_{\nu}(0)=1$.
By analytic continuation the determinant
of the linear system \eqref{eq:ab1}--\eqref{eq:ab2} is also equal to the same value.

Therefore the solution reads
\[
\begin{cases}
A_\nu'(x)=- \pi s_\nu(x-X_\nu) \left[ A_\nu(x) U_\nu(x-X_\nu)
+
B_\nu(x) V_\nu(x-X_\nu) \right] V_\nu(x-X_\nu),
\\[3pt]
B_\nu'(x)= \phantom{-} \pi s_\nu(x-X_\nu) \left[ A_\nu(x) U_\nu(x-X_\nu)
+
B_\nu(x) V_\nu(x-X_\nu) \right] U_\nu(x-X_\nu).
\end{cases}
\]
After integration, it yields the representation formula \eqref{eq:c}.
The boundedness of the kernel comes from the properties of
$s_\nu$, $U_\nu$ and $V_\nu$ in Lemma \ref{lemma:bound}.
It completes the proof.
\end{proof}

Under the same conditions, one can define the integral operator ${\mathcal K}_\nu$
\be \label{eq:knuknu}
\mathcal K_\nu(f)(x)= \pi\int_{x^*}^x \mathcal M_\nu (t) f(t) dt,
\ee
so that
the integral equation \eqref{eq:c} reads
$
\bC_\nu-
\mathcal K_\nu (\bC_\nu)
=\bC_\nu(x^*)$, with a constant vector right hand side.
Classically for this Volterra second type integral equation (see~\cite{Trico}), the solution is expressed with
the resolvent kernel
\be \label{eq:qnuqnu}
{\mathfrak Q}_\nu [f]=\sum_{n\geq 1}{\mathcal K}_\nu^{(n)}(f),
\ee
under the form
\begin{multline} \label{eq:duha1}
\bC_\nu(x)= \bC_\nu(x^*)+{\mathfrak Q}_\nu [\bC_\nu(x^*)] (x) =
\Big(\mathrm{Id}+{\mathfrak Q}_\nu \Big) [\bC_\nu](x^*),
\\
\forall x\in B(x_h,\Lambda_8)\cap \mathbb R,
\end{multline}
where $\bC_\nu(x^*)$ is the initial condition.

Foreseeing the limit process, the following results state the limit properties of $\rho_\nu$
and $\mathcal M_\nu$, and of other useful quantities, as $\nu$ approaches zero.

\begin{Lemma} \label{lm:lim1}
Under Assumptions \ref{ass:UnA}, \ref{ass:SDT} and \ref{ass:eps}, consider for all $\vert\nu\vert\leq \Lambda_7$
the stretching function $\rho_\nu$ described in Lemma \ref{lm:ronu}.
One has that $\lim_{\nu\to 0^\pm}\rho_\nu(z)=
z \sigma_0(z)$ on $B(x_h, \Lambda_7)$, where $\sigma_0$ is the solution of Equation \eqref{eq:sigma} such that
$\sigma_0(0)=1$. This limit function will naturally be called $\rho_0$.

Then the
$\log
$ and square root terms have limits which depend on the sign of $\nu$: we~set
\begin{equation}
\label{eq:rhomupm}
\begin{split}
\log {\rho_0(x-x_h)}^\pm &= \lim_{\nu \to 0^\pm}\log {\rho_\nu(x-X_\nu)}\\
&=
\begin{cases}
\log\left({{ \rho_0(x-x_h)}}\right), & \mbox{for }x_h< x, \\
\log\left({{- \rho_0(x-x_h)}}{}\right) \mp i \pi \mathfrak s, & \mbox{for } x<x_h,
\end{cases}
\end{split}
\end{equation}
where
$
\mathfrak s= \sign\left(\sfrac{i \partial_\nu \epsilon_{11}^\nu(x)}{ \partial_x \epsilon_{11}^\nu(x) } \right)_{\left| x=x_h,\nu=0 \right.}$, and
\be \label{eq:sqrhomupm}
\sqrt{ \rho_0(x-x_h)}^\pm=
\lim_{\nu \to 0^\pm}\sqrt{ \rho_\nu(x-X_\nu)}=
\begin{cases}
\sqrt{ \rho_0(x-x_h)}, & \mbox{for }x_h< x, \\[2pt]
\mp i \mathfrak s \sqrt{- \rho_0(x-x_h)}, & \mbox{for } x<x_h.
\end{cases}
\ee
\end{Lemma}
\begin{proof}
The claim on $\rho_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)\in (-\pi, \pi]$. And the limit of the square root term comes from the principal value of the complex square root: $\sqrt z =\sqrt{ |z|}\,e^{i\sfrac{\theta(z)}{2}}$, where $\theta(z)\in (-\pi, \pi]$ is the argument of~$z$.
\end{proof}

\begin{Lemma} \label{lm:lim11}
Under Assumptions \ref{ass:UnA}, \ref{ass:SDT} and \ref{ass:eps},
the matrix $\mathcal M_\nu$ defined in \eqref{eq:defM} has a limit in the following sense.
For all $\Lambda_9\in(0, \Lambda_8)$ (where $\Lambda_8$ is defined in Lemma \ref{lemma:bound}), one has that
\begin{itemize}
\item $\dpl\mathcal M_\nu(x)\underset{\nu \to 0^+}{\To{C^0[x_h-\Lambda_9,x_h+ \Lambda_9]}}\mathcal M^+(x) $
\item $\dpl\mathcal M_\nu(x)\underset{\nu \to 0^-}{\To{C^0[x_h-\Lambda_9,x_h+ \Lambda_9]}} \mathcal M^-(x) $

\smallskip
\item The continuous kernels $\mathcal M^+$ and $\mathcal M^-$ are such that:
$\mathcal M^+-\mathcal M^-=0$ vanishes on $[x_h,x_h+\Lambda_9]$,
while $\mathcal M^+-\mathcal M^-\neq 0$ is non-identically zero on $[x_h-\Lambda_9,x_h]$.
\end{itemize}
Consider (for simplicity) a Cauchy data $C_\nu(x^*)$ which admits a finite (given) limit in~$\nu$ in the
representation \eqref{eq:duha1} of the solution of integral equation.
Then the functions~$A_\nu$ and $B_\nu$ also converge in $C^0[x_h-\Lambda_9, x_h+\Lambda_9] $ to continuous limit functions $A_0^\pm$ and~$B_0^\pm$.
\end{Lemma}

\begin{proof}
The continuous limits of the kernel $\mathcal M_\nu$ are given by application
of Lemma~\ref{lm:lim1} since the kernel entries only depend on the functions $U_\nu$ and $V_\nu$ which are locally bounded, see \eqref{eq:bound_uv}, and continuous.
The limit of functions $A_\nu$ and $B_\nu$ then directly stem from the limit of
\eqref{eq:duha1}, where the kernel defined by \eqref{eq:defM}--\eqref{eq:knuknu}--\eqref{eq:qnuqnu}
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^\pm$ and $B_0^\pm$ are continuous.
It ends the proof.
\end{proof}

\pagebreak[2]
\begin{df} \label{df:W13}
Recall that $s=-\sign\left((\epsilon_{11}^0)'(x_h) \right)$. We define three additional functions as follows:
\[
\begin{cases}
\dpl
W_1^\nu(x)= \sqrt{ -s c_\nu(x) \frac{\rho_\nu(x-X_\nu)}{\rho_\nu'(x-X_\nu)} }\;J_0
\bigl(\lambda_\nu \sqrt{\rho_\nu (x-X_\nu) } \bigr),
\\[12pt]
\dpl
W_2^\nu(x) = \sqrt{-s c_\nu(x) \frac{\rho_\nu(x-X_\nu)}{\rho_\nu'(x-X_\nu)} }\\
\dpl
\hspace*{2cm}
\Bigl(T_0
\bigl(\lambda_\nu \sqrt{\rho_\nu (x-X_\nu) }
\bigr)+\frac1\pi \log(\rho_\nu (x-X_\nu))J_0
\bigl(\lambda_\nu \sqrt{\rho_\nu (x-X_\nu) }
\bigr)
\Bigr), \\[7pt]
\dpl
W_3^\nu(x) = \sqrt{- s c_\nu(x) \frac{\rho_\nu(x-X_\nu)}{\rho_\nu'(x-X_\nu)} }\;T_0
\bigl(\lambda_\nu \sqrt{\rho_\nu (x-X_\nu) } \bigr).
\end{cases}
\]
Since $c_\nu(x)\rho_\nu(x-X_\nu)$ is bounded, the functions $W_1^\nu$ and $W_3^\nu$ are uniformly bounded.

On the other hand, the function
\[
W_2^\nu(x)= W_3^\nu(x) + \frac1\pi \log(\rho_\nu (x-X_\nu)) W_1^\nu(x)
\]
is not bounded, due the $\log$ term.
\end{df}

The analyticity of these functions is described below.

\begin{Lemma}\label{lm:lim2}
Under Assumptions \ref{ass:UnA}, \ref{ass:SDT} and \ref{ass:eps}, consider that $0<\vert\nu\vert < \Lambda_9$ so that
the stretching function $\rho_\nu$ inherits of the properties described in Lemma \ref{lm:ronu}.
One has
\begin{multline}
\label{eq:byby2}
\lim_{\nu \to 0^\pm}
W_1^\nu(x)
=
W_1(x)
:=\sqrt{ \frac{d^0(x) \rho_0(x-x_h)}{s \epsilon_{11}^0(x) \rho_0'(x-x_h)} }\;J_0(\lambda_0 \sqrt{\rho_0(x-x_h)}), \\ \forall x\in B(x_h,\Lambda_9)\cap\mathbb R,
\end{multline}
where the limit function $W_1$
is analytic on $B(x_h,\Lambda_9)\cap\mathbb R$. Similarly
\begin{multline*}
\lim_{\nu \to 0^\pm}
W_3^\nu(x)
=
W_3(x)
:=\sqrt{ \frac{d^0(x) \rho_0(x-x_h)}{s \epsilon_{11}^0(x) \rho_0'(x-x_h)} }\;T_0(\lambda_0 \sqrt{\rho_0(x-x_h)}), \\ \forall x\in B(x_h,\Lambda_9)\cap\mathbb R,
\end{multline*}
and $W_3$ is analytic on $B(x_h,\Lambda_9)\cap\mathbb R$.
\end{Lemma}

\begin{proof}
Take the limit in the series \eqref{eq:j0} that defines $J_0$ is straightforward.
Only even terms
appear in this series, so that
$J_0(\lambda_0 \sqrt{z\sigma_0(z)})$ is analytic.
For small $x-x_h$ one uses \eqref{eq:cM} to show that the weight
\[
\sqrt{-s c_\nu(x)\, \frac{\rho_\nu(x-X_\nu)}{\rho_\nu'(x-X_\nu)} }
\]
is close to
\[
\sqrt{ \frac{d^0(x) \rho_0(x-x_h)}{s\epsilon_{11}^0(x) \rho_0'(x-x_h)} }.
\]
It is away from the branch cut of the square root since
\be \label{eq:byby}
\frac{d^0(x) \rho_0(x-x_h)}{s\epsilon_{11}^0(x) \rho_0'(x-x_h)}\approx \frac{ d^0(x_h)}{ s (\epsilon_{11}^0)'(x_h) }
= \frac{ - |\epsilon_{12}(x_h)|^2 }{ s (\epsilon_{11}^0)'(x_h) }
= \frac{ |\epsilon_{12}(x_h)|^2 }{ \left| (\epsilon_{11}^0)'(x_h) \right| } >0
\ee
by \eqref{eq:ddd} and the value of $s$ recalled in Definition \ref{df:W13}.
So $W_1$ is analytic in the neighborhood of $x_h$.
Mutatis mutandis the proof is the same for the function $W_3^\nu$.
\end{proof}

Considering the logarithmic part in Definition
\eqref{eq:y0} of $Y_0$, the
situation is a little more involved concerning $W_2^\nu$ since
the limits depend on the sign of $\nu$.

\begin{Lemma} \label{lm:lim3}
Under Assumptions \ref{ass:UnA}, \ref{ass:SDT} and \ref{ass:eps}, consider for all $0< \vert\nu\vert < \Lambda_9$
the stretching function $\rho_\nu$ described in Lemma \ref{lm:ronu}.
One has
$$
\lim_{\nu \to 0^\pm}
W_2^\nu(x)
=
W_2^\pm(x), \quad \forall x\in B(x_h,\Lambda_9)\cap\mathbb R,
$$
where the limit functions
are
\begin{multline}
\label{eq:w_2pm}
W_{2}^\pm=
\sqrt{ \frac{d^0(x) \rho_0(x-x_h)}{s\epsilon_{11}^0(x) \rho_0'(x-x_h)} }\\
{}\cdot\Bigl(T_0\bigl(\lambda_0 \sqrt{\rho_0(x-x_h)}\bigr)
+ \frac1\pi
\log\left({ \rho_0}(x-x_h)
\right)^\pm
J_0
\bigl(\lambda_0 \sqrt{ \rho_0(x-x_h)} \bigr)
\Bigl),
\end{multline}
and $s$ defined in Definition \ref{df:W13}.
The limit functions $W_2^+$ and $W_2^-$ are equal on the interval
$(x_h,x_h+\Lambda_9)$, but are different on the interval $(x_h-\Lambda_9,x_h)$, where they satisfy the jump relation
\be \label{eq:jump}
W_2^+-W_2^-= - 2i \mathfrak s W_1.
\ee
\end{Lemma}

\begin{proof}
Formula \eqref{eq:w_2pm} is an immediate consequence of the definition of $W_3^\nu$.
The second relation comes from \eqref{eq:rhomupm}:
we note that a similar relation has been proved in~\cite{jmpa}.
\end{proof}

\subsection{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)$.

\begin{prop}\label{prop:PU}
Assumptions \ref{ass:UnA}, \ref{ass:SDT} and \ref{ass:eps} hold.
Consider for all $\vert\nu\vert\leq \Lambda_9$
the stretching function $\rho_\nu$ described in Lemma \ref{lm:ronu}. Consider for all $\nu$ a given reference point $x^* \in \mathbb R$ close to
$x_h$ in the ball of analyticity $B(x_h,\Lambda_9)$, as well as a Cauchy data $C_\nu(x^*)=\begin{smallpmatrix}A_\nu(x^*)\cr
B_\nu(x^*)\end{smallpmatrix}$ given through
\be\label{Cauchydata}
C_\nu(x^*)= \pi \begin{pmatrix}
V_\nu'(x^* -X_\nu) & -V_\nu(x^*-X_\nu) \\[3pt]
-U_\nu'(x^*-X_\nu) & U_\nu(x^* -X_\nu)
\end{pmatrix}
G_\nu(x^*)
\begin{pmatrix}
E^\nu_2(x^*)\\[3pt]
W^\nu(x^*)
\end{pmatrix},
\ee
where the matrix $G_\nu(x^*)$ is given by
$$
G_\nu(x)\!=\!
\begin{pmatrix}
0 &\dfrac{\sqrt{x-X_\nu}}{ \sqrt{-sc_\nu(x)(x-X_\nu)} }
\\[15pt]
-
\dfrac{\sqrt{-sc_\nu(x)(x-X_\nu)}}{\sqrt{x-X_\nu}}
&
\dfrac{d}{dx} \biggl(
\dfrac{\sqrt{x-X_\nu}}{ \sqrt{-sc_\nu(x)(x-X_\nu)} } \biggr)
-
\dfrac{a_\nu(x)\sqrt{x-X_\nu}}{ \sqrt{-sc_\nu(x)(x-X_\nu)} }
\end{pmatrix}.
$$
When $\nu$ goes to $0$, this Cauchy data has a finite limit which can be any value in~$\mathbb C^2$. Then the unique solution $(E^\nu_1,E^\nu_2,W^\nu)$ of the regularized system \eqref{eq:sysinit} on $B(x_h,\Lambda_9)\cap\mathbb R$ is given by
\be \label{eq:we2e1}
\begin{cases}
\dpl
W^{\nu}
=
A_\nu W_1^\nu
+
B_\nu W_3^\nu
+ B_\nu\frac1 \pi \log {\rho_\nu(\cdot-X_\nu)}{ } W_1^\nu,
\\[10pt]
\dpl
E_2^\nu
=
A_\nu \Bigl[ \frac{ i \theta \epsilon_{12}^\nu}{d^\nu} W_1^\nu- \frac{\epsilon_{11}^\nu}{d^\nu} \frac d {dx} W_1^\nu
\Bigr]
+B_\nu \Bigl[ \frac{ i \theta \epsilon_{12}^\nu}{d^\nu} W_3^\nu
- \frac{\epsilon_{11}^\nu}{d^\nu} \frac d {dx} W_3^\nu
\Bigr]
\\[8pt]
\dpl
\hfill+
B_\nu \frac1\pi
\log {\rho_\nu(\cdot-X_\nu)}{ }
\Bigl[ \frac{i\theta \epsilon_{12}^\nu}{d^\nu} W_1^\nu
- \frac{\epsilon_{11}^\nu}{d^\nu} \frac d {dx} W_1^\nu
\Bigr]\\[8pt]
\dpl
\hfill-B_\nu
\frac{
1}\pi
\frac{\rho_\nu'(\cdot-X_\nu)
}{\rho_\nu(\cdot-X_\nu)
} \frac{\epsilon_{11}^\nu}{d^\nu} W_1^\nu,
\\[10pt]
\dpl
E_1^\nu
=
A_\nu \Bigl[ \frac{i\theta \epsilon_{11}^\nu}{d^\nu} W_1^\nu
+
\frac{\epsilon_{12}^\nu}{d^\nu} \frac d {dx} W_1^\nu \Bigr]
+
B_\nu \Bigl[ \frac{i\theta \epsilon_{11}^\nu}{d^\nu} W_3^\nu
+
\frac{\epsilon_{12}^\nu}{d^\nu} \frac d {dx} W_3^\nu \Bigr]
\\[8pt]
\dpl
\hfill+
B_\nu \frac1 \pi \log {\rho_\nu(\cdot-X_\nu)}{ }
\Bigl[ \frac{i\theta \epsilon_{11}^\nu}{d^\nu} W_1^\nu
+
\frac{\epsilon_{12}^\nu}{d^\nu} \frac d {dx} W_1^\nu \Bigr]\\[8pt]
\dpl
\hfill+B_\nu\frac{
1}\pi
\frac{\rho_\nu'(\cdot-X_\nu)
}{\rho_\nu(\cdot-X_\nu)
}
\frac{\epsilon_{12}^\nu(x)}{d^\nu}
W_1^\nu.
\end{cases}
\ee
\end{prop}

\begin{proof}
The expression of
\[
W^\nu=\sqrt{-sc_\nu} \left(A_\nu U_\nu(\cdot-X_\nu)
+
B_\nu V_\nu(\cdot-X_\nu)\right)
\]
is immediate from Lemma \ref{lem:Cnu},
since it gives the general solution
\eqref{eq:yw1} of the corresponding second order equation. It yields the first identity of \eqref{eq:we2e1}.

Concerning the second identity of \eqref{eq:we2e1}, one can start from
\[
E_2^\nu= \frac{a_\nu }{c_\nu }\, W^\nu+ \frac1{c_\nu}\, \frac d {dx} W^\nu
\]
by means of \eqref{eq:fond1} and \eqref{eq:matm}.
From \eqref{eq:cM} one has $\sfrac1{c_\nu}=-\sfrac{\epsilon_{11}^\nu}{d^\nu}$
and $\sfrac{a_\nu}{c_\nu}=\sfrac{i\theta \epsilon_{12}^\nu}{d^\nu}$.
So one can write
\be \label{eq:e2e2}
E_2^\nu= \frac{i\theta \epsilon_{12}^\nu}{d^\nu}\, W^\nu
-\frac{\epsilon_{11}^\nu}{d^\nu}\, \frac d {dx} W^\nu.
\ee
To compute $ \frac d {dx} W^\nu$ we directly differentiate all terms in the expression
already obtained for $W^\nu$ (first line of \eqref{eq:we2e1}).
We observe that the derivatives of $A_\nu$ and $B_\nu$ vanish since $A_\nu ' U_\nu(\cdot-x_h) +B'_\nu V_\nu(\cdot-x_h)=0 $, see \eqref{eq:ab1}. It yields
\begin{multline}
\label{eq:dwdw}
\frac d {dx} W^\nu=
A_\nu \frac d {dx} W_1^\nu +
B_\nu \left(\frac d {dx} W_3^\nu + \frac1 \pi \log {\rho_\nu(\cdot-X_\nu)}{ } \frac d {dx} W_1^\nu \right)\\
+
B_\nu\,
\frac{
1}\pi\,
\frac{\rho_\nu'(\cdot-X_\nu)
}{\rho_\nu(\cdot-X_\nu)
}\,W_1^\nu.
\end{multline}
It is then sufficient to plug this expression into \eqref{eq:e2e2}, and to reorganize the sum to get
the $E_2^\nu$ formula in the second line of \eqref{eq:we2e1}.

From \eqref{eq:sysinit}, summing the first equation multiplied by $\epsilon_{11}^{\nu}$ and the second equation by $-\epsilon_{12}^{\nu}$, and from Assumption \ref{ass:SDT},
one gets
\be \label{eq:eeespec2}
E_{1}^\nu=\frac{i\theta \epsilon_{11}^\nu}{d^\nu}\,
W^\nu
+\frac{\epsilon_{12}^\nu }{d^\nu}\, \frac d {dx} W^\nu.
\ee
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^\nu$ and $\frac d {dx} W^\nu$
in \eqref{eq:eeespec2}
to obtain the last part of the claim.
The Cauchy data $C_\nu(x^*)$ is equal to $\begin{smallpmatrix}A_\nu(x^*)\cr
B_\nu(x^*)\end{smallpmatrix}$, where $A_\nu$ and $B_\nu$ are known through \eqref{eq:yw1}--\eqref{eq:ywprime}. This gives
$$\begin{pmatrix}A_\nu(x^*)\cr
B_\nu(x^*)\end{pmatrix}= \pi 
\begin{pmatrix}
V_\nu'(x^* -X_\nu) & -V_\nu(x^*-X_\nu) \\
-U_\nu'(x^*-X_\nu) & U_\nu(x^* -X_\nu)
\end{pmatrix}
\begin{pmatrix}y(x^*)\cr
y'(x^*)\end{pmatrix}.
$$
This linear system can be inverted to obtain $y(x^*),y'(x^*)$.
Then, use Definition \eqref{rescaled} for $y$ as well as the equation on $W^{\nu}, E^{\nu}_2$ given in Definition \ref{def:M} to get $y(x^*), y'(x^*)$ in terms of $W^{\nu}(x^*), E^{\nu}_2(x^*)$, which leads to \eqref{Cauchydata}.
Since $A_\nu(x^*),B_\nu(x^*)$
is arbitrary, the Cauchy data is any vector in $\mathbb C^2$.
\end{proof}

From Formulas \eqref{eq:we2e1} it is clear that
the convergence with respect to $\nu$ is not the same for $E_2^\nu$ and $W_\nu$ on the one hand,
and for $E_1^\nu$ on the other hand.
Indeed the most singular term in $E_2^\nu$ and $W_\nu$ is the logarithm
$\log {\rho_\nu(\cdot-X_\nu)}$: the last term in $E_2^\nu$ is non-singular since
$\sfrac{ \epsilon_{11}^\nu} {\rho_\nu(\cdot-X_\nu)
} $ is the ratio of two terms vanishing at order one at $X_\nu$---so is non-singular.

So $E_2^\nu$ and $W_\nu$
are bounded in $L^p_{\rm loc}(x_h-\Lambda_9,x_h+\Lambda_9)$ for $p<\infty$
and pass to the limit pointwise except at the singularity $x_h=\lim_{\nu\to 0} X_\nu$.
This result was already obtained for the slab geometry
with a completely different technique in \cite[Prop.\,5.14]{jmpa} and is generalized here.
So we state without detail the representation formulas for the limits.

\begin{prop} \label{propo:flim}
Suppose Assumptions \ref{ass:UnA}, \ref{ass:SDT} and \ref{ass:eps} hold.
Consider for all $\vert\nu\vert\leq \Lambda_9$
the stretching function $\rho_\nu$ described in Lemma \ref{lm:ronu}, and the Cauchy data
$C_\nu(x^*)$ (given by \eqref{Cauchydata}) which goes to a finite limit when $\nu$ goes to $0$.
Consider the unique solution $(E^\nu_1,E^\nu_2,W^\nu)$ of the regularized system \eqref{eq:sysinit} on $B(x_h,\Lambda_9)$.

Then for all $1\leq p<\infty$
the $L^p-$limits as $\nu \to 0^\pm$
of $W^\nu$ and $E_2^\nu$ are
\be \label{eq:lim5}
W^\pm=A_0^\pm W_1+B_0^\pm(x)W_3
+ B_0^\pm\,\frac1 \pi\, \log {\rho_0(\cdot-x_h)^\pm}{ } W_1,
\ee
and
\begin{equation}
\label{eq:lim10}
\begin{split}
E_2^\pm=
A_0^\pm \Bigl[ \frac{ i \theta \epsilon_{12}^0}{d^0}\, W_1- \frac{\epsilon_{11}^0}{d^0}\, \frac d {dx} W_1
\Bigr]
&+B_0^\pm \Bigl[ \frac{ i \theta \epsilon_{12}^0}{d^0}\, W_3
- \frac{\epsilon_{11}^0}{d^0}\, \frac d {dx} W_3
\Bigr]
\\
&+
B_0^\pm \frac1\pi
\log {\rho_0(\cdot-x_h)^\pm}{ }
\Bigl[ \frac{i\theta \epsilon_{12}^0}{d^0}\, W_1
- \frac{\epsilon_{11}^0}{d^0}\, \frac d {dx} W_1
\Bigr]\\
&-B_0^\pm\,
\frac{
1}\pi\,
\frac{\rho_0'(\cdot-x_h)
}{\rho_0(\cdot-x_h)
} \,\frac{\epsilon_{11}^0}{d^0}\, W_1,
\end{split}
\end{equation}
where $\log {\rho_0(\cdot-x_h)}{ } ^\pm
$ is defined in \eqref{eq:rhomupm}.
The limits also hold in \hbox{$
\mathbb R\!-\!\{x_h\}\!\cap\!B(x_h,\Lambda_9)$} pointwise.
\end{prop}

The situation is different for the $E_1^\nu$ component of \eqref{eq:we2e1} since
the division by $\rho_\nu(\cdot-X_\nu)$ in the last term
yields a singularity as $\nu$ approaches zero because
\[
W_1(x_h)=\sqrt{ \frac{d^0(x_h) }{s \epsilon_{11}'(x_h)} } \neq0.
\]
This behavior is the same that already demonstrated in \cite{jmpa} using a singular
integral equation technique. The difference is that we have now
an explicit representation of this singular behavior for small $\nu$.
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 \cite{jmpa} and
all calculations are now evident starting from the representation
formula \eqref{eq:we2e1} we state the result without details of the proof.
The notations are that
$
\mathfrak s= \sign\left(\sfrac{i \partial_\nu \epsilon_{11}^\nu(x)}{ \partial_x \epsilon_{11}^\nu(x) } \right)_{\left| x=x_h,\nu=0 \right.}$, $D$ is the Dirac mass and
the principal value is defined by
\[
\langle \PV (\alpha),\phi\rangle= \lim_{\tau\to 0^+} \int_{x\not \in [x_h-\tau, x_h +\tau]}
\alpha(x) \phi(x)dx.
\]

\begin{theor} \label{theore:fin0}
Suppose Assumptions \ref{ass:UnA}, \ref{ass:SDT} and \ref{ass:eps} hold.
Consider for all $\vert\nu\vert\leq \Lambda_9$
the stretching function $\rho_\nu$ described in Lemma \ref{lm:ronu}, and $C_\nu(x^*)$ (given by \eqref{Cauchydata}) which goes to a finite limit when $\nu\to 0$.
Consider the unique solution $(E^\nu_1,E^\nu_2,W^\nu)$ of the regularized system \eqref{eq:sysinit} on $B(x_h,\Lambda_9)$.
Then the limit of $E_1^\nu$, in $\mathcal{D}'(B(x_h,\Lambda_9))$,
\begin{align*}
E_1^\pm=
A_0^\pm \Bigl[ \frac{i\theta \epsilon_{11}^0}{d^0}\;W_1
+
\frac{\epsilon_{12}^0}{d^0}\; \frac d {dx} W_1 \Bigr]
&+
B_0^\pm \Bigl[ \frac{i\theta \epsilon_{11}^0}{d^0}\;W_3
+
\frac{\epsilon_{12}^0}{d^0}\; \frac d {dx} W_3 \Bigr]\\
&+
B_0^\pm\; \frac1 \pi \;\log {\rho_0(\cdot-x_h)^\pm}{ }
\Bigl[ \frac{i\theta \epsilon_{11}^0}{d^0}\;W_1
+
\frac{\epsilon_{12}^0}{d^0}\; \frac d {dx} W_1 \Bigr]
\\
&{}\pm i \mathfrak s\; B_0^\pm\;\frac{\epsilon_{12}^0(x_h)}{d^0}\;W_1(x_h)\, D(\cdot-x_h)\\
&+\PV \Bigl(B_0^\pm\;\frac{
1}\pi\;
\frac{\rho_0'(\cdot-x_h)
}{\rho_0(\cdot-x_h)
}\;
\frac{\epsilon_{12}^0(x)}{d^0}\;
W_1\Bigr).
\end{align*}
\end{theor}

Of course the right hand side terms on the first line converge also
in the space \hbox{$L^p_\mathrm{loc}(x_h-\Lambda_9,x_h+\Lambda_9)$} as in the previous proposition, and the pointwise limit holds
away from the hybrid singularity $x_h$.

\subsection{Resonant heating term}

The heating was defined in Formula \eqref{eq:enbal}.
For $a<b$ in $B(x_h,\Lambda_9)\cap\mathbb R^*$ and
$\nu\in (-\Lambda_9,\Lambda_9)$
we note
\begin{equation} \label{eq:heahea}
\begin{split}
Q^\nu(a,b)&=\Im \int_a^b \left(
\epsilon^\nu(x) \bE^\nu(x), { \bE^\nu(x) } \right)dx\\
&=-\Im\bigl(W^\nu(b) \overline{ E_2^\nu(b) } \bigr)
+
\Im \bigl(W^\nu(a) \overline{E_2^\nu(a) } \bigr),
\end{split}
\end{equation}
where $(E^\nu_1,E^\nu_2,W^\nu)$ is the regularized solution system \eqref{eq:sysinit} on $B(x_h,\Lambda_9)$.
The limit value for $\nu=0$ is the resonant heating and can be characterized as follows.

\begin{theor}[Resonant heating] \label{theore:fin}
Suppose Assumptions \ref{ass:UnA}, \ref{ass:SDT} and \ref{ass:eps} hold.
Consider for all $\vert\nu\vert\leq \Lambda_9$
the stretching function $\rho_\nu$ described in Lemma \ref{lm:ronu}, and $C_\nu(x^*)$ given by \eqref{Cauchydata} which goes to a finite limit when $\nu$ goes to $0$.

For
$a,b\in B(x_h,\Lambda_9)\cap \mathbb R\backslash\{x_h \}$,
one can pass to the limit with respect to $\nu$:
\begin{enumeratei}
\item
If $(a-x_h)(b- x_h)>0$, then $\lim_{\nu \to 0} Q^\nu(a,b)=0$.
\item
If $a<x_h<b$, then
\be \label{eq:qhqh}
\displaystyle
\lim_{\nu \to 0^\pm}Q^\nu(a,b)
= \pm \frac1\pi |B_0^\pm(x_h)|^2 \left| \epsilon_{12}^0(x_h) \right| ^2 \sign \bigl(\Im \left(\partial_\nu \epsilon_{11}^0 (x_h)\right) \bigr).
\ee
\end{enumeratei}
\end{theor}

\begin{proof}
There are many possibilities to compute this limit value from
the previous representation formulas.

\subsubsection*{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 $\Im \left(\epsilon^0(x)\bX,
{\bX} \right) =0 $ for any complex vector $\bX$ since $ \epsilon^0(x)$
is Hermitian.
\subsubsection*{Second case}
Our proof of the result relies on three arguments which
are a decom\-position of the solution between a regular part and a singular part,
an adapted decom\-position of the resonant heating $Q^\nu(a,b)$
where some terms naturally tend to zero,
and a careful
analysis of the remaining singular part.

\subsubsection*{Step 1: decomposition of the solution}
Let us decompose the electric field in two parts, a regular part $\bR^\nu$ and a singular part $\bS^\nu$. The regular part $\bR^\nu$ is regular enough, for example bounded
in $L^2(-\Lambda_9,\Lambda_9)$ uniformly with respect to $\nu$ and with a
pointwise limit for $x\neq x_h$. The singular part $\bS^\nu$ is such that $(x-X_\nu)\bS^\nu$ is bounded in $L^{\infty}(-\Lambda_9,\Lambda_9)$.
One can choose
$$
\bS^\nu(x) = \Bigl(B_\nu\;\frac{
1}\pi\;
\frac{\rho_\nu'(\cdot-X_\nu)
}{\rho_\nu(\cdot-X_\nu)
}\;
\frac{\epsilon_{12}^\nu}{d^\nu}\;
W_1^\nu, 0 \Bigl).
$$
Consider
on the one hand
\eqref{eq:dnu2:rab}--\eqref{eq:dnu2.2} and \eqref{eq:xnuxnu}
which yield an estimate of the shift of $X_\nu$ in the complex plane, and
on the other hand of the representation formulas
for~$E_2^\nu$ and $E_1^\nu$. Then
one has the estimates in $L^2$
\be \label{eq:rnusnu}
\left\| \bR_\nu \right\|_{L^2(a,b)} \leq C
\quad\text{and}\quad \left\| \bS_\nu \right\|_{L^2(a,b)} \leq \frac{C}{\sqrt{ |\nu |}},
\ee
for some constant $C>0$ independent of $\nu$.
The
regular part is estimated using the integrability in space of the logarithm function.
The singular part is estimated thanks
to the basic bound
$
\int_a^b \left| \sfrac1{x+i\nu} \right|^2dx \leq \sfrac{\pi}{| \nu|}$.

\subsubsection*{Step 2: decomposition of $Q(a,b)$}
Let us introduce the Hermitian matrix
\be\label{def:Dnu}
D^\nu = \frac{ \epsilon^\nu - \left(\epsilon^{\nu}\right)^* }{2i\nu} =(D^\nu)^*.
\ee
This Hermitian matrix is uniformly bounded in view of all assumptions, that is
\be \label{eq:rnusnu2}
\left\| D^\nu \right\|_{L^\infty(a,b)}\leq K,
\ee
for a constant uniform with respect to $\nu$.
The decomposition writes
\begin{align*}
Q(a,b)&=\Im \int_a^b \left(
\epsilon^\nu(x) \bE^\nu(x), { \bE^\nu(x) } \right) dx\\
&=\nu \Re \int_a^b \left(
D^\nu(x) \left(\bR^\nu(x)+
\bS^\nu(x)
\right), { \bR^\nu(x)+
\bS^\nu(x)
} \right) dx\\
&=\nu \Re \int_a^b \left(
D^\nu(x) \bR^\nu(x), { \bR^\nu(x)} \right) dx
+
2 \nu \Re \int_a^b \left(
D^\nu(x) \bR^\nu(x), { \bS^\nu(x)} \right) dx\\
&\hspace*{6cm}+
\nu \Re \int_a^b \left(
D^\nu(x) \bS^\nu(x), { \bS^\nu(x)} \right) dx\\
&=\nu \Re \int_a^b \left(
D^\nu(x) \bR^\nu(x), { \bR^\nu(x)} \right) dx
+
2 \nu \Re \int_a^b \left(
D^\nu(x) \bR^\nu(x), { \bS^\nu(x)} \right) dx\\
&\hspace*{6cm}+
\Im \int_a^b \left(
\epsilon^\nu(x) \bS^\nu(x), { \bS^\nu(x)} \right) dx.
\end{align*}
\subsubsection*{Step 3: estimate of all terms}
The two first terms of this expression tend to zero with~$\nu$ in view
of \eqref{eq:rnusnu} and \eqref{eq:rnusnu2}.
Let us study the remaining term, which is
\begin{equation} \label{eq:gnu:1}
\begin{split}
G_\nu&= \Im \int_a^b\left(
\epsilon^\nu(x) \bS^\nu(x), { \bS^\nu(x) } \right)\\
&=\Im
\int_a^b \epsilon_{11}^\nu(x)
\Bigl| B_\nu(x)\;\frac{
1}\pi\;
\frac{\rho_\nu'(x-X_\nu) }{\rho_\nu(x-X_\nu)
}\;
\frac{\epsilon_{12}^\nu(x)}{d^\nu(x)}\;
W_1^\nu (x) \Bigr|^2 dx.
\end{split}
\end{equation}
Define
\[
z_\nu(x)=\frac{\epsilon_{11}^\nu(x) }{x-X_\nu}\quad\text{and}\quad K_{\nu}(x)=(x-X_\nu)B_\nu(x)\;\frac{
1}\pi\;
\frac{\rho_\nu'(x-X_\nu) }{\rho_\nu(x-X_\nu)
}\;
\frac{\epsilon_{12}^\nu(x)}{d^\nu(x)}\;
W_1^\nu (x).
\]
These two terms are regular and have a natural limit as $\nu\to 0^\pm$. One has
\begin{equation} \label{eq:gnu:2}
G_\nu=\Im \int_a^bz_{\nu}(x)\vert K_{\nu}(x)\vert^2\;\frac{x-X_\nu}{\vert x-X_{\nu}\vert^2}\;dx.
\end{equation}
Introduce
\[
M_\nu=\Im \int_a^bz_{\nu}(X_\nu)\vert K_{\nu}(X_\nu)\vert^2\;\frac{x-X_\nu}{\vert x-X_{\nu}\vert^2}\;dx.
\]
One can check that there exists a regular function $H_{\nu}(x)$ such that
\[
z_{\nu}(x)\vert K_{\nu}(x)\vert^2-z_{\nu}(X_\nu)\vert K_{\nu}(X_\nu)\vert^2=(x-X_\nu)H_\nu(x).
\]
Note that $M_\nu(x)\in \mathbb R$. We thus deduce that
\[
G_\nu-M_\nu=\Im \int_a^bH_\nu(x)\;\frac{(x-X_\nu)^2}{\vert x-X_\nu\vert^2}\;dx.
\]
As the limit exists in $L^1$, we use the dominated convergence theorem to deduce that
\begin{equation} \label{eq:gnu:3}
\lim_{\nu=0^\pm}(G_\nu-M_\nu)= \Im \int_a^bH_0(x)dx=0.
\end{equation}
Therefore it is sufficient to pass to the limit
for $M_\nu$. Use
\[
\frac{x-X_\nu}{\vert x-X_{\nu}\vert^2}=\frac{x-\Re X_\nu-i\Im X_{\nu}}{(x-\Re X_\nu)^2+(\Im X_{\nu})^2}.
\]
The change of variable $t=\sfrac{x-\Re X_\nu}{\Im X_{\nu}}$ yields
$$
M_\nu=\Im \biggl(z_{\nu}(X_{\nu})\vert K_{\nu}(X_{\nu})
\vert^2\int_{\psfrac{a-\Re X_{\nu}}{\Im X_{\nu}}}^{\psfrac{b-\Re X_{\nu}}{\Im X_{\nu}}}\frac{t-i}{t^2+1}\;dt
\biggl).
$$
After integration
\be \label{eq:ccll1}
\begin{split}
M_\nu&=\Im \left(z_{\nu}(X_{\nu})\vert K_{\nu}(X_{\nu})\vert^2[\mathfrak A^{\nu}-i\mathfrak B^{\nu}] \right)\\
&=
\Im \left(z_{\nu}(X_{\nu})\vert K_{\nu}(X_{\nu}) \vert^2\mathfrak A^{\nu}\right)
-
\Im \left(z_{\nu}(X_{\nu})\vert K_{\nu}(X_{\nu})\vert^2 i\mathfrak B^{\nu} \right),
\end{split}
\ee
where
\begin{align*}
\mathfrak A^{\nu}&=\frac12\ln\frac{(b-\Re X_{\nu})^2+(\Im
X_{\nu})^2}{(a-\Re X_{\nu})^2+(\Im X_{\nu})^2},\\
\mathfrak B^{\nu}&=\tan^{-1}\Bigl(\frac{b-\Re X_{\nu}}{\Im X_{\nu}} \Bigr)-\tan^{-1} \Bigl(\frac{a-\Re X_{\nu}}{\Im X_{\nu}} \Bigr).
\end{align*}
We claim that all terms have a natural limit in \eqref{eq:ccll1}. One has
$$
\lim_{\nu=0^\pm}z_\nu(X_\nu)= \partial_x\epsilon_{11}^0(x_h)\in \mathbb R.
$$
Moreover,
$
\mathfrak A^{\nu}
\to \frac12\ln\sfrac{(b-x_h)^2}{(a-x_h)^2}\in \mathbb R$.
So
$$
\lim_{\nu=0^\pm }\left(\Im \left(z_{\nu}(X_{\nu})\vert K_{\nu}(X_{\nu}) \vert^2\mathfrak A^{\nu}\right) \right)=0.
$$
When $a<x_h<b$, the limit of $\psfrac{a-\Re X_{\nu}}{\Im X_\nu}$ is $-\infty$ when $\Im X_\nu$ goes to $0$ by positive values, and the limit of $\psfrac{a-\Re X_{\nu}}{\Im X_\nu}$ is $+\infty$ when $\Im X_\nu$ goes to $0$ by negative values.
The sign of the imaginary part $ \sign(\Im X_{\nu})$ is given by
\eqref{eq:xnuxnu}
$$
\sign(\Im X_{\nu})= - \sign \biggl(\Im \Bigl(
\frac{\partial_\nu \epsilon_{11}^0(x_h) }{\partial_x \epsilon_{11}^0(x_h)}
\Bigr)
\biggr) \sign (\nu).
$$
We thus deduce
\[
\lim_{\nu=0^\pm} \mathfrak B^{\nu} = \mp \pi \sign \biggl(\Im \Bigl(
\frac{\partial_\nu \epsilon_{11}^0(x_h) }{\partial_x \epsilon_{11}^0(x_h)}
\Bigr)
\biggr).
\]
Hence
\begin{multline*}
\lim_{\nu=0^\pm }\Im \left(z_{\nu}(X_{\nu})\left| K_{\nu}(X_{\nu}) \right|^2i \mathfrak B^{\nu}\right)\\
=
\partial_x\epsilon_{11}^0(x_h) \Bigl |
B_0^{\pm}(x_h)W_1(x_h)\;\frac{\epsilon_{12}^0(x_h)}{\pi d^0(x_h)}
\Bigr|^2
\biggl(
\mp \pi \sign \biggl(\Im \Bigl(
\frac{\partial_\nu \epsilon_{11}^0(x_h) }{\partial_x \epsilon_{11}^0(x_h)}
\Bigr)
\biggr)
\biggr).
\end{multline*}
Using $d^0(x_h)=\epsilon_{12}^0(x_h)^2 $ thanks
to $\epsilon_{12}^0(x_h)=-\epsilon_{21}^0(x_h)$, one obtains that
$$
\lim_{\nu=0^\pm } M_\nu= \pm \frac1\pi \partial_x\epsilon_{11}^0(x_h)
\Bigl|
\frac{ B_0^{\pm}(x_h) W_1(x_h)}
{\epsilon_{12}^0(x_h) }
\Bigr|^2
\sign \biggl(\Im \Bigl(
\frac{\partial_\nu \epsilon_{11}^0(x_h) }{\partial_x \epsilon_{11}^0(x_h)}
\Bigr)
\biggl).
$$
One more simplification is
$
\left| W_1(x_h) \right|^2 = \sfrac{ \left| \epsilon_{12}^0(x_h) \right| ^2}{
\left| \partial_x \epsilon_{11}^0(x_h) \right|}
$ as a consequence of \eqref{eq:byby2}--\eqref{eq:byby}.
One gets
\begin{equation} \label{eq:gnu:4}
\lim_{\nu=0^\pm} M_\nu=\pm
\frac1\pi\, |B_0^\pm(x_h)|^2 \left| \epsilon_{12}^0(x_h) \right| ^2 \sign \left(\Im \left(\partial_\nu \epsilon_{11}^0 (x_h)\right) \right).
\end{equation}
The claim \eqref{eq:qhqh} is obtained as a consequence of \eqref{eq:gnu:1}--\eqref{eq:gnu:2}--\eqref{eq:gnu:3}--\eqref{eq:gnu:4}.
The proof is ended.
\end{proof}

An interesting and evident corollary is
the following.

\begin{cor} \label{cor:last}
Consider any dielectric tensor satisfying \eqref{eq:epsilonmu} and Assumption \ref{ass:SDT}, not necessarily
with the exact form
\eqref{eq:hyp1} but with the linearized form \eqref{eq:princip}--\eqref{eq:dede22}.
Assume $\epsilon^0$ is a smooth Hermitian matrix, and assume $x_h=0$ is a hybrid resonance.
Define the local dissipation tensor
$$
\mathbb D= -i \left(\partial_\nu \epsilon \right)^0(x_h) =
\begin{pmatrix}
d_{1} & d_{2} \\
-d_{2} & d_{1}
\end{pmatrix}
\quad\text{with $d_1\neq 0$ and $d_2\in i \mathbb R$.}
$$

Then the value of the resonant heating is independent of $\mathbb D$.
The same for the pointwise limits of the electric field, the magnetic field and the numerical value
of the resonant
heating.
\end{cor}

\begin{proof}
It is sufficient to realize that the major assumption that was made, that is
\eqref{eq:dnu2}--\eqref{eq:dnu2.2}, is independent of the exact value of $d_1\neq 0$, and that
the other coefficients in $\mathbb D$ do not appear in the proof of the previous theorem. This is also clear in view of the value of the resonant heating \eqref{eq:qhqh}.
\end{proof}

An interesting consequence is that the initial dielectric
tensor \eqref{eq:hyp1} can be replaced by the linear approximation
$\epsilon^0+i \nu \mathbb D$ with $\mathbb D$ defined in \eqref{eq:dede}.
If one is interested only in the limit value, which is a reasonable assumption
for fusion plasmas where $\nu\approx 10^{-7}$ may be encountered, this is a
valid assumption.

One notices that the positivity condition $d_1^2>|d_2|^2$ that stems
from the physics \eqref{eq:dede22} 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
\be \label{eq:ddspec}
\mathbb D^{\rm art}=
\begin{pmatrix}
d_{1} & d_{2} \\
-d_{2} & d_{3}
\end{pmatrix}
\quad\text{with $d_1\neq 0$, $d_2\in i\mathbb R$ and $d_3\in \mathbb R$.}
\ee
Since only $d_1$ really matters in $\epsilon_{11}^\nu$ which is the only
quantity in the denominators visible in the matrix $\bM_\nu$ defined in \eqref{eq:mm}, the results
of this work are still valid.
We did not elaborate on such a matrix in this work since the condition
$d_3\neq d_1$ is non-physical and of less interest, nevertheless it must be stressed that the numerical results in next section
with tensor $\mathbb D_2$ validate this statement.

\section{Numerical illustrations} \label{sec:num}
\subsection{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 \cite{lima:phd}.
The code solves the system under the form
\eqref{eq:fond1} 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 $\epsilon^0$ for $\nu=0$.
The three dielectric tensors are defined by three different dissipation tensors.
The first one denoted as ${\mathbb D}_0$ corresponds to the
exact formula \eqref{eq:hyp1}.
The second one
$\mathbb D_1 = \begin{smallpmatrix} 1 &0 \\0 & 1\end{smallpmatrix}$
corresponds to the linear approximation
$\epsilon^0+i \nu \mathbb D_1$. One can notice that $\mathbb D_1=\mathbb I$
corresponds to the limiting absorption principle that was studied in \cite{jmpa}.
The last
$\mathbb D_2=\mathbb D_2^* = \begin{smallpmatrix} 1 &-2 i \\2 i & 5\end{smallpmatrix}$
is a Hermitian and positive matrix as is required on the physical basis (see \eqref{eq:dede} for example).
These three methods satisfy the main hypothesis
of Corollary \ref{cor:last} that is
$
\left(\mathbb D_{i} \right)_{11}>0$ for $i=0,1,2$.
In each case, the matrix $\bM_\nu(x,\theta)$ is constructed and
the problem solved numerically in the interval $\left[ -5,5 \right]$
with arbitrary Cauchy data $E_2^\nu(5)=1$ and $W^\nu(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.

\begin{table}[h]
\begin{center}
\begin{tabular}{ | l || cccccc |}
\hline
$\nu $& $.5$ & $.5\times 10^{-1} $ & $.5\times 10^{-2}$ &
$.5\times 10^{-3}$ & $.5\times 10^{-4}$ & $.5\times 10^{-5}$ \\
\hline
\hline
${\mathbb D}_0$ & 9898 & 19484 & 18724 & 18623 & 18604 & 18590
\\
$\mathbb D_1$ & 1534268 & 36976 & 20272 & 18771 &18613 & 18588 \\
$\mathbb D_2$ & 109836517 & 89603 & 22761 & 19001 & 18637 & 18591 \\
\hline
\end{tabular}
\end{center}
\smallskip
\caption{Values of the heating calculated with Formula
\eqref{eq:heahea} with different $\nu$ and different dissipation tensors.
As predicted by the theory, the limit resonant heating in the last column is independent
of the dissipation tensor.}
\label{tab:1}
\end{table}

\begin{figure}
\begin{center}
\includegraphics[width=5cm]{./despres-et-al_figs/Xmode1DBesselD_1nu0Re}\hspace*{1cm}\includegraphics[width=5cm]{./despres-et-al_figs/Xmode1DBesselD_2nu0Re}
\end{center}
\begin{center}
\includegraphics[width=5cm]{./despres-et-al_figs/Xmode1DBesselD_1nu1Re}\hspace*{1cm}\includegraphics[width=5cm]{./despres-et-al_figs/Xmode1DBesselD_2nu1Re}
\end{center}
\begin{center}
\includegraphics[width=5cm]{./despres-et-al_figs/Xmode1DBesselD_1nu2Re}\hspace*{1cm}\includegraphics[width=5cm]{./despres-et-al_figs/Xmode1DBesselD_2nu2Re}
\end{center}
\begin{center}
\includegraphics[width=5cm]{./despres-et-al_figs/Xmode1DBesselD_1nu3Re}\hspace*{1cm}\includegraphics[width=5cm]{./despres-et-al_figs/Xmode1DBesselD_2nu3Re}
\end{center}
\caption{Real part of $V^\nu=E_2^\nu$,
$W^\nu$ and $U^\nu=E_1^\nu$. Parameters are as follows:
$\mathbb D_1$ on the left, $\mathbb D_2$ on the right; from top to bottom
$\nu =.5,.5\cdot 10^{-1},.5\cdot 10^{-2},.5\cdot 10^{-3}$.
One observes the convergence to the same limit.
}
\label{fig:res1}
\end{figure}

\begin{figure}
\begin{center}
\includegraphics[width=5cm]{./despres-et-al_figs/Xmode1DBesselD_1nu0Im}\hspace*{1cm}\includegraphics[width=5cm]{./despres-et-al_figs/Xmode1DBesselD_2nu0Im}
\end{center}
\begin{center}
\includegraphics[width=5cm]{./despres-et-al_figs/Xmode1DBesselD_1nu1Im}\hspace*{1cm}\includegraphics[width=5cm]{./despres-et-al_figs/Xmode1DBesselD_2nu1Im}
\end{center}
\begin{center}
\includegraphics[width=5cm]{./despres-et-al_figs/Xmode1DBesselD_1nu2Im}\hspace*{1cm}\includegraphics[width=5cm]{./despres-et-al_figs/Xmode1DBesselD_2nu2Im}
\end{center}
\begin{center}
\includegraphics[width=5cm]{./despres-et-al_figs/Xmode1DBesselD_1nu3Im}\hspace*{1cm}\includegraphics[width=5cm]{./despres-et-al_figs/Xmode1DBesselD_2nu3Im}
\end{center}
\caption{Imaginary part of $V^\nu=E_2^\nu$,
$W^\nu$ and $U^\nu=E_1^\nu$. Parameters are as follows:
$\mathbb D_1$ on the left, $\mathbb D_2$ on the right;
from top to bottom $\nu =.5,.5\cdot 10^{-1},.5\cdot 10^{-2},.5\cdot 10^{-3}$.
One observes the convergence to the same limit.
}
\label{fig:res2}
\end{figure}

\begin{figure}
\begin{center}
\includegraphics[width=5cm]{./despres-et-al_figs/Xmode1DBesselDtnu0Re}\hspace*{1cm}\includegraphics[width=5cm]{./despres-et-al_figs/Xmode1DBesselDtnu0Im}
\end{center}
\begin{center}
\includegraphics[width=5cm]{./despres-et-al_figs/Xmode1DBesselDtnu1Re}\hspace*{1cm}\includegraphics[width=5cm]{./despres-et-al_figs/Xmode1DBesselDtnu1Im}
\end{center}
\begin{center}
\includegraphics[width=5cm]{./despres-et-al_figs/Xmode1DBesselDtnu2Re}
\hspace*{1cm}\includegraphics[width=5cm]{./despres-et-al_figs/Xmode1DBesselDtnu2Im}
\end{center}
\begin{center}
\includegraphics[width=5cm]{./despres-et-al_figs/Xmode1DBesselDtnu3Re}
\hspace*{1cm}\includegraphics[width=5cm]{./despres-et-al_figs/Xmode1DBesselDtnu3Im}
\end{center}
\caption{Real part of $V^\nu=E_2^\nu$,
$W^\nu$ and $U^\nu=E_1^\nu$ on the left, imaginary part on the right.
The diffusion tensor is the same and equal to ${\mathbb D}_0$.
Regularization parameter from top to bottom $\nu =.5,.5\cdot\nobreak 10^{-1},.5\cdot 10^{-2},.5\cdot 10^{-3}$.
For the smallest value of $\nu$, the results are almost the same
as in Figures \ref{fig:res1} and \ref{fig:res2} computed with the linear models
$\mathbb D_1$ and $\mathbb D_2$.}
\label{fig:res3}
\end{figure}

\vspace*{-\baselineskip}%
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^\nu(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 $\nu>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 $\nu$ goes to the $0$ of the heating term appears clearly in Table \ref{tab:1}.
In Figure \ref{fig:res1} we display the real part of $E_1^\nu,E_2^\nu,W^\nu$
for four different values of $\nu$, and for the linear models with
$\mathbb D_1$ and $\mathbb D_2$.
We observe the numerical convergence: even if the results for $\nu=0.5$
are very different, the functions converge numerically
to the same limit for $\nu=0.5\times 10^{-3}$.
Figure \ref{fig:res2}
confirms this analysis for the imaginary part of the electromagnetic field.
Finally Figure \ref{fig:res3} shows the same results but computed with the
``real'' dissipation tensor. One clearly observes the same limit
as in Figures~\ref{fig:res1} and~\ref{fig:res2}, perhaps with
smaller oscillations
for larger values of $\nu$.

\subsection{A short conclusion}

All our numerical results confirm the singular solution predicted by the theoretical results.
It allows us now to have a comprehensive understanding
of the local structure of the hybrid singular solutions.
The connection of the local structure near the singularity
to the condition at infinity, this problem was addressed in \cite{jmpa}, can also be addressed starting with the tools
of the present work and
using the methods developed in \cite{LWZ}.

\appendix

\section{Proof of Lemma \ref{lemme:rana}}
\label{annex}

\begin{proof}
From Definition \ref{df:rmu}, the function $R_\nu(x)$ is a priori a Laurent series in $x$, with bounded coefficients in $\nu$. Hence we will show that
there exists $\Lambda_6>0$ such that
$R_\nu(x)$ is uniformly bounded, \ie bounded for $x\in B(x_h,\Lambda_6)$ uniformly for $\nu\in B(0,\Lambda_6)$.
And such a $\Lambda_6$ can be replaced by $\Lambda_3$ in case this one was smaller. This will prove that it is
an analytic function with a uniform radius of convergence.

The function $R_\nu$ is the product of $x-X_\nu$ times a contribution
which is the sum of four terms.
The first term is $ a_\nu^2 + b_\nu c_\nu=-\mathfrak D_\nu=c_\nu+\theta^2$, with $\mathfrak D_\nu$ defined in \eqref{eq:detdet}, which is at most $O\left((\epsilon_{11}^{\nu})^{-1}\right)$. The product by $x-X_\nu$
is bounded and admits a finite limit for $x\to X_\nu$.
The second term depends on $\sfrac{a_\nu}{ c_\nu} $ which is locally bounded, so does
not yield any difficulty either: its contribution to $R_\nu$ is bounded.
The third term can be written as
$$
\frac{\sqrt{x-X_{\nu}}}{\sqrt{-sc_{\nu}(x)(x-X_{\nu})}}\;
\frac {d^2}{dx^2}\biggl(
\frac{\sqrt{-sc_{\nu}(x)(x-X_{\nu})}}{\sqrt{x-X_{\nu}}}
\biggl) =
- \frac12 \Bigl(\frac{c'_{\nu}(x)}{c_{\nu}(x)} \Bigr)'+\frac14 \Bigl(\frac{c'_{\nu}(x)}{c_{\nu}(x)}
\Bigr)^2.
$$
Considering the bounded function $k_{\nu}(x)=c_{\nu}(x)(x-X_\nu)$, such that
$$k_{\nu}(x)=-\frac{d^{\nu}(x)}{\epsilon_{11}^{\nu}(x)}\;(x-X_\nu)=-\frac{d^{\nu}(x)}{\int_0^1\partial_x\epsilon_{11}^{\nu}(X_\nu+s(x-X_\nu))ds}.$$
Note that $k_{\nu}$ is analytic in $B(x_h,\Lambda_5)$, whit $\Lambda_5$ defined in Lemma \ref{lemma:3lam}, and that
\[
k_{\nu}(X_\nu)=-\frac{d^{\nu}(X_\nu)}{\partial_x\epsilon_{11}^{\nu}(X_\nu)}=-\frac{(\epsilon_{12}^{\nu}(X_\nu))^2}{\partial_x\epsilon_{11}^{\nu}(X_\nu)},
\]
hence $k_\nu$ is analytic and non-vanishing in a ball $B(x_h,\Lambda_6)$. In addition, one checks that $\sfrac{k'_\nu}{k_{\nu}}$ is analytic in $B(x_h,\Lambda_6)$ hence bounded and its derivative is also bounded. Hence one has
\begin{align*}
\frac{\sqrt{x-X_{\nu}}}{\sqrt{-sc_{\nu}(x)(x-X_{\nu})}}&\;
\frac {d^2}{dx^2}\biggl(
\frac{\sqrt{-sc_{\nu}(x)(x-X_{\nu})}}{\sqrt{x-X_{\nu}}}
\biggr)\\
&\hspace*{-5mm}=
- \frac12 \left(
\frac{k'_{\nu}(x)}{k_{\nu}(x)}
\right)'-\frac12 \frac{1}{(x-X_\nu)^2}+\frac14
\left(\frac{k'_{\nu}(x)}{k_{\nu}(x)}-\frac{1}{x-X_\nu}\right)^2\\
&\hspace*{-5mm}= - \frac14 \frac{1}{(x-X_\nu)^2}-\frac12\frac{k'_{\nu}(x)}{k_{\nu}(x)}\frac{1}{x-X_\nu}
- \frac12 \left(
\frac{k'_{\nu}(x)}{k_{\nu}(x)}\right)'+\frac14
\left(
\frac{k'_{\nu}(x)}{k_{\nu}(x)}\right)^2.
\end{align*}
It yields that
\[
\frac{\sqrt{x-X_{\nu}}}{\sqrt{-sc_{\nu}(x)(x-X_{\nu})}}\;
\frac {d^2}{dx^2}\biggl(
\frac{\sqrt{-sc_{\nu}(x)(x-X_{\nu})}}{\sqrt{x-X_{\nu}}}
\biggr)
=
-\frac1{ 4(x-X_\nu)^2 }
+ \frac{ \tau_\nu(x)}{ x-X_\nu },
\]
where the function $\tau_\nu$ is uniformly bounded in the sense that
$
|\tau_\nu(x)| \leq C $, \ie bounded for $x\in B(x_h,\Lambda_6)$ uniformly for $\nu\in B(0,\Lambda_6)$.
It shows that
the most singular part of the two last terms cancel
in the definition of $R_\nu$.
The exact value of $R_\nu(X_\nu)$ may be computed as follows
using $d^\nu$ defined in \eqref{eq:dnunu}.
From $\sfrac{a_{\nu}}{c_{\nu}}=i\theta \sfrac{\epsilon_{12}^{\nu}}{d^{\nu}}$, one has
by elimination of $\epsilon_{22}^\nu=\epsilon_{11}^\nu$
$$
\Bigl(\frac{a_{\nu}}{c_{\nu}} \Bigr)'(X_\nu)=-i\theta\;\frac{(\epsilon_{12}^{\nu})'}{(\epsilon_{12}^{\nu})^2}(X_{\nu}).
$$
To obtain an explicit expression for the term $R_\nu(X_\nu)$ we check that
$$
\frac{k'_\nu}{k_{\nu}}(X_\nu)=2\frac{(\epsilon_{12}^{\nu})'}{\epsilon_{12}^{\nu}}(X_\nu)-\frac12 \frac{(\epsilon_{11}^{\nu})''}{(\epsilon_{11}^{\nu})'}(X_\nu).
$$
It completes the proof since
\begin{align*}
R_\nu(X_\nu)&= k_{\nu}(X_{\nu})\biggl(1-\Bigl(\frac{a_{\nu}}{c_{\nu}} \Bigr)'\biggr)(X_\nu)-\frac12\frac{k'_\nu}{k_{\nu}}(X_\nu)\\
&=
\left[- \frac{i\theta(\epsilon_{12}^\nu)'+(\epsilon_{12}^{\nu})^2}{(\epsilon_{11}^{\nu})'}-\frac{(\epsilon_{12}^{\nu})'}{\epsilon_{12}^{\nu}}+\frac14 \frac{(\epsilon_{11}^{\nu})''}{(\epsilon_{11}^{\nu})'}\right](X_\nu).\qedhere
\end{align*}
\end{proof}

\section{Analytic solutions}

In this section, we intend to construct analytic solutions to the cold plasma model solution. Analytic solutions are of interest because they can be used as reference solutions for testing numerical codes. The solutions we construct here are easily computable, using special functions. They improve model solutions of the literature (Chen-White, \cite{CW}, \cite{chen:lipton} for example) in the sense that they have a richer behavior in the neighborhood of the points studied (namely the hybrid resonance mode).

\subsection{A first family without singularity}

The first case of interest is such that the function $s$ in equation \eqref{eq:g} is identically $0$.
It greatly simplifies the algebra.
In this case $\rho(x)=\spfrac{x}{1+\ell_0^{-1}x}$, $\ell_0$ being an arbitrary length, does not depend on $\nu$.
The matrix $\bM_\nu$ shall depend on~$\nu$ through an arbitrary choice of $X_{\nu}$ depending on~$\nu$, (such as $X_{\nu}=-i\nu$).
Consider
\begin{align*}
c_{\nu}(x)&=-
\frac{\rho'(x-X_{\nu})}{\rho (x-X_{\nu}) }
=-\frac{1}{(x-X_\nu)(1+\ell_0^{-1}(x-X_\nu))},\\
a_{\nu}(x)&=-c_\nu(x),\\
b_{\nu}(x)&=-c_\nu(x)-\frac14\frac{\lambda_{\nu}^2}{(1+\ell_0^{-1}(x-X_\nu))^2},
\end{align*}
where $\lambda$, $\ell_0$ are arbitrary constants.
For the simplicity of notations, the sign in Definition \ref{df:W13} is assume to be positive, that is $s=1$.
We define $U^0_{\nu}$ and $V^0_{\nu}$ the functions
\begin{equation} \label{eq:logi1}
\begin{cases}
U^0_\nu(z) = \sqrt{\dfrac{\rho(z) }{\rho'(z)}}\;
J_0
\bigl(\lambda \sqrt{ \rho(z)} \bigr),\\[15pt]
V^0_\nu(z) = \sqrt{\dfrac{\rho(z) }{\rho'(z) }}
\left[Y_0
\bigl(\lambda \sqrt{ \rho(z)} \bigl) -\dfrac2\pi \log\left(\sfrac{\lambda}2 \right)\;
J_0
\bigl(\lambda \sqrt{ \rho(z)} \bigr) \right].
\end{cases}
\end{equation}
For $\lambda=\lambda_{\nu}$, the functions $U^0_{\nu}$ and $V^0_{\nu}$ correspond to
the functions $U_{\nu}$ and $V_{\nu}$ in Definition \ref{df:UnuVnu}.

\begin{Lemma}
The family of solutions $W^\nu$ and $E_2^\nu$ of \eqref{eq:fond1}, with the above coefficients of $\bM_\nu$ defined in \eqref{eq:matm} can be written as
\begin{align*}
W^\nu(x)&=(-c_{\nu}(x))^{\sfrac12}(AU^0_{\nu}(x-X_\nu)+BV^0_{\nu}(x-X_\nu)),\\
E_2^\nu(x)&=A\Bigl[c_\nu(x)^{-1}\bigl(((-c_{\nu}(x))^{\sfrac12}U^0_{\nu}(x-X_\nu))'+ a_\nu(x)(-c_{\nu}(x))^{\sfrac12}U^0_{\nu}(x-X_\nu) \bigl) \Bigr]
\cr
&\hspace*{.5cm}+\ B \Bigl[c_\nu(x)^{-1} \bigl(((-c_{\nu}(x))^{\sfrac12}V^0_{\nu}(x-X_\nu))'+ a_\nu(x)(-c_{\nu}(x))^{\sfrac12}V^0_{\nu}(x-X_\nu) \bigl) \Bigr],
\end{align*}
where $A$ and $B$ are two given arbitrary constants.
The solutions are non singular.
\end{Lemma}

\begin{proof}
Recall that $U^0_{\nu}$ and $V^0_{\nu}$ solve
\[
\frac{d^2}{dx^2}U^0_{\nu}=\Bigl(-\Bigl(\frac{\rho'_\nu}{2\rho_\nu} \Bigr)^2-\frac{(\rho'_\nu)^2\lambda^2}{4\rho_{\nu}}+s_{\nu} \Bigr)U^0_{\nu}
\]
with $\rho_{\nu}:=\rho$ and $s_{\nu}:=0$. One checks that $(\rho_\nu)'(x)= \sfrac{1}{(1+\ell_0^{-1}(x-X_{\nu}))^2}$, hence $(\rho_{\nu}'(x))^{-\sfrac12}$ is a polynomial of degree 1, hence $s_{\nu}=0$ as announced.

Consider
\[
W^\nu(x)=(-c_{\nu}(x))^{\sfrac12}(AU^0_{\nu}(x-X_\nu)+BV^0_{\nu}(x-X_\nu)).
\]
$E_2^\nu$ is obtained through
\[
(W^\nu)'=c_\nu E_2^\nu-a_\nu W^\nu= c_\nu E_2^\nu+c_\nu W^\nu.
\]
As one has
\[
\Bigl(\frac{1}{c_\nu}(-c_\nu)^{\sfrac12} \Bigr)'+\frac{1}{c_\nu} \Bigl((-c_\nu)^{\sfrac12} \Bigr)'=0,
\]
one gets
\begin{align*}
(E_2^\nu)'-a_\nu E_2^{\nu}&= (E_2^\nu)'+c_\nu E_2^\nu\\
&= \Bigl(\frac{1}{c_\nu}((-c_\nu)^{\sfrac12})''\Bigr)(AU^0_\nu+BV^0_\nu)+\frac{1}{c_\nu}(-c_\nu)^{\sfrac12}(AU^0_\nu+BV^0_\nu)''\cr
&= \Bigl(\frac{1}{c_\nu}((-c_\nu)^{\sfrac12})''+\frac{1}{c_\nu}(-c_\nu)^{\sfrac12} \Bigl(-\frac14c_\nu^2-\frac{\lambda^2\rho_\nu}{4}c_\nu^2 \Bigl)\Bigl)(AU^0_\nu+BV^0_\nu).
\end{align*}
It is enough to call $b_\nu$ the coefficient such that $(E_2^{\nu})'-a_{\nu}E_2^{\nu}=b_{\nu}W^\nu$ using
\[
(AU^0_{\nu}+BV^0_{\nu})''=\Bigl(-\frac{(\rho'_\nu)^2}{4\rho_\nu^2}-\frac{\lambda^2(\rho'_\nu)^2}{4\rho_\nu} \Bigr)(AU^0_{\nu}+BV^0_{\nu}).
\]
The solutions are non singular since the $\log$ contribution in $V_\nu^0$ is counterbalanced
by the different square roots in front of it.
\end{proof}

\subsection{A second family with singularity}

The second case corresponds to a solution. It writes as a combination of Bessel functions
with $\rho(x)= x(1+\ell_0^{-2}x^2)$ independent of $\nu$. One may choose in this case
\[
c_{\nu}(x)=K\;\frac{\rho'(x-X_\nu)}{\rho(x-X_\nu)}\quad\text{and}\quad b_{\nu}(x)=c_{\nu}(x) \Bigl(\frac{a_{\nu}(x)}{c_\nu(x)} \Bigr)^2-\Bigl(\frac{a_{\nu}(x)}{c_\nu(x)} \Bigr)'-\frac{1}{K}\rho\rho'(x-X_\nu).
\]
It yields
\[
c_\nu(x)=K \Bigl[\frac{1}{x-X_\nu}+\frac{2(x-X_\nu)}{\ell_0^2+(x-X_\nu)^2} \Bigr],\quad
a_\nu(x)=C_0c_\nu(x)
\]
and
\begin{align*}
b_\nu(x)&=-C_0^2c_\nu(x)-\frac{1}{K}\rho'(x-X_\nu)\rho(x-X_\nu)\\
&=-C_0^2c_\nu(x)-\frac{1}{K}(x-X_\nu)(1+\ell_0^{-2}(x-X_\nu)^2)(1+3\ell_0^{-2}(x-X_\nu)^2).
\end{align*}

With a convenient
choice of $X_\nu\not \in \mathbb R$, one has

\begin{Lemma}
The solutions of
$$\begin{pmatrix}E_2^{\nu}\\[7pt]
W^\nu\end{pmatrix}'=\begin{pmatrix}C_0c_\nu&-C_0^2c_{\nu}-\dfrac{1}{K}\rho\rho'(x-X_\nu)\\[7pt]
c_{\nu}(x)&-C_0c_{\nu}(x)\end{pmatrix}\begin{pmatrix}E_2^{\nu}\\[7pt]
W^\nu\end{pmatrix},
$$
with
\[
c_{\nu}(x)=K\;\frac{\rho'}{\rho}(x-X_\nu)\quad\text{and}\quad\rho(z)=z(1+\ell_0^{-2}z^2),
\]
are
\begin{align*}
W^\nu(x)&=AJ_0(\rho(x-X_\nu))+BY_0(\rho(x-X_\nu)), \\
E_2^\nu(x)&=K^{-1}\rho(x-X_\nu)(AJ'_0(\rho((x-X_\nu)) +BY'_0(\rho((x-X_\nu)))\\
&\hspace*{5.45cm}+C_0(AJ_0(\rho(x))+BY_0(\rho(x))),
\end{align*}
where $A$ and $B$ are two arbitrary positive constants.
The solutions are singular for $B\neq 0$.
\end{Lemma}

\begin{proof}
The proof of this Lemma is constructive, namely we will deduce from the fact that $W^\nu= AJ_0(\rho (x-X_\nu))+BY_0(\rho
(x-X_\nu))$. We drop all the indices ${\nu}$.

The two equations are (using $a=C_0c$, where $C_0$ is a given positive constant) $W'=cE_2-C_0cW$, $E_2'=C_0cE_2+bW$. We replace $E_2$ by $\sfrac{W'}{c}+C_0W$ in the second equation, to obtain
$(\sfrac{W'}{c})'=(C_0^2c+b)W$.
Note that $W=J_0(\rho)$ is a solution of this equation, hence
$$
(C_0^2c+b)J_0(\rho)= \Bigl(\frac{\rho'}{c}J'_0(\rho) \Bigr)'= \Bigl(\frac{\rho'}{c} \Bigr)'J'_0(\rho)+\frac{(\rho')^2}{c}\;J''_0(\rho).
$$
Our aim is to remove all the terms containing $J'_0(\rho)$, using in particular
\[
J''_0(\rho)=-\frac{1}{\rho}\;J'_0(\rho)-J_0(\rho).
\]
Hence one deduces
$$(C_0^2c+b)J_0(\rho)= \Bigl(\Bigl(\frac{\rho'}{c} \Bigr)'-\frac{(\rho')^2}{\rho c} \Bigr)J'_0(\rho)-\frac{(\rho')^2}{c}\;J_0(\rho).$$
The equation that must be satisfied by $c$ in order to have $J_0(\rho)$ as a solution is $(\sfrac{\rho'}{c})'=\sfrac{(\rho')^2}{\rho c}$, which is equivalent to $\sfrac{\rho'}{\rho c}$ is constant, namely
$
c=K\sfrac{\rho'}{\rho}
$.
The coefficient $b$ can thus be identified thanks to
\[
b=-C_0^2c-\frac{(\rho')^2}{c}=-C_0^2K\;\frac{\rho'}{\rho}-K^{-1}\rho\rho',
\]
which is the given coefficient $b$ in the matrix $\bM_{\nu}$.
The solutions are singular because the $\log$ term in $Y_0$ is nude.
This ends the
proof.
\end{proof}

\section{Multi species}

Our aim here is to show that our method can be used
to study the resonances in the case of multi-species.

For simplicity, let us consider the case of the electron-ion system.
After elementary calculations based on the fundamental equations
that one may find in \cite{branbilla} for the cold plasma model (from which the author deduces expressions (17.6) and (17.7) for the effective dielectric tensor), one obtains the generalization
of
\eqref{eq:matm}
\begin{equation} \label{eq:AZ12}
\bM_\nu(x,\theta)=
\frac{1}{D_\nu(x)}
\begin{pmatrix}
A_\nu(x,\theta)&B_\nu(x,\theta) \cr
C_\nu (x)&-A_\nu(x,\theta)
\end{pmatrix}.
\end{equation}
Here $D_\nu$ satisfies
\begin{itemize}
\item $D_0(x)$ vanishes at points $x_{hi}$ such that $A_0(x_{hi}), B_0(x_{hi}), C_0(x_{h_i})$ do not vanish and $(A_0^2+B_0C_0)(x_{hi})=0$
\item $i\partial_{\nu}D_\nu(x_{hi})\vert_{\nu=0}\in \mathbb R \not=0$
\end{itemize}
which are the unique conditions needed for the theoretical set-up of our main result.
This coefficient $D_\nu$ is analogous to $\epsilon_{11}^\nu$.

Note that the method described above does not need expressions (17.6) and (17.7) of \cite{branbilla}, these expressions classically derive the effective dielectric tensor as the sum of the effective dielectric tensors of each species, hence introducing artificial poles in the dielectric tensor at all the cyclotron resonances. The formulation we introduce here avoids introducing these poles.

We describe a two species case, one is electrons and one is ions with $Z$ protons and just give the main ideas.
To obtain \eqref{eq:AZ12}, one can start from the
Maxwell equations
$$\nabla \wedge E=i\omega B,
\quad
c^2\nabla \wedge B=-i\omega E+\epsilon_0^{-1} j.$$
Thanks to the electro-neutrality, the ions density is $Z^{-1}N_e(x)$, $N_e(x)$ being the electrons density, hence the current $j$ is given by
\[
j=-eN_e(x)v_e+eZ(Z^{-1}N_e(x)) v,
\]
$v$ being the ion velocity, $v_e$ being the electron velocity, hence one has
\[
j=eN_e(x)(v-v_e).
\]
The external magnetic field is $\bB_0= B_0(x) b$, $b$ is a unit vector.
The so-called fluid equations are given through Newton's law
$$m_e(-i\omega +\nu)v_e=- e(E+B_0(x)v_e\wedge b),
\quad
m_i(-i\omega +\nu)v= Ze(E+B_0(x)v\wedge b)
$$
from which one deduces the equation for the electric current thanks to $v_e=v-\frac{1}{eN_e(x)}j$:
$$
(-i\omega +\nu)j=\epsilon_0\omega_p^2(1+p)(E+B_0(x)v\wedge b)-\frac{eB_0(x)}{m_e}j\wedge b,
\qquad p=\frac{Zm_e}{m_i}\in (0,1).
$$
Introduce the ion cyclotron frequency and the electron cyclotron frequency
\[
\omega_c^i(x)=\frac{ZeB_0(x)}{m_i},\quad \omega_c^e(x)=\frac{eB_0(x)}{m_e}.
\]
The traditional approach is to calculate $v$ in terms of $E$ (which introduces a singularity at $x$ such that $\omega_c^i(x)=\omega$), then to deduce $j$ (or $v_e$) in terms of $E$ (which introduces a singularity at $x$ such that $\omega_c^e(x)=\omega$) 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 $\nu$ goes to $0$. Instead of using this approach, we seek the TE field $(E_2(x), B_3(x):=W(x))e^{i\theta y}$. Let $(j_1,j_2)$ be the components of the electric current transverse to $b$. The TE field being solution of
$$
\begin{cases}
c^2W'=i\omega E_2-\epsilon_0^{-1} j_2,\\
E_2'=i\theta E_1+i\omega W,
\end{cases}
$$
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
$$
\begin{cases}
(-i\omega +\nu)j_2-\omega_c^ej_1+\epsilon_0\omega_p^2(1+p)B_0v_1=\epsilon_0\omega_p^2(1+p)E_2,\cr
\omega_c^ej_2+(-i\omega+\nu)j_1-\epsilon_0\omega_p^2(1+p)B_0v_2-\epsilon_0\omega_p^2(1+p)E_1=0,\\[5pt]
(-i\omega+\nu)v_1-\omega_c^iv_2-\dfrac{Ze}{m_i}E_1=0,\\[7pt]
\omega_c^iv_1+(-i\omega+\nu)v_2=\dfrac{Ze}{m_i}E_2,\cr
\epsilon_0^{-1} j_1-i\omega E_1=c^2i\theta W,
\end{cases}
$$
where $(v_1,v_2)$ are the components of the ion velocity transverse to $b$.
Treating this system globally leads to a determinant $D_\nu$:
$$D_{\nu}(x)=-i\omega ((\omega_c^i)^2+(i\omega-\nu)^2)((\omega_c^e)^2+(i\omega-\nu)^2)+(-i\omega+\nu)((i\omega-\nu)^2+\omega_c^i\omega_c^e)\omega_p^2(1+p).$$
For the purpose of finding the hybrid mode associated with this system, one checks that
$$
-\frac{D_0(x)}{i\omega}=(p^2(\omega_c^e)^2-\omega^2)((\omega_c^e)^2-\omega^2)+(1+p)\omega_p^2(p(\omega_c^e)^2-\omega^2).
$$
Note that: for $\omega=\omega_c^i$, $-\sfrac{D_0(x)}{i\omega}$ is equal to $(1+p)\omega_p^2\omega_c^i(\omega_c^e-\omega_c^i)>0$; while for $\omega=\omega_c^e$, $-\sfrac{D_0(x)}{i\omega}$ is equal to $(1+p)\omega_p^2\omega_c^e(\omega_c^i-\omega_c^e)<0$, which has
an opposite sign.
If one denotes $X=(\sfrac{\omega}{\omega_c^e})^2$, we note that
$$-\frac{D_0(x)}{i\omega (\omega_c^e)^4}=(p^2-X)(1-X)+(1+p)\,\frac{(\omega_p(x))^2}{(\omega_c^e)^4}\,(p-X),$$
hence the hybrid resonances, in this case, correspond to $X$ solution of
$$
X^2-\Bigl(p^2+1+(p+1)\frac{\omega_p^2}{(\omega_c^e)^2} \Bigl)X+p \Bigl(p+(1+p)\frac{\omega_p^2}{(\omega_c^e)^2} \Bigr)=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 $\omega>0$, this gives rise to two hybrid resonances $X_-^h$ and $X_+^h$ such that $\omega_c^e(X_-^h)=\omega X_-^{\sfrac12}$ and $\omega_c^e(X_+^h)=\omega X_+^{\sfrac12}$. They give rise to the behavior described in this paper.
Other formulas can be sought for.

\backmatter
\bibliographystyle{jepplain}
\bibliography{despres-et-al}
\end{document}