Model of Propagation of VLF Beams in the Waveguide Earth-Ionosphere. Principles of Tensor Impedance Method in 1 Multilayered Gyrotropic Waveguides

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


Introduction
The results of 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 waveguide "Earth-Ionosphere" (WGEI) are presented.The amplitude and phase of the VLF wave propagates in the earth-ionosphere waveguide (WGEI) can change, and these changes may be observable using ground-based and/or satellite detectors.This reflects the variations in the ionospheric electrodynamic characteristics (complex dielectric permittivity) and the influences on the ionosphere, for example, "from above" by chain Sun -Solar Wind -Magnetosphere -Ionosphere (Patra et al., 2011;Koskinen, 2011;Boudjada et al., 2012;Wu et al., 2016;Yiğit 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/Earth, such as cyclones and hurricanes (Nina et al., 2017;Rozhnoi et al., 2014;Chou et al., 2015) as well as from earthquakes (Hayakawa, 2015;Surkov and Hayakawa, 2014;Sanchez-Dulcet et al., 2015) and tsunamis.From inside the ionosphere the strong thunderstorms, lightning discharges, and terrestrial gamma-ray flashes or sprite streamers (Cummer et al., 1998;Qin et al., 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 the atmospheric physics and space plasma physics with astrophysics and high-energy physics.The corresponding "intersection area" for these four disciplines includes cosmic rays and very popular now 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 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 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 pairs "transmitted-receiver" separated by a distance of a couple of thousand km (Biagi et al., 2011;2015).The World Wide Lightning Location Network is the one more international facility 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 to the VLF plasma waves, in particular whistlers (Artemyev et al., 2013(Artemyev et al., , 2015;;Agapitov et al., 2014;Agapitov 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, 2010;Kim and Kim, 2016), our model provides an optimal balance between the analytical and numerical approaches.It combines analytical-numerical approaches basing on matrix sweep method (Samarskii, 2001).As a result, this model allows to obtain analytically the tensor impedance and, at the same time, provides high effectiveness and stability of the modeling.(ii) In distinction to the full-wave finite difference time domain models (Chevalier and Inan, 2006;Marshall et al., 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/magnetosphere.Ltym hf,jnf.C jndtnjv htwtyptynfv b ghb 'njv enjxyz.B ntrcn (iii) In distinction to the models considered in (Kuzichev and Shklyar, 2010;Kuzichev et al., 2018;Lehtinen and Inan, 2009;Lehtinen and Inan, 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, mode approach, and finally, the beam method.This combined approach, finally, creates the possibility to interpret adequately a data of both ground-based and satellite detection on 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 metal-dielectric, 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 waveguide structures metal-semiconductordielectric, 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 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 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 details in the Subsection 3.3.Our new approach, i.e a new tensor impedance method for modeling propagation of electromagnetic beams (TIMEB), includes a set of very attractive features for practical purposes.These features are: (i) the surface impedance characterizes cover layer of finite thickness, and this impedance is expressed analytically; (ii) the method allows an effective modelling of 3D beam propagating in the gyrotropic waveguiding structure; (iii) finally, 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 modelled.
Our approach was targeting properly 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 inclusion also ingoing waves; (ii) to consider a plasma-like system placed into the external magnetic field, such as the LAIM system (Grimalsky et al., 1999 a, 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 conditions on the upper boundary of the covering layer are derived.
Based on this and absence of ingoing waves, the leakage modes above the upper boundary of the structure (in other words, 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, 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, absent in the ideal metal planar waveguide without gyrotropy and anisotropy.We will search, how such an effect depends on the carrier frequency of the beam, propagating in the WGEI, inclination of the geomagnetic field and perturbations in the electron concentration, which could vary under the influences of the powerful enough sources placed "below", "above" and/or "inside" the ionosphere.
In Sect. 2 formulation of the problem is presented.In Sect. 3 the algorithm is discussed including the determination of the VLF waves/beams radiation conditions into the upper ionosphere/magnetosphere at the upper boundary, placed in the F-region at the altitude (250-400) km.The effective tensor impedance boundary conditions at the upper boundary (~ 85 km) of the effective earth-ionosphere waveguide and the 3D 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 (EM) waves with frequencies f = (10 -100) kHz can propagate along the Earth's surface for long distances >1000 km.The Earth's surface of a high conductivity z = 0 ( z is vertical coordinate) and the ionosphere F-layer 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 0 < z < 60 km, the approximately isotropic ionosphere D-layer 60 km < z < 75 km, and the anisotropic E-and F-layers of the ionosphere, due to the geomagnetic field 0 H  , added by the boundary conditions at the Earth's surface and at the F-layer.In

Algorithm
The boundary conditions and calculations of impedance and beam propagation in the WGEI is considered in this Secion.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 to the upper ionosphere/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's 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 L z ~ 250 -400 km to the lower ionosphere L z ~ 70 -90 km.All components of the monochromatic EM field is 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 , 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'YZ' where OZ' axis is aligned along the geomagnetic field 0 H  :   (Spiegel, 1959) dependent on angle  .For the medium with a scalar conductivity σ , e.g.lower ionosphere or atmosphere, the effective permittivity in (1) reduces to ε = 1 -4πiσ/ω.

The equations for the EM field and upper boundary conditions
The case of the VLF waveguide modes k x is slightly complex and should be calculated accounting for boundary conditions at the Earth's surface and upper surface of the effective WGEI.The EM field depends on the horizontal coordinate x as exp(-k x x).Taking into account k x ≤ k 0 (k 0 = ω/c), in simulations of the VLF beam propagation, we assume k x = k 0 .
Therefore, Maxwell's equations are: Here, 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: The expressions for the horizontal components of the electric field E x , E y are: 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 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 z = L max ~ (250 -400) km, where the conditions of the radiation are formulated.The dispersion equation connected the wave numbers and the frequency of the outgoing waves is obtained from Eqs.
(3), where , ~exp( ) x y z H ik z   , while the derivatives like ∂β 11 /∂z and the inhomogeneity of the media are neglected: (1 ) ( ( 1) )(1 ) ) Thus, generally Eq. ( 5) determined the wave numbers for the outgoing waves is 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 Eqs. ( 3), where the dependence on x is neglected: |k 1,2 |>> k 0 .This approximation is valid within F-layer where the first outgoing wave corresponds to the whistlers with small dissipation, the second one is the highly dissipating slow wave.To formulate the boundary conditions for Eqs.
(3a, b) at z ≥ L max , the EM field components can be presented as: In the relations ( 6) (3) are simplified in the approximation described above: The solution of Eqs. ( 7) is searched for as: , ~exp( ) 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 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 boundary conditions are rewritten in terms of H x , H y only: The relations ( 12) are the upper boundary conditions of the radiation for the boundary z=L max ~ 250 -400 km.These conditions will be transformed/recalculated using the analytical numerical recurrent procedure into equivalent impedance boundary conditions at z = L z ~ 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 0 x k  , we receive the dispersion equation using Eqs.( 5), ( 8), ( 9), in the form: x k and ' z k are the transverse and longitudinal components of wave number relative to geomagnetic field.For F-region of the ionosphere, where e He      , Eq. ( 13) reduces to the standard form of the In a special case of the waves propagating along geomagnetic field, ' 0 x k  , for the propagating whistler waves, well-known dispersion dependence is and Sagdeev, 1979).For the formulated problem we can reasonably assume Therefore eq. ( 13) is reduced to 4 , and then, similarly to the relations ( 12), the boundary conditions can be presented, in terms of the tangential components of electric field as: Conditions ( 12) or ( 14) are the conditions of radiation (absence of ingoing waves) formulated at the upper boundary z=L max and suitable for the determination of the energy of the wave leaking from the WGEI into the upper ionosphere/magnetosphere.Note, the equations ( 12), ( 14) expressed the boundary conditions of the radiation (more accurately speaking, an absence of incoming waves, what is the consequence to the causality principle), are obtained as a result of limitation by the small parameter 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 finite 0 x k k .

Equivalent Tensor Impedance Boundary Conditions
The tensor impedance at the upper boundary of the effective WGEI z=L z (see Fig. 1), is obtained by the conditions of radiation ( 12) or ( 14), recalculated from the level z = L z ~ 80 -90 km, placed in the F region of the ionosphere, at z=L max ~ 250 -400 km.
The main idea of the effective tensor impedance method is the unification of analytical and numerical approaches and 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 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.Eqs.(3) jointly with the boundary conditions ( 12), have been solved by finite differences.The derivatives in Eqs.
(3) are approximated as In Eq. ( 15) Here h is the discretization step along OZ axis; N is the total number of nodes.At each step j the difference approximations of Eqs.
(3) take the form: цhere , 1, 2,...,1, , Due to the complexity of expressions for the matrix coefficients in Eq. ( 16) we have shown them in Appendix.The set of the matrix Eq. ( 16) has been solved by factorization method also known as an elimination/matrix sweep method (see Samarskii, 2001).It can be written as: 1 ˆ, ,...,1 This method is a variant of the Gauss elimination method for the matrix 3-diagonal set of the Eq. ( 16).The value of ˆN b is obtained from the boundary conditions (12) as: Then the matrices ˆj b 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 for ˆj b follow from ( 16), (17) as: Therefore, for (17), we obtain The derivatives in Eqs. ( 4) have been approximated as: and similar equation can be obtained for ( ) ; The equivalent tensor impedance is obtained using a two-step procedure.( 1 ( 1) , 1 1 The proposed method of the transfer of the boundary conditions from the ionosphere F-layer L max = 250 -400 km into the lower part of the E-layer 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 the 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, gives 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 , H y in Eqs.
(3) seems natural, as well as the transfer (17a), because the impedance conditions are the expressions of the electric E x , E y through these magnetic components H x , H y at the upper boundary of the VLF waveguide 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 upper layer with large impedance value to lower altitude layer with 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.1.This determines the choice of the upper boundary z = L z for the effective WGEI.Due to small impedance, EM waves incident from below on this boundary are reflected effectively back.Therefore, the region 0 z z L   indeed can be presented as an effective WGEI.This waveguide includes not only lower boundary at L ISO ~ 65 -75 km with rather small losses, but also thin dissipative and anisotropic/gyrotropic layer 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 summarised as follows: (i) in contrast to invariant imbedding method currently proposed method can be used for direct recalculation of tensor impedance as it determined analytically, see Eqs. ( 22).
(ii) for the media without non-locality, proposed method does not require to solve integral equation(s).
(iii) the proposed method does not require forward and reflected waves.The conditions for the radiation at the upper boundary z = L max (see Eqs. ( 12)) are determined through the total field components H x,y , which simplify the overall calculations.
(iv) the overall calculation procedure is very effective and computationally stable.Note, that even for the very lowloss 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's use the transverse components of electric E y , and magnetic H y fields to derive equations for the slow varying amplitudes A(x,y,z), 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 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(x,y,z).Using ( 21) and ( 22), the boundary conditions are determined at the height z = L z for the slowly varying amplitudes A(x,y,z), B(x,y,z) of the transverse components E y , H y .From Maxwell's equations the components E x and H x through E y , H y in the method of beams: where  21) and ( 24), the boundary conditions for A, B can be defined as: The evolution equations for the slowly varying amplitudes A(x,y,z), 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 , ~exp( ).
x y

E H i t ik x ik y
     , Maxwell's equations are: Here In Eq. ( 27 The equations for E y , H y obtained from the Maxwell equations are: After substitution of ( 27) for E x , E z into Eqs.( 28), the coupled equations for E y , H y can be derived.The follow expansion should be used: Then, according to (Weiland and Wilhelmsson, 1977): The expansions should be until the quadratic terms of k y and the linear terms of δk x .As a result, parabolic equations (Levy 2000) for the slow varying amplitudes A and B are derived.In the lower ionosphere/atmosphere, where the effective permittivity reduces to a scalar ε(ω,z), they are independent: ( 1) In Eq. (30b), ; ; Eqs. (30b) are reduced to Eqs. (30a) when the effective permittivity is scalar.At the Earth's surface z = 0, the impedance conditions are reduced, accounting for that the medium is isotropic and the conductivity of the Earth is finite, to the form: Here σ E ~ 10 8 s -1 is the Earth's conductivity.The boundary conditions (31a) at the Earth's surface, where 022 Eqs. ( 30), combined with the boundary conditions (25) at the upper boundary of the VLF waveguide z = L z , and with the boundary conditions at the Earth's surface (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 solution of Eqs. ( 30), ( 25), (31b) are chosen in the form The size of the computing region along 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 OY axis is quite small, i.e. the terms ∂ 2 A/dy 2 , ∂ 2 B/dy 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: ˆ, 0 In the terms ˆy L C  , the derivatives with respect to y are included, whereas all other terms are included into ˆz L C  .Then the following fractional steps have been applied, the first one is along y, the second one is along z: The region of simulation is 0 < x < Lx = 1000 -2000 km, 0 < y < L y = (2000 -3000) km, 0 < z < L z = 80 -90 km.The numerical scheme (34) is absolutely stable.Here h x is the step along OX axis, x p = p h x , p = 0, 1, 2,….This step has been chosen from the conditions of the simulation results independence on the diminishing h x .On the simulation at each step along OX axis, the correction on the Earth's curvature has been inserted in adiabatic manner applying the rotation of the local coordinate frame XOZ.Because the step along x is small h x ~ 1 km << L i , this correction of the C  results in the multiplier exp(-ik 0 •δx), where δx = z•(h x /R E ), R E >> L i is the Earth's radius (see Fig. 2 and the capture to this figure).At the distances x ≤ 1000 km, the simulation results do not depend on the insertion of this correction, whereas at higher distances some 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 modes excitations by a current source such as dipole-like 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 will be possible to demonstrate that this structure has indeed good enough waveguiding properties.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 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) (12) can be used as the boundary conditions for the VLF beam radiated to the upper ionosphere/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 a beam propagation, even accounting for the wave dispersion in magnetized plasma.The proper effective boundary condition, similarly to (Rapoport et al. 2014) would allow to make relatively accurate matching between the regions, described by the full wave electromagnetic approach with Maxwell's equations and complex geometrical optics (FWEM-CGO approach).All of these materials are not included in this paper, but will be delivered in the two future papers.The first paper will be dedicated to the modeling VLF waves propagating in 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 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 TIMEB method.It was shown that more than five lowest modes of the WGEI are strongly localized in the atmosphere-lower ionosphere.Their longitudinal wavenumbers are close to the corresponding wavenumbers of EM-waves in the atmosphere.This fact convinced 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 , ε h in the coordinate frame associated with the geomagnetic field 0 H  are given in Fig. 3.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 o .The VLF frequency is ω = 10 5 s -1 , f = ω/2π ≈ 15.9 kHz.The Earth's surface is assumed as ideally conductive at the level Z = 0.The values of EM-field are given in absolute units.The magnetic field is measured in Oersteds (Oe), or Gauss (Gs), 1 Gs = 10 -4 T, whereas the electric field is also in Gs, 1 Gs = 300 V/cm.
Note that in the absolute (Gaussian) units the magnitudes of the magnetic field component  The wave beams are localized within the WGEI, 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 75 km < z < 85 km (Fig. 4b, d).These transformations depend on the permittivity component values of the ionosphere at the altitude 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 (32), we suppose that when entering the system (at X = 0), only one of the two polarization modes is excited, namely, the TM mode, i.e. at X = 0, 0; 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  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.(curves 3 in Fig. 5)

Influence of the parameters of 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, 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 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 WGEI are small, |Z0αβ| ≤ 0.1 (see, for example, Table 1).( 2) Respectively, the carrier modes of the VLF beam are close to the modes of the ideal metallized planar waveguide.These modes are subdivided into the sets of uncoupled (Ex,Hy,Ez) and (Hx,Ey,Hz) modes.The detailed search of the propagation of the separate eigenmodes of the WGEI is not a goal of this paper, and respectively, will be the subject of the separate paper.
(3) Because we have adopted for the initial beam(s) the input boundary conditions in the form (32) (with characterizes the mode coupling and corresponding transformation of the polarization at the distance 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 (Biagi et al. 2015).Other In Fig. 8a-c  ) in the considered frequency range (Fig. 8 c).The total losses increase monotonically with increasing frequency and depend weakly on the value of  (Fig. 8 d).
To model the effect of increasing and decreasing the electron concentration e n in the lower ionosphere on the polarization parameter, we have used the following parameterization for the e n change  The change in the concentration in the lower ionosphere causes rather nontrivial effect on the parameter of the polarization transformation |E y /H y |, Fig. 9 a-c.Note that either increase or decrease in the ionosphere plasma concentration have been reported as a result of seismogenic phenomena, tsunamis, particle precipitation in the ionosphere due to waveparticle interaction in the radiation belts (Pulinets et al. 2005;Shinagawa et al. 2013;Arnoldy et al. 1989;Glukhov et al. 1992;Tolstoy et al. 1986) etc.The changes in the |E y /H y | due to increase or decrease in electron concentration vary by absolute values from dozens to thousands percent, as it is seen from the comparison between Figs. 9b, 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.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 on the polarization transformation parameter minimum near the E region of the ionosphere (lines 1, 2 in Fig. 10c).The decrease of surface impedance Z 0 at the lower boundary Earth-atmosphere 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. 10 c).

Discussion
The observations presented in (Rozhnoi et al., 2015), shows 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, 1974) is also taken into account.It is assumed that the mesospheric current has only the Z-component and is positive, which means that it is directed vertically downward, as is 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. 12 b) is directed opposite to those excited by the corresponding mesospheric current (curve 1, Fig. 12 b).As a result, in the presence of both mesospheric and a seismogenic surface current, the total mesospheric electric field (curve 3, Fig. 12 b) is smaller in absolute value than in the presence of only a mesospheric current (curve 1, Fig. 12 b).It has been shown by Rapoport et al., (2006)   decrease and increase, respectively, due to the appearance of a seismogenic surface electric current, in addition to the mesospheric current (curve 3, Fig. 12, b).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 Figs. 13 a, 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. 12 b), in addition to the existing mesospheric current source (curve 2 in Fig. 12 b).This seismogenic increase 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 inside the WGEI (see Figs. 4-7).Nevertheless, the TIMEB method described by Eqs. ( 15)-( 19), (22-24), ( 27), ( 30) and allows to determine all field components in the range of altitudes max 0 z L   , where max 300 km L  .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 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 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 x=600 km from the beam input (Figs.4b, c), the few lowest modes of the WGEI along z and y coordinates still persist.At distance x=1000 km (Fig. 4d, c; Fig. 6e, f; and Fig. 7a, b), only the main mode persists in the z direction.Note, the described field structure correspond to real WGEI with losses.The gyrotropy and anisotropy causes the volume effects and surface impedance, in distinction to the ideal planar metallized 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 Is., VA and Siple Station in Antarctica (Kintner et al., 1983;Brittain et al., 1983;Siefring and Kelly, 1991;Arnoldy and Kintner, 1989), where singleaxis 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 J. W. Bonnell, UC Berkeley, NASA Grant 80NSSC18K0782).The VIPER sounding rocket campaign is consist of a summer nighttime launch during quiet magnetosphere conditions from Wallops Flight Facility, VA, collecting data through the D, E, and F regions of the ionosphere with a payload carrying the following instrumentation: 2D E-and 3D B-field waveforms, DC-1 kHz; 3D ELF to VLF waveforms, 100 Hz to 50 kHz; 1D wideband E-field measurement 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 3D EM field measurement, DC through VLF, and relevant plasma and neutral particle measurements at mid-latitudes through the radiation fields of (1) an existing VLF transmitter (the VLF transmitter Cutler with call sign NAA, the very low frequency (VLF) shore radio station at Cutler, Maine, USA, which transmits, at a frequency of 24 kHz an input power of up to 1.8 megawatts, see Fig. 11) and ( 2) naturally-occurring lightning transients through and above the leaky upper boundary of the WGEI supported by a vigorous theory and modeling effort in order to explore the vertical and horizontal profile of the observed 3D 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 B-field profile.

