<html><head><meta name="color-scheme" content="light dark"></head><body><pre style="word-wrap: break-word; white-space: pre-wrap;">\makeatletter
\@ifundefined{HCode}
{
\documentclass[screen,CRPHYS,Unicode,biblatex,published]{cedram}
\addbibresource{crphys20250054.bib}
\newenvironment{noXML}{}{}
\newenvironment{Table}{\begin{table}}{\end{table}}
\def\thead{\noalign{\relax}\hline}
\def\endthead{\noalign{\relax}\hline}
\def\tabnote#1{\vskip4pt\parbox{.8\linewidth}{#1}}
\def\tsup#1{$^{{#1}}$}
\def\tsub#1{$_{{#1}}$}
\def\sfrac#1#2{{#1}/{#2}}
\def\sbfrac#1#2{({#1})/({#2})}
\def\sofrac#1#2{({#1}/{#2})}
\def\stfrac#1#2{({#1})/{#2}}
\def\stofrac#1#2{(({#1})/{#2})}
\usepackage{amssymb}
\RequirePackage{etoolbox}
\usepackage{amsmath}
\csdef{Seqnsplit}{\\}
\def\botline{\\\hline}
\def\tbody{\noalign{\relax}\hline}
\def\jobid{crphys20250054}
\def\rmpi{\uppi}
\allowdisplaybreaks
%\graphicspath{{/tmp/\jobid_figs/web/}}
\graphicspath{{./figures/}}
\newcounter{runlevel}
\let\MakeYrStrItalic\relax
\def\refinput#1{}
\def\0{\phantom{0}}
\def\back#1{}
\DOI{10.5802/crphys.238}
\datereceived{2024-09-13}
\daterevised{2024-11-15}
\datererevised{2025-01-20}
\dateaccepted{2025-02-04}
\ItHasTeXPublished
\skip\footins 15pt
\def\myRef{\ref}
\def\ndash{\text{--}}
\def\mathbi#1{\text{\textit{\textbf{#1}}}}
}
{\documentclass[crphys]{article}
\def\CDRdoi{10.5802/crphys.238}
\let\newline\break
\def\sfrac#1#2{{#1}/{#2}}
\def\sofrac#1#2{({#1}/{#2})}
\makeatletter
\def\selectlanguage#1{} 
\def\href#1#2{\url[#1]{#2}}
\def\sfrac#1#2{{#1}/{#2}}
\def\sbfrac#1#2{({#1})/({#2})}
\def\sofrac#1#2{({#1}/{#2})}
\def\stfrac#1#2{({#1})/{#2}}
\def\stofrac#1#2{(({#1})/{#2})}
}
\makeatother

\usepackage{upgreek}

\dateposted{2025-03-18}
\begin{document}

\begin{noXML}

\CDRsetmeta{articletype}{research-article}

\title{Ray tracing of very low frequency waves produced by active
experiments or lightning events at low Earth orbit}

\alttitle{Lancer de rayons d'ondes \`{a} tr\`{e}s basse fr\'{e}quence
produites par des exp\'{e}riences actives ou des \'{e}clairs \`{a}
orbite terrestre basse}

