Journal topic
Ann. Geophys., 38, 207–230, 2020
https://doi.org/10.5194/angeo-38-207-2020

Special issue: Solar magnetism from interior to corona and beyond

Ann. Geophys., 38, 207–230, 2020
https://doi.org/10.5194/angeo-38-207-2020

Regular paper 10 Feb 2020

Regular paper | 10 Feb 2020

# Model of the propagation of very low-frequency beams in the Earth–ionosphere waveguide: principles of the tensor impedance method in multi-layered gyrotropic waveguides

Model of the propagation of very low-frequency beams in the Earth–ionosphere waveguide: principles of the tensor impedance method in multi-layered gyrotropic waveguides
Yuriy Rapoport1,6, Vladimir Grimalsky2, Viktor Fedun3, Oleksiy Agapitov1,4, John Bonnell4, Asen Grytsai1, Gennadi Milinevsky1,5, Alex Liashchuk6, Alexander Rozhnoi7, Maria Solovieva7, and Andrey Gulin1 Yuriy Rapoport et al.
• 1Faculty of Physics, Taras Shevchenko National University of Kyiv, Kyiv, Ukraine
• 2IICBA, CIICAp, Autonomous University of the State of Morelos (UAEM), Cuernavaca, Morelos, Mexico
• 3Department of Automatic Control and Systems Engineering, The University of Sheffield, Sheffield, UK
• 4Space Science Laboratory, University of California, Berkeley, Berkeley, California, USA
• 5International Center of Future Science, College of Physics, Jilin University, Changchun, China
• 6National Space Facilities Control and Test Center, State Space Agency of Ukraine, Kyiv, Ukraine
• 7Shmidt Institute of Physics of the Earth, Russian Academy of Sciences, Moscow, Russia

Correspondence: Yuriy Rapoport (yuriy.rapoport@gmail.com)

Abstract

The modeling of very low-frequency (VLF) electromagnetic (EM) beam propagation in the Earth–ionosphere waveguide (WGEI) is considered. A new tensor impedance method for modeling the propagation of electromagnetic beams in a multi-layered and inhomogeneous waveguide is presented. The waveguide is assumed to possess the gyrotropy and inhomogeneity with a thick cover layer placed above the waveguide. The influence of geomagnetic field inclination and carrier beam frequency on the characteristics of the polarization transformation in the Earth–ionosphere waveguide is determined. The new method for modeling the propagation of electromagnetic beams allows us to study the (i) propagation of the very low-frequency modes in the Earth–ionosphere waveguide and, in perspective, their excitation by the typical Earth–ionosphere waveguide sources, such as radio wave transmitters and lightning discharges, and (ii) leakage of Earth–ionosphere waveguide waves into the upper ionosphere and magnetosphere. The proposed approach can be applied to the variety of problems related to the analysis of the propagation of electromagnetic waves in layered gyrotropic and anisotropic active media in a wide frequency range, e.g., from the Earth–ionosphere waveguide to the optical waveband, for artificial signal propagation such as metamaterial microwave or optical waveguides.

1 Introduction

The results of the analytical and numerical study of very low-frequency (VLF) electromagnetic (EM) wave/beam propagation in the lithosphere–atmosphere–ionosphere–magnetosphere system (LAIM), in particular in the Earth–ionosphere waveguide (WGEI), are presented. The amplitude and phase of the VLF wave propagates in the Earth–ionosphere waveguide can change, and these changes may be observable using ground-based and/or satellite detectors. This reflects the variations in ionospheric electrodynamic characteristics (complex dielectric permittivity) and the influences on the ionosphere, for example, “from above” by the Sun–solar-wind–magnetosphere–ionosphere chain . Then the influence on the ionosphere “from below” comes from the most powerful meteorological, seismogenic and other sources in the lower atmosphere and lithosphere and Earth, such as cyclones and hurricanes as well as from earthquakes and tsunamis. From inside the ionosphere, strong thunderstorms, lightning discharges, and terrestrial gamma-ray flashes or sprite streamers influence the ionospheric electrodynamic characteristics as well. Note that the VLF signals are very important for the merging of atmospheric physics and space plasma physics with astrophysics and high-energy physics. The corresponding “intersection area” for these four disciplines includes cosmic rays and currently very popular objects of investigation – high-altitude discharges (sprites), anomalous X-ray bursts and powerful gamma-ray bursts. The key phenomena for the occurrence of all of these objects is the appearance of the runaway avalanche in the presence of high-energy seed electrons. In the atmosphere, there are cosmic-ray secondary electrons . Consequently, these phenomena are intensified during the air shower generating by cosmic particles . The runaway breakdown and lightning discharges including high-altitude ones can cause radio emission both in the high-frequency (HF) range, which could be observed using the Low-Frequency Array (LOFAR) radio telescope network facility and other radio telescopes , and in the VLF range . The corresponding experimental research includes the measurements of the VLF characteristics by the international measurement system of the “transmitted-receiver” pairs separated by a distance of a couple thousand kilometers . The World Wide Lightning Location Network is one of the international facilities for VLF measurements during thunderstorms with lightning discharges . Intensification of magnetospheric research, wave processes, particle distribution and wave–particle interaction in the magnetosphere including radiation belts leads to the great interest in VLF plasma waves, in particular whistlers .

The differences of the proposed model for the simulation of VLF waves in the WGEI from others can be summarized in three main points. (i) In distinction to the impedance invariant imbedding model , our model provides an optimal balance between the analytical and numerical approaches. It combines analytical and numerical approaches based on the matrix sweep method . As a result, this model allows for analytically obtaining the tensor impedance and, at the same time, provides high effectiveness and stability for modeling. (ii) In distinction to the full-wave finite-difference time domain models , our method provides the physically clear lower and upper boundary conditions, in particular physically justified upper boundary conditions corresponding to the radiation of the waves propagation in the WGEI to the upper ionosphere and magnetosphere. This allows for the determination of the leakage modes and the interpretation not only of ground-based but also satellite measurements of the VLF beam characteristics. (iii) In distinction to the models considered in , , based on the mode presentations and made in the frequency domain, we use the combined approach. This approach includes the radiation condition at the altitudes of the F region, equivalent impedance conditions in the lower E region and at the lower boundary of the WGEI, the mode approach, and finally, the beam method. This combined approach, finally, creates the possibility to adequately interpret data of both ground-based and satellite detection of the VLF EM wave/beam propagating in the WGEI and those, which experienced a leakage from the WGEI into the upper ionosphere and magnetosphere. Some other details on the distinctions from the previously published models are given below in Sect. 3.

The methods of effective boundary conditions such as effective impedance conditions are well known and can be used, in particular, for the layered metal-dielectric, metamaterial and gyrotropic active layered and waveguiding media of different types including plasma-like solid state and space plasma (Wait1996) media. The plasma wave processes in the metal–semiconductor–dielectric waveguide structures, placed into the external magnetic field, were widely investigated from radio to optical-frequency ranges. Corresponding waves are applied in modern plasmonics and in the non-destructive testing of semiconductor interfaces. It is interesting to realize the resonant interactions of volume and surface electromagnetic waves in these structures, so the simulations of the wave spectrum are important. To describe such complex layered structures, it is very convenient and effective to use the impedance approach . As a rule, impedance boundary conditions are used when the layer covering waveguide is thin . One of the known exclusions is the impedance invariant imbedding model. The difference between our new method and that model is already mentioned above and is explained in more detail in Sect. 3.3. Our new approach, i.e., a new tensor impedance method for modeling the propagation of electromagnetic beams (TIMEB), includes a set of very attractive features for practical purposes. These features are (i) that the surface impedance characterizes cover layer of a finite thickness, and this impedance is expressed analytically; (ii) that the method allows for an effective modeling of 3-D beam propagating in the gyrotropic waveguiding structure; (iii) finally that if the considered waveguide can be modified by any external influence such as bias magnetic or electric fields, or by any extra wave or energy beams (such as acoustic or quasistatic fields, etc.), the corresponding modification of the characteristics (phase and amplitude) of the VLF beam propagating in the waveguide structure can be modeled.

Our approach was properly employed and is suitable for the further development which will allow to solve also the following problems: (i) the problem of the excitation of the waveguide by the waves incident on the considered structure from above could be solved as well with the slight modification of the presented model, with the inclusion also ingoing waves; (ii) the consideration of a plasma-like system placed into the external magnetic field, such as the LAIM system or dielectric-magnetized semiconductor structure. The electromagnetic waves radiated outside the waveguiding structure, such as helicons or whistlers (Wait1996), and the waveguide modes could be considered altogether. An adequate boundary radiation condition on the upper boundary of the covering layer is derived. Based on this and an absence of ingoing waves, the leakage modes above the upper boundary of the structure (in other words, the upper boundary of covering layer) will be searched with the further development of the model delivered in the present paper. Namely, it will be possible to investigate the process of the leakage of electromagnetic waves from the open waveguide. Then their transformation into magnetized plasma waves and propagating along magnetic field lines, and the following excitation of the waveguiding modes by the waves incident on the system from external space (Walker1976) can be modeled as a whole. Combining with the proper measurements of the phases and amplitudes of the electromagnetic waves, propagating in the waveguiding structures and leakage waves, the model can be used for searching and even monitoring the external influences on the layered gyrotropic active artificial or natural media, for example, microwave or optical waveguides or the system LAIM and the Earth–ionosphere waveguide, respectively.

An important effect of the gyrotropy and anisotropy is the corresponding transformation of the field polarization during the propagation in the WGEI, which is absent in the ideal metal planar waveguide without gyrotropy and anisotropy. We will determine how such an effect depends on the carrier frequency of the beam, propagating in the WGEI, and the inclination of the geomagnetic field and perturbations in the electron concentration, which could vary under the influences of the sources powerful enough placed “below”, “above” and/or “inside” the ionosphere.

In Sect. 2 the formulation of the problem is presented. In Sect. 3 the algorithm is discussed including the determination of the VLF wave/beam radiation conditions into the upper ionosphere and magnetosphere at the upper boundary, placed in the F region at 250–400 km altitude. The effective tensor impedance boundary conditions at the upper boundary (∼85km) of the effective Earth–ionosphere waveguide and the 3-D model TIMEB of the propagation of the VLF beam in the WGEI are discussed as well. The issues regarding the VLF beam leakage regimes are considered only very briefly, since the relevant details will be presented in the following articles. In Sect. 4 the results of the numerical modeling are presented. In Sect. 5 the discussion is presented, including an example of the qualitative comparison between the results of our theory and an experiment including the future rocket experiment on the measurements of the characteristics of VLF signal radiated from the artificial VLF transmitter, which is propagating in the WGEI and penetrating into the upper ionosphere.

2 Formulation of the problem

The VLF electromagnetic waves with frequencies of f=(10–100) kHz can propagate along the Earth's surface for long distances >1000km. The Earth's surface of a high conductivity of z=0 (z is vertical coordinate) and the ionosphere F layer of z=300km form the VLF waveguide (see Fig. 1). The propagation of the VLF electromagnetic radiation excited by a near-Earth antenna within the WGEI should be described by the full set of Maxwell's equations in the isotropic atmosphere at $\mathrm{0}km, the approximate altitude of the nearly isotropic ionosphere D layer at $\mathrm{60}\phantom{\rule{0.125em}{0ex}}\mathrm{km}, and the anisotropic E and F layers of the ionosphere, due to the geomagnetic field H0, added by the boundary conditions at the Earth's surface and at the F layer. In Fig. 1, θ is the angle between the directions of the vertical axis z and geomagnetic field H0. Note that the angle of θ is complementary to the angle of inclination of the geomagnetic field. The geomagnetic field H0 is directed along z axis and lies in the plane xz, while the planes ${x}^{\prime }{z}^{\prime }$ and xz coincide with each other.

Figure 1The geometry of the anisotropic and gyrotropic waveguide. EM waves propagate in the OX direction. H0 is the external magnetic field. The effective WGEI for EM waves occupies the region $\mathrm{0}. The isotropic medium occupies the region $\mathrm{0}, LISO<Lz. The anisotropic and gyrotropic medium occupies the region ${L}_{\mathrm{ISO}}. The covering layer occupies the region ${L}_{z}. WGEI includes the isotropic region $\mathrm{0} and a part of the anisotropic region ${L}_{z}. It is supposed that the anisotropic region is relatively small part of the WG, $\left({L}_{z}-{L}_{\mathrm{ISO}}\right)/{L}_{z}\sim$ (0.1–0.2). At the upper boundary of the covering layer (z=Lmax) EM radiation to the external region (z>Lmax) is accounted for with the proper boundary conditions. The integration of the equations describing the EM field propagation allows for obtaining effective impedance boundary conditions at the upper boundary of the effective WG (z=Lz). These boundary conditions effectively include all the effects on the wave propagation of the covering layer and the radiation (at z=Lmax) to the external region (z>Lmax).

3 Algorithm

The boundary conditions and calculations of impedance and beam propagation in the WGEI are considered in this section. The other parts of the algorithm, e.g., the reflection of the EM waves from the WGEI effective upper boundary and the leakage of EM waves from the WGEI into the upper ionosphere and magnetosphere, will be presented very briefly as they are the subjects of the next papers.

## 3.1 Direct and inverse tensors characterizing the ionosphere

In the next subsections we derived the formulas describing the transfer of the boundary conditions at the upper boundary (z=Lmax) (Fig. 1), resulting in the tensor impedance conditions at the upper boundary of the effective WGEI (z=Li). Firstly let us describe the tensors characterizing the ionosphere.

The algorithm's main goal is to transfer the EM boundary conditions from the upper ionosphere at the height of Lz at ∼250–400 km to the lower ionosphere at Lz∼70–90 km. All components of the monochromatic EM field are considered to be proportional to exp (iωt). The anisotropic medium is inhomogeneous along OZ axis only and characterized by the permittivity tensor $\stackrel{\mathrm{^}}{\mathit{\epsilon }}\left(\mathit{\omega },z\right)$ or by the inverse tensor $\stackrel{\mathrm{^}}{\mathit{\beta }}\left(\mathit{\omega },z\right)={\stackrel{\mathrm{^}}{\mathit{\epsilon }}}^{-\mathrm{1}}\left(\mathit{\omega },z\right)$: $\mathbit{E}=\stackrel{\mathrm{^}}{\mathit{\beta }}\left(\mathit{\omega },z\right)\cdot \mathbit{D}$, where D is the electric induction. Below the absolute units are used. The expressions for the components of the effective permittivity of the ionosphere are in the coordinate frame ${X}^{\prime }Y{Z}^{\prime }$, where the OZ axis is aligned along the geomagnetic field H0.