Figure 14. Proposed VIPER Trajectory
The aim of this experiment will be the investigation of the VLF beams launched by the near-ground source/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 (TIMEB)" in a multi-layered/inhomogeneous waveguide.The main differences/advantages of the proposed tensor impedance method in comparison with the known method of the impedance recalculating, in particular invariant imbedding methods (Shalashov and Gospodchikov, 2010;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 of integral equation(s), as the invariant imbedding method; (iii) the proposed tensor impedance method does not need the revealing the forward and reflected waves.Moreover, even the conditions of the radiation in Eq. ( 12) at the upper boundary z=L max is 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 the results, which qualitatively are in agreement with the experimental data on the seismogenic increasing losses of VLF waves/beams propagating in the WGEI.
(4) The observable qualitative effect is 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 D-and E-layers of the ionosphere at the altitudes 75 -120 km.
(5) Changes in complex tensors of both volume dielectric permittivity and impedances at the lower and upper boundaries of effective WGEI influence remarkably on the VLF losses in the WGEI.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 rather nontrivial effect on the parameter of 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 percent.In the case of decreasing electron concentration, the main maximum of |Ey/Hy | appears in the lower atmosphere at an altitude of around 20 km.In the case of 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.    ) . 1 (1 ) 1

Fig. 1 ,.
Fig.1,  is the angle between the directions of the vertical axis z and geomagnetic field 0 H 