\author{\firstname{Leo} \lastname{Cerfolli}\CDRorcid{0009-0005-2166-7706}}
\address{CEA, DAM, DIF, F-91297, Arpajon, France}
\address{UPS, CEA, LMCE, 91297 Bruy\`eres-le-Ch\^atel, France}
\email[L. Cerfolli]{leo.cerfolli@cea.fr}

\author{\firstname{Jean-Francois} \lastname{Ripoll}\CDRorcid{0000-0002-1177-519X}\IsCorresp}
\addressSameAs{1}{CEA, DAM, DIF, F-91297, Arpajon, France}
\addressSameAs{2}{UPS, CEA, LMCE, 91297 Bruy\`eres-le-Ch\^atel, France}
\email[J.-F. Ripoll]{jean-francois.ripoll@cea.fr}

\author{\firstname{Richard Bertram} \lastname{Horne}\CDRorcid{0000-0002-0412-6407}}
\address{British Antarctic Survey, Natural Environment Research Council, Cambridge, UK}
\email[R. B. Horne]{rh@bas.ac.uk}

\author{\firstname{Alexandra Michelle}\nobreakauthor\lastname{Wold}\CDRorcid{0009-0002-1298-2390}}
\address{Ann \&amp; H.J. Smead Department of Aerospace Engineering Sciences,
University of Colorado Boulder, CO, USA}
\email[A. M. Wold]{alexandra.wold@colorado.edu}

\author{\firstname{Robert Andrew} \lastname{Marshall}\CDRorcid{0000-0003-0974-1609}}
\addressSameAs{4}{Ann \&amp; H.J. Smead Department of Aerospace Engineering Sciences,
University of Colorado Boulder, CO, USA}
\email[R. A. Marshall]{robert.marshall@colorado.edu}

\keywords{\kwd{Ray tracing}
\kwd{Very low frequency waves}
\kwd{Magnetic field models}
\kwd{Low Earth orbit}}

\altkeywords{\kwd{Ray-tracing}
\kwd{Ondes tr\`{e}s basses fr\'{e}quences}
\kwd{Mod\`{e}les de champs g\'{e}omagn\'{e}tiques}
\kwd{Basses orbites}}

\begin{abstract}
We investigate the propagation in the plasmasphere of Very Low
Frequency (VLF) electromagnetic waves, such as natural
lightning-generated whistler waves and waves produced by active
experiments. An active experiment is an artificial controlled
disturbance of the low orbit space or the ionosphere. The aim is often
to produce electromagnetic waves for removing high-energy particles
(mostly electrons). We study the wave propagation parameters whether
they are geometric, background, or intrinsic parameters, such as the
magnetic field model, the ambient plasma density model of the
plasmasphere, and the wave frequency. All of these parameters cause
different behaviors of propagation, which are discussed in this
article.
\end{abstract}

\begin{altabstract}
Nous explorons la propagation d'ondes \'{e}lectromagn\'{e}tiques de
tr\`{e}s basses fr\'{e}quences dans la plasmasph\`{e}re, telle que des
ondes de type siffleur g\'{e}n\'{e}r\'{e}es par les \'{e}clairs ou des
ondes produites par des exp\'{e}riences actives. Une exp\'{e}rience
active est une perturbation artificielle et contr\^{o}l\'{e}e de
l'espace proche de la Terre ou de l'ionosph\`{e}re souvent dans
l'id\'{e}e de d\'{e}peupler cette r\'{e}gion de particules de hautes
\'{e}nergies (souvent des \'{e}lectrons).
Pour le faire, nous ordonnons les param\`{e}tres de propagation selon
trois classes : les param\`{e}tres de fond, les param\`{e}tres
g\'{e}om\'{e}triques et les param\`{e}tres intrins\`{e}ques. Les
param\`{e}tres de fond regroupent le mod\`{e}le de champ magn\'{e}tique
ainsi que celui de la densit\'{e} plasmasph\'{e}rique. Ils influent via
leurs valeurs aux diff\'{e}rents points de calcul ainsi que leurs
gradients. Le param\`{e}tre intrins\`{e}que ici \'{e}tabli est la
fr\'{e}quence de l'onde. Avec ces param\`{e}tres nous verrons que des
comportements tr\`{e}s diff\'{e}rents sont accessibles.
\end{altabstract}

\thanks{Direction G\'{e}n\'{e}rale de l'Armement (DGA), Agence pour
l'Innovation de D\'{e}fense (AID)}

\shortrunauthors

\maketitle

\end{noXML}

{\vspace*{-6pt}}

\section{Introduction}\label{sec1}

The plasmasphere, named by~\cite{1}, is the toroidal region of cold
electrons (1--5 eV) located between the ionosphere from ${\sim}$90~km of
altitude up to ${\sim}$1000~km~\cite{2} and its external limit, the
plasmapause. The plasmapause is generally defined as the location where
the electron density decreases by a factor~5 on a radial distance of
${\sim}$3000~km~\cite{3}. The plasmapause generally locates at the 
100~\#/cc density level~\cite{4}. Many types of electromagnetic waves
propagate within the plasmasphere and we focus in this article only on
very low frequency (VLF) waves as natural lightning-generated whistler
(LGW) waves as well as waves produced by active experiments. An active
experiment is an artificial controlled disturbance in the low orbit
space or in the ionosphere. The aim is often to produce waves for
removing high-energy particles (mostly electrons) forming radiation
belts~\cite{5}, which are source of hard radiations for orbiting
satellites. As instances of active experiments, the most noticeable are
Barium releases as CRRES-1991~\cite{6} and more recently the
amplification of local waves by the rocket exhausts of hydrazine 
burn~\cite{7}. Here, we study the propagation of VLF waves through a
parametric numerical study involving the HOTRAY  code~\cite{8}, varying
various parameters such as the magnetic field model, the ambient plasma
density model of the plasmasphere, and the wave frequency. 

In Section~\ref{sec2}, we briefly describe the numerical method and
code for wave propagation computation. Section~\ref{sec3} shows the importance
of the magnetic field model through model comparisons. Section~\ref{sec4} is
dedicated to the description and role of the diffusive equilibrium
model driving the electron density in the Earth's plasmasphere.
Section~\ref{sec5} shows the influence of the wave frequency for wave
propagation. Conclusions are gathered in the last section.

\section{Simulation of wave propagation}\label{sec2}

We use the HOTRAY code~\cite{8} to simulate VLF waves. The HOTRAY code
uses the eikonal approximation of geometric optics to simulate the
curved propagation and find the ray  path (\cite{8} and references in
there). Based on Eikonal approximation, the main resolved equations
are given by: 
{\begin{eqnarray}
\frac{\partial \vec {R}}{\partial t}&amp;=&amp;-\frac{\partial D/\partial \vec {k}}{\partial D/\partial \omega }=\vec {v}_{g}\label{eq1} \Seqnsplit
\frac{\partial \vec {k}}{\partial t}&amp;=&amp;\frac{\partial D/\partial \vec {R}}{\partial D/\partial \omega }\label{eq2}
\end{eqnarray}}\unskip

Equation~(\ref{eq1}) solves the propagation in space, where $R$ is the
position vector, $t$ the time, $D$ the dispersion term and $k$ the complex
wave vector. HOTRAY assumes that the ambient plasma is non-relativistic
and that wave amplitudes are small for applying both the linear theory
and linear instability growth and damping. The code works with the
underlying WKB (Wentzel--Kramers--Brillouin), $|1/k (\mathrm{d} k/\mathrm{d} x) |
\ll k$, which means that the change of wavelength over a wavelength
should be small~\cite{9}. This approximation fails when $k \rightarrow
0$ or when $\mathrm{d} k/\mathrm{d} x \rightarrow \infty$, or whenever the wave
approaches either a cutoff or a resonance. This implies the spatial
gradients in the ambient cold plasma are assumed to be small compared
with the wavelength. From a numerical point of view, the WKB
approximation is generally satisfied, as we define the inhomogeneous
plasma by a succession of straight propagation into layers of
homogeneous plasma where density and magnetic field are constant (see
below the case of magneto-reflection).

In HOTRAY, the chosen strategy is to launch given frequency waves so
the frequency of the wave, $f_{0}=\sfrac{\omega _{0}}{2\uppi }$, or its
pulsation $\omega _{0}$, is always fixed. For computing the
amplification/damping, since we fix the pulsation real, the complex
contribution given by the relation dispersion comes from the wave
vector (see exception below for magnetoreflection). In
Equation~(\ref{eq1}) the right-hand side can be identified as the wave
group velocity. Equation~(\ref{eq2}) solves for the evolution of the
wave vector according to the dispersion relation. In all points along
the ray path,~we have to solve the full dispersion relation, given by
{\begin{equation}\label{eq3}
D(\vec {R}, \vec {k}, \omega)=0.
\end{equation}}\unskip
The complete expression of the dispersion relation is given
in~\cite{8}. Equation~(\ref{eq3}) can be written as
$An^{4}+Bn^{2}+C=0$.

Furthermore, the numerator in the right-hand side of
Equation~(\ref{eq2}) can be developed as
{\begin{equation}\label{eq4}
\frac{\partial D}{\partial \vec {R}}=\frac{\partial D}
{\partial \vec {B}}\frac{\partial \vec {B}}{\partial \vec {R}}+
\frac{\partial D}{\partial N}\frac{\partial N}{\partial \vec {R}}+
\frac{\partial D}{\partial \vec {k}}\frac{\partial \vec {k}}
{\partial \vec {R}}
\end{equation}}\unskip

The three terms in the right-hand side of Equation~(\ref{eq4})
highlight the competition between magnetic field gradients, density
gradients, and variations of the wave vector. In the third term,
$\mathrm{d} k/\mathrm{d} R$, the variations in position of the wave
vector account for changes in amplitude and the orientation and norm of
the wave vector. 

The fact that the wave vector is complex implies the wave amplitude can
be damped or amplified. The wave vector direction is defined by two
angles: the wave normal angle defined as the angle between the wave
vector and the local magnetic field and the azimuthal angle defined as
the geographic eastward component of the wave vector. 

Background parameters regrouping the magnetic field and the density
model appear as gradients in Equation~(\ref{eq4}) as well as single-point
values in the dispersion relation. An important condition for wave
propagation is dictated by the lower hybrid (LH) frequency. The LH
frequency is the limit for allowing the bounce of the ray at high
latitude (see more below). Its definition is given by: 
{\begin{equation}\label{eq5}
\omega _{\mathrm{LH}}=\left[\left(\Omega _{i}\Omega _{e}\right)^{-1}
+{\omega }_{pi}^{-2}\right]^{-1/2}
\end{equation}}\unskip
where $\Omega _{i}=\sfrac{q_{i}B}{m_{i}}$ is the cyclotron pulsation,
$q_{i}$ is the charge, $m_{i}$ is the mass, $e$ and $i$ indices refer to
electrons and ions respectively, and $\omega _{p}$ is the plasma pulsation.
$\omega _{\mathrm{LH}}$ The lower hybrid frequency varies as B,  through both
the cyclotron pulsation terms, and as the square root of the electron
density. An increase in $B$ or in $\sqrt n_{e}$ increases the lower
hybrid frequency. This shows another competition between the magnetic
field and the plasma density for each single location of the wave
during its propagation. 

Landau resonant interactions are included in the code and occur when
the electron speed is near the phase velocity of the wave. This creates
an energetic flux between the wave and the electron. If the wave gives
energy to the electron, then the wave is damped. Otherwise, the wave is
amplified. The spatial scale for the Landau instability growth rate is
dictated by each cell size and the local growth rate, $\Gamma _{i}$,
is driven by the imaginary part of the local wave vector $k_{i}$ in the
layer number $i$ by $\Gamma
_{i}(\mathrm{k})=\mathfrak{I}(k_{i})\Delta r_{i}$, where
$\mathfrak{I}$ is the imaginary operator and $r_{i}$ the spatial step.
Expanding the dispersion relation in Taylor series for  $\sfrac{\gamma
}{\omega }\ll 1$, \cite{14} show that
$\mathfrak{I}(k_{i})=-\sbfrac{\gamma
_{i}}{\sfrac{\sofrac{\partial D}{\partial k}}{\sofrac{\partial
D}{\partial \omega }}_{i}}$ where the denominator can be
identified as the group velocity in cell $i$. Near a magnetoreflection,
the parallel group velocity changes sign, passing to zero, while the
perpendicular group velocity, thus the group velocity, does not go to
zero---otherwise the wave would not continue to propagate.

Near zero parallel group velocity can bring numerical errors because
the spatial step becomes numerically too small. In order to avoid these
errors, the growth rate, $\Gamma _{i}({k})$, is rather
computed in the time domain. In that case, the numerical method is to
assume the pulsation is not real anymore but switch to a complex number
written as $\omega _{i}=\omega _{0}+j\gamma _{i}$, with $j$ the imaginary
unit, $\gamma _{i}$~the imaginary part of the pulsation leading to a
local growth rate, $\Gamma _{i}=\gamma _{i}\Delta t_{i}$, and
$\omega _{0}=2\uppi f_{0}$~\cite{8}. The complex pulsation then serves to
compute the damping/growth of the wave amplitude but the wave frequency
is still maintained constant at $f_{0}$. The reflection point is
neither a cut-off where $k$ tends to zero nor a resonance where $k$ becomes
very large. $k$ rather remains finite. However, near magneto-reflection
(with the parallel group velocity passing to zero) the WKB modeling
fails in principle but still gives an approximation accurate enough to
describe the wave propagation at the turning point for the propagation
of a low frequency wave in the Earth's magnetosphere (e.g.\ \cite{10}).

There are three available thermal ions causing Landau resonance, often
physically leading to Landau damping. They correspond to the common
species in the plasmasphere, namely H$^+$, He$^+$, O$^+$. Their density,
$N_{b}$, and proportions are adjusted following \cite{11} such that
their presence adds a warm plasma correction to the cold electron
background. The proportions used here is given in Table~\ref{tab1} and are the
same as in \cite{11}. Distribution functions are given by
{\begin{equation*}
f_{i}=\frac{1}{\uppi ^{3/2}{\alpha }_{\bot i}^{2}\alpha _{\parallel i}}
\mathrm{e}^{-\sofrac{(v_{z}-v_{di})^{2}}{{\alpha }_{\parallel i}^{2}}}
\left[\Delta _{i}\mathrm{e}^{-\sofrac{{v}_{\bot }^{2}}{{\alpha }_{\bot i}^{2}}}+
\frac{1-\Delta _{i}}{1-\beta _{i}}\left(\mathrm{e}^{-\sofrac{{v}_{\bot }^{2}}
{{\alpha }_{\bot i}^{2}}}-\mathrm{e}^{-\sofrac{{v}_{\bot }^{2}}{\beta _{i}
{\alpha }_{\bot i}^{2}}}\right)\right]
\end{equation*}}\unskip
where $i$ refers to the specie, $\Delta $ determines the depth of the
loss cone, $v_{d}$ is the drift velocity along the magnetic field
line, $\alpha _{\bot}$ or $\alpha _{\parallel }$ is the thermal
velocity perpendicular/parallel to the magnetic field given in keV
using $\sfrac{1}{2}m_{i}\alpha ^{2}$, and $\beta $ measures the
width of the loss cone. Here, we assume there is no loss cone
($\Delta =1$) and the drift velocity along the magnetic field line
is 0. Then the distribution function reduces as
{\begin{equation*}
f_{i}=\frac{1}{\uppi ^{3/2}{\alpha }_{\bot i}^{2}\alpha _{\parallel i}}
\mathrm{e}^{-\sofrac{{v}_{z}^{2}}{{\alpha }_{\parallel i}^{2}}}
\mathrm{e}^{-\sofrac{{v}_{\bot }^{2}}{{\alpha }_{\bot i}^{2}}}
\end{equation*}}\unskip
In our simulations, we choose both the perpendicular and parallel thermal
velocities so that they are equals.

\begin{table}[t!]
\caption{\label{tab1}Species and plasma parameters used in this study.
The drift velocity is always zero.}
\begin{tabular}{ccccc}
\thead
\parbox[t]{1.8cm}{\centering Species index} &amp;  Species &amp; 
\parbox[t]{3.8cm}{\centering Thermal velocity  perpendicular to B  (keV)} 
&amp; \parbox[t]{2.8cm}{\centering Parallel thermal velocity (keV)} &amp; 
Density (\#/m$^{-3}$)\vspace*{2pt} \\
\endthead
1 &amp;O$^{+}$ &amp; $10^{-4}$ &amp; $10^{-4}$ &amp; $10^{9}$ \raisebox{.35cm}{}\\
2 &amp; $\mathrm{e}^{-}$ &amp; $10^{-3}$ &amp; $10^{-3}$ &amp; $10^{10}$ \\
3 &amp; He$^{+}$ &amp; $10^{-4}$ &amp; $10^{-4}$ &amp; $10^{9}$ \\
4 &amp; H$^{+}$ &amp; $10^{-4}$ &amp; $10^{-4}$ &amp; $8\times 10^{9}$ \\
5 &amp; H$^{+}$ &amp; $10^{-4}$ &amp; $10^{-4}$ &amp; $2.782\times 10^{7}$ \\
6 &amp; $\mathrm{e}^{-}$ &amp; $2.5\times 10^{-2}$ &amp; $2.5\times 10^{-2}$ &amp; $2.5\times 10^{7}$ \\
7 &amp; $\mathrm{e}^{-}$ &amp; $1.5\times 10^{-1}$ &amp; $1.5\times 10^{-1}$ &amp; $2.5\times 10^{6}$ \\
8 &amp; $\mathrm{e}^{-}$ &amp; $1$ &amp; $1$ &amp; $2.5\times 10^{5}$ \\
9 &amp; $\mathrm{e}^{-}$ &amp; $5$ &amp; $5$ &amp; $5\times 10^{4}$ \\
10 &amp; $\mathrm{e}^{-}$ &amp; $20$ &amp; $20$ &amp; $2\times 10^{4}$ 
\botline
\end{tabular}
\end{table}


\section{Influence of the magnetic field model}\label{sec3}

Active experiments are generally performed at low Earth orbit (from
100~km to 1500~km) and produce electromagnetic waves in a large range
of L-shells according to the experiment's altitude and latitude. The
L-shell is the Earth's radius normalized distance of a given field line
at the magnetic equator. The L-shell discriminates the field lines and
depends on the chosen magnetic field model. At low Earth orbit (LEO),
there are three main models for Earth's magnetic field: the pure dipole
used by many codes for simplicity, as HOTRAY, the eccentric tilted
dipole, which is the most accurate dipolar approach of the Earth's
magnetic field, and the full International Geomagnetic Reference Field
(IGRF) model. Dipolar models have orthogonal metrics in classic space
coordinates (like geographic or magnetic coordinates) which provide
mathematical simplifications. The eccentric tilted dipole is a pure
dipole for which the main axis is not the Earth's rotation axis, since
it is tilted, and the magnetic center does not correspond to the
geographic center, since it is shifted. The IGRF is defined by a
potential as a harmonic development around a tilted dipole field and is
based on a potential defined from ground measurements. The potential
expression is:
{\begin{equation}
\label{eq6}
V(r,\theta,\varphi,t)=R_{e}{\sum }_{n=1}^{N}R^{n+1}
\left[{g}_{n}^{m}\cos (m\varphi)+{h}_{n}^{m}\sin 
(m\varphi)\right]{P}_{n}^{m}(\cos(\theta))
\end{equation}}\unskip
where coefficients noted $h$ and $g$ are established from ground
measurements, $R_{e}$ is the Earth's radius, $R=R_{e}/r$, $r$ is the
radial distance, $\theta $ is the latitude, $\varphi $ is the longitude
and ${P}_{n}^{m}$ are associated Legendre polynomials. 

The magnetic field is then derived from the potential as 
{\begin{equation}\label{eq7}
\vec {B}=-\vec {\nabla }(V)
\end{equation}}\unskip
We use the 13th generation of the IGRF coefficients~\cite{12}. In 
Figure~\ref{fig1}, we plot the magnetic amplitude with respect to
latitude and longitude at 100~km of altitude. This altitude is chosen,
since it is the commonly adopted limit for electron precipitations in
the upper atmosphere. Both the eccentric tilted dipole model and the
IGRF model have the same main structures. First, the ``S shape'' is the
global shape along longitude produced by the variation of the magnetic
equator along longitude. The second is the South Atlantic Anomaly
(SAA), which is a zone of large depression in magnetic amplitude. The
center of this depression is different between the eccentric tilted
dipole and the IGRF models. For the eccentric tilted dipole, the center
is at ${-}$30\textdegree~of latitude and ${-}$10\textdegree~of longitude,
while, for IGRF, it is centered at ${-}$40\textdegree~of latitude and
between ${-}$60\textdegree~and ${-}$30\textdegree~of longitude. This
depression is a consequence of both the tilted axis and the shift of
the magnetic center. The global shape of the SAA is also different for
the two models. In the IGRF model, high order harmonics ($n \geq 3$ in
Equation~(\ref{eq5})) constrain the SAA to fit with ground measurements.

\begin{figure}
\includegraphics{fig01}
\caption{\label{fig1}Magnetic field amplitude projected on Earth at
100~km for (left) a pure dipole field model as used in HOTRAY, (center)
an eccentric tilted dipole field model, and (right) the IGRF 13th
generation model.}
\end{figure}

In Figure~\ref{fig2}, we plot the relative error between the two
dipolar models and the IGRF. Differences for each comparison are
located around the IGRF SAA. The relative error reaches a maximum of
101\% for the dipole and 53\% for the eccentric tilted dipole.
Elsewhere, dipolar models are rather accurate. However, the relative
error does not show topologic differences seen on the field lines (not
discussed here). Field lines end up different for each model, with
implication on the propagation as field lines act as a wave guide. Note
that, for very low frequency waves, the plasmapause also acts as a wave
guide \cite{13}. The consequence is that ray paths can
change. In the following sections, the HOTRAY code will be used with
the dipole field in order to provide first results for future
comparisons with the IGRF model, currently being integrated in HOTRAY.

\begin{figure}
\includegraphics{fig02}
\caption{\label{fig2}Magnetic field amplitude relative error at 100~km,
(left) for a pure dipole field model compared with the IGRF model and
(right) for an eccentric tilted dipole filed model compared with the
IGRF model. Most of the errors are concentrated in the SAA region.}
\end{figure}

\section{Influence of the plasmaspheric density model}\label{sec4}

The HOTRAY code computes ray-path evolving in an ambient medium filled
of cold ions and electrons. Their density is given by the diffusive
equilibrium model, which equates the electron density with the ionic
densities $N_{e}=N_{i}$. The diffusive equilibrium further decomposes
$N_{i}$ in: 
{\begin{equation}\label{eq8}
N_{i}(R)=N_{b}N_{\mathrm{de}}(R){N_{li}}(R)N_{pl}(R)
\end{equation}}\unskip
where the density $N_{b}$ is the number of particles (\#) per meter
cube on a given point of reference at an altitude $R_{b}$. $N_{b}$ is a
scale density which scales the density profile, defined at all radii,
$R$, of the whole plasmasphere and will be studied below. The
dimensionless density  $N_{\mathrm{de}}$ is the profile of diffusive
equilibrium along the radial direction, which depends on temperature
$T_{\mathrm{DE}}$, and gravity, $g(R_{b})$, for each plasma component, as
follows
{\begin{equation}\label{eq9}
N_{\mathrm{de}}(R)^{2}={\sum }_{i=1}^{n}\eta _{i}\exp 
\left(-\frac{G}{H_{i}}\right)
\end{equation}}\unskip
with 
{\begin{equation}\label{eq10}
H_{i}=\frac{k_{\mathrm{B}}T_{\mathrm{DE}}}{M_{i}m_{p}g(R_{b})}
\end{equation}}\unskip

In Equations~(\ref{eq9})--(\ref{eq10}) giving the profile of the
diffuse equilibrium, each exponential term is weighted by the
proportion of the plasma component ions, ${\eta}_i$, which is a
function of the radial distance. There are three available species: H$^+$,
He$^+$, and O$^+$. This highlights the importance of the characteristic
length of thermal diffusion, $H_{i}$, prescribed for each component of
the ambient plasma. The associated thermal velocity is given by
$v_{\mathrm{th}}\approx \sqrt{} (3k_{\mathrm{B}}T_{\mathrm{DE}}/m)$. G is the
geopotential height, $k_{\mathrm{B}}$ is the Boltzmann constant, $M_{i}$ is the
mass number for the ion labeled $i$, and $m_{p}$ the proton mass. 

The dimensionless density $N_{{li}}$ connects with the bottom of
the ionosphere through a decreasing exponential function defined from
the chosen altitude of the bottom ionosphere $R_{0}$ (in km), up to its
height $H $ (in km) as follows
{\begin{equation}\label{eq11}
N_{li}=1-\exp \left[-(R-R_{0})^{2}/H^{2}\right]
\end{equation}}\unskip
The dimensionless density profiler $N_{pl}$ forms the plasmapause by
introducing a latitude dependence through the L-shell value, $L$, and
curves the density isolines with the plasmapause shape, at position
$L_{p} $ expressed in L-shell.
{\begin{equation}\label{eq12}
N_{pl}=\exp^{-\sfrac{(L-L_{p})^{2}}{{H}_{p}^{2}}}+
\left(1-\exp^{-\sfrac{(L-L_{p})^{2}}{{H}_{p}^{2}}}\right)
\left(\frac{R_{c}}{R}\right)^{a}+
\left(1-\left(\frac{R_{c}}{R}\right)^{a}\right)
\left(\exp^{-\sfrac{(R-R_{c})^{2}}{{H}_{s}^{2}}}\right)
\end{equation}}\unskip
$H_{p}$ is the half width of the plasmapause boundary in $L$ value.
$R_{c}$ is the geocentric distance in km to the level at which the
density outside the plasmapause field line is equal to the density
inside. $H_{s}$ is the scale height of radial density decrease.
Finally, $a$ determines how the density falls off outside the
plasmapause.

The following subsections are dedicated to the influence of the
diffusive equilibrium temperature, directly influencing the diffusion
of the electron density from Earth's atmosphere and the absolute value
of the electron density in HOTRAY simulations of VLF waves. 


\subsection{Influence of the diffusive equilibrium
temperature}\label{sec41}

For each component, the characteristic length of thermal diffusion is
proportional to the ratio of temperature with gravity so that, for a
given position, as gravity is fixed, the scale length becomes directly
proportional to temperature. In Figure~\ref{fig3}, we vary the
diffusive equilibrium temperature from $T=1000$ to 2000~K and show the
variation of the electron density. For these simulations, the scale
density is fixed at $n_{e}=2.7\times 10^{9}~\mathrm{m}^{-3}$ and the
plasmapause at $L_{p}=4.5$. The temperature radially diffuses density
at higher radial distance and follows the plasmapause. The diffusive
scale length increases with temperature so that, as the temperature
increases, more plasma components can reach higher altitudes driven by
their thermal velocity.

\begin{figure}
\includegraphics{fig03}
\caption{\label{fig3}
Variation of the diffusive equilibrium electron density with
temperature for (left) 1000~K, (center) 1500~K, and (right) 2000~K.}
\end{figure}

Figure~\ref{fig4} shows four rays launched at a latitude of
${-}$40\textdegree~and 400~km of altitude into the plasmasphere computed
using the diffusive equilibrium at temperatures varying from 1000 to
4000~K. For all rays, we find the ray path is the same, except that it
is shortened as the temperature increases. This is due to the shape of
density isolines, which remains the same while the value of the
isodensity is different according to the temperature. The more the
density value is large, the more matter the wave passes through, which
causes the waves to be more damped by the Landau damping effect. Thus,
as the temperature increases, waves vanish earlier on the ray path.

\begin{figure}
\includegraphics{fig04}
\caption{\label{fig4} Ray paths of a 5~kHz whistler wave for different
temperatures, from $T=1000$ to 4000~K. The ray path is the same but
shortens as the temperature increases.}
\end{figure}

\subsection{Influence of the density scaling}\label{sec42}

In Equation~(\ref{eq8}) the density parameter, $N_{b}$, scales all
values of density in the plasmasphere leading to different behaviors
for a given wave. In Figure~\ref{fig5}, we propagate a 5~kHz whistler
wave for different values of $N_{b}$. For low density, the magnetic
field dominates the propagation, allowing the wave to reach the
atmosphere on the conjugate point of the launch location and to vanish
in the upper atmosphere without bouncing (e.g.\ red line in
Figure~\ref{fig5}). The magnetic field acts as a wave guide for the
propagation. At high latitude, when density is high enough compared
with the magnetic field, the lower hybrid frequency can exceed the wave
frequency and permits the lower hybrid resonance to occur. In that
case, the wave resonance cone is reduced to zero as the surface of
refractive index becomes closed. A normal propagation, i.e.\ a
perpendicular propagation to the local magnetic field, is thus
theoretically allowed so that the wave's normal angle can reach
90\textdegree. As the wave propagates perpendicularly, its wave normal
angle will now exceed 90\textdegree, so that the wave is reflected
back. A bounce has happened. Several of these magnetic bounces can
occur in a common propagation. This effect is combined with the
previously described effect of the density effect. The more the
encountered density is high, the more the wave damping is strong, so
that the ray path is shortened along the ray path.

\begin{figure}
\includegraphics{fig05}
\caption{\label{fig5} Ray paths of a 5~kHz whistler wave for different
scale density at 2000~K. With higher density at fixed magnetic
amplitude magnetoreflections are allowed.}
\end{figure}

\section{Influence of wave parameters: the frequency}\label{sec5}


Here, we study the importance of the wave frequency within the whistler
range, as this parameter influences the escape point from the
ionosphere, the bounce at high latitude, and the region reached by the
wave. Using an electron temperature of 2000~K and an electron density
of $n_{e}=2.7\times 10^{9}~\mathrm{m}^{-3}$, we propagate
several waves at frequencies ranging from 1~kHz to 15~kHz within the
plasmasphere. Results are presented in Figure~\ref{fig6}. The 15~kHz
ray only follows the field line because the lower hybrid frequency
never comes near the wave frequency so that magnetoreflections are
never allowed and the ray continues to follow the field line. This ray
reaches the Earth's surface before its wave normal angle can reach the
90\textdegree~critical value for allowing propagation in the opposite
direction. Lower frequency rays, from 1 to 10~kHz, are conversely able
to bounce between the two hemispheres. The lowest frequency waves at 1
and 2~kHz have bounced trajectories that are outward directed, i.e.\ 
each mirror point increases in altitude. However, higher frequency
waves of 5 and 10~kHz have an inward propagation after the first
bounce, i.e.\ each mirror point's altitude successively decreases. This
is due to a competition between the magnetic field shape and the
density.

\begin{figure}
\includegraphics{fig06}
\caption{\label{fig6}Ray path for several wave frequencies in the
whistler range, from 1~kHz to 10~kHz, for an ambient electron plasma at
2000~K. The 15~kHz wave does not bounce, since the lower hybrid
frequency is never reached. Inward and outward propagation is allowed
according to the frequency values.}
\end{figure}

\section{Conclusions}\label{secVI}

In this article, we study the main parameters of the propagation of
very low frequency waves in the plasmasphere using the British
Antarctic Survey (BAS) HOTRAY code. Comparing various magnetic field
models, we find that the relative error between a pure dipole field
model and the IGRF model reaches a maximum of 101\% in the SAA region
at 100~km. The relative error reduces to 53\% for the eccentric tilted
dipole field model. An ongoing task is to implement the IGRF model into
HOTRAY. The electron density of the plasmasphere model has two main
parameters, the temperature of the diffusive equilibrium and the scale
density. The temperature diffuses the density to higher L-shells. We
find it does not change the ray path but shortens it as temperature
increases due to stronger Landau damping. The scale density defines
isodensity lines, which can allow (or not) magnetoreflections according
to the value of the local lower hybrid frequency. Finally, the
frequency determines whether or not the propagating wave bounces. As
the bounce is allowed (or not), the propagation remains outward for
low-frequency waves (${&lt;}$5~kHz in our simulation). For higher
frequencies, the propagation evolves to an earthward propagation as
soon as the first bounce has occurred. 

\section*{Declaration of interests}
The authors do not work for, advise, own shares in, or receive funds
from any organisation that could benefit from this article, and have
declared no affiliation other than their research organisations.

\section*{Funding}

J-FR and LC thank the Direction G\'{e}n\'{e}rale de l'Armement (DGA)
and the Agence pour l'Innovation de D\'{e}fense (AID) for funding the
ASTRID project ``PACTE- ESPACE''.

\back{}

\printbibliography
\refinput{crphys20250054-reference.tex}

\end{document}
</pre></body></html>