$\begin{array}{}\text{(1)}& \begin{array}{rl}{\stackrel{\mathrm{^}}{\mathit{\epsilon }}}^{\prime }& =\left(\begin{array}{ccc}{\mathit{\epsilon }}_{\mathrm{1}}& {\mathit{\epsilon }}_{h}& \mathrm{0}\\ -{\mathit{\epsilon }}_{h}& {\mathit{\epsilon }}_{\mathrm{1}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& {\mathit{\epsilon }}_{\mathrm{3}}\end{array}\right),\\ {\mathit{\epsilon }}_{\mathrm{1}}& =\mathrm{1}-\frac{{\mathit{\omega }}_{\mathrm{pe}}^{\mathrm{2}}\cdot \left(\mathit{\omega }-i{\mathit{\nu }}_{\mathrm{e}}\right)}{\left(\left(\mathit{\omega }-i{\mathit{\nu }}_{\mathrm{e}}{\right)}^{\mathrm{2}}-{\mathit{\omega }}_{H\mathrm{e}}^{\mathrm{2}}\right)\cdot \mathit{\omega }}\\ & -\frac{{\mathit{\omega }}_{\mathrm{pi}}^{\mathrm{2}}\cdot \left(\mathit{\omega }-i{\mathit{\nu }}_{\mathrm{i}}\right)}{\left(\left(\mathit{\omega }-i{\mathit{\nu }}_{\mathrm{i}}{\right)}^{\mathrm{2}}-{\mathit{\omega }}_{H\mathrm{i}}^{\mathrm{2}}\right)\cdot \mathit{\omega }},\\ {\mathit{\epsilon }}_{h}& \equiv ig;\\ g& =-\frac{{\mathit{\omega }}_{\mathrm{pe}}^{\mathrm{2}}\cdot {\mathit{\omega }}_{H\mathrm{e}}}{\left(\left(\mathit{\omega }-i{\mathit{\nu }}_{\mathrm{e}}{\right)}^{\mathrm{2}}-{\mathit{\omega }}_{H\mathrm{e}}^{\mathrm{2}}\right)\cdot \mathit{\omega }}\\ & +\frac{{\mathit{\omega }}_{\mathrm{pi}}^{\mathrm{2}}\cdot {\mathit{\omega }}_{H\mathrm{i}}}{\left(\left(\mathit{\omega }-i{\mathit{\nu }}_{\mathrm{i}}{\right)}^{\mathrm{2}}-{\mathit{\omega }}_{H\mathrm{i}}^{\mathrm{2}}\right)\cdot \mathit{\omega }},\\ {\mathit{\epsilon }}_{\mathrm{3}}& =\mathrm{1}-\frac{{\mathit{\omega }}_{\mathrm{pe}}^{\mathrm{2}}}{\left(\mathit{\omega }-i{\mathit{\nu }}_{\mathrm{e}}\right)\cdot \mathit{\omega }}-\frac{{\mathit{\omega }}_{\mathrm{pi}}^{\mathrm{2}}}{\left(\mathit{\omega }-i{\mathit{\nu }}_{\mathrm{i}}\right)\cdot \mathit{\omega }};\\ {\mathit{\omega }}_{\mathrm{pe}}^{\mathrm{2}}& =\frac{\mathrm{4}\mathit{\pi }{e}^{\mathrm{2}}n}{{m}_{\mathrm{e}}},\\ {\mathit{\omega }}_{\mathrm{pi}}^{\mathrm{2}}& =\frac{\mathrm{4}\mathit{\pi }{e}^{\mathrm{2}}n}{{m}_{\mathrm{i}}},\\ {\mathit{\omega }}_{H\mathrm{e}}& =\frac{e{H}_{\mathrm{0}}}{{m}_{\mathrm{e}}c},\\ {\mathit{\omega }}_{H\mathrm{i}}& =\frac{e{H}_{\mathrm{0}}}{{m}_{\mathrm{i}}c},\end{array}\end{array}$

where ωpe, ωpi, ωHe and ωHi are the plasma and cyclotron frequencies for electrons and ions, respectively, me, mi, νe and νi are the masses and collision frequencies for electrons and ions, respectively, and n is the electron concentration. The approximations of the three-component plasma-like ionosphere (including one electron component, one effective ion and one effective neutral component) and quasi-neutrality are accepted. The expressions for the components of the permittivity tensor $\stackrel{\mathrm{^}}{\mathit{\epsilon }}\left(\mathit{\omega }$, z) are obtained from Eq. (1) by means of multiplication with the standard rotation matrices (Spiegel1959) dependent on angle θ. For the medium with a scalar conductivity σ, e.g., lower ionosphere or atmosphere, the effective permittivity in Eq. (1) reduces to $\mathit{\epsilon }=\mathrm{1}-\mathrm{4}\mathit{\pi }i\mathit{\sigma }/\mathit{\omega }$.

## 3.2 The equations for the EM field and upper boundary conditions

In the case of the VLF waveguide modes, the longitudinal wave number kx is slightly complex and should be calculated accounting for boundary conditions at the Earth's surface and the upper surface of the effective WGEI. The EM field depends on the horizontal coordinate x as exp (−ikxx). Taking into account kxk0 $\left({k}_{\mathrm{0}}=\mathit{\omega }/c\right)$, in simulations of the VLF beam propagation, it is possible to put kx=k0. Therefore, Maxwell's equations are

$\begin{array}{}\text{(2)}& \begin{array}{rl}& -\frac{\partial {H}_{y}}{\partial z}=i{k}_{\mathrm{0}}{D}_{x},\phantom{\rule{0.25em}{0ex}}\frac{\partial {H}_{x}}{\partial z}+i{k}_{x}{H}_{z}=i{k}_{\mathrm{0}}{D}_{y},\\ & -i{k}_{x}{H}_{y}=i{k}_{\mathrm{0}}{D}_{z}\\ & -\frac{\partial {E}_{y}}{\partial z}=-i{k}_{\mathrm{0}}{H}_{x},\phantom{\rule{0.25em}{0ex}}\frac{\partial {E}_{x}}{\partial z}+i{k}_{x}{E}_{z}=-i{k}_{\mathrm{0}}{H}_{y},\\ & -i{k}_{x}{E}_{y}=-i{k}_{\mathrm{0}}{H}_{z},\end{array}\end{array}$

where ${E}_{x}={\mathit{\beta }}_{\mathrm{11}}{D}_{x}+{\mathit{\beta }}_{\mathrm{12}}{D}_{y}+{\mathit{\beta }}_{\mathrm{13}}{D}_{z}$, etc. All components of the EM field can be represented through the horizontal components of the magnetic field Hx and Hy. The equations for these components are given below.

$\begin{array}{}\text{(3a)}& \begin{array}{rl}& \frac{\partial }{\partial z}\left(\frac{{\mathit{\beta }}_{\mathrm{22}}}{\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}}\frac{\partial {H}_{x}}{\partial z}\right)-\frac{\partial }{\partial z}\left(\frac{{\mathit{\beta }}_{\mathrm{21}}}{\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}}\frac{\partial {H}_{y}}{\partial z}\right)\\ & -i{k}_{x}\frac{\partial }{\partial z}\left(\frac{{\mathit{\beta }}_{\mathrm{23}}}{\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}}{H}_{y}\right)+{k}_{\mathrm{0}}^{\mathrm{2}}{H}_{x}=\mathrm{0}\end{array}\text{(3b)}& \begin{array}{rl}& \frac{\partial }{\partial z}\left(\left({\mathit{\beta }}_{\mathrm{11}}+\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}\frac{{\mathit{\beta }}_{\mathrm{12}}\cdot {\mathit{\beta }}_{\mathrm{21}}}{\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}}\right)\frac{\partial {H}_{y}}{\partial z}\right)-\frac{\partial }{\partial z}\left(\frac{{\mathit{\beta }}_{\mathrm{12}}}{\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}}\frac{\partial {H}_{x}}{\partial z}\right)\\ & +i{k}_{x}\frac{\partial }{\partial z}\left(\left({\mathit{\beta }}_{\mathrm{13}}+\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}\frac{{\mathit{\beta }}_{\mathrm{12}}\cdot {\mathit{\beta }}_{\mathrm{23}}}{\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}}\right){H}_{y}\right)\\ & +i{k}_{x}\left({\mathit{\beta }}_{\mathrm{31}}+\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}\frac{{\mathit{\beta }}_{\mathrm{32}}\cdot {\mathit{\beta }}_{\mathrm{21}}}{\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}}\right)\frac{\partial {H}_{y}}{\partial z}-\\ & -i{k}_{x}\frac{{\mathit{\beta }}_{\mathrm{32}}}{\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}}\frac{\partial {H}_{x}}{\partial z}\\ & +{k}_{\mathrm{0}}^{\mathrm{2}}\left(\mathrm{1}-{\mathit{\beta }}_{\mathrm{33}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}-\frac{{k}_{x}^{\mathrm{4}}}{{k}_{\mathrm{0}}^{\mathrm{4}}}\frac{{\mathit{\beta }}_{\mathrm{23}}\cdot {\mathit{\beta }}_{\mathrm{32}}}{\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}}\right){H}_{y}=\mathrm{0}\end{array}\end{array}$

The expressions for the horizontal components of the electric field Ex, Ey are given below.

$\begin{array}{}\text{(4)}& \begin{array}{rl}{E}_{x}& =\frac{i}{{k}_{\mathrm{0}}}\left(\left({\mathit{\beta }}_{\mathrm{11}}+\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}\frac{{\mathit{\beta }}_{\mathrm{12}}\cdot {\mathit{\beta }}_{\mathrm{21}}}{\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}}\right)\frac{\partial {H}_{y}}{\partial z}-\frac{{\mathit{\beta }}_{\mathrm{12}}}{\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}}\frac{\partial {H}_{y}}{\partial z}\right)\\ & -\frac{{k}_{x}}{{k}_{\mathrm{0}}}\left({\mathit{\beta }}_{\mathrm{13}}+\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}\frac{{\mathit{\beta }}_{\mathrm{12}}\cdot {\mathit{\beta }}_{\mathrm{23}}}{\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}}\right){H}_{y}\\ {E}_{y}& =\frac{i}{{k}_{\mathrm{0}}}\left(-\frac{{\mathit{\beta }}_{\mathrm{22}}}{\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}}\frac{\partial {H}_{x}}{\partial z}+\frac{{\mathit{\beta }}_{\mathrm{21}}}{\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}}\frac{\partial {H}_{y}}{\partial z}\right)\\ & -\frac{{k}_{x}}{{k}_{\mathrm{0}}}\frac{{\mathit{\beta }}_{\mathrm{23}}}{\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}}{H}_{y}\end{array}\end{array}$

In the region zLmax the upper ionosphere is assumed to be weakly inhomogeneous, and the geometric optics approximation is valid in the VLF range there. We should note that due to high inhomogeneity of the ionosphere in the vertical direction within the E layer (i.e., at the upper boundary of the effective VLF WGEI) such an approximation is not applicable. These conditions determine the choice of the upper boundary of z=Lmax at  250–400 km, where the conditions of the radiation are formulated. The dispersion equation which connected the wave numbers and the frequency of the outgoing waves is obtained from Eq. (3a, 3b), where ${H}_{x,y}\sim \mathrm{exp}\left(-i{k}_{z}z\right)$, while the derivatives like $\partial {\mathit{\beta }}_{\mathrm{11}}/\partial z$ and the inhomogeneity of the media are neglected.

$\begin{array}{}\text{(5)}& \begin{array}{rl}& \left({\mathit{\beta }}_{\mathrm{22}}{k}_{z}^{\mathrm{2}}-{k}_{\mathrm{0}}^{\mathrm{2}}\left(\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}\right)\right)\cdot \left(\left({\mathit{\beta }}_{\mathrm{11}}\left(\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}\right)\right\right\\ & +\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}{\mathit{\beta }}_{\mathrm{12}}\cdot {\mathit{\beta }}_{\mathrm{21}}){k}_{z}^{\mathrm{2}}\\ & +\left({\mathit{\beta }}_{\mathrm{13}}+{\mathit{\beta }}_{\mathrm{31}}\right)\left(\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}\right)+\\ & \frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}\left({\mathit{\beta }}_{\mathrm{12}}\cdot {\mathit{\beta }}_{\mathrm{23}}+{\mathit{\beta }}_{\mathrm{32}}\cdot {\mathit{\beta }}_{\mathrm{21}}\right){k}_{x}{k}_{z}\\ & -{k}_{\mathrm{0}}^{\mathrm{2}}\left(\left(\mathrm{1}-{\mathit{\beta }}_{\mathrm{33}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}\right)\left(\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}\right)-\frac{{k}_{x}^{\mathrm{4}}}{{k}_{\mathrm{0}}^{\mathrm{4}}}{\mathit{\beta }}_{\mathrm{23}}\cdot {\mathit{\beta }}_{\mathrm{32}}\right))-\\ & -\left({\mathit{\beta }}_{\mathrm{21}}{k}_{z}^{\mathrm{2}}+{\mathit{\beta }}_{\mathrm{23}}{k}_{x}{k}_{z}\right)\cdot \left({\mathit{\beta }}_{\mathrm{12}}{k}_{z}^{\mathrm{2}}-{\mathit{\beta }}_{\mathrm{32}}{k}_{x}{k}_{z}\right)=\mathrm{0}\end{array}\end{array}$

Thus, generally Eq. (5) determined the wave numbers for the outgoing waves to be of the fourth order (Wait1996). The boundary conditions at the upper boundary z=Lmax within the ionosphere F layer are the absence of the ingoing waves, i.e., the outgoing radiated (leakage) waves are present only. Two roots should be selected that possess the negative imaginary parts $\mathrm{Im}\left({k}_{z\mathrm{1},z\mathrm{2}}\right)<\mathrm{0}$; i.e., the outgoing waves dissipate upwards. However, in the case of VLF waves, some simplification can be used. Namely, the expressions for the wave numbers k1,2 are obtained from Eq. (3a, 3b), where the dependence on x is neglected: $|{k}_{\mathrm{1},\mathrm{2}}|\gg {k}_{\mathrm{0}}$. This approximation is valid within the F layer where the first outgoing wave corresponds to the whistlers with a small dissipation; the second one is the highly dissipating slow wave. To formulate the boundary conditions for Eq. (3a, 3b) at zLmax, the EM field components can be presented as

$\begin{array}{}\text{(6)}& \begin{array}{rl}{H}_{x}& ={A}_{\mathrm{1}}{e}^{-i{k}_{z\mathrm{1}}\stackrel{\mathrm{̃}}{z}}+{\mathit{\alpha }}_{\mathrm{2}}{A}_{\mathrm{2}}{e}^{-i{k}_{z\mathrm{2}}\stackrel{\mathrm{̃}}{z}},\\ {H}_{y}& ={\mathit{\alpha }}_{\mathrm{1}}{A}_{\mathrm{1}}{e}^{-i{k}_{z\mathrm{1}}\stackrel{\mathrm{̃}}{z}}+{A}_{\mathrm{2}}{e}^{-i{k}_{z\mathrm{2}}\stackrel{\mathrm{̃}}{z}}.\end{array}\end{array}$

In Eq. (6), $\stackrel{\mathrm{̃}}{z}=z-{L}_{z}$ . Equation (3a, 3b) are simplified in the approximation described above.

$\begin{array}{}\text{(7)}& \begin{array}{rl}{\mathit{\beta }}_{\mathrm{22}}\frac{{\partial }^{\mathrm{2}}{H}_{x}}{\partial {z}^{\mathrm{2}}}& -{\mathit{\beta }}_{\mathrm{21}}\frac{{\partial }^{\mathrm{2}}{H}_{y}}{\partial {z}^{\mathrm{2}}}+{k}_{\mathrm{0}}^{\mathrm{2}}{H}_{x}=\mathrm{0},\\ {\mathit{\beta }}_{\mathrm{11}}\frac{{\partial }^{\mathrm{2}}{H}_{y}}{\partial {z}^{\mathrm{2}}}& -{\mathit{\beta }}_{\mathrm{12}}\frac{{\partial }^{\mathrm{2}}{H}_{x}}{\partial {z}^{\mathrm{2}}}+{k}_{\mathrm{0}}^{\mathrm{2}}{H}_{y}=\mathrm{0}\end{array}\end{array}$

The solution of Eq. (7) is searched for as ${H}_{x,y}\sim \mathrm{exp}\left(-i{k}_{z}\stackrel{\mathrm{̃}}{z}\right)$. The following equation has been obtained to get the wave numbers kz1,z2 from Eq. (7):

$\begin{array}{}\text{(8)}& \begin{array}{rl}& {\mathit{\kappa }}^{\mathrm{4}}-\left({\mathit{\beta }}_{\mathrm{22}}+{\mathit{\beta }}_{\mathrm{11}}\right){\mathit{\kappa }}^{\mathrm{2}}+{\mathit{\beta }}_{\mathrm{11}}{\mathit{\beta }}_{\mathrm{22}}-{\mathit{\beta }}_{\mathrm{12}}{\mathit{\beta }}_{\mathrm{21}}=\mathrm{0},\\ & {\mathit{\kappa }}^{\mathrm{2}}=\frac{{k}_{\mathrm{0}}^{\mathrm{2}}}{{k}_{z}^{\mathrm{2}}}.\end{array}\end{array}$

Therefore, from Eq. (8) follows,

$\begin{array}{}\text{(9)}& \begin{array}{rl}{\mathit{\kappa }}_{\mathrm{1},\mathrm{2}}^{\mathrm{2}}& =\frac{{\mathit{\beta }}_{\mathrm{11}}+{\mathit{\beta }}_{\mathrm{22}}}{\mathrm{2}}±{\left({\left(\frac{{\mathit{\beta }}_{\mathrm{11}}+{\mathit{\beta }}_{\mathrm{22}}}{\mathrm{2}}\right)}^{\mathrm{2}}+{\mathit{\beta }}_{\mathrm{12}}{\mathit{\beta }}_{\mathrm{21}}\right)}^{\mathrm{1}/\mathrm{2}};\\ {\mathit{\alpha }}_{\mathrm{1}}& =\frac{{\mathit{\beta }}_{\mathrm{22}}-{\mathit{\kappa }}_{\mathrm{1}}^{\mathrm{2}}}{{\mathit{\beta }}_{\mathrm{21}}}=\frac{{\mathit{\beta }}_{\mathrm{12}}}{{\mathit{\beta }}_{\mathrm{11}}-{\mathit{\kappa }}_{\mathrm{1}}^{\mathrm{2}}};\\ {\mathit{\alpha }}_{\mathrm{2}}& =\frac{{\mathit{\beta }}_{\mathrm{11}}-{\mathit{\kappa }}_{\mathrm{2}}^{\mathrm{2}}}{{\mathit{\beta }}_{\mathrm{12}}}=\frac{{\mathit{\beta }}_{\mathrm{21}}}{{\mathit{\beta }}_{\mathrm{22}}-{\mathit{\kappa }}_{\mathrm{2}}^{\mathrm{2}}};\\ {k}_{z\mathrm{1},z\mathrm{2}}^{\mathrm{2}}& =\frac{{k}_{\mathrm{0}}^{\mathrm{2}}}{{\mathit{\kappa }}_{\mathrm{1},\mathrm{2}}^{\mathrm{2}}}.\end{array}\end{array}$

The signs of kz1,z2 have been chosen from the condition $\mathrm{Im}\left({k}_{z\mathrm{1},z\mathrm{2}}\right)<\mathrm{0}$. From Eq. (5) at the upper boundary of z=Lmax, the following relations are valid:

$\begin{array}{}\text{(10)}& {H}_{x}={A}_{\mathrm{1}}+{\mathit{\alpha }}_{\mathrm{2}}{A}_{\mathrm{2}},\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}{H}_{y}={\mathit{\alpha }}_{\mathrm{1}}{A}_{\mathrm{1}}+{A}_{\mathrm{2}}.\end{array}$

From Eq. (10) one can get

$\begin{array}{}\text{(11)}& \begin{array}{rl}{A}_{\mathrm{1}}& ={\mathrm{\Delta }}^{-\mathrm{1}}\left({H}_{x}-{\mathit{\alpha }}_{\mathrm{2}}{H}_{y}\right);\phantom{\rule{0.25em}{0ex}}{A}_{\mathrm{2}}={\mathrm{\Delta }}^{-\mathrm{1}}\left({H}_{y}-{\mathit{\alpha }}_{\mathrm{1}}{H}_{x}\right);\\ \mathrm{\Delta }& =\mathrm{1}-{\mathit{\alpha }}_{\mathrm{1}}{\mathit{\alpha }}_{\mathrm{2}}.\end{array}\end{array}$

Thus, it is possible to exclude the amplitudes of the outgoing waves A1,2 from Eqs. (9). As a result, at z=Lmax the boundary conditions are rewritten in terms of Hx and Hy only.

$\begin{array}{}\text{(12)}& \begin{array}{rl}\frac{\partial {H}_{x}}{\partial z}& =-i\left({k}_{z\mathrm{1}}{A}_{\mathrm{1}}+{k}_{z\mathrm{2}}{\mathit{\alpha }}_{\mathrm{2}}{A}_{\mathrm{2}}\right)\\ & =-\frac{i}{\mathrm{\Delta }}\left(\left({k}_{z\mathrm{1}}-{\mathit{\alpha }}_{\mathrm{1}}{\mathit{\alpha }}_{\mathrm{2}}{k}_{z\mathrm{2}}\right){H}_{x}+{\mathit{\alpha }}_{\mathrm{2}}\left({k}_{z\mathrm{2}}-{k}_{z\mathrm{1}}\right){H}_{y}\right)\\ \frac{\partial {H}_{y}}{\partial z}& =-i\left({k}_{z\mathrm{1}}{\mathit{\alpha }}_{\mathrm{1}}{A}_{\mathrm{1}}+{k}_{z\mathrm{2}}{A}_{\mathrm{2}}\right)\\ & =-\frac{i}{\mathrm{\Delta }}\left(\left({k}_{z\mathrm{2}}-{\mathit{\alpha }}_{\mathrm{1}}{\mathit{\alpha }}_{\mathrm{2}}{k}_{z\mathrm{1}}\right){H}_{y}+{\mathit{\alpha }}_{\mathrm{1}}\left({k}_{z\mathrm{1}}-{k}_{z\mathrm{2}}\right){H}_{x}\right)\end{array}\end{array}$

The relations in Eq. (12) are the upper boundary conditions of the radiation for the boundary z=Lmax at ∼250–400 km. These conditions will be transformed and recalculated using the analytical numerical recurrent procedure into equivalent impedance boundary conditions at z=Lz at ∼70–90 km.

Note that in the “whistler–VLF approximation” is valid at frequencies ∼10kHz for the F region of the ionosphere. In this approximation and kx≈0, we receive the dispersion equation using Eqs. (5), (8), (9), in the form

$\begin{array}{}\text{(13)}& {{k}_{z}^{\prime }}^{\mathrm{2}}{k}^{\mathrm{2}}={k}_{\mathrm{0}}^{\mathrm{2}}{g}^{\mathrm{2}},\end{array}$

