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

. 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 inﬂuence of geomagnetic ﬁeld


Introduction
The results of the analytical and numerical study of very low-frequency (VLF) electromagnetic (EM) wave/beam propagation in the lithosphere-atmosphere-ionospheremagnetosphere system (LAIM), in particular in the Earthionosphere waveguide (WGEI), are presented. The amplitude and phase of the VLF wave propagates in the Earthionosphere 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 (Patra et al., 2011;Koskinen, 2011;Boudjada et al., 2012;Wu et al., 2016;Yigit et al., 2016). 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 (Nina et al., 2017;Rozhnoi et al., 2014;Chou et al., 2017) as well as from earthquakes (Hayakawa, 2015;Surkov and Hayakawa, 2014;Sanchez-Dulcet et al., 2015) and tsunamis. From inside the ionosphere, strong thunderstorms, lightning discharges, and terrestrial gamma-ray flashes or sprite streamers (Cummer et al., 1998;Qin et al., Published by Copernicus Publications on behalf of the European Geosciences Union. 208 Y. Rapoport et al.: Tensor impedance method 2012;Dwyer, 2012;Dwyer and Uman, 2014;Cummer et al., 2014;Mezentsev et al., 2018) 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 (Gurevich and Zybin, 2001). Consequently, these phenomena are intensified during the air shower generating by cosmic particles (Gurevich and Zybin, 2001;Gurevich et al., 2009). 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 (Buitink et al., 2014;Scholten et al., 2017;Hare et al., 2018), and in the VLF range (Gurevich and Zybin, 2001). 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 (Biagi et al., 2011. The World Wide Lightning Location Network is one of the international facilities for VLF measurements during thunderstorms with lightning discharges (Lu et al., 2019). 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 (Artemyev et al., 2013(Artemyev et al., , 2015Agapitov et al., 2014Agapitov et al., , 2018. 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 (Shalashov and Gospodchikov, 2011;Kim and Kim, 2016), our model provides an optimal balance between the analytical and numerical approaches. It combines analytical and numerical approaches based on the matrix sweep method (Samarskii, 2001). 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 finitedifference time domain models (Chevalier and Inan, 2006;Marshall and Wallace, 2017;Yu et al., 2012;Azadifar et al., 2017), 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 consid-ered in Kuzichev and Shklyar (2010), Kuzichev et al. (2018), Inan (2009, 2008) 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 (Tretyakov, 2003;Senior and Volakis, 1995;Kurushin and Nefedov, 1983) are well known and can be used, in particular, for the layered metaldielectric, metamaterial and gyrotropic active layered and waveguiding media of different types (Tretyakov, 2003;Senior and Volakis, 1995;Kurushin and Nefedov, 1983;Collin, 2001;Wait, 1996) including plasma-like solid state (Ruibys and Tolutis, 1983) and space plasma (Wait, 1996) media. The plasma wave processes in the metal-semiconductordielectric waveguide structures, placed into the external magnetic field, were widely investigated (Ruibys and Tolutis, 1983;Maier, 2007;Tarkhanyan and Uzunoglu, 2006) 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 (Tretyakov, 2003;Senior and Volakis, 1995;Kurushin and Nefedov, 1983). As a rule, impedance boundary conditions are used when the layer covering waveguide is thin (Senior and Volakis, 1995;Kurushin and Nefedov, 1983). 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 Y. Rapoport et al.: Tensor impedance method 209 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 (Grimalsky et al., 1999a, b) or dielectric-magnetized semiconductor structure. The electromagnetic waves radiated outside the waveguiding structure, such as helicons (Ruibys and Tolutis, 1983) or whistlers (Wait, 1996), 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 (Walker, 1976) 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 (∼ 85 km) 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.

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 > 1000 km. The Earth's surface of a high conductivity of z = 0 (z is vertical coordinate) and the ionosphere F layer of z = 300 km 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 0 < z < 60 km, the approximate altitude of the nearly isotropic ionosphere D layer at 60 km < z < 75 km, and the anisotropic E and F layers of the ionosphere, due to the geomagnetic field H 0 , 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 H 0 . Note that the angle of θ is complementary to the angle of inclination of the geomagnetic field. The geomagnetic field H 0 is directed along z ′ axis and lies in the plane xz, while the planes x ′ z ′ and xz coincide with each other.

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.

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 = L max ) ( Fig. 1), resulting in the tensor impedance conditions at the upper boundary of the effective WGEI (z = L i ). 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 L z at ∼ 250-400 km to the lower ionosphere at L z ∼ 70-90 km. All Figure 1. The geometry of the anisotropic and gyrotropic waveguide. EM waves propagate in the OX direction. H 0 is the external magnetic field. The effective WGEI for EM waves occupies the region 0 < z < L z . The isotropic medium occupies the region 0 < z < L ISO , L ISO < L z . The anisotropic and gyrotropic medium occupies the region L ISO < z < L max . The covering layer occupies the region L z < z < L max . WGEI includes the isotropic region 0 < z < L ISO and a part of the anisotropic region L z < z < L max . It is supposed that the anisotropic region is relatively small part of the WG, (L z − L ISO )/L z ∼ (0.1-0.2). At the upper boundary of the covering layer (z = L max ) EM radiation to the external region (z > L max ) 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 = L z ). These boundary conditions effectively include all the effects on the wave propagation of the covering layer and the radiation (at z = L max ) to the external region (z > L max ). 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ε(ω, z) or by the inverse tensorβ(ω, z) = ε −1 (ω, z): E =β(ω, z)·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 ′ Y Z ′ , where the OZ ′ axis is aligned along the geomagnetic field H 0 .ε where ω pe , ω pi , ω H e and ω H i are the plasma and cyclotron frequencies for electrons and ions, respectively, m e , m i , ν 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 plasmalike ionosphere (including one electron component, one effective ion and one effective neutral component) and quasineutrality are accepted. The expressions for the components of the permittivity tensorε(ω, z) are obtained from Eq. (1) by means of multiplication with the standard rotation matrices (Spiegel, 1959) dependent on angle θ. For the medium with a scalar conductivity σ , e.g., lower ionosphere or atmosphere, the effective permittivity in Eq. (1) reduces to ε = 1 − 4π iσ/ω.

The equations for the EM field and upper boundary conditions
In the case of the VLF waveguide modes, the longitudinal wave number k x 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(−ik x x). Taking into account k x ≈ k 0 (k 0 = ω/c), in simulations of the VLF beam propagation, it is possible to put k x = k 0 . Therefore, where E x = β 11 D x + β 12 D y + β 13 D z , etc. All components of the EM field can be represented through the horizontal components of the magnetic field H x and H y . The equations for these components are given below.
The expressions for the horizontal components of the electric field E x , E y are given below.
In the region z ≥ L max 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 = L max 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 ∼ exp(−ik z z), while the derivatives like ∂β 11 /∂z and the inhomogeneity of the media are neglected.
Thus, generally Eq. (5) determined the wave numbers for the outgoing waves to be of the fourth order (Wait, 1996). The boundary conditions at the upper boundary z = L max 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 Im(k z1,z2 ) < 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 k 1,2 are obtained from Eq. (3a, 3b), where the dependence on x is neglected: |k 1,2 | ≫ k 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 z ≥ L max , the EM field components can be presented as In Eq. (6),z = z−L z . Equation (3a, 3b) are simplified in the approximation described above.
The solution of Eq. (7) is searched for as H x,y ∼ exp(−ik zz ).
The following equation has been obtained to get the wave numbers k z1,z2 from Eq. (7): Therefore, from Eq. (8) follows, The signs of k z1,z2 have been chosen from the condition Im(k z1,z2 ) < 0. From Eq. (5) at the upper boundary of z = L max , the following relations are valid: From Eq. (10) one can get Thus, it is possible to exclude the amplitudes of the outgoing waves A 1,2 from Eqs. (9). As a result, at z = L max the bound-ary conditions are rewritten in terms of H x and H y only.
The relations in Eq. (12) are the upper boundary conditions of the radiation for the boundary z = L max at ∼ 250-400 km. These conditions will be transformed and recalculated using the analytical numerical recurrent procedure into equivalent impedance boundary conditions at z = L z at ∼ 70-90 km. Note that in the "whistler-VLF approximation" is valid at frequencies ∼ 10 kHz for the F region of the ionosphere. In this approximation and k x ≈ 0, we receive the dispersion equation using Eqs. (5), (8), (9), in the form where and k ′ x and k ′ z are the transverse and longitudinal components of the wave number relative to geomagnetic field. For the F region of the ionosphere, where ν e ≪ ω ≪ ω H e , Eq. (13) reduces to the standard form of the whistler dispersion equation |k ′ z |k = k 0 |g|, g ≈ −ω 2 pe /(ωω H e ) and ω = c 2 k|k ′ z |(ω H e /ω 2 pe ). In a special case of the waves propagating along geomagnetic field, k ′ x = 0, for the propagating whistler waves, the well-known dispersion dependence is ω = c 2 k ′ z 2 (ω H e /ω 2 pe ) (Artsimovich and Sagdeev, 1979). For the formulated problem we can reasonably assume k x ≈ 0. Therefore Eq. (13) is reduced to k 4 z cos 2 θ = k 4 0 g 2 . As a result, we get k z1 = √ g/ cos θk 0 and k z2 = −i √ g/ cos θk 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 Conditions in Eq. (12) or Eq. (14) are the conditions of radiation (absence of ingoing waves) formulated at the upper boundary at z = L max 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 | → 0 in Eq. (5). In spite of the disappearance of the dependence of these boundary conditions explicitly on k x , the dependence of the characteristics of the wave propagation process on k x , 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 k x ∼ k 0 .

Equivalent tensor impedance boundary conditions
The tensor impedance at the upper boundary of the effective WGEI of z = L z (see Fig. 1) is obtained by the conditions of radiation in Eqs. (12) or (14), recalculated from the level of z = L max ∼ 250-400 km, placed in the F region of the ionosphere, to the level of z = L z ∼ 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., Tretyakov, 2003;Senior and Volakis, 1995;Kurushin and Nefedov, 1983;Alperovich and Fedorov, 2007). There is one known exception, namely the invariant imbedding impedance method (Shalashov and Gospodchikov, 2011;Kim and Kim, 2016). 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 In Eq. (15) z j +1/2 = h · (j + 0.5). In Eq. (10) the approxima- 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 where 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 Samarskii, 2001). It can be writ-ten as This method is a variant of the Gauss elimination method for the matrix three-diagonal set of Eq. (16). The value ofb N is obtained from the boundary conditions (12) aŝ Then the matricesb j have been computed sequentially down to the desired value of z = L z = h · N z , where the impedance boundary conditions are assumed to be applied. At each step the expression forb j follows from Eqs. (16), (17a) and (17b) as and a similar equation can be obtained for ∂H y ∂z N z . Note that as a result of this discretization, only the values at the grid level N z are included in the numerical approximation of the derivatives ∂H x,y /∂z at z = L z . We determine tensor impedanceẐ at z = L z at ∼ 85 km. 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 N z and corresponding to this altitude).
The equivalent tensor impedance is obtained using a two-step procedure.
(1) We obtain the matrixb j using Eq. (3a, 3b) with the boundary conditions of Eq. (12) and the procedure of Eqs. (17) The proposed method of the transfer of the boundary conditions from the ionosphere F layer at L max = 250-400 km into the lower part of the E layer at L z = 80-90 km is stable and easily realizable in comparison with some alternative approaches based on the invariant imbedding methods (Shalashov and Gospodchikov, 2011;Kim and Kim, 2016). 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 H x and H y in Eq. (3a, 3b) seems natural, as well as the transfer of Eq. (17a), because the impedance conditions are the expressions of the electric E x and E y through the magnetic components H x and H y 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 = L max to z = L z , 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 0αβ | ≤ 0.1. This determines the choice of the upper boundary at z = L z for the effective WGEI. Due to the small impedance, EM waves incident from below on this boundary are reflected effectively back. Therefore, the region 0 ≤ z ≤ L z indeed can be presented as an effective WGEI. This waveguide includes not only lower boundary at L ISO 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 (Shalashov and Gospodchikov, 2011;Kim and Kim, 2016), 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 = L max (see Eq. 12) are determined through the total field components of H x,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.

Propagation of electromagnetic waves in the gyrotropic waveguide and the TIMEB method
Let us use the transverse components of the electric E y and magnetic H y fields to derive equations for the slow varying amplitudes A(x, y, z) and B(x, y, z) of the VLF beams. These components can be represented as Here we assumed k x = k 0 to reflect beam propagation in the WGEI with the main part in the atmosphere and lower iono-sphere (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 (Guglielmi and Pokhotelov, 1996) of the ionosphere causes, altogether with the impedance boundary condition, the proper z dependence of B(x, y, z). Using Eqs. (21) and (22), the boundary conditions are determined at the height of z = L z for the slowly varying amplitudes A(x, y, z) and B(x, y, z) of the transverse components E y and H y . As it follows from Maxwell's equations, the components E x and H x through E y and H y in the method of beams have the form where γ 12 = −1 0 (ε 13 ε 32 −ε 12 ε 33 ),β 13 = −1 0 ε 13 andβ 33 = −1 0 ε 33 ; 0 = ε 11 ε 33 − ε 13 ε 31 . From Eqs. (21) and (24), the boundary conditions for A and B can be defined as The evolution equations for the slowly varying amplitudes A(x, y, z) and B(x, y, z) 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 E, H ∼ exp(iωt −ik x x− ik y y), Maxwell's equations are Here D x = ε 11 E x +ε 12 E y +ε 13 E z . From Eq. (21), the equations for E x and E z through E y andH y are In Eq. (27) equations for E y and H y obtained from the Maxwell equations are After the substitution of Eq. (27) for E x and E z into Eqs. (28), the coupled equations for E y and H y can be derived. The following expansion should be used: k x = k 0 + δk x , |δk x | ≪ k 0 ; also |k y | ≪ k 0 . Then, according to Weiland and Wilhelmsson (1977), The expansions should be used until the quadratic terms in k y and the linear terms in δk x . As a result, parabolic equations (Levy, 2000) 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.
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 where σ E ∼ 10 8 s −1 is the Earth's conductivity. The boundary conditions (31a) at the Earth's surface, where Z 022 = Z 021 ≡ Z 0E , Z 012 = Z 021 = 0, β 33 = ε −1 (z = 0), γ 12 = 0 andβ 13 = 0 can be rewritten as In the relations of Eq. (32), z 1 , z 0 , y 0 and B 0 are the vertical position of maximum value, the vertical and transverse characteristic dimensions of the spatial distribution and the maximum value of H y , 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, L y ∼ 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 ∂ 2 A/∂y 2 and ∂ 2 B/∂y 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 (Samarskii, 2001). Namely, the problem has been approximated by the finite differences In the terms ofL y C, the derivatives with respect to y are included, whereas all other terms are included inL z C. Then the following fractional steps have been applied; the first one is along y, and the second one is along z.
Here h x is the step along the OX axis, and x p = ph x and p = 0, 1, 2, . . .. This step has been chosen from the conditions of the simulation result independence of the diminishing h x .
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 h x and ∼ 1 km ≪ L z , this correction of the C results in the multiplier exp(−ik 0 · δx), where δx = z · (h x /R E ) and R E ≫ L z 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.

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 dipolelike VLF radio source and lightning discharge. Then, the reflection of the waves incident on the upper boundary (z = L z ) 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 (Rapoport et al., 2014) would be quite suitable for modeling beam propagation, even accounting for the wave dispersion in magnetized plasma. The proper effective boundary condition, similar to Rapoport et al. (2014), would allow for making a rel-atively 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.

Modeling results
The dependencies of the permittivity components ε 1 , ε 3 and ε h in the coordinate frame associated with the geomagnetic field H 0 are given in Fig. 3. The parameters of the ionosphere used for modeling are taken from Al'pert (1972), Alperovich and Fedorov (2007), Kelley (2009), Schunk andNagy (2010), andJursa (1985). 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 ω = 10 5 s −1 and f = ω/2π ≈ 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) (1 Gs = 10 −4 T), whereas the electric field is also in Gs, 1 Gs = 300 V 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 = 75 km. 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 surfaceionosphere waveguide. The computed components of the tensor impedance at z = 85 km are Z 011 = 0.087 + i0.097, Z 021 = 0.085 + i0.063, Z 012 = −0.083 − i0.094 and Z 022 = 0.093 + i0.098. So, a condition of |Z 0αβ | ≤ 0.15 is satisfied there, which is necessary for the applicability of the boundary conditions in Eq. (25). The maximum value of the H y component is 0.1 Oe = 10 −5 T in Fig. 4a for the ini-tial VLF beam at x = 0. This corresponds to the value of the E z component of 0.1 Gs = 30 V cm −1 . At the distance of x = 1000 km the magnitudes of the magnetic field H y are of about 3 × 10 −5 Oe = 3 nT, whereas the electric field E y is 3 × 10 −6 Gs ≈ 1 mV cm −1 .
The wave beams are localized within the WGEI at 0 < z < 75 km, 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 75 km < z < 85 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, H y = 0 and E y = 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 H y = 0 and E y = 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 E y , is also excited. This effect is illustrated in Fig. 4b and d.
The magnitude of the E y 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. 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 H y 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.
The distributions of |E y | and |H y | on z at x = 1000 km in the center of the waveguide, y = 1500 km, 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 = 120 km af-fects the excitation of the E y component within the waveguide rather weakly.

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 = L z (Fig. 1). For the corresponding transformation of the field polarization determination, we introduce the characteristic polarization relation |E y /H y |(z; y = L y /2; x = x 0 ), taken at the central plane of the beam (y = L y ) at the characteristic distance (x = x 0 ) 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ǫ is different from the isotropicÎ 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 0αβ | ≤ 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 H y = 0, E y = 0) for the initial beam(s), the above-mentioned value |E y /H y |(z; y = L y /2; x = x 0 ) characterizes the mode coupling and corresponding transformation of the polarization at the distance of x 0 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 x 0 = 1000 km, 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  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 3 × 10 −6 Gs ≈ 1 mV cm −1 ; for the magnetic field, the maximum value (e) is 3 × 10 −5 Gs ≈ 3 nT. At the altitudes z < 75 km, |E z | ≈ |H y |. ω = 1.0 × 10 5 s −1 and θ = 45 • .
of |H ymax (x = x 0 )/H ymax (x = 0)|. 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.
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 (ω = 1.14 × 10 5 s −1 ), the larger maximum of the polarization parameter corresponds to the intermediate value of the angle θ = 45 • (Fig. 8b); for the lower frequency (ω = 0.86×10 5 s −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 (ω = 1.14 × 10 5 s −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 n e in the lower ionosphere on the polar-  ization parameter, we have used the following parameterization for the n e change n e = n e (z) − n 0e (z) of the electron concentration, where n 0e (z) is the unperturbed altitude distribution of the electron concentration.
n e (z) = n 0e (z) (z); In Eq. (35), z 12 ≡ z 2 − z 1 ; z is the effective width of the electron concentration perturbation altitude distribution. The perturbation n e is concentrated in the range of altitudes z 1 ≤ z ≤ z 2 and is equal to zero outside this region, n e (z 1 ) = n e (z 2 ) = 0, while (z 1 ) = (z 2 ) = 0. 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 concentra-tion have been reported as a result of seismogenic phenomena, tsunamis, particle precipitation in the ionosphere due to wave-particle interaction in the radiation belts (Pulinets and Boyarchuk, 2005;Shinagawa et al., 2013;Arnoldy and Kintner, 1989;Glukhov et al., 1992;Tolstoy et al., 1986), 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 ω = 1.14 × 10 5 s −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).
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(Z 011 ) and Im(Z 011 ) values in general decrease with an increasing angle of θ (see Fig. 10a and b). The average values of Re(Z 011 ) at θ = 5, 30, 45 and 60 • (lines 1-4 in Fig. 10a) and Im(Z 011 ) at θ = 45 • and θ = 60 • (curves 3 and 4 in Fig. 10b) increase with an increasing frequency in the frequency range (0.86-1.14) ×10 5 s −1 . The average Im(Z 011 ) value at θ = 5 and 30 • changes in the frequency range (0.86-1.14) ×10 5 s −1 non-monotonically with the maximum at (1-1.1) ×10 5 s −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 Z 0 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 ∼ 75 km (Fig. 10c).

Discussion
The observations presented in Rozhnoi et al. (2015) show a possibility for seismogenic increasing losses of VLF waves in the WGEI (Fig. 11; see details in Rozhnoi et al., 2015). We discuss the qualitative correspondence of our results to these experimental data.
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 (Rapoport et al., 2006), the presence of the mesospheric current source, which followed from the observations (Martynenko et al., 2001;Meek et al., 2004;Bragin et al., 1974) is also taken into account. It is assumed that the mesospheric current only has the z com-ponent 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 Rapoport et al. (2006) that the decrement of losses |k ′′ | for VLF waves in the WGEI is proportional to |k ′′ | ∼ |ε ′′ | ∼ N e /ν e . N e and ν e decrease and increase, respectively, due to the appearance of a Figure 8. Characteristics of the polarization transformation parameter |E y /H y | (a-c) and the effective coefficient of the total losses at the distance of x 0 = 1000 km 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 ω = 0.86 × 10 5 s −1 and ω = 1.14 × 10 5 s −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 × 10 5 , 1 × 10 5 and 1.14 × 10 5 s −1 , respectively, and θ = 45 • . Figure 9. (a) Decreased and increased electron concentration (line 2, red) and (line 3, blue) correspond to f 00 = −1.25 and f 00 = 250, respectively, relative to the reference concentration (line 1, black) with parametrization conditions (see Eq. 35) z 1 = 50 km, z 2 = 90 km and = 20 km. 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 0.86 × 10 5 s −1 , 1.0 × 10 5 s −1 and 1.14 × 10 5 s −1 , respectively. Angle θ is equal to 45 • . 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 in- Figure 10. (a, b) Frequency dependencies of the real (a) and imaginary (b) parts of the effective tensor impedance Z 011 component at the upper boundary (z = L z , 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 0.86 × 10 5 s −1 and angle θ = 45 • for the isotropic surface impedance Z 0 at the lower surface of the WGEI equal to 10 −4 . Earth conductivity σ is equal to 10 9 s −1 for line 1 and Z 0 = 10 −2 (σ = 10 7 s −1 ) for line 2. crease in losses corresponds qualitatively to the results presented in Rozhnoi et al. (2015).
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 in-side the WGEI (see Figs. 4-7). 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 0 ≤ z ≤ L max , where L max = 300 km. The structure and behavior of these eigen- Figure 12. Modification of the ionosphere by the electric field of seismogenic origin based on the theoretical model (Rapoport et al., 2006). (a) Geometry of the model (Rapoport et al., 2006) for the determination of the electric field excited by the seismogenic current source J z (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 (T e /T e0 ), electron concentration (N e /N e0 ) and electron collision frequency (ν e /ν e0 ). Figure 13. Altitude 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 = 1000 km 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 N e 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 N e and small ν e . Lines 1 and 2 in (a) and (b) correspond qualitatively to the lines 1 and 3, respectively, in Fig. 12b, with ω = 1.5 × 10 5 s −1 and θ = 45 • . modes 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 ≤ z ≤ L max , where L z = 85 km 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 = 600 km 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 = 1000 km (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 (Collin, 2001).
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 Siefring and Kelley, 1991;Arnoldy and Kintner, 1989). 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 fixedfrequency 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 neutraldensity 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.
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.

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 (Shalashov and Gospodchikov, 2011;Kim and Kim, 2016), 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 = L max are determined through the total field components H x,y that makes the proposed procedure technically less cumbersome and practically more convenient.
2. The waveguide includes the region for the altitudes 0 < z < 80-90 km. The boundary conditions are the radiation conditions at z = 300 km, 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) ×10 5 s −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(Z 11 ) and Im(Z 11 ), in general, decrease with an increasing angle of θ. Average values of Re(Z 11 ) for θ equal to 5, 30, 45 and 60 • and Im(Z 11 ) cor-respond to θ equal to 45 and 60 • , and these increase with an increasing frequency in the considered frequency range of 0.86-1.14 ×10 5 s −1 . The average value of Im(Z 11 ) corresponds to θ equal to 5 and 30 • and changes in the frequency range of 0.86-1.14 ×10 5 s −1 non-monotonically while having maximum values around the frequency of 1-1.1 ×10 5 s −1 .
(iv) The value of finite impedance at the lower Earthatmosphere 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 ≤ z ≤ L max , where L z = 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 higherorder 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 = 600 km from the beam input, the few lowest modes of the WGEI along z coordinates still survived. Further, at x = 1000 km, 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-solarwind-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). α (0) N = 1 + ih z (k 1 − α 1 α 2 k 2 ); ih z α 2 (k 2 − k 1 ) ih z α 1 (k 1 − k 2 ); 1 + ih z (k 2 − α 1 α 2 k 1 ) ,α 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 (Al'pert, 1972;Alperovich and Fedorov, 2007;Kelley, 2009;Schunk and Nagy, 2010;Jursa, 1985).
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. The authors declare that they have no conflict of interest.
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