Figure 1 .
Figure 1.The geometry of the anisotropic/gyrotropic waveguide.EM waves propagate in OX direction.0 H 


) We obtain the matrix ˆj b using Eqs.(3a, b) with the boundary conditions (12) and the procedure (17) -(19) described above.(2) Placing the expressions (21) with tensor impedance into the left parts and the derivatives in the form (20) into the right parts of Eqs.(4), the analytical expressions for the components of the tensor impedance are: the presence of gyrotropic layer and the tensor impedance boundary conditions at the upper boundary z = L z of the VLF waveguide, the equations for the slowly varying amplitudes in the general case are:

Figure 2 .
Figure 2. The rotation of the local Cartesian coordinate frame at each step along the Earth's surface h x on a small angle ≈ ≈ Δx/R E , radians, while Δx=h x .The following strong inequalities are valid h x << L z << R E .At the Earth's surface Z = 0.
Figure 3. (a) The vertical dependencies of the modules of components of the permittivity in the frame associated with the geomagnetic field | 1 |, | 3 |, | h |,, curves 1, 2, 3 correspondingly.(b) -(g) The real and imaginary parts of the components  1 ,  3 ,  h , general and detailed views.

Figure 4 .
Figure 4. Part a) is the initial distribution of |H y | at x = 0. Parts b), c) are |E y | and |H y | at x = 600 km.Parts d), e) are |E y | and |H y | at x = 1000 km.For the electric field maximum value (Fig. d) is 3•10 -6 Gs  1 mV/cm, for the magnetic field maximum value (Fig. e) is 3•10-5 Gs  3 nT 10 mV/cm.At the altitudes z < 75 km, |E z |  |H y |.