where ${k}^{\mathrm{2}}={k}_{x}^{\mathrm{2}}+{k}_{z}^{\mathrm{2}}={{k}_{x}^{\prime }}^{\mathrm{2}}+{{k}_{z}^{\prime }}^{\mathrm{2}}$ and ${k}_{x}^{\prime }$ and ${k}_{z}^{\prime }$ are the transverse and longitudinal components of the wave number relative to geomagnetic field. For the F region of the ionosphere, where ${\mathit{\nu }}_{\mathrm{e}}\ll \mathit{\omega }\ll {\mathit{\omega }}_{H\mathrm{e}}$, Eq. (13) reduces to the standard form of the whistler dispersion equation $|{k}_{z}^{\prime }|k={k}_{\mathrm{0}}|g|$, $g\approx -{\mathit{\omega }}_{\mathrm{pe}}^{\mathrm{2}}/\left(\mathit{\omega }{\mathit{\omega }}_{H\mathrm{e}}\right)$ and $\mathit{\omega }={c}^{\mathrm{2}}k|{k}_{z}^{\prime }|\left({\mathit{\omega }}_{H\mathrm{e}}/{\mathit{\omega }}_{\mathrm{pe}}^{\mathrm{2}}\right)$. In a special case of the waves propagating along geomagnetic field, ${k}_{x}^{\prime }=\mathrm{0}$, for the propagating whistler waves, the well-known dispersion dependence is $\mathit{\omega }={c}^{\mathrm{2}}{{k}_{z}^{\prime }}^{\mathrm{2}}\left({\mathit{\omega }}_{H\mathrm{e}}/{\mathit{\omega }}_{\mathrm{pe}}^{\mathrm{2}}\right)$ . For the formulated problem we can reasonably assume kx≈0. Therefore Eq. (13) is reduced to ${k}_{z}^{\mathrm{4}}{\mathrm{cos}}^{\mathrm{2}}\mathit{\theta }={k}_{\mathrm{0}}^{\mathrm{4}}{g}^{\mathrm{2}}$. As a result, we get ${k}_{z\mathrm{1}}=\sqrt{g/\mathrm{cos}\mathit{\theta }}{k}_{\mathrm{0}}$ and ${k}_{z\mathrm{2}}=-i\sqrt{g/\mathrm{cos}\mathit{\theta }}{k}_{\mathrm{0}}$, and then, similar to the relations in Eq. (12), the boundary conditions can be presented in terms of the tangential components of the electric field as

$\begin{array}{}\text{(14)}& \begin{array}{rl}& \frac{\partial \mathbit{U}}{\partial z}+\stackrel{\mathrm{^}}{B}\mathbit{U}=\mathrm{0};\\ & \mathbit{U}=\left(\begin{array}{c}{E}_{x}\\ {E}_{y}\end{array}\right);\\ & \stackrel{\mathrm{^}}{B}=\frac{\mathrm{1}}{\mathrm{2}}\sqrt{\frac{g}{\mathrm{cos}\mathit{\theta }}}{k}_{\mathrm{0}}\left(\begin{array}{cc}\mathrm{1}+i& -\mathrm{1}-i\\ \mathrm{1}+i& \mathrm{1}+i\end{array}\right).\end{array}\end{array}$

Conditions in Eq. (12) or Eq. (14) are the conditions of radiation (absence of ingoing waves) formulated at the upper boundary at z=Lmax and suitable for the determination of the energy of the wave leaking from the WGEI into the upper ionosphere and magnetosphere. Note that Eqs. (12) and (14) expressed the boundary conditions of the radiation (more accurately speaking, an absence of incoming waves, which is the consequence to the causality principle) are obtained as a result of the passage to limit $|{k}_{x}/{k}_{z}|\to \mathrm{0}$ in Eq. (5). In spite of the disappearance of the dependence of these boundary conditions explicitly on kx, the dependence of the characteristics of the wave propagation process on kx, as a whole, is accounted for, and all results are still valid for the description of the wave beam propagation in the WGEI along the horizontal axis x with a finite kxk0.

## 3.3 Equivalent tensor impedance boundary conditions

The tensor impedance at the upper boundary of the effective WGEI of z=Lz (see Fig. 1) is obtained by the conditions of radiation in Eqs. (12) or (14), recalculated from the level of $z={L}_{\mathrm{max}}\sim \mathrm{250}$–400 km, placed in the F region of the ionosphere, to the level of $z={L}_{z}\sim \mathrm{80}$–90 km, placed in the E region.

The main idea of the effective tensor impedance method is the unification of analytical and numerical approaches and the derivation of the proper impedance boundary conditions without “thin-cover-layer” approximation. This approximation is usually used in the effective impedance approaches, applied either for artificial or natural layered gyrotropic structures (see e.g., ). There is one known exception, namely the invariant imbedding impedance method . The comparison of our method with the invariant imbedding impedance method will be presented at the end of this subsection. Equation (3a, 3b) jointly with the boundary conditions of Eq. (12) have been solved by finite differences. The derivatives in Eq. (3a, 3b) are approximated as

$\begin{array}{}\text{(15)}& \begin{array}{rl}\frac{\partial }{\partial z}\left(C\left(z\right)\frac{\partial {H}_{x}}{\partial z}\right)& \approx \frac{\mathrm{1}}{h}\left(C\left({z}_{j+\mathrm{1}/\mathrm{2}}\right)\frac{\left({H}_{x}{\right)}_{j+\mathrm{1}}-\left({H}_{x}{\right)}_{j}}{h}\right\\ & -C\left({z}_{j-\mathrm{1}/\mathrm{2}}\right)\frac{\left({H}_{x}{\right)}_{j}-\left({H}_{x}{\right)}_{j-\mathrm{1}}}{h})\\ \frac{\partial }{\partial z}\left(F\left(z\right){H}_{x}\right)& \approx \frac{\mathrm{1}}{\mathrm{2}h}\left(F\left({z}_{j+\mathrm{1}}\right)\left({H}_{x}{\right)}_{j+\mathrm{1}}\\ & -F\left({z}_{j-\mathrm{1}}\right)\left({H}_{x}{\right)}_{j-\mathrm{1}}\right).\end{array}\end{array}$

In Eq. (15) ${z}_{j+\mathrm{1}/\mathrm{2}}=h\cdot \left(j+\mathrm{0.5}\right)$. In Eq. (10) the approximation is $\partial {H}_{x}/\partial z\approx \left[\left({H}_{x}{\right)}_{N}-\left({H}_{x}{\right)}_{N-\mathrm{1}}\right]/h$. Here h is the discretization step along the OZ axis, and N is the total number of nodes. At each step j the difference approximations of Eq. (3a, 3b) take the form

$\begin{array}{}\text{(16)}& {\stackrel{\mathrm{^}}{\mathit{\alpha }}}_{j}^{\left(-\right)}\cdot {\mathbit{H}}_{j-\mathrm{1}}+{\stackrel{\mathrm{^}}{\mathit{\alpha }}}_{j}^{\left(\mathrm{0}\right)}\cdot {\mathbit{H}}_{j}+{\stackrel{\mathrm{^}}{\mathit{\alpha }}}_{j}^{\left(+\right)}\cdot {\mathbit{H}}_{j+\mathrm{1}}=\mathrm{0},\end{array}$

where ${\mathbit{H}}_{j}=\left(\begin{array}{c}{H}_{x}\\ {h}_{y}\end{array}\right)$, $j=N-\mathrm{1},N-\mathrm{2},\mathrm{\dots },\mathrm{1}$, ${z}_{j}=h\cdot j$ and ${L}_{z}=h\cdot N$. Due to the complexity of expressions for the matrix coefficients in Eq. (16), we have shown them in Appendix A. The set of the matrix Eq. (16) has been solved by the factorization method also known as an elimination and matrix sweep method (see ). It can be written as

$\begin{array}{}\text{(17a)}& {\mathbit{H}}_{j}={\stackrel{\mathrm{^}}{b}}_{j}\cdot {\mathbit{H}}_{j+\mathrm{1}},\phantom{\rule{0.25em}{0ex}}j=N,\mathrm{\dots }\mathrm{1}\text{(17b)}& \begin{array}{rl}& {H}_{xj+\mathrm{1}}={b}_{\mathrm{11}j+\mathrm{1}}{H}_{\mathrm{1}}+{b}_{\mathrm{12}j+\mathrm{1}}{H}_{\mathrm{2}};\\ & {H}_{yj+\mathrm{1}}={b}_{\mathrm{21}j+\mathrm{1}}{H}_{\mathrm{1}}+{b}_{\mathrm{22}j+\mathrm{1}}{H}_{\mathrm{2}};\\ & {H}_{\mathrm{1}}\equiv {H}_{xj};\\ & {H}_{\mathrm{2}}\equiv {H}_{yj}\end{array}.\end{array}$

This method is a variant of the Gauss elimination method for the matrix three-diagonal set of Eq. (16). The value of ${\stackrel{\mathrm{^}}{b}}_{N}$ is obtained from the boundary conditions (12) as

$\begin{array}{}\text{(18)}& {\stackrel{\mathrm{^}}{\mathit{\alpha }}}_{N}^{\left(-\right)}\cdot {\mathbit{H}}_{N-\mathrm{1}}+{\stackrel{\mathrm{^}}{\mathit{\alpha }}}_{N}^{\left(\mathrm{0}\right)}\cdot {\mathbit{H}}_{N}=\mathrm{0}.\end{array}$

Therefore ${\stackrel{\mathrm{^}}{b}}_{N}=-\left({\stackrel{\mathrm{^}}{\mathit{\alpha }}}_{N}^{\left(\mathrm{0}\right)}{\right)}^{-\mathrm{1}}\cdot {\stackrel{\mathrm{^}}{\mathit{\alpha }}}_{N}^{\left(-\right)}$. Then the matrices ${\stackrel{\mathrm{^}}{b}}_{j}$ have been computed sequentially down to the desired value of $z={L}_{z}=h\cdot {N}_{z}$, where the impedance boundary conditions are assumed to be applied. At each step the expression for ${\stackrel{\mathrm{^}}{b}}_{j}$ follows from Eqs. (16), (17a) and (17b) as

$\begin{array}{}\text{(19)}& \left({\stackrel{\mathrm{^}}{\mathit{\alpha }}}_{j}^{\left(\mathrm{0}\right)}+{\stackrel{\mathrm{^}}{\mathit{\alpha }}}_{j}^{\left(+\right)}\cdot {\stackrel{\mathrm{^}}{b}}_{j+\mathrm{1}}\right)\cdot {\mathbit{H}}_{j}=-{\stackrel{\mathrm{^}}{\mathit{\alpha }}}_{j}^{\left(-\right)}\cdot {\mathbit{H}}_{j-\mathrm{1}}=\mathrm{0}.\end{array}$

Therefore, for Eq. (17a, 17b), we obtain ${\stackrel{\mathrm{^}}{b}}_{j}=-\left({\stackrel{\mathrm{^}}{\mathit{\alpha }}}_{j}^{\left(\mathrm{0}\right)}+{\stackrel{\mathrm{^}}{\mathit{\alpha }}}_{j}^{\left(+\right)}\cdot {\stackrel{\mathrm{^}}{b}}_{j+\mathrm{1}}{\right)}^{-\mathrm{1}}\cdot {\stackrel{\mathrm{^}}{\mathit{\alpha }}}_{j}^{\left(-\right)}$. The derivatives in Eq. (4) have been approximated as

$\begin{array}{}\text{(20)}& \begin{array}{rl}{\left(\frac{\partial {H}_{x}}{\partial z}\right)}_{{N}_{z}}& \approx \frac{\left({H}_{x}{\right)}_{{N}_{z}+\mathrm{1}}-\left({H}_{x}{\right)}_{{N}_{z}}}{h}\\ & =\frac{\left({b}_{{N}_{z}+\mathrm{1}\phantom{\rule{0.25em}{0ex}}\mathrm{11}}-\mathrm{1}\right)\cdot \left({H}_{x}{\right)}_{{N}_{z}}+{b}_{{N}_{z}+\mathrm{1}\phantom{\rule{0.25em}{0ex}}\mathrm{12}}\cdot \left({H}_{y}{\right)}_{{N}_{z}}}{h},\end{array}\end{array}$

and a similar equation can be obtained for ${\left(\frac{\partial {H}_{y}}{\partial z}\right)}_{{N}_{z}}$. Note that as a result of this discretization, only the values at the grid level Nz are included in the numerical approximation of the derivatives $\partial {H}_{x,y}/\partial z$ at z=Lz. We determine tensor impedance $\stackrel{\mathrm{^}}{Z}$ at z=Lz at ∼85km. The tensor values are included in the following relations, all of which are corresponded to altitude (in other words, to the grid with the number Nz and corresponding to this altitude).

The equivalent tensor impedance is obtained using a two-step procedure. (1) We obtain the matrix ${\stackrel{\mathrm{^}}{b}}_{j}$ using Eq. (3a, 3b) with the boundary conditions of Eq. (12) and the procedure of Eqs. (17)–(19) described above. (2) Placing the expressions of Eq. (21) with tensor impedance into the left parts and the derivatives $\partial {H}_{x,y}/\partial z$ in Eq. (20) into the right parts of Eq. (4), the analytical expressions for the components of the tensor impedance are

$\begin{array}{}\text{(22)}& \begin{array}{rl}{Z}_{\mathrm{011}}& =-\frac{i}{{k}_{\mathrm{0}}h}\left(\frac{{\mathit{\beta }}_{\mathrm{21}}}{\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}}\cdot {b}_{\mathrm{21}}-\frac{{\mathit{\beta }}_{\mathrm{22}}}{\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}}\cdot \left({b}_{\mathrm{11}}-\mathrm{1}\right)\right),\\ {Z}_{\mathrm{012}}& =-\frac{i}{{k}_{\mathrm{0}}h}\left(\frac{{\mathit{\beta }}_{\mathrm{21}}}{\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}}\frac{\partial {H}_{y}}{\partial z}\cdot \left({b}_{\mathrm{22}}-\mathrm{1}\right)\right\\ & -\frac{{\mathit{\beta }}_{\mathrm{22}}}{\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}}\cdot {b}_{\mathrm{12}}-{k}_{x}h\cdot \frac{{\mathit{\beta }}_{\mathrm{23}}}{\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}}),\\ {Z}_{\mathrm{021}}& =\frac{i}{{k}_{\mathrm{0}}h}\left({\mathit{\beta }}_{\mathrm{11}}+\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}\frac{{\mathit{\beta }}_{\mathrm{12}}\cdot {\mathit{\beta }}_{\mathrm{21}}}{\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}}\cdot {b}_{\mathrm{21}}\right\\ & -\frac{{\mathit{\beta }}_{\mathrm{12}}}{\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}}\cdot \left({b}_{\mathrm{11}}-\mathrm{1}\right)),\\ {Z}_{\mathrm{022}}& =\frac{i}{{k}_{\mathrm{0}}h}\left(\left({\mathit{\beta }}_{\mathrm{11}}+\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}\frac{{\mathit{\beta }}_{\mathrm{12}}\cdot {\mathit{\beta }}_{\mathrm{21}}}{\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}}\right)\cdot \left({b}_{\mathrm{22}}-\mathrm{1}\right)\right\\ & -{k}_{x}h\cdot \left({\mathit{\beta }}_{\mathrm{13}}+\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}\frac{{\mathit{\beta }}_{\mathrm{12}}\cdot {\mathit{\beta }}_{\mathrm{23}}}{\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}}\right)-\frac{{\mathit{\beta }}_{\mathrm{12}}}{\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}}\cdot {b}_{\mathrm{12}}).\end{array}\end{array}$

The proposed method of the transfer of the boundary conditions from the ionosphere F layer at Lmax=250–400 km into the lower part of the E layer at Lz=80–90 km is stable and easily realizable in comparison with some alternative approaches based on the invariant imbedding methods . The stability of our method is due to the stability of the Gauss elimination method when the coefficients at the matrix central diagonal are dominating. The last is valid for the ionosphere with electromagnetic losses where the absolute values of the permittivity tensor are large. The application of the proposed matrix sweep method in the media without losses may require the use of the Gauss method with the choice of the maximum element to ensure stability. However, as our simulations (not presented here) demonstrated, for the electromagnetic problems in the frequency domain, the simple Gauss elimination and the choice of the maximal element give the same results. The accumulation of errors may occur in evolutionary problems in the time domain when the Gauss method should be applied sequentially many times. The use of the independent functions Hx and Hy in Eq. (3a, 3b) seems natural, as well as the transfer of Eq. (17a), because the impedance conditions are the expressions of the electric Ex and Ey through the magnetic components Hx and Hy at the upper boundary of the VLF waveguide at 80–90 km. The naturally chosen direction of the recalculation of the upper boundary conditions from z=Lmax to z=Lz, i.e., from the upper layer with a large impedance value to the lower-altitude layer with a relatively small impedance value, provides, at the same time, the stability of the simulation procedure. The obtained components of the tensor impedance are small, $|{Z}_{\mathrm{0}\mathit{\alpha }\mathit{\beta }}|\le \mathrm{0.1}$. This determines the choice of the upper boundary at z=Lz for the effective WGEI. Due to the small impedance, EM waves incident from below on this boundary are reflected effectively back. Therefore, the region $\mathrm{0}\le z\le {L}_{z}$ indeed can be presented as an effective WGEI. This waveguide includes not only lower boundary at LISO at ∼65–75 km with a rather small losses, but it also includes thin dissipative and anisotropic and gyrotropic layers between 75 and 85–90 km.

Finally, the main differences and advantages of the proposed tensor impedance method from other methods for impedance recalculating, and in particular invariant imbedding methods , can be summarized as follows:

1. In contrast to the invariant imbedding method, the currently proposed method can be used for the direct recalculation of tensor impedance, as is determined analytically (see Eq. 22).

2. For the media without non-locality, the proposed method does not require the solving of one or more integral equations.

3. The proposed method does not require forward and reflected waves. The conditions for the radiation at the upper boundary of z=Lmax (see Eq. 12) are determined through the total field components of Hx,y, which simplify the overall calculations.

4. The overall calculation procedure is very effective and computationally stable. Note that even for the very low-loss systems, the required level of stability can be achieved with modification based on the choice of the maximal element for matrix inversion.

## 3.4 Propagation of electromagnetic waves in the gyrotropic waveguide and the TIMEB method

Let us use the transverse components of the electric Ey and magnetic Hy fields to derive equations for the slow varying amplitudes $A\left(x,y,z\right)$ and $B\left(x,y,z\right)$ of the VLF beams. These components can be represented as

$\begin{array}{}\text{(23)}& \begin{array}{rl}{E}_{y}& =\frac{\mathrm{1}}{\mathrm{2}}A\left(x,y,z\right)\cdot {e}^{i\mathit{\omega }t-i{k}_{\mathrm{0}}x}+c.c.,\\ {H}_{y}& =\frac{\mathrm{1}}{\mathrm{2}}B\left(x,y,z\right)\cdot {e}^{i\mathit{\omega }t-i{k}_{\mathrm{0}}x}+c.c.\end{array}\end{array}$

Here we assumed kx=k0 to reflect beam propagation in the WGEI with the main part in the atmosphere and lower ionosphere (D region), which are similar to free space by its electromagnetic parameters. The presence of a thin anisotropic and dissipative layer belonging to the E region of the ionosphere causes, altogether with the impedance boundary condition, the proper z dependence of $B\left(x,y,z\right)$. Using Eqs. (21) and (22), the boundary conditions are determined at the height of z=Lz for the slowly varying amplitudes $A\left(x,y,z\right)$ and $B\left(x,y,z\right)$ of the transverse components Ey and Hy. As it follows from Maxwell's equations, the components Ex and Hx through Ey and Hy in the method of beams have the form

$\begin{array}{}\text{(24)}& {H}_{x}\approx -\frac{i}{{k}_{\mathrm{0}}}\frac{\partial {E}_{y}}{\partial z},\phantom{\rule{0.25em}{0ex}}{E}_{x}\approx {\mathit{\gamma }}_{\mathrm{12}}{E}_{y}+i\frac{{\stackrel{\mathrm{̃}}{\mathit{\beta }}}_{\mathrm{33}}}{{k}_{\mathrm{0}}}\frac{\partial {H}_{y}}{\partial z}+{\stackrel{\mathrm{̃}}{\mathit{\beta }}}_{\mathrm{13}}{H}_{y},\end{array}$

where ${\mathit{\gamma }}_{\mathrm{12}}={\mathrm{\Delta }}_{\mathrm{0}}^{-\mathrm{1}}\left({\mathit{\epsilon }}_{\mathrm{13}}{\mathit{\epsilon }}_{\mathrm{32}}-{\mathit{\epsilon }}_{\mathrm{12}}{\mathit{\epsilon }}_{\mathrm{33}}\right)$, ${\stackrel{\mathrm{̃}}{\mathit{\beta }}}_{\mathrm{13}}={\mathrm{\Delta }}_{\mathrm{0}}^{-\mathrm{1}}{\mathit{\epsilon }}_{\mathrm{13}}$ and ${\stackrel{\mathrm{̃}}{\mathit{\beta }}}_{\mathrm{33}}={\mathrm{\Delta }}_{\mathrm{0}}^{-\mathrm{1}}{\mathit{\epsilon }}_{\mathrm{33}}$; ${\mathrm{\Delta }}_{\mathrm{0}}={\mathit{\epsilon }}_{\mathrm{11}}{\mathit{\epsilon }}_{\mathrm{33}}-{\mathit{\epsilon }}_{\mathrm{13}}{\mathit{\epsilon }}_{\mathrm{31}}$. From Eqs. (21) and (24), the boundary conditions for A and B can be defined as

$\begin{array}{}\text{(25)}& \begin{array}{rl}& A-\frac{i}{{k}_{\mathrm{0}}}{Z}_{\mathrm{011}}\cdot \frac{\partial A}{\partial z}+{Z}_{\mathrm{012}}\cdot B\approx \mathrm{0},\\ & {\mathit{\gamma }}_{\mathrm{12}}\cdot A+\frac{i}{{k}_{\mathrm{0}}}{Z}_{\mathrm{021}}\cdot \frac{\partial A}{\partial z}+\left({\stackrel{\mathrm{̃}}{\mathit{\beta }}}_{\mathrm{13}}-{Z}_{\mathrm{022}}\right)\\ & \phantom{\rule{0.25em}{0ex}}\cdot B+\frac{i}{{k}_{\mathrm{0}}}{\stackrel{\mathrm{̃}}{\mathit{\beta }}}_{\mathrm{33}}\cdot \frac{\partial B}{\partial z}\approx \mathrm{0}.\end{array}\end{array}$

The evolution equations for the slowly varying amplitudes $A\left(x,y,z\right)$ and $B\left(x,y,z\right)$ of the VLF beams are derived. The monochromatic beams are considered when the frequency ω is fixed and the amplitudes do not depend on time t. Looking for the solutions for the EM field as $\mathbit{E},\mathbit{H}\sim \mathrm{exp}\left(i\mathit{\omega }t-i{k}_{x}x-i{k}_{y}y\right)$, Maxwell's equations are

$\begin{array}{}\text{(26)}& \begin{array}{rl}& -i{k}_{y}{H}_{z}-\frac{\partial {H}_{y}}{\partial z}=i{k}_{\mathrm{0}}{D}_{x},\\ & \frac{\partial {H}_{x}}{\partial z}+i{k}_{x}{H}_{z}=i{k}_{\mathrm{0}}{D}_{y},\\ & -i{k}_{x}{H}_{y}+i{k}_{y}{H}_{x}=i{k}_{\mathrm{0}}{D}_{z}\\ & -i{k}_{y}{E}_{z}-\frac{\partial {E}_{y}}{\partial z}=-i{k}_{\mathrm{0}}{H}_{x},\\ & \frac{\partial {E}_{x}}{\partial z}+i{k}_{x}{E}_{z}=-i{k}_{\mathrm{0}}{H}_{y},\\ & -i{k}_{x}{E}_{y}+i{k}_{y}{E}_{x}=-i{k}_{\mathrm{0}}{H}_{z}.\end{array}\end{array}$

Here ${D}_{x}={\mathit{\epsilon }}_{\mathrm{11}}{E}_{x}+{\mathit{\epsilon }}_{\mathrm{12}}{E}_{y}+{\mathit{\epsilon }}_{\mathrm{13}}{E}_{z}$. From Eq. (21), the equations for Ex and Ez through EyandHy are

$\begin{array}{}\text{(27)}& \begin{array}{rl}{E}_{x}& =\frac{\mathrm{1}}{{\mathrm{\Delta }}_{y}}\left\{\left[{\mathit{\epsilon }}_{\mathrm{13}}{\mathit{\epsilon }}_{\mathrm{32}}-\left({\mathit{\epsilon }}_{\mathrm{12}}+\frac{{k}_{x}{k}_{y}}{{k}_{\mathrm{0}}^{\mathrm{2}}}\right)\cdot \left({\mathit{\epsilon }}_{\mathrm{33}}-\frac{{k}_{y}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}\right)\right]{E}_{y}\right\\\ & +\frac{i}{{k}_{\mathrm{0}}}\left({\mathit{\epsilon }}_{\mathrm{33}}-\frac{{k}_{y}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}\right)\frac{\partial {H}_{y}}{\partial z}+\frac{{k}_{x}}{{k}_{\mathrm{0}}}{\mathit{\epsilon }}_{\mathrm{13}}\cdot {H}_{y}+i\frac{{k}_{y}}{{k}_{\mathrm{0}}^{\mathrm{2}}}{\mathit{\epsilon }}_{\mathrm{13}}\frac{\partial {E}_{y}}{\partial z}}\\ {E}_{z}& =\frac{\mathrm{1}}{{\mathrm{\Delta }}_{y}}\left\{\left[{\mathit{\epsilon }}_{\mathrm{31}}\left({\mathit{\epsilon }}_{\mathrm{12}}+\frac{{k}_{x}{k}_{y}}{{k}_{\mathrm{0}}^{\mathrm{2}}}\right)-{\mathit{\epsilon }}_{\mathrm{32}}\left({\mathit{\epsilon }}_{\mathrm{11}}-\frac{{k}_{y}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}\right)\right]{E}_{y}\right\\\ & -\frac{i}{{k}_{\mathrm{0}}}{\mathit{\epsilon }}_{\mathrm{31}}\frac{\partial {H}_{y}}{\partial z}-\frac{{k}_{x}}{{k}_{\mathrm{0}}}\left({\mathit{\epsilon }}_{\mathrm{11}}-\frac{{k}_{y}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}\right){H}_{y}\\ & -i\frac{{k}_{y}}{{k}_{\mathrm{0}}^{\mathrm{2}}}\left({\mathit{\epsilon }}_{\mathrm{11}}-\frac{{k}_{y}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}\right)\frac{\partial {E}_{y}}{\partial z}}.\end{array}\end{array}$

In Eq. (27), ${\mathrm{\Delta }}_{y}\equiv \left({\mathit{\epsilon }}_{\mathrm{11}}-\frac{{k}_{y}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}\right)\cdot \left({\mathit{\epsilon }}_{\mathrm{33}}-\frac{{k}_{y}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}\right)-{\mathit{\epsilon }}_{\mathrm{31}}\cdot {\mathit{\epsilon }}_{\mathrm{13}}$. The equations for Ey and Hy obtained from the Maxwell equations are

$\begin{array}{}\text{(28)}& \begin{array}{rl}\left(\frac{{\partial }^{\mathrm{2}}}{\partial {z}^{\mathrm{2}}}\right& -{k}_{x}^{\mathrm{2}}-{k}_{y}^{\mathrm{2}}){E}_{y}+i{k}_{y}\left(\frac{\partial {E}_{z}}{\partial z}-i{k}_{x}{E}_{x}-i{k}_{y}{E}_{y}\right)\\ & +{k}_{\mathrm{0}}^{\mathrm{2}}{D}_{y}=\mathrm{0};\phantom{\rule{0.25em}{0ex}}-i{k}_{\mathrm{0}}\frac{\partial {E}_{x}}{\partial z}+{k}_{x}{k}_{\mathrm{0}}{E}_{z}+{k}_{\mathrm{0}}^{\mathrm{2}}{H}_{y}=\mathrm{0}.\end{array}\end{array}$

After the substitution of Eq. (27) for Ex and Ez into Eqs. (28), the coupled equations for Ey and Hy can be derived. The following expansion should be used: ${k}_{x}={k}_{\mathrm{0}}+\mathit{\delta }{k}_{x}$, $|\mathit{\delta }{k}_{x}|\ll {k}_{\mathrm{0}}$; also $|{k}_{y}|\ll {k}_{\mathrm{0}}$. Then, according to ,

$\begin{array}{}\text{(29)}& -i\cdot \mathit{\delta }{k}_{x}\to \frac{\partial }{\partial x},-i\cdot {k}_{y}\to \frac{\partial }{\partial y}.\end{array}$

The expansions should be used until the quadratic terms in ky and the linear terms in δkx. As a result, parabolic equations (Levy2000) for the slowly varying amplitudes A and B are derived. In the lower ionosphere and atmosphere, where the effective permittivity reduces to a scalar ε(ω,z), they are independent.

$\begin{array}{}\text{(30a)}& \begin{array}{rl}& \frac{\partial A}{\partial x}+\frac{i}{\mathrm{2}{k}_{\mathrm{0}}}\left(\frac{{\partial }^{\mathrm{2}}A}{\partial {y}^{\mathrm{2}}}+\frac{{\partial }^{\mathrm{2}}A}{\partial {z}^{\mathrm{2}}}\right)+\frac{i{k}_{\mathrm{0}}}{\mathrm{2}}\cdot \left(\mathit{\epsilon }-\mathrm{1}\right)A=\mathrm{0}\\ & \frac{\partial B}{\partial x}+\frac{i}{\mathrm{2}{k}_{\mathrm{0}}}\left(\frac{\mathrm{1}}{\mathit{\beta }}\frac{\partial }{\partial z}\left(\mathit{\beta }\frac{\partial B}{\partial z}\right)+\frac{{\partial }^{\mathrm{2}}B}{\partial {y}^{\mathrm{2}}}\right)+\frac{i{k}_{\mathrm{0}}}{\mathrm{2}}\cdot \left(\mathit{\epsilon }-\mathrm{1}\right)B=\mathrm{0},\end{array}\end{array}$

where $\mathit{\beta }\equiv {\mathit{\epsilon }}^{-\mathrm{1}}$. Accounting for the presence of the gyrotropic layer and the tensor impedance boundary conditions at the upper boundary of z=Lz of the VLF waveguide, the equations for the slowly varying amplitudes in the general case are

$\begin{array}{}\text{(30b)}& \begin{array}{rl}& \frac{\partial A}{\partial x}+\frac{i}{\mathrm{2}{k}_{\mathrm{0}}}\left(\frac{{\partial }^{\mathrm{2}}A}{\partial {y}^{\mathrm{2}}}+\frac{{\partial }^{\mathrm{2}}A}{\partial {z}^{\mathrm{2}}}\right)+\frac{i{k}_{\mathrm{0}}}{\mathrm{2}}\cdot \left({\stackrel{\mathrm{̃}}{\mathit{\epsilon }}}_{\mathrm{22}}-\mathrm{1}\right)\cdot A\\ & \phantom{\rule{0.25em}{0ex}}+\frac{{\mathit{\gamma }}_{\mathrm{21}}}{\mathrm{2}}\frac{\partial B}{\partial z}+\frac{i{k}_{\mathrm{0}}}{\mathrm{2}}\cdot {\mathit{\gamma }}_{\mathrm{23}}B=\mathrm{0}\\ & \phantom{\rule{0.25em}{0ex}}\frac{\partial B}{\partial x}+\frac{i}{\mathrm{2}{k}_{\mathrm{0}}}\left(\frac{\mathrm{1}}{{\stackrel{\mathrm{̃}}{\mathit{\beta }}}_{\mathrm{11}}}\frac{\partial }{\partial z}\left({\stackrel{\mathrm{̃}}{\mathit{\beta }}}_{\mathrm{33}}\frac{\partial B}{\partial z}\right)+\frac{{\partial }^{\mathrm{2}}B}{\partial {y}^{\mathrm{2}}}\right)\\ & \phantom{\rule{0.25em}{0ex}}+\frac{i}{\mathrm{2}{\stackrel{\mathrm{̃}}{\mathit{\beta }}}_{\mathrm{11}}}\frac{\partial }{\partial z}\left({\mathit{\gamma }}_{\mathrm{12}}A\right)+\frac{\mathrm{1}}{\mathrm{2}{\stackrel{\mathrm{̃}}{\mathit{\beta }}}_{\mathrm{11}}}\frac{\partial }{\partial z}\left({\stackrel{\mathrm{̃}}{\mathit{\beta }}}_{\mathrm{13}}B\right)+\frac{i{k}_{\mathrm{0}}}{\mathrm{2}{\stackrel{\mathrm{̃}}{\mathit{\beta }}}_{\mathrm{11}}}{\mathit{\gamma }}_{\mathrm{32}}A+\\ & \phantom{\rule{0.25em}{0ex}}+\frac{{\stackrel{\mathrm{̃}}{\mathit{\beta }}}_{\mathrm{31}}}{\mathrm{2}{\stackrel{\mathrm{̃}}{\mathit{\beta }}}_{\mathrm{11}}}\frac{\partial B}{\partial z}+\frac{i{k}_{\mathrm{0}}}{\mathrm{2}}\cdot \left(\frac{\mathrm{1}}{{\stackrel{\mathrm{̃}}{\mathit{\beta }}}_{\mathrm{11}}}-\mathrm{1}\right)\cdot B=\mathrm{0}.\end{array}\end{array}$

In Eq. (30b),

$\begin{array}{rl}& {\mathit{\gamma }}_{\mathrm{12}}\equiv \frac{{\mathit{\epsilon }}_{\mathrm{13}}\cdot {\mathit{\epsilon }}_{\mathrm{32}}-{\mathit{\epsilon }}_{\mathrm{12}}\cdot {\mathit{\epsilon }}_{\mathrm{33}}}{\mathrm{\Delta }},\\ & {\mathit{\gamma }}_{\mathrm{21}}\equiv \frac{{\mathit{\epsilon }}_{\mathrm{23}}\cdot {\mathit{\epsilon }}_{\mathrm{31}}-{\mathit{\epsilon }}_{\mathrm{21}}\cdot {\mathit{\epsilon }}_{\mathrm{33}}}{\mathrm{\Delta }},\\ & {\mathit{\gamma }}_{\mathrm{23}}\equiv \frac{{\mathit{\epsilon }}_{\mathrm{21}}\cdot {\mathit{\epsilon }}_{\mathrm{13}}-{\mathit{\epsilon }}_{\mathrm{23}}\cdot {\mathit{\epsilon }}_{\mathrm{11}}}{\mathrm{\Delta }},\\ & {\mathit{\gamma }}_{\mathrm{32}}\equiv \frac{{\mathit{\epsilon }}_{\mathrm{31}}\cdot {\mathit{\epsilon }}_{\mathrm{12}}-{\mathit{\epsilon }}_{\mathrm{32}}\cdot {\mathit{\epsilon }}_{\mathrm{11}}}{\mathrm{\Delta }},\\ & {\stackrel{\mathrm{̃}}{\mathit{\beta }}}_{\mathrm{11}}\equiv \frac{{\mathit{\epsilon }}_{\mathrm{11}}}{\mathrm{\Delta }},\\ & {\stackrel{\mathrm{̃}}{\mathit{\beta }}}_{\mathrm{13}}\equiv \frac{{\mathit{\epsilon }}_{\mathrm{13}}}{\mathrm{\Delta }},\\ & {\stackrel{\mathrm{̃}}{\mathit{\beta }}}_{\mathrm{31}}\equiv \frac{{\mathit{\epsilon }}_{\mathrm{31}}}{\mathrm{\Delta }},\\ & {\stackrel{\mathrm{̃}}{\mathit{\beta }}}_{\mathrm{33}}\equiv \frac{{\mathit{\epsilon }}_{\mathrm{33}}}{\mathrm{\Delta }},\\ & \mathrm{\Delta }\equiv {\mathit{\epsilon }}_{\mathrm{11}}\cdot {\mathit{\epsilon }}_{\mathrm{33}}-{\mathit{\epsilon }}_{\mathrm{13}}\cdot {\mathit{\epsilon }}_{\mathrm{31}};\\ & {\stackrel{\mathrm{̃}}{\mathit{\epsilon }}}_{\mathrm{22}}\equiv {\mathit{\epsilon }}_{\mathrm{22}}+\frac{{\mathit{\epsilon }}_{\mathrm{21}}\left({\mathit{\epsilon }}_{\mathrm{13}}{\mathit{\epsilon }}_{\mathrm{32}}-{\mathit{\epsilon }}_{\mathrm{12}}{\mathit{\epsilon }}_{\mathrm{33}}\right)+{\mathit{\epsilon }}_{\mathrm{23}}\left({\mathit{\epsilon }}_{\mathrm{31}}{\mathit{\epsilon }}_{\mathrm{12}}-{\mathit{\epsilon }}_{\mathrm{32}}{\mathit{\epsilon }}_{\mathrm{11}}\right)}{\mathrm{\Delta }}.\end{array}$

Equation (30b) is reduced to Eq. (30a) when the effective permittivity is scalar. At the Earth's surface of z=0, the impedance conditions are reduced, accounting for the medium being isotropic and the conductivity of the Earth being finite, to the form

$\begin{array}{}\text{(31a)}& {E}_{y}={Z}_{\mathrm{0}E}{H}_{x},\phantom{\rule{0.25em}{0ex}}{E}_{x}=-{Z}_{\mathrm{0}E}{H}_{y},\phantom{\rule{0.25em}{0ex}}{Z}_{\mathrm{0}E}\equiv {\left(\frac{i\mathit{\omega }}{\mathrm{4}\mathit{\pi }{\mathit{\sigma }}_{E}}\right)}^{\mathrm{1}/\mathrm{2}},\end{array}$

where ${\mathit{\sigma }}_{E}\sim {\mathrm{10}}^{\mathrm{8}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$ is the Earth's conductivity. The boundary conditions (31a) at the Earth's surface, where ${Z}_{\mathrm{022}}={Z}_{\mathrm{021}}\equiv {Z}_{\mathrm{0}E}$, ${Z}_{\mathrm{012}}={Z}_{\mathrm{021}}=\mathrm{0}$, ${\mathit{\beta }}_{\mathrm{33}}={\mathit{\epsilon }}^{-\mathrm{1}}\left(z=\mathrm{0}\right)$, γ12=0 and ${\stackrel{\mathrm{̃}}{\mathit{\beta }}}_{\mathrm{13}}=\mathrm{0}$ can be rewritten as

$\begin{array}{}\text{(31b)}& \begin{array}{rl}& {E}_{y}+\frac{i}{{k}_{\mathrm{0}}}{Z}_{\mathrm{0}E}\frac{\partial {E}_{y}}{\partial z}=\mathrm{0},\\ & \frac{i}{\mathit{\epsilon }\left(z=\mathrm{0}\right){k}_{\mathrm{0}}}\frac{\partial {H}_{y}}{\partial z}+{Z}_{\mathrm{0}E}{H}_{y}=\mathrm{0}.\end{array}\end{array}$

Equation (30a, 30b), combined with the boundary conditions of Eq. (25) at the upper boundary of the VLF waveguide of z=Lz, and with the boundary conditions at the Earth's surface (Eq. 31b), are used to simulate the VLF wave propagation. The surface impedance of the Earth has been calculated from the Earth's conductivity (see Eq. 31a). The initial conditions to the solution of Eqs. (30a, 30b), (25) and (31b) are chosen in the form

$\begin{array}{}\text{(32)}& \begin{array}{rl}& A\left(x=\mathrm{0},y,z\right)=\mathrm{0},\\ & B\left(x=\mathrm{0},y,z\right)={B}_{\mathrm{0}}\mathrm{exp}\left(-\left(\left(y-\mathrm{0.5}{L}_{y}\right)/{y}_{\mathrm{0}}{\right)}^{\mathrm{2}n}\right)\\ & \phantom{\rule{0.25em}{0ex}}\cdot \mathrm{exp}\left(-\left(\left(z-{z}_{\mathrm{1}}\right)/{z}_{\mathrm{0}}{\right)}^{\mathrm{2}n}\right),\\ & n=\mathrm{2}.\end{array}\end{array}$

In the relations of Eq. (32), z1, z0, y0 and B0 are the vertical position of maximum value, the vertical and transverse characteristic dimensions of the spatial distribution and the maximum value of Hy, respectively, at the input of the system, x=0. The size of the computing region along the OY axis is, by the order of value, Ly∼2000 km. Because the gyrotropic layer is relatively thin and is placed at the upper part of the VLF waveguide, the beams are excited near the Earth's surface. The wave diffraction in this gyrotropic layer along the OY axis is quite small; i.e., the terms ${\partial }^{\mathrm{2}}A/\partial {y}^{\mathrm{2}}$ and ${\partial }^{\mathrm{2}}B/\partial {y}^{\mathrm{2}}$ are small there as well. Contrary to this, the wave diffraction is very important in the atmosphere in the lower part of the VLF waveguide near the Earth's surface. To solve the problem of the beam propagation, the method of splitting with respect to physical factors has been applied . Namely, the problem has been approximated by the finite differences

$\begin{array}{}\text{(33)}& \mathbit{C}\equiv \left(\begin{array}{c}A\\ B\end{array}\right),\phantom{\rule{1em}{0ex}}\frac{\partial \mathbit{C}}{\partial x}+{\stackrel{\mathrm{^}}{L}}_{y}\mathbit{C}+{\stackrel{\mathrm{^}}{L}}_{z}\mathbit{C}=\mathrm{0}.\end{array}$

In the terms of $\stackrel{\mathrm{^}}{{L}_{y}}\mathbit{C}$, the derivatives with respect to y are included, whereas all other terms are included in $\stackrel{\mathrm{^}}{{L}_{z}}\mathbit{C}$. Then the following fractional steps have been applied; the first one is along y, and the second one is along z.

$\begin{array}{}\text{(34)}& \begin{array}{rl}& \frac{{\mathbit{C}}^{p+\mathrm{1}/\mathrm{2}}-{\mathbit{C}}^{p}}{{h}_{x}}+{\stackrel{\mathrm{^}}{L}}_{y}{\mathbit{C}}^{p+\mathrm{1}/\mathrm{2}}=\mathrm{0},\\ & \frac{{\mathbit{C}}^{p+\mathrm{1}}-{\mathbit{C}}^{p+\mathrm{1}/\mathrm{2}}}{{h}_{x}}+{\stackrel{\mathrm{^}}{L}}_{z}{\mathbit{C}}^{p+\mathrm{1}}=\mathrm{0}\end{array}\end{array}$

The region of simulation is $\mathrm{0}–2000 km, $\mathrm{0} 2000–3000 km and $\mathrm{0}–90 km. The numerical scheme of Eq. (34) is absolutely stable. Here hx is the step along the OX axis, and xp=phx and $p=\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{\dots }$. This step has been chosen from the conditions of the simulation result independence of the diminishing hx.

Figure 2The rotation of the local Cartesian coordinate frame at each step along the Earth's surface hx on a small angle $\mathit{\delta }\mathit{\phi }\approx \mathrm{\Delta }x/{R}_{E}$ in radians, while Δx=hx. The following strong inequalities are valid for ${h}_{x}\ll {L}_{z}\ll {R}_{\mathrm{E}}$. At the Earth's surface z=0.

On the simulation at each step along the OX axis, the correction on the Earth's curvature has been inserted in an adiabatic manner applying the rotation of the local coordinate frame XOZ. Because the step along x is small hx and $\sim \mathrm{1}\phantom{\rule{0.125em}{0ex}}\mathrm{km}\ll {L}_{z}$, this correction of the C results in the multiplier $\mathrm{exp}\left(-i{k}_{\mathrm{0}}\cdot \mathit{\delta }x\right)$, where $\mathit{\delta }x=z\cdot \left({h}_{x}/{R}_{\mathrm{E}}\right)$ and RELz is the Earth's radius (see Fig. 2 and the caption to this figure). At the distances of x≤1000 km, the simulation results do not depend on the insertion of this correction, whereas at higher distances a quantitative difference occurs: the VLF beam propagates more closely to the upper boundary of the waveguide.

## 3.5 VLF waveguide modes and reflection from the VLF waveguide upper effective boundary

In general, our model needs the consideration of the waveguide mode excitations by a current source such as a dipole-like VLF radio source and lightning discharge. Then, the reflection of the waves incident on the upper boundary (z=Lz) of the effective WGEI can be considered. There it will be possible to demonstrate that this structure indeed has waveguiding properties that are good enough. Then, in the model described in the present paper, the VLF beam is postulated already on the input of the system. To understand how such a beam is excited by the, say, dipole antenna near the lower boundary of z=0 of the WGEI, the formation of the beam structure based on the mode presentation should be searched. Then the conditions of the radiation (absence of ingoing waves; Eq. 12) can be used as the boundary conditions for the VLF beam radiated to the upper ionosphere and magnetosphere. Due to a relatively large scale of the inhomogeneity in this region, the complex geometrical optics would be quite suitable for modeling beam propagation, even accounting for the wave dispersion in magnetized plasma. The proper effective boundary condition, similar to , would allow for making a relatively accurate match between the regions, described by the full-wave electromagnetic approach with Maxwell's equations and the complex geometrical optics (FWEM-CGO approach). All of these materials are not included in this paper, but they will be delivered in the two future papers. The first paper will be dedicated to the modeling of VLF waves propagating in the WGEI based on the field expansion as a set of eigenmodes of the waveguide (the mode presentation approach). The second paper will deal with the leakage of the VLF beam from the WGEI into the upper ionosphere and magnetosphere and the VLF beam propagation in these media. Here we describe only one result, which concerns the mode excitation in the WGEI, because this result is principally important for the justification of the TIMEB method. It was shown that more than the five lowest modes of the WGEI are strongly localized in the atmosphere or lower ionosphere. Their longitudinal wave numbers are close to the corresponding wave numbers of EM waves in the atmosphere. This fact demonstrates that the TIMEB method can be applied to the propagation of the VLF electromagnetic waves in the WGEI.

4 Modeling results

The dependencies of the permittivity components ε1, ε3 and εh in the coordinate frame associated with the geomagnetic field H0 are given in Fig. 3. The parameters of the ionosphere used for modeling are taken from , , , , and . The typical results of simulations are presented in Fig. 4. The parameters of the ionosphere correspond to Fig. 3. The angle θ (Fig. 1) is 45. The VLF frequency is ω=105s−1 and $f=\mathit{\omega }/\mathrm{2}\mathit{\pi }\approx \mathrm{15.9}$kHz. The Earth's surface is assumed to as ideally conductive at the level of z=0. The values of the EM field are given in absolute units. The magnetic field is measured in oersted (Oe) or gauss (Gs) ($\mathrm{1}\phantom{\rule{0.125em}{0ex}}\mathrm{Gs}={\mathrm{10}}^{-\mathrm{4}}$T), whereas the electric field is also in Gs, 1 Gs=300V cm−1. Note that in the absolute (Gaussian) units of the magnitudes of the magnetic field component $|{H}_{y}|$ are the same as ones of the electric field component $|{E}_{z}|$ in the atmosphere region where the permittivity is ε≈1. The correspondence between the absolute units and practical SI units is given in the Fig. 4 caption.

It is seen that the absolute values of the permittivity components increase sharply above z=75km. The behavior of the permittivity components is step-like, as seen from Fig. 3a. Therefore, the results of simulations are tolerant to the choice of the upper wall position of the Earth's surface–ionosphere waveguide. The computed components of the tensor impedance at z=85km are ${Z}_{\mathrm{011}}=\mathrm{0.087}+i\mathrm{0.097}$, ${Z}_{\mathrm{021}}=\mathrm{0.085}+i\mathrm{0.063}$, ${Z}_{\mathrm{012}}=-\mathrm{0.083}-i\mathrm{0.094}$ and ${Z}_{\mathrm{022}}=\mathrm{0.093}+i\mathrm{0.098}$. So, a condition of $|{Z}_{\mathrm{0}\mathit{\alpha }\mathit{\beta }}|\le \mathrm{0.15}$ is satisfied there, which is necessary for the applicability of the boundary conditions in Eq. (25). The maximum value of the Hy component is $\mathrm{0.1}\phantom{\rule{0.125em}{0ex}}\mathrm{Oe}={\mathrm{10}}^{-\mathrm{5}}$T in Fig. 4a for the initial VLF beam at x=0. This corresponds to the value of the Ez component of $\mathrm{0.1}\phantom{\rule{0.125em}{0ex}}\mathrm{Gs}=\mathrm{30}\phantom{\rule{0.125em}{0ex}}\mathrm{V}\phantom{\rule{0.125em}{0ex}}{\mathrm{cm}}^{-\mathrm{1}}$. At the distance of x=1000km the magnitudes of the magnetic field Hy are of about $\mathrm{3}×{\mathrm{10}}^{-\mathrm{5}}\phantom{\rule{0.125em}{0ex}}\mathrm{Oe}=\mathrm{3}\phantom{\rule{0.125em}{0ex}}\mathrm{nT}$, whereas the electric field Ey is $\mathrm{3}×{\mathrm{10}}^{-\mathrm{6}}\phantom{\rule{0.125em}{0ex}}\mathrm{Gs}\approx \mathrm{1}\phantom{\rule{0.125em}{0ex}}\mathrm{mV}\phantom{\rule{0.125em}{0ex}}{\mathrm{cm}}^{-\mathrm{1}}$.

Figure 3(a) The vertical dependencies of the modules of components of the permittivity in the frame associated with the geomagnetic field $|{\mathit{\epsilon }}_{\mathrm{1}}|$, $|{\mathit{\epsilon }}_{\mathrm{3}}|$ and $|{\mathit{\epsilon }}_{h}|$ with curves 1, 2 and 3, correspondingly. (b)(g) The real (corresponding lines with the values denoted by one prime) and imaginary parts (corresponding lines with the values denoted by two primes) of the components ε1, ε3 and εh in general and detailed views.

Figure 4Panel (a) is the initial distribution of $|{H}_{y}|$ at x=0. Panels (b) and (c) are $|{E}_{y}|$ and $|{H}_{y}|$ at x=600 km. Panels (d) and (e) are $|{E}_{y}|$ and $|{H}_{y}|$ at x=1000 km. For the electric field, the maximum value (d) is $\mathrm{3}×{\mathrm{10}}^{-\mathrm{6}}\phantom{\rule{0.125em}{0ex}}\mathrm{Gs}\approx \mathrm{1}\phantom{\rule{0.125em}{0ex}}\mathrm{mV}\phantom{\rule{0.125em}{0ex}}{\mathrm{cm}}^{-\mathrm{1}}$; for the magnetic field, the maximum value (e) is $\mathrm{3}×{\mathrm{10}}^{-\mathrm{5}}\phantom{\rule{0.125em}{0ex}}\mathrm{Gs}\approx \mathrm{3}\phantom{\rule{0.125em}{0ex}}\mathrm{nT}$. At the altitudes $z<\mathrm{75}\phantom{\rule{0.125em}{0ex}}\mathrm{km},|{E}_{z}|\approx |{H}_{y}|$. $\mathit{\omega }=\mathrm{1.0}×{\mathrm{10}}^{\mathrm{5}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$ and θ=45.

The wave beams are localized within the WGEI at $\mathrm{0}, mainly in the regions with the isotropic permittivity (see Fig. 4b–e). The mutual transformations of the beams of different polarizations occur near the waveguide upper boundary due to the anisotropy of the ionosphere within the thin layer at $\mathrm{75}\phantom{\rule{0.125em}{0ex}}\mathrm{km} (Fig. 4b, d). These transformations depend on the permittivity component values of the ionosphere at the altitude of z>80 km and on the components of the tensor impedance. Therefore, the measurements of the phase and amplitude modulations of different EM components near the Earth's surface can provide information on the properties of the lower and middle ionosphere.

In accordance with boundary condition in Eq. (32), we suppose that when entering the system at X=0, only one of the two polarization modes is excited, namely, the transverse magnetic (TM) mode, i.e., at x=0, Hy≠0 and Ey=0 (Fig. 4a). Upon further propagation of the beam with such boundary conditions at the entrance to the system in a homogeneous isotropic waveguide, the property of the electromagnetic field described by the relation Hy≠0 and Ey=0 will remain valid. The qualitative effect due to the presence of gyrotropy (a) in a thin bulk layer near the upper boundary of WGEI and (b) in the upper boundary condition with complex gyrotropic and anisotropic impedance is as follows: during beam propagation in the WGEI, the transverse electric (TE) polarization mode with the corresponding field components, including Ey, is also excited. This effect is illustrated in Fig. 4b and d.

The magnitude of the Ey component depends on the values of the electron concentration at the altitudes z=75–100 km. In Fig. 5a and b the different dependencies of the electron concentration n(z) are shown (see the solid, dashed and dotted lines 1, 2 and 3, respectively). The corresponding dependencies of the component absolute values of the permittivity are shown in Fig. 4c and d.

Figure 5Different profiles of the electron concentration n(z) used in simulations: solid, dashed and dotted lines correspond to undisturbed, decreased and increased concentrations, respectively. Shown are the (a) detailed view and (b) general view. Panels (c) and (d) show the permittivity $|{\mathit{\epsilon }}_{\mathrm{3}}|$ and εh modules.

The distributions of $|{E}_{y}|$ and $|{H}_{y}|$ at x=1000 km are given in Fig. 6. Results in Fig. 6a and b correspond to the solid (1) curve n(z) in Fig. 5; Fig. 6c and d correspond to the dashed (2) curve; Fig. 6e and f correspond to the dotted (3) curve in Fig. 5. The initial Hy beam is the same and is given in Fig. 4a. The values of the tensor impedance for these three cases are presented in Table 1.

Table 1The values of tensor impedance components corresponding to the data shown in Fig. 5.

Figure 6Panels (a), (c) and (e) are dependencies of $|{E}_{y}|$. Panels (b), (d) and (f) are dependencies of $|{H}_{y}|$ at x=1000km, with $\mathit{\omega }=\mathrm{1.0}×{\mathrm{10}}^{\mathrm{5}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$ and θ=45. The initial beams are the same as in Fig. 4a. Panels (a) and (b) correspond to the solid (1) curves in Fig. 5. Panels (c) and (d) are for the dashed (2) curves. Panels (e) and (f) correspond to the dotted (3) curves there.

Figure 7The dependencies of EM components at altitude z in the center of the waveguide at y=1500km for the different profiles of the electron concentration. The solid (1), dashed (2) and dotted (3) curves correspond to the different profiles of the electron concentration in Fig. 5. Panels (a) and (b) are the same kinds of curves, with $\mathit{\omega }=\mathrm{1.0}×{\mathrm{10}}^{\mathrm{5}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$ and θ=45.

The distributions of $|{E}_{y}|$ and $|{H}_{y}|$ on z at x=1000km in the center of the waveguide, y=1500km, are given in Fig. 7. These simulations show that the change in the complex tensors of both volume dielectric permittivity and impedance at the lower and upper boundaries of the effective WGEI influence the VLF losses remarkably. The modulation of the electron concentration at the altitudes above z=120km affects the excitation of the Ey component within the waveguide rather weakly.

5 Influence of the parameters of the WGEI on the polarization transformation and losses of the propagating VLF waves

An important effect of the gyrotropy and anisotropy is the corresponding transformation of the field polarization during the propagation in the WGEI, which is absent in the ideal metal planar waveguide without gyrotropy and anisotropy. We will show that this effect is quite sensitive to the carrier frequency of the beam, propagating in the WGEI, and the inclination of the geomagnetic field and perturbations in the electron concentration, which can vary under the influences of the powerful sources placed below, above and inside the ionosphere. In the real WGEI, the anisotropy and gyrotropy are connected with the volume effect and effective surface tensor impedances at the lower and upper surfaces of the effective WGEI, where z=0 and z=Lz (Fig. 1). For the corresponding transformation of the field polarization determination, we introduce the characteristic polarization relation $|{E}_{y}/{H}_{y}|\left(z;y={L}_{y}/\mathrm{2};x={x}_{\mathrm{0}}\right)$, taken at the central plane of the beam (y=Ly) at the characteristic distance (x=x0) from the beam input and VLF transmitter. The choice of the characteristic polarization parameter ($|{E}_{y}/{H}_{y}|$) and its dependence on the vertical coordinate (z) is justified by conditions (1)–(3). (1) The WGEI is similar to the ideal planar metalized waveguide because, first, the tensor $\stackrel{\mathrm{^}}{\mathit{ϵ}}$ is different from the isotropic $\stackrel{\mathrm{^}}{I}$ only in the relatively small upper part of the WGEI in the altitude range from 75–80 to 85 km (see Fig. 1). Second, both the Earth and ionosphere conductivity are quite high, and corresponding impedances are quite low. In particular, the elements of the effective tensor impedance at the upper boundary of the WGEI are small, $|{Z}_{\mathrm{0}\mathit{\alpha }\mathit{\beta }}|\le \mathrm{0.1}$ (see, for example, Table 1). (2) Respectively, the carrier modes of the VLF beam are close to the modes of the ideal metalized planar waveguide. These uncoupled modes are subdivided into sets of (${E}_{x},{H}_{y},{E}_{z}$) and (${H}_{x},{E}_{y},{H}_{z}$). The detailed search of the propagation of the separate eigenmodes of the WGEI is not a goal of this paper and will be the subject of the separate paper. (3) Because we have adopted the input boundary conditions in Eq. (32) (with Hy≠0, Ey=0) for the initial beam(s), the above-mentioned value $|{E}_{y}/{H}_{y}|\left(z;y={L}_{y}/\mathrm{2};x={x}_{\mathrm{0}}\right)$ characterizes the mode coupling and corresponding transformation of the polarization at the distance of x0 from the beam input due to the presence of the volume and surface gyrotropy and anisotropy in the real WGEI. The results presented below are obtained for x0=1000km, that is, by the order of value, a typical distance, for example, between the VLF transmitter and receiver of the European VLF/LF radio network . Another parameter characterizing the propagation of the beam in the WGEI is the effective total loss parameter of $|{H}_{y\mathrm{max}}\left(x={x}_{\mathrm{0}}\right)/{H}_{y\mathrm{max}}\left(x=\mathrm{0}\right)|$. Note that this parameter characterizes both dissipative and diffraction losses. The last are connected with beam spreading in the transverse (y) direction during the propagation in the WGEI.

In Fig. 8 the polarization and loss characteristic dependencies on both the carrier beam frequency and the angle θ between the geomagnetic field and the vertical directions (see Fig. 1) are shown.

Figure 8Characteristics of the polarization transformation parameter $|{E}_{y}/{H}_{y}|$ (a–c) and the effective coefficient of the total losses at the distance of x0=1000km from the beam input (d). Corresponding altitude dependence of the electron concentration is shown in line 1 of Fig. 5b. Panels (a), (b) and (d) show dependences of the polarization parameter (a, b) and total losses (d) on the vertical coordinate and frequency for different θ angles, respectively. Black (1), red (2), green (3) and blue (4) curves in panels (a), (b) and (d) correspond to 5, 30, 45 and 60, respectively. Panels (a) and (b) correspond to the frequencies $\mathit{\omega }=\mathrm{0.86}×{\mathrm{10}}^{\mathrm{5}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$ and $\mathit{\omega }=\mathrm{1.14}×{\mathrm{10}}^{\mathrm{5}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$, respectively. Panel (c) shows the dependence of the polarization parameter on the vertical coordinate for the different frequencies. Black (1), red (2) and green (3) lines correspond to the frequencies 0.86×105, 1×105 and $\mathrm{1.14}×{\mathrm{10}}^{\mathrm{5}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$, respectively, and θ=45.

In Fig. 8a–c the altitude dependence of the polarization parameter $|{E}_{y}/{H}_{y}|$ exhibits two main maxima in the WGEI. The first one lies in the gyrotropic region above 70 km, while the second one in the isotropic region of the WGEI. As seen from Fig. 8a and b, the value of the (larger) second maximum increases, while the value of the first maximum decreases, and its position shifts to the lower altitudes with increasing frequency. At the higher frequency ($\mathit{\omega }=\mathrm{1.14}×{\mathrm{10}}^{\mathrm{5}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$), the larger maximum of the polarization parameter corresponds to the intermediate value of the angle θ=45 (Fig. 8b); for the lower frequency ($\mathit{\omega }=\mathrm{0.86}×{\mathrm{10}}^{\mathrm{5}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$), the largest value of the first (higher) maximum corresponds to the almost vertical direction of the geomagnetic field (θ=5; Fig. 8a). For the intermediate value of the angle (θ=45), the largest value of the main maximum corresponds to the higher frequency ($\mathit{\omega }=\mathrm{1.14}×{\mathrm{10}}^{\mathrm{5}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$) in the considered frequency range (Fig. 8c). The total losses increase monotonically with increasing frequency and depend weakly on the value of θ (Fig. 8d).

To model the effect of increasing and decreasing the electron concentration ne in the lower ionosphere on the polarization parameter, we have used the following parameterization for the ne change $\mathrm{\Delta }{n}_{\mathrm{e}}={n}_{\mathrm{e}}\left(z\right)-{n}_{\mathrm{0}\mathrm{e}}\left(z\right)$ of the electron concentration, where n0e(z) is the unperturbed altitude distribution of the electron concentration.

$\begin{array}{}\text{(35)}& \begin{array}{rl}& \mathrm{\Delta }{n}_{\mathrm{e}}\left(z\right)={n}_{\mathrm{0}\mathrm{e}}\left(z\right)\mathrm{\Phi }\left(z\right);\\ & \mathrm{\Phi }\left(z\right)=\left[F\left(z\right)\right]-\frac{\left(z-{z}_{\mathrm{2}}{\right)}^{\mathrm{2}}}{\mathrm{\Delta }{z}_{\mathrm{12}}^{\mathrm{2}}}\left[F\left({z}_{\mathrm{1}}\right)\right]-\frac{\left(z-{z}_{\mathrm{1}}{\right)}^{\mathrm{2}}}{\mathrm{\Delta }{z}_{\mathrm{12}}^{\mathrm{2}}}\left[F\left({z}_{\mathrm{2}}\right)\right];\\ & F\left(z\right)={f}_{\mathrm{00}}\cdot {\mathrm{cosh}}^{-\mathrm{2}}\left\{\left[z-\left(\frac{{z}_{\mathrm{1}}+{z}_{\mathrm{2}}}{\mathrm{2}}\right)\right]/\mathrm{\Delta }z\right\}\end{array}\end{array}$

In Eq. (35), $\mathrm{\Delta }{z}_{\mathrm{12}}\equiv {z}_{\mathrm{2}}-{z}_{\mathrm{1}}$; Δz is the effective width of the electron concentration perturbation altitude distribution. The perturbation Δne is concentrated in the range of altitudes ${z}_{\mathrm{1}}\le z\le {z}_{\mathrm{2}}$ and is equal to zero outside this region, $\mathrm{\Delta }{n}_{\mathrm{e}}\left({z}_{\mathrm{1}}\right)=\mathrm{\Delta }{n}_{\mathrm{e}}\left({z}_{\mathrm{2}}\right)=\mathrm{0}$, while $\mathrm{\Phi }\left({z}_{\mathrm{1}}\right)=\mathrm{\Phi }\left({z}_{\mathrm{2}}\right)=\mathrm{0}$.

Figure 9(a) Decreased and increased electron concentration (line 2, red) and (line 3, blue) correspond to ${f}_{\mathrm{00}}=-\mathrm{1.25}$ and f00=250, respectively, relative to the reference concentration (line 1, black) with parametrization conditions (see Eq. 35) z1=50km, z2=90km and Δ=20km. Panels (b) and (c) are the polarization parameter $|{E}_{y}/{H}_{y}|$ altitude distribution for decreased and increased electron concentration, respectively. In (b) and (c) lines 1, 2 and 3 correspond to ω values $\mathrm{0.86}×{\mathrm{10}}^{\mathrm{5}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$, $\mathrm{1.0}×{\mathrm{10}}^{\mathrm{5}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$ and $\mathrm{1.14}×{\mathrm{10}}^{\mathrm{5}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$, respectively. Angle θ is equal to 45.

The change in the concentration in the lower ionosphere causes a rather nontrivial effect on the parameter of the polarization transformation $|{E}_{y}/{H}_{y}|$ (Fig. 9a–c). Note that either an increase or a decrease in the ionosphere plasma concentration have been reported as a result of seismogenic phenomena, tsunamis, particle precipitation in the ionosphere due to wave–particle interaction in the radiation belts , etc. The changes in the $|{E}_{y}/{H}_{y}|$ due to an increase or decrease in electron concentration vary by absolute values from dozens to thousands of a percent. This can be seen from the comparison between Fig. 9b and c (lines 3) and Fig. 8c (line 3), which corresponds to the unperturbed distribution of the ionospheric electron concentration (see also lines 1 in Figs. 5b and 9a). It is even more interesting that in the case of decreasing (Fig. 9a, curve 2) electron concentration, the main maximum of $|{E}_{y}/{H}_{y}|$ appears in the lower atmosphere (at the altitude around 20 km, Fig. 9b, curve 3, which corresponds to $\mathit{\omega }=\mathrm{1.14}×{\mathrm{10}}^{\mathrm{5}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$). In the case of increasing electron concentration (Fig. 9a, curve 3), the main maximum of $|{E}_{y}/{H}_{y}|$ appears near the E region of the ionosphere (at the altitude around 77 km, Fig. 9c). The secondary maximum, which is placed, in the absence of the perturbation of the electron concentration, in the lower atmosphere (Fig. 8c, curves 2 and 3) or mesosphere and ionosphere D region ((Fig. 8c, curve 1) practically disappears or just is not seen in the present scale in the case under consideration (Fig. 9c, curves 1–3).

Figure 10(a, b) Frequency dependencies of the real (a) and imaginary (b) parts of the effective tensor impedance Z011 component at the upper boundary (z=Lz, see Fig. 1) of the WGEI. Lines 1 (black), 2 (red), 3 (blue) and 4 (green) correspond to θ=5, 30, 45 and 60, respectively. Panel (c) shows the polarization parameter $|{E}_{y}/{H}_{y}|$ altitude dependency at the frequency $\mathrm{0.86}×{\mathrm{10}}^{\mathrm{5}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$ and angle θ=45 for the isotropic surface impedance Z0 at the lower surface of the WGEI equal to 10−4. Earth conductivity σ is equal to 109 s−1 for line 1 and ${Z}_{\mathrm{0}}={\mathrm{10}}^{-\mathrm{2}}$ ($\mathit{\sigma }={\mathrm{10}}^{\mathrm{7}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$) for line 2.

In Fig. 10, the real (a) and imaginary (b) parts of the surface impedance at the upper boundary of the WGEI have a quasi-periodical character with the amplitude of oscillations occurring around the effective average values (not shown explicitly in Fig. 10a and b), which decreases with an increasing angle of θ. The average Re(Z011) and Im(Z011) values in general decrease with an increasing angle of θ (see Fig. 10a and b). The average values of Re(Z011) at θ=5, 30, 45 and 60 (lines 1–4 in Fig. 10a) and Im(Z011) at θ=45 and θ=60 (curves 3 and 4 in Fig. 10b) increase with an increasing frequency in the frequency range (0.86–1.14) $×{\mathrm{10}}^{\mathrm{5}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$. The average Im(Z011) value at θ=5 and 30 changes in the frequency range (0.86–1.14) $×{\mathrm{10}}^{\mathrm{5}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$ non-monotonically with the maximum at (1–1.1) $×{\mathrm{10}}^{\mathrm{5}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$. The value of finite impedance at the lower Earth–atmosphere boundary of the WGEI influences the polarization transformation parameter minimum near the E region of the ionosphere (lines 1 and 2 in Fig. 10c). The decrease of surface impedance Z0 at the lower Earth–atmosphere boundary of the WGEI by two orders of magnitude produces the 100 % increase of the corresponding $|{E}_{y}/{H}_{y}|$ minimum at z∼75km (Fig. 10c).

6 Discussion

The observations presented in show a possibility for seismogenic increasing losses of VLF waves in the WGEI (Fig. 11; see details in ). We discuss the qualitative correspondence of our results to these experimental data.

Figure 11Averaged residual VLF and LF signals from ground-based observations at the wave paths of JJY-Moshiri, JJI-Kamchatka, JJY-Kamchatka, Australia-Kamchatka and Hawaii-Kamchatka. Horizontal dotted lines show the 2σ level. The color-filled zones highlight values exceeding the −2σ level. In panel b Dst (disturbance storm time index) variations and earthquakes magnitude values are shown (from , see their Fig. 1 but not including the Detection of Electro-Magnetic Emissions Transmitted from Earthquake Regions (DEMETER) data; the work of , is licensed under a Creative Commons Attribution 4.0 International License (CC BY 4.0)). See other details in .

Figure 12Modification of the ionosphere by the electric field of seismogenic origin based on the theoretical model . (a) Geometry of the model for the determination of the electric field excited by the seismogenic current source Jz(x,y) and penetrated into the ionosphere with isotropic (I) and anisotropic (II) regions of the “atmosphere–ionosphere” system. (b) Electric field in the mesosphere in the presence of the seismogenic current sources only in the mesosphere (1), in the lower atmosphere (2), and both in the mesosphere and in the lower atmosphere (3). (c) Relative perturbations caused by the seismogenic electric field, normalized on the corresponding steady-state values in the absence of the perturbing electric field, denoted by the index “0”, electron temperature (TeTe0), electron concentration (NeNe0) and electron collision frequency (νeνe0).

Figure 13Altitude distributions of the normalized tangential (y) electric (a) and magnetic (b) VLF beam field components in the central plane of the transverse beam distribution (y=0) at the distance of x=1000km from the input of the system. Line 1 in (a, b) corresponds to the presence of the mesospheric electric current source only with a relatively small value of Ne and a large νe. Line 2 corresponds to the presence of both mesospheric and near-ground seismogenic electric current sources with a relatively large value of Ne and small νe. Lines 1 and 2 in (a) and (b) correspond qualitatively to the lines 1 and 3, respectively, in Fig. 12b, with $\mathit{\omega }=\mathrm{1.5}×{\mathrm{10}}^{\mathrm{5}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$ and θ=45.

The modification of the ionosphere due to electric field excited by the near-ground seismogenic current source has been taken into account. In the model , the presence of the mesospheric current source, which followed from the observations is also taken into account. It is assumed that the mesospheric current only has the z component and is positive, which means that it is directed vertically downward, as in the fair-weather current (curve 1, Fig. 12b). Then suppose that the surface seismogenic current is directed in the same way as the mesospheric current. We first consider the case when the mesospheric current is zero and only the corresponding seismogenic current is present near the Earth. The corresponding mesospheric electric field under the condition of a given potential difference between the Earth and the ionosphere (curve 2, Fig. 12b) is directed opposite to those excited by the corresponding mesospheric current (curve 1, Fig. 12b). As a result, in the presence of both mesospheric and a seismogenic surface current, the total mesospheric electric field (curve 3, Fig. 12b) is smaller in absolute value than in the presence of only a mesospheric current (curve 1, Fig. 12b). It has been shown by that the decrement of losses $|{k}^{\prime \prime }|$ for VLF waves in the WGEI is proportional to $|{k}^{\prime \prime }|\sim |{\mathit{\epsilon }}^{\prime \prime }|\sim {N}_{\mathrm{e}}/{\mathit{\nu }}_{\mathrm{e}}$. Ne and νe decrease and increase, respectively, due to the appearance of a seismogenic surface electric current, in addition to the mesospheric current (curve 3, Fig. 12b). As a result, the losses increase compared with the case when the seismogenic current is absent and the electric field has a larger absolute value (curve 1, Fig. 12). The increase in losses in the VLF beam, shown in Fig. 13 (compare curves 2 and 1 in Fig. 13a and b), corresponds to an increase in losses with an increase in the absolute value of the imaginary part of the dielectric constant when a near-surface seismogenic current source appears (curve 3 in Fig. 12b) in addition to the existing mesospheric current source (curve 2 in Fig. 12b). This seismogenic increase in losses corresponds qualitatively to the results presented in .

The TIMEB is a new method of modeling characteristics of the WGEI. The results of beam propagation in WGEI modeling presented above include the range of altitudes inside the WGEI (see Figs. 47). Nevertheless, the TIMEB method described by Eqs. (15)–(19), (22)–(24), (27), (30a) and (30b) and allows for the determination of all field components in the range of altitudes of $\mathrm{0}\le z\le {L}_{\mathrm{max}}$, where Lmax=300km. The structure and behavior of these eigenmodes in the WGEI and leakage waves will be a subject of separate papers. We present here only the final qualitative result of the simulations. In the range ${L}_{z}\le z\le {L}_{\mathrm{max}}$, where Lz=85km is the upper boundary of the effective WGEI, all field components (1) are at least one order of altitude less than the corresponding maximal field value in the WGEI and (2) have the oscillating character along the z coordinate and describe the modes leaking from the WGEI.

Let us make a note also on the dependences of the field components in the WGEI on the vertical coordinate (z). The initial distribution of the electromagnetic field with altitude z (Fig. 4a) is determined by the boundary conditions of the beam (see Eq. 32). The field component includes higher eigenmodes of the WGEI. The higher-order modes experienced quite large losses and practically disappear after beam propagation on a 1000 km distance. This determines the change in altitude (z) and transverse (y) distributions of the beam field during propagation along the WGEI. In particular, at the distance of x=600km from the beam input (Figs. 4b, c), the few lowest modes of the WGEI along the z and y coordinates still persist. At distance x=1000km (Figs. 4d, c; 6e, f; and 7a, b), only the main mode persists in the z direction. Note that the described field structure correspond to the real WGEI with losses. The gyrotropy and anisotropy causes the volume effects and surface impedance, in distinction to the ideal planar metalized waveguide with isotropic filling (Collin2001).

The closest approach of the direct investigation of the VLF electromagnetic field profile in the Earth–ionosphere waveguide was a series of sounding rocket campaigns at mid and high latitudes at Wallops Island, Virginia, USA, and Siple Station in Antarctica . Here single-axis E field and three-axis B field antennas, supplemented in some cases with in situ plasma density measurements, were used to detect the far-field fixed-frequency VLF signals radiated by US Navy and Stanford ground transmitters.

The most comprehensive study of the WGEI will be provided by the ongoing NASA VIPER (VLF Trans-Ionospheric Propagation Experiment Rocket) project (PI John Bonnell, University of California, Berkeley; NASA grant no. 80NSSC18K0782). The VIPER sounding rocket campaign consists of a summer nighttime launch during quiet magnetosphere conditions from Wallops Flight Facility, Virginia, USA, collecting data through the D, E and F regions of the ionosphere with a payload carrying the following instrumentation: 2-D E and 3-D H field waveforms (DC-1 kHz); 3-D waveforms ranging from an ELF (extremely low frequency) to VLF (100 Hz to 50 kHz); 1-D wideband E field measurements of plasma and upper hybrid lines (100 kHz to 4 MHz); and Langmuir probe plasma density and ion gauge neutral-density measurements at a sampling rate of at least tens of Hz. The VIPER project will fly a fully 3-D EM field measurement, direct current (DC) through VLF, and relevant plasma and neutral particle measurements at mid latitudes through the radiation fields of (1) an existing VLF transmitter (the very low-frequency shore radio station with the call sign NAA at Cutler, Maine, USA, which transmits at a frequency of 24 kHz and an input power of up to 1.8 MW; see Fig. 11) and (2) naturally occurring lightning transients through and above the leaky upper boundary of the WGEI. This is supported by a vigorous theory and modeling effort in order to explore the vertical and horizontal profile of the observed 3-D electric and magnetic radiated fields of the VLF transmitter and the profile related to the observed plasma and neutral densities. The VLF wave's reflection, absorption and transmission processes as a function of altitude will be searched making use of the data on the vertical VLF E and H field profile.

Figure 14Proposed VIPER trajectory.

The aim of this experiment will be the investigation of the VLF beams launched by the near-ground source and VLF transmitter with the known parameters and propagating both in the WGEI and leaking from WGEI into the upper ionosphere. Characteristics of these beams will be compared with the theory proposed in the present paper and the theory on leakage of the VLF beams from WGEI, which we will present in the next papers.

7 Conclusions
1. We have developed the new and highly effective robust method of tensor impedance for the VLF electromagnetic beam propagation in the inhomogeneous waveguiding media: the “tensor impedance method for modeling propagation of electromagnetic beams” in a multi-layered and inhomogeneous waveguide. The main differences and advantages of the proposed tensor impedance method in comparison with the known method of impedance recalculating, in particular invariant imbedding methods , are the following: (i) our method is a direct method of the recalculation of tensor impedance, and the corresponding tensor impedance is determined analytically (see Eq. 22); (ii) our method applied for the media without non-locality does not need a solution for integral equation(s), as in the invariant imbedding method; and (iii) the proposed tensor impedance method does not need to reveal the forward and reflected waves. Moreover, even the conditions of the radiation in Eq. (12) at the upper boundary z=Lmax are determined through the total field components Hx,y that makes the proposed procedure technically less cumbersome and practically more convenient.

2. The waveguide includes the region for the altitudes $\mathrm{0}–90 km. The boundary conditions are the radiation conditions at z=300km, which can be recalculated to the lower altitudes as the tensor relations between the tangential components of the EM field. In another words, the tensor impedance conditions have been used at z=80–90 km.

3. The application of this method jointly with the previous results of the modification of the ionosphere by seismogenic electric field gives results which qualitatively are in agreement with the experimental data on the seismogenic increasing losses of VLF wave/beam propagation in the WGEI.

4. The observable qualitative effect is the mutual transformation of different polarizations of the electromagnetic field occur during the propagation. This transformation of the polarization depends on the electron concentration, i.e., the conductivity, of the D and E layers of the ionosphere at altitudes of 75–120 km.

5. Changes in complex tensors of both volume dielectric permittivity and impedances at the lower and upper boundaries of effective the WGEI influence the VLF losses in the WGEI remarkably.

6. An influence is demonstrated on the parameters characterizing the propagation of the VLF beam in the WGEI, in particular, the parameter of the transformation polarization $|{E}_{y}/{H}_{y}|$ and tensor impedance at the upper boundary of the effective WGEI, the carrier beam frequency, and the inclination of the geomagnetic field and the perturbations in the altitude distribution of the electron concentration in the lower ionosphere.

• (i)

The altitude dependence of the polarization parameter $|{E}_{y}/{H}_{y}|$ has two main maxima in the WGEI: the higher maximum is in the gyrotropic region above 70 km, while the other is in the isotropic region of the WGEI. The value of the (larger) second maximum increases, while the value of the first maximum decreases, and its position shifts to the lower altitudes with increasing frequency. In the frequency range of ω= (0.86–1.14) ×105s−1, at the higher frequency, the larger maximum polarization parameter corresponds to the intermediate value of the angle θ=45; for the lower frequency, the largest value of the first (higher) maximum corresponds to the nearly vertical direction of the geomagnetic field. The total losses increase monotonically with increasing frequency and depend weakly on the value of θ (Fig. 1).

• (ii)

The change in the concentration in the lower ionosphere causes a rather nontrivial effect on the parameter of the polarization transformation $|{E}_{y}/{H}_{y}|$. This effect does include the increase and decrease of the maximum value of the polarization transformation parameter $|{E}_{y}/{H}_{y}|$. The corresponding change of this parameter has large values from dozens to thousands of a percent. In the case of a decreasing electron concentration, the main maximum of $|{E}_{y}/{H}_{y}|$ appears in the lower atmosphere at an altitude of around 20 km. In the case of an increasing electron concentration, the main maximum of $|{E}_{y}/{H}_{y}|$ appears near the E region of the ionosphere (at the altitude around 77 km), while the secondary maximum practically disappears.

• (iii)

The real and imaginary parts of the surface impedance at the upper boundary of the WGEI have a quasi-periodical character with the amplitude of “oscillations” occurring around some effective average value decreases with an increasing angle of θ. Corresponding average values of Re(Z11) and Im(Z11), in general, decrease with an increasing angle of θ. Average values of Re(Z11) for θ equal to 5, 30, 45 and 60 and Im(Z11) correspond to θ equal to 45 and 60, and these increase with an increasing frequency in the considered frequency range of 0.86–1.14 $×{\mathrm{10}}^{\mathrm{5}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$. The average value of Im(Z11) corresponds to θ equal to 5 and 30 and changes in the frequency range of 0.86–1.14 $×{\mathrm{10}}^{\mathrm{5}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$ non-monotonically while having maximum values around the frequency of 1–1.1 $×{\mathrm{10}}^{\mathrm{5}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$.

• (iv)

The value of finite impedance at the lower Earth–atmosphere boundary of the WGEI quite observably influences the polarization transformation parameter minimum near the E region of the ionosphere. The decrease of surface impedance Z at the lower Earth–atmosphere boundary of the WGEI in two orders causes the increase of the corresponding minimum value of $|{E}_{y}/{H}_{y}|$ in ∼100 %.

7. In the range ${L}_{z}\le z\le {L}_{\mathrm{max}}$, where Lz=85 km is the upper boundary of the effective WGEI, all field components (a) are at least one order of altitude less than the corresponding maximal value in the WGEI and (b) have the oscillating character (along the z coordinate), which describes the modes leaking from the WGEI. The detail consideration of the electromagnetic waves leaking from the WGEI will be presented in a separate paper. The initial distribution of the electromagnetic field with z (vertical direction) is determined by the initial conditions on the beam. This field includes higher eigenmodes of the WGEI. The higher-order modes, in distinction to the lower ones, have quite large losses and practically disappear after a beam propagation for 1000 km distance. This circumstance determines the change in altitude (z) distribution of the field of the beam during its propagation along the WGEI. In particular, at the distance of x=600km from the beam input, the few lowest modes of the WGEI along z coordinates still survived. Further, at x=1000km, practically, only the main mode in the z direction remains. This fact is reflected in a minimum number of oscillations of the beam field components along z at a given value of x.

8. The proposed propagation of VLF electromagnetic beams in the WGEI model and results will be useful to explore the characteristics of these waves as an effective instrument for diagnostics of the influences on the ionosphere from above in the Sun–solar-wind–magnetosphere–ionosphere system; from below from the most powerful meteorological, seismogenic and other sources in the lower atmosphere and lithosphere and Earth, such as hurricanes, earthquakes and tsunamis; from inside the ionosphere by strong thunderstorms with lightning discharges; and even from far space by gamma flashes and cosmic rays events.

Appendix A: The matrix coefficients included in Eq. (16)

Here the expressions of the matrix coefficients are presented that are used in the matrix factorization to compute the tensor impedance (see Eq. 16).

$\begin{array}{}\text{(A1)}& \begin{array}{rl}& {\stackrel{\mathrm{^}}{\mathit{\alpha }}}_{N}^{\left(\mathrm{0}\right)}=\left(\begin{array}{cc}\mathrm{1}+\frac{i{h}_{z}}{\mathrm{\Delta }}\left({k}_{\mathrm{1}}-{\mathit{\alpha }}_{\mathrm{1}}{\mathit{\alpha }}_{\mathrm{2}}{k}_{\mathrm{2}}\right);& \frac{i{h}_{z}}{\mathrm{\Delta }}{\mathit{\alpha }}_{\mathrm{2}}\left({k}_{\mathrm{2}}-{k}_{\mathrm{1}}\right)\\ \frac{i{h}_{z}}{\mathrm{\Delta }}{\mathit{\alpha }}_{\mathrm{1}}\left({k}_{\mathrm{1}}-{k}_{\mathrm{2}}\right);& \mathrm{1}+\frac{i{h}_{z}}{\mathrm{\Delta }}\left({k}_{\mathrm{2}}-{\mathit{\alpha }}_{\mathrm{1}}{\mathit{\alpha }}_{\mathrm{2}}{k}_{\mathrm{1}}\right)\end{array}\right),\phantom{\rule{1em}{0ex}}{\stackrel{\mathrm{^}}{\mathit{\alpha }}}_{N}^{\left(-\right)}=\left(\begin{array}{cc}-\mathrm{1};& \mathrm{0}\\ \mathrm{0};& -\mathrm{1}\end{array}\right);\phantom{\rule{1em}{0ex}}\mathrm{\Delta }\equiv \mathrm{1}-{\mathit{\alpha }}_{\mathrm{1}}{\mathit{\alpha }}_{\mathrm{2}};\\ & {\stackrel{\mathrm{^}}{\mathit{\alpha }}}_{j}^{\left(-\right)}=\left(\begin{array}{cc}\left(\frac{{\mathit{\beta }}_{\mathrm{22}}}{\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}}{\right)}_{j-\mathrm{1}/\mathrm{2}};& -\left(\frac{{\mathit{\beta }}_{\mathrm{21}}}{\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}}{\right)}_{j-\mathrm{1}/\mathrm{2}}+\frac{i{k}_{x}{h}_{z}}{\mathrm{2}}\left(\frac{{\mathit{\beta }}_{\mathrm{23}}}{\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}}{\right)}_{j-\mathrm{1}}\\ -\left(\frac{{\mathit{\beta }}_{\mathrm{12}}}{\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}}{\right)}_{j-\mathrm{1}/\mathrm{2}}+\frac{i{k}_{x}{h}_{z}}{\mathrm{2}}\left(\frac{{\mathit{\beta }}_{\mathrm{32}}}{\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}}{\right)}_{j-\mathrm{1}};& \left({\mathit{\beta }}_{\mathrm{11}}+\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}\frac{{\mathit{\beta }}_{\mathrm{12}}\cdot {\mathit{\beta }}_{\mathrm{21}}}{\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}}{\right)}_{j-\mathrm{1}/\mathrm{2}}-\\ & -\frac{i{k}_{x}{h}_{z}}{\mathrm{2}}\left({\mathit{\beta }}_{\mathrm{13}}+\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}\frac{{\mathit{\beta }}_{\mathrm{12}}\cdot {\mathit{\beta }}_{\mathrm{23}}}{\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}}{\right)}_{j-\mathrm{1}}-\frac{i{k}_{x}{h}_{z}}{\mathrm{2}}\left({\mathit{\beta }}_{\mathrm{31}}+\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}\frac{{\mathit{\beta }}_{\mathrm{32}}\cdot {\mathit{\beta }}_{\mathrm{21}}}{\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}}{\right)}_{j}\end{array}\right)\\ & {\stackrel{\mathrm{^}}{\mathit{\alpha }}}_{j}^{\left(+\right)}=\left(\begin{array}{cc}\left(\frac{{\mathit{\beta }}_{\mathrm{22}}}{\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}}{\right)}_{j+\mathrm{1}/\mathrm{2}};& -\left(\frac{{\mathit{\beta }}_{\mathrm{21}}}{\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}}{\right)}_{j+\mathrm{1}/\mathrm{2}}-\frac{i{k}_{x}{h}_{z}}{\mathrm{2}}\left(\frac{{\mathit{\beta }}_{\mathrm{23}}}{\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}}{\right)}_{j+\mathrm{1}}\\ -\left(\frac{{\mathit{\beta }}_{\mathrm{12}}}{\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}}{\right)}_{j+\mathrm{1}/\mathrm{2}}-\frac{i{k}_{x}{h}_{z}}{\mathrm{2}}\left(\frac{{\mathit{\beta }}_{\mathrm{32}}}{\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}}{\right)}_{j+\mathrm{1}};& \left({\mathit{\beta }}_{\mathrm{11}}+\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}\frac{{\mathit{\beta }}_{\mathrm{12}}\cdot {\mathit{\beta }}_{\mathrm{21}}}{\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}}{\right)}_{j+\mathrm{1}/\mathrm{2}}+\\ & +\frac{i{k}_{x}{h}_{z}}{\mathrm{2}}\left({\mathit{\beta }}_{\mathrm{13}}+\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}\frac{{\mathit{\beta }}_{\mathrm{12}}\cdot {\mathit{\beta }}_{\mathrm{23}}}{\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}}{\right)}_{j+\mathrm{1}}+\frac{i{k}_{x}{h}_{z}}{\mathrm{2}}\left({\mathit{\beta }}_{\mathrm{31}}+\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}\frac{{\mathit{\beta }}_{\mathrm{32}}\cdot {\mathit{\beta }}_{\mathrm{21}}}{\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}}{\right)}_{j}\end{array}\right)\\ & {\stackrel{\mathrm{^}}{\mathit{\alpha }}}_{j}^{\left(\mathrm{0}\right)}=\left(\begin{array}{cc}-\left(\frac{{\mathit{\beta }}_{\mathrm{22}}}{\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}}{\right)}_{j-\mathrm{1}/\mathrm{2}}-\left(\frac{{\mathit{\beta }}_{\mathrm{22}}}{\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}}{\right)}_{j+\mathrm{1}/\mathrm{2}}+{k}_{\mathrm{0}}^{\mathrm{2}}{h}_{z}^{\mathrm{2}};& \left(\frac{{\mathit{\beta }}_{\mathrm{21}}}{\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}}{\right)}_{j-\mathrm{1}/\mathrm{2}}+\left(\frac{{\mathit{\beta }}_{\mathrm{21}}}{\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}}{\right)}_{j+\mathrm{1}/\mathrm{2}}\\ \left(\frac{{\mathit{\beta }}_{\mathrm{12}}}{\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}}{\right)}_{j-\mathrm{1}/\mathrm{2}}+\left(\frac{{\mathit{\beta }}_{\mathrm{12}}}{\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}}{\right)}_{j+\mathrm{1}/\mathrm{2}};& -\left({\mathit{\beta }}_{\mathrm{11}}+\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}\frac{{\mathit{\beta }}_{\mathrm{12}}\cdot {\mathit{\beta }}_{\mathrm{21}}}{\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}}{\right)}_{j-\mathrm{1}/\mathrm{2}}-\left({\mathit{\beta }}_{\mathrm{11}}+\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}\frac{{\mathit{\beta }}_{\mathrm{12}}\cdot {\mathit{\beta }}_{\mathrm{21}}}{\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}}{\right)}_{j+\mathrm{1}/\mathrm{2}}+\\ & +{k}_{\mathrm{0}}^{\mathrm{2}}{h}_{z}^{\mathrm{2}}\cdot \left(\mathrm{1}-{\mathit{\beta }}_{\mathrm{33}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}-\frac{{k}_{x}^{\mathrm{4}}}{{k}_{\mathrm{0}}^{\mathrm{4}}}\frac{{\mathit{\beta }}_{\mathrm{23}}\cdot {\mathit{\beta }}_{\mathrm{32}}}{\mathrm{1}-{\mathit{\beta }}_{\mathrm{22}}\frac{{k}_{x}^{\mathrm{2}}}{{k}_{\mathrm{0}}^{\mathrm{2}}}}{\right)}_{j}\end{array}\right)\end{array}\end{array}$
Data availability
Data availability.

The VLF–LF data (Fig. 11) are property of the Shmidt Institute of Physics of the Earth and the University of Sheffield groups, and they are not publicly accessible. According to an agreement between all participants, we cannot make the data openly accessible. Data can be provided under commercial conditions via direct request to rozhnoi@ifz.ru. The ionospheric data used for modeling the electrodynamics characteristics of the VLF waves in the ionosphere are shown in part in Fig. 5 (namely, the altitude distribution of the electron concentration). The other data necessary for the determination of the components of tensor of dielectric permittivity and the electrodynamics modeling in the accepted simple approximation of the three-component plasma-like ionosphere (including electron, one-effective-ion and one-effective-neutral components) and quasi-neutrality are mentioned in Sect. 3.1. The corresponding ionospheric data have been taken from well-known published handbooks referred in the paper .

Author contributions
Author contributions.

YR and VG proposed the idea and concept of the paper, made analytical calculations, wrote the initial version of the paper, and took part in the revision of the paper. VG developed the code, and YR took part in its verification. YR, VG, AsG and AnG performed the numerical modeling. YR and AsG administered the project. VG, YR, AsG and AnG prepared the figures. AL and AsG took part in the writing of the initial version of the paper and the preparation of the revised version. VF, AR and MS took part in the analysis of the data and preparation of the data on VLF propagation in the ionosphere. OA, JB and GM took part in developing the concept of the paper and writing the initial version of the paper. OA, GM, VF and AsG contributed to the preparation of the revised version of the paper. All participants took part in the analysis of the results.

Competing interests
Competing interests.

The authors declare that they have no conflict of interest.

Special issue statement
Special issue statement.

This article is part of the special issue “Solar magnetism from interior to corona and beyond”. It is a result of the Dynamic Sun II: Solar Magnetism from Interior to Corona, Siem Reap, Angkor Wat, Cambodia, 12–16 February 2018.

Financial support
Financial support.

The work of Oleksiy Agapitov and John Bonnell was supported by NASA (grant no. 80NSSC18k0782). The work of Oleksiy Agapitov was partially supported by the NASA LWS program (grant no. 80NSSC20K0218). Viktor Fedun was financially supported by the Royal Society International Exchanges Scheme Standard Programme (grants nos. IE170301, R1/191114) and the Natural Environment Research Council (NERC) (grant no. NE/P017061/1).

Review statement
Review statement.

This paper was edited by Dalia Buresova and reviewed by Pier Francesco Biagi and one anonymous referee.

References

Agapitov, O. V., Artemyev, A. V., Mourenas, D., Kasahara, Y., and Krasnoselskikh, V.: Inner belt and slot region electron lifetimes and energization rates based on AKEBONO statistics of whistler waves, J. Geophys. Res.-Space, 119, 2876–2893, https://doi.org/10.1002/2014JA019886, 2014. a

Agapitov, O. V., Mourenas, D., Artemyev, A. V., Mozer, F. S., Hospodarsky, G., Bonnell, J., and Krasnoselskikh, V: Synthetic Empirical Chorus Wave Model From Combined Van Allen Probes and Cluster Statistics, J. Geophys. Res.-Space, 123, 2017JA024843, https://doi.org/10.1002/2017JA024843, 2018. a

Alperovich, L. S. and Fedorov, E. N.: Hydromagnetic Waves in the Magnetosphere and the Ionosphere, Springer, New York, 426 pp., 2007. a, b, c

Al'pert, Y. L: Propagation of Electromagnetic Waves and the Ionosphere, Nauka Publ., Moscow, 563 pp., 1972 (in Russian). a, b

Arnoldy, R. L. and Kintner, P. M.: Rocket observations of the precipitation of electrons by ground VLF transmitters, J. Geophys. Res.-Space, 94, 6825–6832, https://doi.org/10.1029/JA094iA06p06825, 1989. a, b

Artemyev, A. V., Agapitov, O. V., Mourenas, D., Krasnoselskikh, V., and Zelenyi, L. M.: Storm-induced energization of radiation belt electrons: Effect of wave obliquity, Geophys. Res. Lett., 40, 4138–4143, https://doi.org/10.1002/grl.50837, 2013. a

Artemyev, A., Agapitov, O., Mourenas, D., Krasnoselskikh, V., and Mozer, F.: Wave energy budget analysis in the Earth's radiation belts uncovers a missing energy, Nat. Commun., 6, 7143, https://doi.org/10.1038/ncomms8143, 2015. a

Artsimovich, L. A. and Sagdeev, R. Z: Physics of plasma for physicists, ATOMIZDAT Publ. House, Moscow, 320 pp., 1979 (in Russian). a

Azadifar, M., Li, D., Rachidi, F., Rubinstein, M., Diendorfer, G., Schulz, W., Pichler, H., Rakov, V. A., Paolone, M., and Pavanello, D.: Analysis of lightning-ionosphere interaction using simultaneous records of source current and 380-km distant electric field, J. Atmos. Sol.-Terr. Phy., 159, 48–56, https://doi.org/10.1016/j.jastp.2017.05.010, 2017. a

Biagi, P. F., Maggipinto, T., Righetti, F., Loiacono, D., Schiavulli, L., Ligonzo, T., Ermini, A., Moldovan, I. A., Moldovan, A. S., Buyuksarac, A., Silva, H. G., Bezzeghoud, M., and Contadakis, M. E.: The European VLF/LF radio network to search for earthquake precursors: setting up and natural/man-made disturbances, Nat. Hazards Earth Syst. Sci., 11, 333–341, https://doi.org/10.5194/nhess-11-333-2011, 2011. a

Biagi, P. F., Maggipinto, T., and Ermini, A.: The European VLF/LF radio network: current status, Acta Geod. Geophys., 50, 109–120, https://doi.org/10.1007/s40328-014-0089-x, 2015. a, b

Boudjada, M. Y., Schwingenschuh, K., Al-Haddad, E., Parrot, M., Galopeau, P. H. M., Besser, B., Stangl, G., and Voller, W.: Effects of solar and geomagnetic activities on the sub-ionospheric very low frequency transmitter signals received by the DEMETER micro-satellite, Ann. Geophys., 55, 1, https://doi.org/10.4401/ag-5463, 2012. a

Bragin, Yu. A., Tyutin Alexander, A., Kocheev, A. A., and Tyutin, A. A.: Direct measurements of the vertical electric field of the atmosphere up to 80 km, Space Res., 12, 279–281, 1974. a

Brittain, R., Kintner, P. M., Kelley, M. C., Siren, J. C., and Carpenter, D. L.: Standing wave patterns in VLF Hiss, J. Geophys. Res.-Space, 88, 7059–7064, https://doi.org/10.1029/JA088iA09p07059, 1983. a

Buitink, S., Corstanje, A., Enriquez, J. E., Falcke, H., Hörandel ,J. R., Huege, T., Nelles, A., Rachen, J. P., Schellart, P., Scholten, O., ter Veen, S., Thoudam, S., and Trinh, T. N. G.: A method for high precision reconstruction of air shower Xmax using two-dimensional radio intensity profiles, arXiv:1408.7001v2 [astro-ph.IM], 1 September 2014. a

Chevalier, M. W. and Inan, U. S.: A Technique for Efficiently Modeling Long-Path Propagation for Use in Both FDFD and FDTD, IEEE Ant. Wireless Propag. Lett., 5, 525–528, https://doi.org/10.1109/LAWP.2006.887551, 2006. a

Chou, M.-Y., Lin, C. C. H., Yue, J., Chang, L. C., Tsai, H.-F., and Chen, C.-H.: Medium-scale traveling ionospheric disturbances triggered by Super Typhoon Nepartak, Geophys. Res. Lett., 44, 7569–7577, https://doi.org/10.1002/2017gl073961, 2017. a

Collin, R. E.: Foundations for Microwave Engineering, John Willey & Sons, Inc., New York, 2001. a, b

Cummer, S. A., Inan, U. S., Bell, T. F., and Barrington-Leigh, C. P.: ELF radiation produced by electrical currents in sprites, Geophys. Res. Lett., 25, 1281–1284, https://doi.org/10.1029/98GL50937, 1998. a

Cummer, S. A., Briggs, M. S., Dwyer, J. R., Xiong, S., Connaughton, V., Fishman, G. J., Lu, G., Lyu, F., and Solanki, R.: The source altitude, electric current, and intrinsic brightness of terrestrial gamma ray flashes, Geophys. Res. Lett., 41, 8586–8593, https://doi.org/10.1002/2014GL062196, 2014. a

Dwyer, J. R.: The relativistic feedback discharge model of terrestrial gamma ray flashes, J. Geophys. Res., 117, A02308, https://doi.org/10.1029/2011JA017160, 2012. a

Dwyer, J. R. and Uman, M. A.: The physics of lightning, Phys. Rep., 534, 147–241, https://doi.org/10.1016/j.physrep.2013.09.004, 2014. a

Glukhov, V. S., Pasko, V. P., and Inan, U. S.: Relaxation of transient lower ionospheric disturbances caused by lightning-whistler-induced electron precipitation bursts, J. Geophys. Res., 97, 16971, https://doi.org/10.1029/92ja01596, 1992. a

Grimalsky, V. V., Kremenetsky, I. A., and Rapoport, Y. G.: Excitation of electromagnetic waves in the lithosphere and their penetration into ionosphere and magnetosphere, J. Atmos. Electr., 19, 101–117, 1999a. a

Grimalsky, V. V., Kremenetsky, I. A., and Rapoport, Y. G.: Excitation of EMW in the lithosphere and propagation into magnetosphere, in: Atmospheric and Ionospheric Electromagnetic Phenomena Associated with Earthquakes, edited by: Hayakawa, M., TERRAPUB, Tokyo, 777–787, 1999b. a

Guglielmi A. V. and Pokhotelov, O. A.: Geoelectromagnetic waves, IOP Publ., Bristol, 402 pp., 1996. a

Gurevich, A. V. and Zybin, K. P.: Runaway breakdown and electric discharges in thunderstorms, Physics-Uspekhi, 44, 1119–1140, https://doi.org/10.1070/PU2001v044n11ABEH000939, 2001. a, b, c

Gurevich, A. V., Karashtin, A. N., Ryabov, V. A., and Chubenko, A. L. S.: Nonlinear phenomena in the ionospheric plasma.Effects of cosmic rays and runaway breakdown on thunderstorm discharges, Physics-Uspekhi, 52, 735–745, https://doi.org/10.3367/UFNe.0179.200907h.0779, 2009. a

Hare, B. M., Scholten, O., Bonardi, A., Buitink, S., Corstanje, A., Ebert, U., Falcke, H., Hörandel, J. R., Leijnse, H., Mitra, P., Mulrey, K., Nelles, A., Rachen, J. P., Rossetto, L., Rutjes, C., Schellart, P., Thoudam, S., Trinh, T. N. G., ter Veen, S., and Winchen, T.: LOFAR Lightning Imaging: Mapping Lightning With Nanosecond Precision, J. Geophys. Res.-Atmos., 123, 2861–2876, https://doi.org/10.1002/2017JD028132, 2018. a

Hayakawa, M.: Earthquake prediction with radio techniques, Wiley, Singapore, 294 pp., 2015. a

Jursa, A. S. (Ed.): Handbook of Geophysics and the Space Environment Air Force Geophyics Lab., Hanscom Air Force Base, Mass.: Air Force Geophysics Laboratory, Air Force Systems Command, United States Air Force, 1042 pp., 1985. a, b

Kelley, M. C.: The Earth's Ionosphere Plasma Physics and Electrodynamics, Acad. Press, Amsterdam, 556 pp., 2009. a, b

Kim, S. and Kim, K.: Invariant imbedding theory of wave propagation in arbitrarily inhomogeneous stratified bi-isotropic media, J. Opt., 18, 065605, https://doi.org/10.1088/2040-8978/18/6/065605, 2016. a, b, c, d, e

Kintner, P. M., Brittain, R., Kelley, M. C., Carpenter, D. L., and Rycroft, M. J.: In situ measurements of transionospheric VLF wave injection, J. Geophys. Res.-Space, 88, 7065–7073, https://doi.org/10.1029/JA088iA09p07065, 1983. a

Koskinen, H. E. J.: Physics of Space Storms. From the Solar Surface to the Earth, Springer-Verlag, Berlin–Heidelberg, 419 pp., 2011. a

Kurushin, E. P. and Nefedov, E. I.: Electrodynamics of anisotropic waveguiding structures, Moscow: Nauka (Science), 1983 (in Russian). a, b, c, d, e

Kuzichev, I. V. and Shklyar, D. R.: On full-wave solution for VLF waves in the near-Earth space, J. Atmos. Solar-Terr. Phy., 72 , 1044–1056, https://doi.org/10.1016/j.jastp.2018.07.002, 2010. a

Kuzichev, I.V., Vasko, I. Yu., Malykhin, A. Yu., and Soto-Chavez, A. R.: On the ionospheric propagation of VLF waves generated by currents in the lower ionosphere, J. Atmos. Solar-Terr. Phy., 179, 138–148, https://doi.org/10.1016/j.jastp.2018.07.002, 2018. a

Lehtinen, N. G. and Inan, U. S.: Radiation of ELF/VLF waves by harmonically varying currents into a stratified ionosphere with application to radiation by a modulated electrojet, J. Geophys. Res., 113, A06301, https://doi.org/10.1029/2007JA012911, 2008. a

Lehtinen, N. G. and Inan, U. S.: Full-wave modeling of transionospheric propagation of VLF waves, Geophys. Res. Lett., 36, L03104, https://doi.org/10.1029/2008GL036535, 2009. a

Levy, M.: Parabolic equation methods for electromagnetic wave propagation, The Inst. of Electrical Eng., Padstow, Cornwall, 336 pp., 2000. a

Lu, G., Zhang, H., Cummer, S. A., Wang, Y., Lyu, F., Briggs, M., Xiong, S., and Chen, A.: A comparative study on the lightning sferics associated with terrestrial gamma-ray flashes observed in Americas and Asia, J. Atmos. Solar-Terr. Phy. 183, 67–75, https://doi.org/10.1016/j.jastp.2019.01.001, 2019. a

Maier, S. A.: Plasmonics: Fundamentals and Applications, Springer, NY, 234 pp., 2007. a

Marshall, R. A. and Wallace, T.: Finite-Difference Modeling of Very-Low-Frequency Propagation in the Earth-Ionosphere Waveguide, IEEE Trans. Ant. Propag., 65, 7185–7197, https://doi.org/10.1109/TAP.2017.2758392, 2017. a

Martynenko, S. I., Rozumenko, V. T., and Tyrnov, O. F.: New possibilities for mesospheric electricity diagnostics, Adv. Space Res., 27, 1127–1132, https://doi.org/10.1016/S0273-1177(01)00208-3, 2001. a

Meek, C. E., Manson, A. H., Martynenko, S. I., Rozumenko, V. T., and Tyrnov, O. F.: Remote sensing of mesospheric electric fields using MF radars, J. Atmos. Solar-Terr. Phy., 66, 881–890, https://doi.org/10.1016/j.jastp.2004.02.002, 2004. a

Mezentsev, A., Lehtinen, N., Ostgaard, N., Perez-Invernon, F. J., and Cummer, S. A.: Spectral characteristics of VLF sferics associated with RHESSI TGFs, J. Geophys. Res.-Atmos., 123, 139–159, https://doi.org/10.1002/2017JD027624, 2018. a

Nina, A., Radovanović, M., Milovanović, B., Kovačević, A., Bajčetić, J., and Popović, L. C.: Low ionospheric reactions on tropical depressions prior hurricanes, Adv. Space Res., 60, 1866–1877, https://doi.org/10.1016/j.asr.2017.05.024, 2017. a

Patra, S., Spencer, E. Horton, W., and Sojka, J.: Study of Dst/ring current recovery times using the WINDMI model, J. Geophys. Res., 116, A02212, https://doi.org/10.1029/2010JA015824, 2011. a

Pulinets, S. and Boyarchuk, K.: Ionospheric Precursors of Earthquakes, Springer, Berlin, 315 pp., 2005. a

Qin, J., Celestin, S., and Pasko, V. P.: Low frequency electromagnetic radiation from sprite streamers, Geophys. Res. Lett., 39, L22803, https://doi.org/10.1029/2012GL053991, 2012. a

Rapoport, Y. G., Gotynyan, O. E., Ivchenko, V. N., Hayakawa, M., Grimalsky, V. V., Koshevaya, S. V., and Juarez, D.: Modeling electrostatic – photochemistry seismoionospheric coupling in the presence of external currents, Phys. Chem. Earth, 31, 437–446, https://doi.org/10.1016/j.pce.2006.02.010, 2006. a, b, c, d

Rapoport, Y. G., Boardman, A. D., Grimalsky, V. V., Ivchenko, V. M., and Kalinich, N: Strong nonlinear focusing of light in nonlinearly controlled electromagnetic active metamaterial field concentrators, J. Opt.-UK, 16, 0552029–0552038, https://doi.org/10.1088/2040-8978/16/5/055202, 2014. a, b

Rozhnoi, A., Solovieva, M., Levin, B., Hayakawa, M., and Fedun, V.: Meteorological effects in the lower ionosphere as based on VLF/LF signal observations, Nat. Hazards Earth Syst. Sci., 14, 2671–2679, https://doi.org/10.5194/nhess-14-2671-2014, 2014. a

Rozhnoi, A., Solovieva, M., Parrot, M., Hayakawa, M., Biagi, P.-F., Schwingenschuh, K., and Fedun, V.: VLF/LF signal studies of the ionospheric response to strong seismic activity in the Far Eastern region combining the DEMETER and ground-based observations, Phys. Chem. Earth, 85–86, 141–149, https://doi.org/10.1016/j.pce.2015.02.005, 2015. a, b, c, d, e, f

Ruibys, G. and Tolutis, R: Nonreciprocal HF signal transmission by surface helicon, Electron. Lett., 19, 273, https://doi.org/10.1049/el:19830191, 1983. a, b, c

Samarskii, A. A.: The Theory of Difference Schemes, Marcel Dekker, NY, 2001. a, b, c

Sanchez-Dulcet, F., Rodriguez-Bouza M., Silva, H. G., Herraiz, M., Bezzeghoud, M., and Biagi, P. F.: Analysis of observations backing up the existence of VLF and ionospheric TEC anomalies before the Mw6.1 earthquake in Greece, January 26, 2014, Phys. Chem. Earth, 85–86, 150–166, https://doi.org/10.1016/j.pce.2015.07.002, 2015. a

Scholten, O., Bonardi, A., Buitink, S., Corstanje, A., Ebert, U., Falcke, H., Hörandel, J., Mitra, P., Mulrey, K., Nelles, A., Rachen, J., Rossetto, L., Rutjes, C., Schellart, P., Thoudam, S., Trinh, G., ter Veen, S., and Winchen, T.: Precision study of radio emission from air showers at LOFAR, EPJ Web of Conferences, 136, 02012, https://doi.org/10.1051/epjconf/201713602012, 2017. a

Schunk, R. W. and Nagy, A. F.: Ionospheres, Physics, Plasma Physics, and Chemistry, CUP, Cambridge, 628 pp., 2010. a, b

Senior, T. B. A. and Volakis, J. L.: Approximate boundary conditions in electromagnetics, Institution of Electrical Engineers, London, 1995. a, b, c, d, e

Shalashov, A. G. and Gospodchikov, E. D.: Impedance technique for modeling electromagnetic wave propagation in anisotropic and gyrotropic media, Physics-Uspekhi, 54, 145–165, https://doi.org/10.3367/UFNe.0181.201102c.0151, 2011. a, b, c, d, e

Shinagawa, H., Tsugawa, T., Matsumura, M., Iyemori, T., Saito, A., Maruyama, T., and Otsuka, Y.: Two-dimensional simulation of ionospheric variations in the vicinity of the epicenter of the Tohoku-oki earthquake on 11 March 2011, Geophys. Res. Lett., 40, 5009–5013, https://doi.org/10.1002/2013gl057627, 2013. a

Siefring, C. L. and Kelley, M. C.: Analysis of standing wave patterns in VLF transmitter signals: Effects of sporadic E layers and in situ measurements of low electron densities, J. Geophys. Res.-Space, 96, 17813–17826, https://doi.org/10.1029/91JA00615, 1991. a

Spiegel, M. R.: Theory and problems of vector analysis and an introduction to tensor analysis, Shnaum Publ., NY, 223 pp., 1959. a

Surkov, V. and Hayakawa, M.: Ultra and Extremely Low Frequency Electromagnetic Fields, Springer, Tokyo, 2014. a

Tarkhanyan, R. H. and Uzunoglu, N. K., Radiowaves and Polaritons in Anisotropic Media, Wiley-VCH, Weinheim, 210 pp., 2006. a

Tolstoy, A., Rosenberg, T. J., Inan, U. S., and Carpenter, D. L.: Model predictions of subionospheric VLF signal perturbations resulting from localized, electron precipitation-induced ionization enhancement regions, J. Geophys. Res., 91, 13473, https://doi.org/10.1029/ja091ia12p13473, 1986. a

Tretyakov, S.: Analytical Modeling in Applied Electromagnetics, Artech House, Boston, 2003. a, b, c, d

Wait, J. R.: Electromagnetic waves in stratified media, IEEE Press and Oxford Univ. Press, New York, 1996. a, b, c, d

Walker, A. D. M.: The Theory of Whistler Propagation, Rev. Geophys. Space Phys., 14, 629–638, 1976. a

Weiland, J. and Wilhelmsson, H.: Coherent Non-Linear Interaction of Waves in Plasmas, Pergamon, London, 1977. a

Wu, C.-C., Liou, K., Lepping, R. P., Hutting, L., Plunkett, S., Howard, R. A., and Socker, D.: The first super geomagnetic storm of solar cycle 24: “The St. Patrick's day event (17 March 2015)”, Earth Planets Space, 68, 151, https://doi.org/10.1186/s40623-016-0525-y, 2016. a

Yiğit, E., Koucká Knížová, P., Georgieva, K., and Ward, W.: A review of vertical coupling in the Atmosphere–Ionosphere system: Effects of waves, sudden stratospheric warmings, space weather, and of solar activity, Solar-Terr. Phy., 141, 1–12, https://doi.org/10.1016/j.jastp.2016.02.011, 2016. a

Yu, Y., Niu, J., and Simpson, J. J.: A 3-D Global Earth-Ionosphere FDTD Model Including an Anisotropic Magnetized Plasma Ionosphere, IEEE Trans. Ant. Propag., 60, 3246–3256, https://doi.org/10.1109/TAP.2012.2196937, 2012. a