Figure 5 .
Fig. 5a, b the different dependencies of the electron concentration n(z) are shown (see solid (1), dash (2) and dot (3) lines).The corresponding dependencies of the component absolute values of the permittivity are shown in Figs.4c and 4d.

Figure 7 .
Figure 7.The dependencies of EM components on the altitude z in the center of the waveguide, y = 1500 km, for the different profiles of the electron concentration.The solid (1), dash (2), and dot (3) curves correspond to the different profiles of the electron concentration in Fig. 5, a), b), the same kinds of curves; the central plane of the beam (y=L y ) at the characteristic distance (x=x 0 ) from the beam input/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 metallized waveguide because, first, the tensor   is different from the isotropic I  only in the relatively small

Figure 8 .
Figure 8. Characteristics of the polarization transformation parameter |E y /H y | (a-c) and effective coefficient of the total losses at the distance x 0 =1000 km from the beam input ( d ); corresponding altitude dependence of the electron concentration is shown in Fig. 5b line (1); (a, b) and ( d ) -dependences of the polarization parameter (a, b) and total losses (d) on the vertical coordinate for different angles  , respectively; black (1), red (2), green (3) and blue (4) curves in Figs.a, b and d correspond to  equal to 5 ,30 ,45    and 60  , respectively; (a) and (b) correspond to frequencies Fig. a, 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 ( Figure 9. (a) Decreased and increased electron concentration (line 2, red color) and (line 3, blue) correspond to 1.25 f   and 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   

Figure 11 .
Figure 11.Averaged residual VLF/LF signals from ground-based observations at the wave paths: JJY-Moshiri, JJI-Kamchatka, JJY-Kamchatka, NWC-Kamchatka, and NPM-Kamchatka.Horizontal dotted lines show 2σ level.The color filled zones highlight values exceeding the -2σ level.In panel below Dst variations and earthquakes magnitude values are shown (from Rozhnoi et al., 2015, but not including the DEMETER data).See other details in (Rozhnoi et al., 2015).

Figure 12 .Figure 13 .
Figure 12.Modification of the ionosphere by electric field of seismogenic origin.(a) Geometry of the model (Rapoport et al. 2006) for the determination of the electric field excited by seismogenic current source ( , ) z J x y and penetrated into the ionosphere with isotropic (I ) and anisotropic (II ) regions of the system "atmosphere-ionosphere".(b) Electric field in the mesosphere in presence of the seismogenic current sources only in the mesosphere (1); in the lower atmosphere (2); both in the mesosphere and in the lower atmosphere (3).( c ) Relative perturbations caused by seismogenic electric field, normalized on the corresponding steady- boundary of the effective WGEI, all field components are (1) at least one order of altitude less than the corresponding maximal field value in the WGEI and (2) field components have the oscillating character along z coordinate and describe the modes, leaking from the WGEI.

( 6 )
An influence is demonstrated on the parameters characterizing the propagation of the VLF beam in the WGEI, in particular, on the parameter of the transformation polarization |E y /H y | and tensor impedance at the upper boundary of the effective WGEI, of the carrier beam frequency, 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 that as a result of this discretization, only the values at the grid Eqs. (21), the equations for E x , E z through E y , H y are: