ULF Activity in the Earth Environment: Penetration of Electric Field from the Near-Ground Source to the Ionosphere under Different Conﬁgurations of the Geomagnetic Field

: The problem with the penetration of electric ﬁelds from atmospheric near-Earth electric current sources to the ionosphere is investigated both within the dynamic simulations of the Maxwell equations in the frequency domain and within the simpliﬁed quasi-electrostatic approach. Two cases of the geomagnetic ﬁeld lines are considered. The ﬁrst case is the penetration of the geomagnetic ﬁeld lines deeply into the magnetosphere (open ﬁeld lines), whereas the second one is the return of these lines into the Earth’s surface (closed ﬁeld lines). The proper boundary conditions are formulated. It is demonstrated that in the case of the open ﬁeld lines the results of the dynamic simulations differ essentially from the quasi-electrostatic approach, which is not valid there. In the case of the closed ﬁeld lines, the results of simulations are practically the same both within the dynamic approach and within the quasi-electrostatic one. From realistic values of the densities of atmospheric electric currents ~0.1 µ A/m 2 , the values of the electric ﬁelds within the ionosphere F- layer may reach about 1–10 mV/m. Author Contributions: Conceptualization, methodology, V.Y., Y.R., V.G.,V.I.; software, data curation V.G., A.G., A.F.; validation, Y.R., A.G., S.P.; investigation, Y.R., V.I.; writing—original draft preparation, Y.R., V.G., A.G., V.I.; writing—review and editing, V.Y., S.P., A.F., V.K.; discussion, interpretation of physical results and conclusions, all authors; visualization (preparation of ﬁgures), Y.R., S.P., V.K.; supervision, V.Y.,Y.R.; project factorization method, the correct approximation of the boundary conditions for different cases. The description of the computer codes is presented in the text of the paper in Appendices B and C, as well as in proper references [48–50].


Introduction
The problem of the lithosphere (Earth)-atmosphere-ionosphere-magnetosphere (LEAIM) coupling has been investigated for a long time  with special attention to the preparation of natural hazards and their ionospheric effects. In particular, this concerns earthquakes (EQs), tropical cyclones, and volcanic eruptions [1][2][3].
For several decades, experimental results have been accumulated and recently generalized, including targeted observations over a number of years, which unambiguously indicate the presence of seismogenic disturbances in TEC, VLF, low frequency (LF) and ultra-low frequency (ULF) ranges of electromagnetic field, surface temperature, outgoing infrared (IR) radiation, chemical potential in the lower atmosphere, ionization anomalies, etc. [2,[23][24][25][26][27]. As it was mentioned above, a set of mechanisms of coupling, in particular seismo-ionospheric one, in the LEAIM system has been proposed; in particular electrostatic [28,29], electro-photochemical [5,30], electromagnetic [15,16], mechanical (via Acoustic-Gravity Waves (AGW)) [31] and mixed mechanical-electrical [1,15,24,25,27,32,33] mechanisms have been discussed; at the same time, most likely, that in any of the above options, the observed physical phenomena are provided, accounting for the achieving by the corresponding perturbations the amplitudes, significantly exceeding the noise level, by synergetic processes in active open nonlinear LEAIM system [6,9,15,18,33,34]. Only a synergetic approach allows one to adequately, and therefore necessarily self-consistently, determine the electric field in the lower atmosphere. Moreover, for this, it is necessary to take into account a variety of physical processes, the development of instabilities and their nonlinear limitations [9] in different regions and at different heights in the LEAIM system, as well as boundary conditions, including in the lower atmosphere [15,34]. In this work, as in [15], we do not apply such a self-consistent synergetic approach and focus on the questions formulated below about comparing the dynamic and quasi-electrostatic approaches, and for situations with different configurations of the geomagnetic field. In this work, as in [15], instead of a self-consistent consideration of the motion and dynamics of charged and uncharged components of a plasma-like medium, an external current in the lower atmosphere is specified. Hereinafter, an external current is understood as a current that characterizes a given source, which is not determined by a self-consistent calculation based on the properties of the medium in which the excitations are investigated, but "assigned" and in this sense is "external" in relation to the given system, in which perturbations are sought [35]. The consequence of this is the inevitable rejection of any claims for an adequate determination of the components of the electric field in the region where the external current is specified, i.e., in the lower atmosphere. To emphasize this, in the graphs of the vertical distribution of the electric field components given in Section 3 (Figures 5,7,and 8), an empty space is left in the region of the lower atmosphere. The external current only simulates the real current generated in the lower atmosphere by all the physical processes taking place in a self-consistent manner, but at the same time, it makes it possible to determine the field penetrating the upper ionosphere/magnetosphere.
In the present paper, we are concentrated on electromagnetic/quasi-static channel of coupling, in particular seismogenic, in the system LEAIM. We consider the electromagnetic/ quasi-static perturbations in the ionosphere as a sensitive indicator of the coupling in the system mentioned above. It is interesting to note the electric fields of what orders of magnitude are observed in different regions of the atmosphere-ionosphere and at what altitudes. In the papers [29,36,37] the measurements of the vertical electric field variations with amplitudes of order (1-few) kV/m on the ground level in the region of the preparation of strong earthquakes have been reported. In [7,38] the observations of ionospheric inhomogeneities using MF-HF (Medium Frequency-High Frequency) radar have been made; the results correspond to the presence of an effective source of a mesospheric current or an electric field with amplitude (1-few V/m) [7,38,39]. It is interesting that, as it follows from the results of the numerical modeling, electric fields with amplitude of the same order (1-few V/m) in the mesosphere correspond to the presence of the (seismogenic) vertical electric field on the lithosphere-atmosphere boundary of order (1-few) kV/m [4,36,37,39]. AIP (Advanced Ionospheric Probe) on board of FORMOSAT-5 [40] provided measurements, in particular, of ion velocities. The amplitude of a seismogenic electric field in the ionosphere has been evaluated as a value of order 1 mV/m. Seismogenic plasma bubbles observed In Equation (1), → j ext is the density of the external current,ε(ω, z) is the dielectric permittivity of the medium. The density of the external current and the field components depend on the frequency and all coordinates x, y, z. It is assumed that the properties of the ionosphere and the lower magnetosphere depend only on the vertical coordinate associated with the Earth's surface, while in the upper magnetosphere they depend on the coordinate associated with the direction of the geomagnetic field. In Equation (1) ω ω ν ω ω ν ε ω ν ω ω ω ν ω ω ω ω ε ε ω ν ω ω ν ω ω ω ω ω ω ν ω ω ω ν ω ω The following notations are used in relations (3):  Figure 4). Hereafter, for convenience, we omit the designation of functional dependency on the height and frequency for components of dielectric permittivity tensor. The parameters of the atmosphere and the ionosphere are assumed as unchanged during the observation time. The characteristic time of essential change of the atmosphere and ionosphere parameters is of about several (3)(4)(5) hours, respectively. Frequency domain presentation is used in our paper approximately for the evaluations of the effects of electric field penetration into the ionosphere, by the order of value. We restrict ourselves by the consideration of the penetration of the electromagnetic In the coordinate frame X'Y'Z' associated with the geomagnetic field → H 0 the dielectric permittivity tensor has the form [51,52]: where The following notations are used in relations (3): In (3), (4) M i , m e , n e0 , ω pe,pi , ω He,Hi , ν e,i , H 0 are the masses of the positive ion and electron, the concentration of electrons, the plasma frequencies of electrons and positive ions, the cyclotron frequencies of electrons and positive ions, the collision frequencies for electrons and ions, and the value of the geomagnetic field, respectively. These values depend on the height i.e., on the local vertical coordinate z. Corresponding dependencies used in modelling are shown in Section 3 ( Figure 4). Hereafter, for convenience, we omit the designation of functional dependency on the height and frequency for components of dielectric permittivity tensor. The parameters of the atmosphere and the ionosphere are assumed as unchanged during the observation time. The characteristic time of essential change of the atmosphere and ionosphere parameters is of about several (3)(4)(5) hours, respectively. Frequency domain presentation is used in our paper approximately for the evaluations of the effects of electric field penetration into the ionosphere, by the order of value. We restrict ourselves by the consideration of the penetration of the electromagnetic (EM) fields from given external current sources and under given ionosphere parameters. For the other parameters, the results concerning the penetration of the electric field into the ionosphere would change, respectively.
In the coordinate frame associated with the Earth's surface, the expression for the dielectric permittivity tensor is obtained, using (2) by rotating the coordinate frame in the XOZ plane:ε Atmosphere 2021, 12, 801 5 of 31 In relation (5), ε ij ; i, j = 1, 2, 3 are the components of the dielectric permittivity of the ionospheric plasma.
When deriving the relations for the atmosphere and the lower part of the ionosphere, the Cartesian coordinate frame is used, where the OZ axis is directed vertically upwards ( Figure 1). When considering the propagation of electromagnetic fields in the upper atmosphere and magnetosphere, a set of local Cartesian coordinates associated with geomagnetic field lines within the magnetosphere is used. The approximation is that the coordinate frame is considered to be locally Cartesian, with environmental parameters that change slowly depending on the height or distance along the geomagnetic field line. The local orientation of the OZ axis is directed along the geomagnetic field line. This approximation, which is actually adiabatic, can be used due to the high anisotropy of the conductivity of the ionosphere and magnetosphere. In this case, the ULF magnetohydrodynamic (MHD) waves mainly propagate along the lines of the geomagnetic field.
Also, for wave perturbations in the ULF range, an attempt has been done to investigate the penetration of an electric field within the framework of a quasi-electrostatic approximation by means of an electric potential ϕ ∼ exp(iωt).
The fundamental difference between quasi-stationary equations (6) and the corresponding electrostatic approximation equations, which were mostly used earlier, including some of our previous works, is the presence in the last equation of system (6), of the second quasi-static, term in the left part, proportional to the frequency ω. This will allow us to investigate for the first time the possibility of the above-mentioned boundary transition between the corresponding results for electric fields using dynamic and quasi-static approximations. The following periodic zero boundary conditions for the domains along x and y are used for the numerical solution of the equations: 0 < x < L x , 0 < y < L y . Accordingly, the fast Fourier transform (FFT) is used, and the corresponding modes are represented as follows: E, H, ϕ∼ exp(−i( k x x + k y y)).
In (7), k x,y are discrete wavenumbers of the corresponding Fourier modes, and for E, H, ϕ the amplitudes of the corresponding Fourier modes are meant, and the indices denoting the corresponding modes are omitted at the Fourier amplitudes. Then the system of two equations is obtained similarly to [15] for the column → F ≡ (E x , E y ) T . Here the superscript "T" means transposition, the Fourier components E x,y depend on the local vertical coordinate z. This system is written in the matrix form as follows: The matricesÂ 1 ,B 01 ,B 02 ,B F ,B 2J z the vectors → B 1J z , → e F , and the values B J tg = − 4πik 0 c , ε 33t included in dynamic Equation (8) are presented in Appendix A (see Equations (A12)-(A17)) along with other details of the derivation of this equation. The dependencies of the matrix and vector coefficients on ω, z are omitted in Equation (8).
It is assumed in the numerical modeling that the source of external current in the lower atmosphere has only vertical component ( Atmosphere 2021, 12, 801 6 of 31 In relation (9), j 0 , r 0 , z 1,0 are the amplitude of the vertical component of the external current, the characteristic scale of the external current in the horizontal plane, the altitude position of the external current maximum, the characteristic scale of the external current in the vertical plane, respectively. Denote L x,y the distance between the imaginary side walls (Figure 1), which determines the periodicity of the system in the horizontal plane. Note that the boundary condition with a zero value of the tangential electric field on the imaginary sidewalls correspond to the case of a strongly localized value of the current density j z , and the condition must be satisfied We illustrate two physically different situations mentioned above in Figure 2. The first physical situation is with geomagnetic lines which go to the magnetosphere and do not return to the Earth (open geomagnetic field lines, Figure 2a). Boundary conditions are formulated in the magnetosphere at z = L z ; these are the conditions of radiation and the absence of waves propagating in the reverse direction to the ionosphere from the magnetosphere. These boundary conditions are equivalent to a system of linear relations between dE x /dz, dE y /dz, and E x , E y analogously to [15]: (see relation (A18)). Nevertheless, we are interested in the comparison between the results based on the dynamic and quasi-static approaches for the determination of the electric field penetration through the system LEAIM. So, it would be quite reasonable to consider the boundary conditions (11), (12) using the more exact approach and accounting for the existence of effectively at least two small parameters in the theory describing ULF waves in the upper ionosphere and magnetosphere. One of these parameters describes anisotropy and is small in the F region of ionosphere and magnetosphere and is equal to Note that the matrixm included in the left side of Equation (12) consists of elements m ij ; i, j = 1, 2 that are included in the left sides of Equation (11).
The coefficients m ij , i, j = 1, 2 included in system (11) (or, (12)) are presented in Appendix A, see relations (A38). In Appendix A, the derivation of the expressions for m ij , i, j = 1, 2 mentioned above has been realized based on the approximation ε 3 → ∞ (see relation (A18)). Nevertheless, we are interested in the comparison between the results based on the dynamic and quasi-static approaches for the determination of the electric field penetration through the system LEAIM. So, it would be quite reasonable to consider the boundary conditions (11), (12) using the more exact approach and accounting for the existence of effectively at least two small parameters in the theory describing ULF waves in the upper ionosphere and magnetosphere. One of these parameters describes anisotropy and is small in the F region of ionosphere and magnetosphere and is equal to where In (14), L ⊥ is the characteristic size of the transverse distribution of the (horizontal) electric field components; the determination of k t corresponds also to one presented after Equation (A10) in Appendix A. The other parameter is |ε 1 |(k 2 0 /k 2 t ), where k 0 = ω/c, and is proportional to ω, when ω is relatively small (see relation (3)).
As it follows from the accurate analysis of the Equations (A25)-(A27), the dispersion relation for Alfven waves (AWs) takes the form: In Equation (15), Accounting for (13), one can get, using (15), (16), more accurate expression than (A29) for AW mode: The dispersion relation for fast magnetosonic waves (FMSW) mode has nevertheless the same form as the one presented in Equation (A30): Consider two cases, corresponding to the different ratios between the two parameters mentioned above. Namely, the first of them describes the anisotropy (see relation (13)). The second one is |ε 1 |(k 2 0 /k 2 t ) and is small and proportional to ω for relatively small frequency, as mentioned above.
(1) In the first case, suppose that both parameters are small, and the condition is satisfied. Under the condition (19), Equations (17) and (18) reduce to Accounting for (20) and applying the approach similar to the one used in Appendix A for the derivation of the relations (A38), one can get, by the order of value, the estimation In relation (21), ∆z e f f is the penetration depth of the MHD waves penetrating into the higher ionosphere/magnetosphere from some altitude z where the (horizontal) electric field distribution has the typical transverse size of order L t . (2) In the second case, suppose that, instead of (19), the second parameter |ε 1 |(k 2 0 /k 2 t ) is of order 1, the following condition is satisfied: In this case, the relations (17), (18) reduce in accordance with the relations (A29), (A30), to the form Let us make the estimation for z = 200 km, ω ∼ 2.5 · 10 −2 s −1 , |ε 1 | ∼ 10 7 , L t ∼ 4000 km (see Figures 4e and 5a,b, respectively, shown below in the Section 3). Using the relations (23), (A38), one can get The second physical situation is when geomagnetic lines return to the Earth's surface (closed geomagnetic field lines, Figure 2b). In this situation, the boundary conditions are E x = E y = 0 on the Earth's surface (or ϕ = 0, for a quasi-electrostatic approximation; in the case of open geomagnetic field lines in this approximation, the same boundary condition applies to the magnetosphere). Such boundary conditions on the Earth's surface correspond to the infinite conductivity of the Earth. The consideration of the finite conductivity of the solid Earth weakly changes the corresponding results, as our simulations have demonstrated.
For the solution of set of Equations (8) the matrix run method, or matrix factorization method [48][49][50], is used. The corresponding boundary condition (11) (or (12)), as well as zero conditions for xand y-components of the electric field on the side walls (see Figure 1) and the lower boundary z = 0 are taken into account. This method is described in Appendix B. When the transverse sizes of the simulation region are quite large, the results of simulations are the same both for the zero boundary conditions at the sidewalls and for the periodic boundary conditions. Because in the upper ionosphere and the magnetosphere the electromagnetic field propagates along the geomagnetic field lines, which are curved, the adiabatic rotation of the locally rectangular coordinate frames (c.f.) has been applied, Figure 3. Namely, near the Earth's surface OZ axis of c.f. is aligned vertically, whereas OX, OY ones are tangential to the Earth's surface. Starting from the altitude z = z 2 ≈ 300 km, the adiabatic rotation of OZ, OX axes takes place, and at the altitude z = z 3 ≈ 600 km OZ axis becomes aligned along the geomagnetic field line. The adiabaticity is due to the small rotation angle at each simulation step along OZ axis h z < 1 km. In the case of open geomagnetic field lines, the OZ axis of the local c.f. preserves this orientation within the magnetosphere. In the case of the closed geomagnetic field lines, the inverse procedure has been realized near the return point, so in the return point at the Earth's surface, OZ axis is aligned vertically. It has been checked that the simulation results depend weakly on the positions of the points z 2 , z 3 , axis of the local c.f. preserves this orientation within the magnetosphere. In the case of the closed geomagnetic field lines, the inverse procedure has been realized near the return point, so in the return point at the Earth's surface, OZ axis is aligned vertically. It has been checked that the simulation results depend weakly on the positions of the points z2, z3, z4, z5.  The necessity of the application of the dynamic simulations for the problem of the penetration of the electric fields into the ionosphere is due to the following factors.
1. The atmosphere, ionosphere, and magnetosphere jointly with the current sources are principally dynamic. Especially, the increase of the electric fields in the atmosphere lasts for several minutes. The properties of the ionosphere change for several hours. The conductivity of the air near the Earth's surface is low: σ < 10 −3 s −1 , so it is necessary to take into account the influence of the nonzero frequency there, if ω > ω 00 ; ω 00 ∼ 3 · 10 −3 s −1 .
2. When the boundary condition for the electric potential is in the deep magnetosphere L z = 25,000 km, the deep magnetosphere itself is a dynamical system due to the Earth's rotation around the Sun. During 1 hour of rotation, the Earth moves >100,000 km in space, so the configuration of the geomagnetic field changes.
3. In the set of papers, in particular in [45], it was shown that there are definite restrictions for the applicability of the electrostatic approximation. Further transition to the quasi-static model [46] confirms the importance of accounting for the non-stationary character of the processes in the near-Earth space. In [46], it was stated that the cases of open and close (see Figure 2 in that paper) geomagnetic field lines should be treated differently. As it is shown in Appendix C for the open geomagnetic field lines, the upper boundary conditions (relations (A53)) have the structure similar to those of the upper boundary conditions [47].
In papers [17,45,47] there were estimations of limited validity of the electrostatic approximation. Then a lot of authors used a more correct term 'the quasi-statics' instead of 'the electrostatics', for instance [46] and references there.
In the paper [46], there is a clear understanding that under different configurations of the geomagnetic field lines the approaches to simulations on the penetration of the electric field with the quasi-statics should be different, especially with respect to the boundary condition from above, see for instance Figure 2 in [46].

The Results of Simulations
The main goal is to estimate the penetration of the electric field components to the Eand F-layers of the ionosphere under different geometries of the geomagnetic field lines.
The values of the parameters for simulations have been selected due to the considered problem as follows. As for the height z = z 1 , at which the maximum of the external current is placed, we consider two examples. In the first of them, the current maximum is somewhat raised directly above the surface region. This study is of interest from several points of view. First, the external current in the model adopted in this work qualitatively corresponds to the effective sources associated with synergetic processes. From the point of view of the prospects of the above-mentioned synergetic approach, one should take into account the presence of a number of different instabilities at different heights in the "lithosphereionosphere" system, including the lower atmosphere, as well as its surface layer [42,43], the mesosphere [39], E [53] and F [1,15,33] regions of the ionosphere. Therefore, it is of interest to consider the case when the external current is raised directly above the surface layer. It should be mentioned that powerful surface sources causing a noticeable response in the ionosphere include also tropical cyclones and typhoons [1,15,33,[53][54][55][56][57]. The kinetic energy carried by tropical cyclones is comparable to the energy of powerful earthquakes, so the tropical cyclones are among the most destructive large-scale atmospheric formations on our planet [54,56,57]. The impact of a tropical cyclone on the ionosphere is carried out by means of a mixed mechanism, including enhanced ionization, electrophysical, hydrodynamic, meteorological, and other processes. As a result of these processes, the self-organized powerful atmospheric formation is formed, namely, an atmospheric vortex. Such a vortex is a known characteristic feature of tropical cyclone. This formation has a horizontal scale of the order of hundreds to thousands of kilometers [54]. Such a formation includes structures with strong convection and clouds, so-called convective towers [58], where both forming AGWs and intense currents and fields arising due to charge separation are clearly manifested, in particular, in lightning discharges. Moreover, such electrically saturated convective structures reach heights of 16 km [59,60]. The resulting electric fields, together with AGW, affect the ionosphere, causing noticeable disturbances of the total electron content (TEC) [58], variations in the penetration frequency foF2 [61], variations in the characteristics of VLF waves propagating in the waveguide Earth-Ionosphere and other phenomena, some of which are considered now, in the processes of tropical cyclone/hurricane formation, as the corresponding precursors [62].
In accordance with the pointed above features, for the first part of simulations, it has been adopted that the electric current source is localized in the lower atmosphere, z 1 = 14 km, above the Earth's surface. The maximum value of the density of electric current j 0 = 0.1 µA/m 2 , [54,56], the transverse and vertical scales r 0 = 500 km [54,56], and z 0 = 7 km [59,60] are adequate parameters, characterizing the external current source for the evaluations by the order of value.
After the calculation for a source in the lower atmosphere with an elevated maximum (z 1 = 14 km), a number of calculations will also be performed for different positions of the maximum of the external current source along with the height in the region 0 ≤ z 1 ≤ z 01 , z 01 = 14 km. Among these calculations, some will be done for the vertical coordinate of the current source maximum, which, altogether with other parameters, would qualitatively correspond to the seismogenic source. Also, calculations based on the dynamic and quasielectrostatic approaches will be performed for some frequencies in the ULF range, with the justification for the choice of these values.
The parameters of the atmosphere, ionosphere, and magnetosphere used in the simulations for different values of height are presented in Figure 4. The parameters of the ionosphere used for modeling are taken from [18,51,52,[63][64][65]. The simulations are tolerant to changes in the parameters and mainly yield qualitative results.
correspond to the seismogenic source. Also, calculations based on the dynamic and quasielectrostatic approaches will be performed for some frequencies in the ULF range, with the justification for the choice of these values.
The parameters of the atmosphere, ionosphere, and magnetosphere used in the simulations for different values of height are presented in Figure 4. The parameters of the ionosphere used for modeling are taken from [18,51,52,[63][64][65]. The simulations are tolerant to changes in the parameters and mainly yield qualitative results.  The choice of the frequencies of the external current source located in the lower atmosphere is due to the following considerations.
(1) We compare calculations based on dynamic and quasi-electrostatic models of electric field penetration into the ionosphere and magnetosphere. In particular, we are interested in the lower frequency limit. The frequency ω = 2.5 · 10 −2 s −1 chosen for the calculations determines the lower frequency limit that we can practically reach at the moment in the calculations on the basis of the dynamic model. This limit is determined by the fact that decreasing model frequency two times requires a fourfold increase in Random Access Memory (RAM), during the calculations based on the dynamic approach. Despite this, we qualitatively investigate the penetration of the field into the ionosphere and magnetosphere also at frequencies ω < 2.5 · 10 −2 s −1 for the quasi-electrostatic approach.
As shown below, the application of the quasi-electrostatic approach is adequate for the case of closed geomagnetic field lines. Moreover, the goals of this work include clarifying the qualitative nature of the dependences of the field penetration on some parameters of the source, including its frequency. The calculation results in the dynamic and quasielectrostatic approximations are compared in the region of minimum frequencies ω ∼ 2.5 · 10 −2 s −1 (namely, for ω = 2.5 · 10 −2 s −1 and ω = 5 · 10 −2 s −1 ), in which our computational capabilities still allow us to apply not only the quasi-electrostatic approximation, which does not have, with the mathematical point of view, frequency restrictions "from below", but also a dynamic approach.
As we will demonstrate in this section, in the aforementioned boundary frequency range, where both dynamic and quasi-static approaches are applicable, they, as expected, with appropriate calculations, give for the maxima of the field components the values of the same order of magnitude. In our qualitative studies of the electric field penetration from the lower atmosphere into the ionosphere and magnetosphere, the accuracy of determining the field components in order of magnitude is quite satisfactory.
(2) Let us take into account the existence of a correspondence, as will be shown below, in order of magnitude and in qualitative behavior depending on the height, between the results obtained in the case of closed geomagnetic field lines for the maximum field value using dynamic and quasi-electrostatic approaches. This allows, by matching the results for these two approaches in the above-designated frequency range, where they are both applicable, to "go down" in frequency on the basis of a quasi-electrostatic approach. Thus, it is possible to study the field penetration from a given (external) current source, with an accuracy of an order of magnitude, in the lower part of the ULF range.
The typical results of simulations are presented in Figures 5-10. The frequency interval is ω = 0.025-0.05 s −1 , where there is a possibility to apply the quasi-electrostatic approximation. The values of the electric field components are in mV/m units. The results of the simulations are of qualitative character and are tolerant to changes in the parameters of the ionosphere. The region of simulations in the locally horizontal plane XOY has been chosen quite large to avoid an influence of boundaries on the electric field distributions, so it should be L x,y ≥ 10,000 km. The initial inclination of the geomagnetic field from the normal to the Earth's surface is θ i = 10 • . The results depend weakly on this parameter when θ i < 30 • .
In the case of open field lines, the results of simulations are presented in Figures 5-7 and 13. Within the dynamic simulations there exists the penetration of the horizontal electric field components to the heights z > 150 km. These dependencies are presented at the heights z ≥ 50 km only because at lower heights the values of the electric field, especially the vertical component E z , cannot be computed correctly in the regions where the external current source j ext is not zero [15]. The used approximation of the given external current sources is too rough there [15]. Within the dynamic approach, the simulated penetration does not depend on the position of the upper boundary condition L z , when L z ≥ 1000 km, see Figure 5. Specifically, the dependencies for the upper boundary condition at L z = 25,000 km coincide with ones for L z = 1500 km, curves 3 in Figure 5. Also, the values of the electric field components depend slightly on the frequencies at ω ≤ 0.1 s −1 . This can be seen from Figure 6 where the spatial distributions of the electric field components are presented at z = 200 km in F-layer; compare  (24), in the magnetosphere at the vertical distances z > 1000 km, the maxima slowly decrease with the growth of z, as our simulations have demonstrated, but the detailed investigation of the penetration of the electromagnetic fields into the magnetosphere is out of scope of this paper.
Contrary to this, within the quasi-electrostatic approximation the simulated penetration is very poor in the ionosphere F-layer, 3 orders smaller than one within the dynamic simulations. The results of the quasi-electrostatic simulations are given in Figure 7. Two types of boundary conditions within the magnetosphere have been applied, namely either ϕ = 0 or E z = 0. The results of the simulations do not depend on the types of the boundary conditions. Note that the upper condition E z = 0 corresponds to the absence of the electric current along the geomagnetic field line [46]. The simulated values of the electric field depend on the position of the upper boundary condition. For the curve 1 the value of L z is L z = 1000 km, for curves 2, 3 they are 1500 km and 25,000 km (Figure 7). At the distance L z = 25,000 km along the geomagnetic field line, i.e., in the deep magnetosphere, the parameters of the magnetosphere are assumed as n e0 = 10 2 cm −3 , M/m e = 10 4 , ω He = 10 3 s −1 , ν e = 10 s −1 , ν i = 10 −4 s −1 , that yields the components of the conductivity σ 1 ≈ 3.1·10 5 + i6.8·10 6 , σ 3 ≈ 2.6·10 9 − i1.1·10 8 , σ h ≈ −1.7·10 6 + i2.9·10 4 . The results of the quasi-electrostatic simulations depend weakly on changing the parameters in the deep magnetosphere.
tion is very poor in the ionosphere F-layer, 3 orders smaller than one within the dynamic simulations. The results of the quasi-electrostatic simulations are given in Figure 7. Two types of boundary conditions within the magnetosphere have been applied, namely either ϕ = 0 or Ez = 0. The results of the simulations do not depend on the types of the boundary conditions. Note that the upper condition Ez = 0 corresponds to the absence of the electric current along the geomagnetic field line [46]. The simulated values of the electric field depend on the position of the upper boundary condition. For the curve 1 the value of Lz is Lz = 1000 km, for curves 2, 3 they are 1500 km and 25,000 km (Figure 7). At the distance Lz = 25,000 km along the geomagnetic field line, i.e., in the deep magnetosphere, the parameters of the magnetosphere are assumed as ne0 = 10 2 cm −3 , M/me = 10 4 , ωHe = 10 3 s −1 , νe = 10 s −1 , νi = 10 −4 s −1 , that yields the components of the conductivity σ1≈ 3.1·10 5 +i6.8·10 6 , σ3 ≈ 2.6·10 9 -i1.1·10 8 , σh ≈ −1.7·10 6 +i2.9·10 4 . The results of the quasi-electrostatic simulations depend weakly on changing the parameters in the deep magnetosphere.  The explanation of the great difference between the results of the dynamic and the quasi-static approximations (compare Figures 5a-c and 7a-c, respectively) is as follows. The quasi-electrostatic approximation can be applied when the following condition of the small retardation is valid: Lz << vA/ω [66], where vA = c/ε1 1/2 ≈ H0/(4πne0M) 1/2 is the Alfven velocity [51,52]. For the rough estimations we accept the condition of the validity of quasi- The explanation of the great difference between the results of the dynamic and the quasi-static approximations (compare Figures 5a-c and 7a-c, respectively) is as follows. The quasi-electrostatic approximation can be applied when the following condition of the small retardation is valid: L z << v A /ω [66], where v A = c/ε 1 1/2 ≈ H 0 /(4πn e0 M) 1/2 is the Alfven velocity [51,52]. For the rough estimations we accept the condition of the validity of quasi-electrostatic approximation in the form: It is of about v A ≈ 400 km/s within the upper ionosphere and v A ≈ 100 km/s in the magnetosphere [51,52]. Thus, it should be L z < 2000 km. Such a condition could not be satisfied for L z = 25,000 km, see curves 3 in Figure 7a In the case of closed geomagnetic field lines, the horizontal components of the electric field penetrate into the ionosphere both within the dynamic simulations and in the quasielectrostatic approximation. The results of simulations are given in Figures 8-12. The distance along the geomagnetic field line is Lz = 2000 km, i.e., it includes two parts, from the Earth's surface to the magnetosphere and back from the upper point in the magnetosphere to the Earth's surface. The values of the electric field have a correspondence in these simulations. The quasi-static approximation is approximately valid, when the total distance along the geomagnetic field line is of about 2000 km between the Earth's surfaces. The quantitative difference is due to both some non-quasi-stationarity and the roughness of the model. In this case, to get more exact results it is necessary to take into account in detail the curvature of the geomagnetic field lines under the quasi-static approximation. Under dynamic simulations, the influence of this curvature is not so essential, as our simulations have been demonstrated.  Figures 7, 8 and 13, corresponds to the rejection of claims for an adequate definition of the electric field in the external current region, introduced instead of an synergetic self-consistent definition of current and field, as explained in detail in the introduction to this work, as well as in [15].
In the case of closed geomagnetic field lines, the horizontal components of the electric field penetrate into the ionosphere both within the dynamic simulations and in the quasielectrostatic approximation. The results of simulations are given in Figures 8-12. The distance along the geomagnetic field line is L z = 2000 km, i.e., it includes two parts, from the Earth's surface to the magnetosphere and back from the upper point in the magnetosphere to the Earth's surface. The values of the electric field have a correspondence in these simulations. The quasi-static approximation is approximately valid, when the total distance along the geomagnetic field line is of about 2000 km between the Earth's surfaces. The quantitative difference is due to both some non-quasi-stationarity and the roughness of the model. In this case, to get more exact results it is necessary to take into account in detail the curvature of the geomagnetic field lines under the quasi-static approximation. Under dynamic simulations, the influence of this curvature is not so essential, as our simulations have been demonstrated.
ulations. The quasi-static approximation is approximately valid, when the total distance along the geomagnetic field line is of about 2000 km between the Earth's surfaces. The quantitative difference is due to both some non-quasi-stationarity and the roughness of the model. In this case, to get more exact results it is necessary to take into account in detail the curvature of the geomagnetic field lines under the quasi-static approximation. Under dynamic simulations, the influence of this curvature is not so essential, as our simulations have been demonstrated.  As it is seen from Figure 9, parts a and d, b and e, c and f, in the case of closed geomagnetic field lines, the differences in corresponding maximum values of electric field components for the frequencies 0.05 s −1 and 0.025 s −1 does not exceed a value of order 25%.
From Figures 5 and 8 it is seen that the dynamically computed values of the electric field in the ionosphere F-layer are a bit higher in the case of closed geomagnetic field lines. This can be explained by a partial reflection of the magnetohydrodynamic waves from the Earth's surface in the place of returning the geomagnetic field lines. As it is seen from Figure 9, parts a and d, b and e, c and f, in the case of closed geomagnetic field lines, the differences in corresponding maximum values of electric field components for the frequencies 0.05 s −1 and 0.025 s −1 does not exceed a value of order 25%.
From Figures 5 and 8 it is seen that the dynamically computed values of the electric field in the ionosphere F-layer are a bit higher in the case of closed geomagnetic field lines. This can be explained by a partial reflection of the magnetohydrodynamic waves from the Earth's surface in the place of returning the geomagnetic field lines.
As it is seen from Figure 9, parts a and d, b and e, c and f, in the case of closed geomagnetic field lines, the differences in corresponding maximum values of electric field components for the frequencies 0.05 s −1 and 0.025 s −1 does not exceed a value of order 25%.
From Figures 5 and 8 it is seen that the dynamically computed values of the electric field in the ionosphere F-layer are a bit higher in the case of closed geomagnetic field lines. This can be explained by a partial reflection of the magnetohydrodynamic waves from the Earth's surface in the place of returning the geomagnetic field lines. Our simulations have been demonstrated that for the quasi-electrostatic simulations when the inclination angle is θi > 10°, the distributions become smeared in the horizontal plane. This numerical effect does not change the results qualitatively; it can be removed accounting for the curvature of the geomagnetic field lines, using a curvilinear coordinate Our simulations have been demonstrated that for the quasi-electrostatic simulations when the inclination angle is θ i > 10 • , the distributions become smeared in the horizontal plane. This numerical effect does not change the results qualitatively; it can be removed accounting for the curvature of the geomagnetic field lines, using a curvilinear coordinate frame. Also, the tuning of the locally horizontal region of simulations in the upper ionosphere and the magnetosphere should be realized. But there is a simpler solution to the pointed above problem [46]. Since the longitudinal conductivity, that is, the component along the geomagnetic field line, is 5-7 orders of magnitude greater than the transverse conductivity in the upper ionosphere and magnetosphere, it is possible to make a connection by short-circuiting two regions of the ionosphere, namely the one from which the geomagnetic field line comes out and the one into which this line enters, i.e., magnetically conjugated region. Thus, the corresponding ionospheric regions located at altitudes up to z = 300 km can be connected. In other words, in the quasi-electrostatic simulations, it is possible to avoid the simulations within the upper ionosphere and the magnetosphere. Of course, this procedure is valid only in quasi-electrostatic simulations in the case of closed geomagnetic field lines. In this simplified approach the results of simulations are of agreement with the dynamic simulations under greater inclination angles θ i > 20 • . The results of simulations at higher inclination angles are given in Figures 11 and 12, the quasi-electrostatic and dynamic simulations correspondingly.
Atmosphere 2021, 12, x FOR PEER REVIEW 17 of 33 frame. Also, the tuning of the locally horizontal region of simulations in the upper ionosphere and the magnetosphere should be realized. But there is a simpler solution to the pointed above problem [46]. Since the longitudinal conductivity, that is, the component along the geomagnetic field line, is 5-7 orders of magnitude greater than the transverse conductivity in the upper ionosphere and magnetosphere, it is possible to make a connection by short-circuiting two regions of the ionosphere, namely the one from which the geomagnetic field line comes out and the one into which this line enters, i.e., magnetically conjugated region. Thus, the corresponding ionospheric regions located at altitudes up to z = 300 km can be connected. In other words, in the quasi-electrostatic simulations, it is possible to avoid the simulations within the upper ionosphere and the magnetosphere. Of course, this procedure is valid only in quasi-electrostatic simulations in the case of closed geomagnetic field lines. In this simplified approach the results of simulations are of agreement with the dynamic simulations under greater inclination angles θi > 20°. The results of simulations at higher inclination angles are given in Figures 11 and 12, the quasi-electrostatic and dynamic simulations correspondingly.  (a) (b) (c) Figure 11. Distributions of the electric field components (in mV/m) in the ionosphere F-layer at the height z = 200 km. The quasi-electrostatic simulations, the case of closed geomagnetic field lines. The frequency is ω = 0.025 s −1 . It is increased inclination angle compared with Figure 8. Part (a) is |Ex(x,y)|, (b) is |Ey(x,y)|, (c) is|Ez(x,y)|. The geomagnetic field inclination is θi = 40°. The part of the higher ionosphere and the magnetosphere above 300 km has been eliminated by means of the short circuit along the geomagnetic field lines.
(a) (b) (c)  From Figures 11 and 12 it is seen that in the case of the closed geomagnetic field lines the correspondence between the values of the penetrated electric field obtained in the dynamic and quasi-static simulations is also valid at higher inclination angles of the geomagnetic field. The maximum values of |E y | component along the axis OY perpendicular to the geomagnetic field are practically the same in Figure 9e, and in Figure 12b, for the dynamic simulations; and in Figure 10b, and Figure 11b, for the quasi-static ones From Figure 6, Figure 9, and Figure 12 it is seen that in the horizontal distributions |E x,y,z |(x,y) there exists the prolonged fine structure of the scales~4000 km. This structure is more expressed at lower frequencies, compared to Figures 6 and 9, parts a The origin of this fine structure is due to propagation of the EM waves in the gyrotropic ionosphere and in the lower magnetosphere, z ≤ 600 km, which are strongly non-uniform in the vertical direction. The strong coupling between outgoing waves and backward ones occurs there. Such interference is more essential at lower frequencies where the longitudinal wave numbers are smaller. Our simulations have demonstrated that at higher frequencies ω > 0.05 s −1 this prolonged structure is practically absent. The prolonged structures are less pronounced at inclination angle θ i = 40 o (Figure 12a-c) than at θ i = 10 • (Figure 9d-f). The details of forming this structure should be a subject of special investigations, especially at lower frequencies ω < 0.025 s −1 .
At lower frequencies ω < 0.025 s −1 the utilized algorithm for the dynamic simulations needs the extremely large amount of the computer memory, namely the decrease of the frequency two times demands the increase of the memory at least four times, i.e., the doubling of the points along both OX and OY axes. The reason is that the field components E x , E y are used as independent functions under the dynamic simulations and the transition to the exactly static case ω → 0 is a non-trivial computation problem. In the quasi-static simulations, the single function ϕ is used at all frequencies and the problem pointed above is absent.
The dependencies of the maximum values of the electric field components on the position of the external current have been simulated. A dynamic approach was used for the cases of open and closed field lines (see the first and second lines in Table 1). It is seen that the maximum values only slightly depend on the position of the maximum external current density when this current lies in the atmosphere, z 1 ≥ z 0 . Some decrease occurs when the current partially cuts of by the Earth's surface. Additionally, in Table 1 (see the third line) there are dependencies of the maximum values of the electric field components for the closed geomagnetic lines obtained under the quasi-electrostatic simulations. Table 1. The dependencies of the maximum values of |E x |, |E y |, |E z |, mV/m, on the position of the excitation current z 1 , the circular frequency is ω = 0.025 s −1 . The inclination angle of the geomagnetic field is θ i = 10 • . There are the same parameters z 0 = 7 km, j 0 = 0.1 µA/m 2 ; x 0 = y 0 = 500 km, besides the data in the last column of the Note that the dependence on the vertical coordinate in the spatial distribution of the external current (9) with the parameters z l = (3-5) km, z 0 = 7 km corresponds qualitatively to the distribution of seismogenic external current in [67]; transverse scale x 0 = y 0 = 500 km and 1000 km correspond to the parameters of powerful earthquakes with magnitudes M = 6.27 and 6.98, respectively [68]. As it is seen from Table 1, electric field penetrating at a given altitude in the ionosphere with the chosen parameters differ less than 1.5 times, for the cases of open and closed geomagnetic field lines, in the simulations based on the dynamic approximation.
The distributions of the electric field components are presented for the case when the external current partially is cut by the Earth's surface, Figure 13. It is seen that the distributions shown in Figure 13a,b are similar to ones in Figure  6d,e, respectively, while curves 1, 2, 3 in Figure 13c are similar to curves 3 in Figure 5a-c, respectively. The difference is in 15% smaller maximum values of |Ex|, |Ey| in Figure 13, due to the cut-off of the external current by the Earth's surface.
In addition to the calculations, presented in Table 1, it was also shown that for the quasi-electrostatic case and closed geomagnetic field lines, in the frequency range 10 −4 s −1 ≤ ω ≤ 0.025 s −1 , the dependence of the maximum values of the electric field components on the frequency, at z = 200 km changes unessentially within 3% difference; this was shown, in particular, for the parameters z1 = 14 km, z0 = 7 km. Note also that, in the frequency range

Discussion
The total vertical distance of the simulation is L z = 1500 km. The geomagnetic field inclination is θ i = 10 • . The parameters are z 1 = 5 km. z 0 = 7 km. In the part (c) the dependencies of the maximum values of the electric field components on z are presented, curve 1 is |E x |, 2 is |E y |, 3 is |E z |.
It is seen that the distributions shown in Figure 13a,b are similar to ones in Figure 6d,e, respectively, while curves 1, 2, 3 in Figure 13c are similar to curves 3 in Figure 5a-c, respectively. The difference is in 15% smaller maximum values of |E x |, |E y | in Figure 13, due to the cut-off of the external current by the Earth's surface.
In addition to the calculations, presented in Table 1, it was also shown that for the quasi-electrostatic case and closed geomagnetic field lines, in the frequency range 10 −4 s −1 ≤ ω ≤ 0.025 s −1 , the dependence of the maximum values of the electric field components on the frequency, at z = 200 km changes unessentially within 3% difference; this was shown, in particular, for the parameters z 1 = 14 km, z 0 = 7 km. Note also that, in the frequency range 10 −4 s −1 < ω < 10 −3 s −1 , the calculations are rather formal and mathematical in nature. The reason is that during a time of order of the corresponding period of oscillations, 2π/ω ∼ 20 h (ω ∼ 10 −4 s −1 ), the changes in the parameters of the atmosphere and ionosphere, ignored in the frequency domain presentation, are already quite essential.

Discussion
The important problem of the penetration of the electric field to the ionosphere from the near-Earth sources is under investigation till now. A standard approach is the utilization of the electrostatic or quasi-electrostatic equation for the electric potential in media with highly anisotropic conductivity supplemented by boundary conditions. At the high conductive solid Earth's surface, the boundary condition is a constant value of the electric potential. Usually, some boundary condition for the electric potential and its derivatives is formulated at the boundary between the upper atmosphere and the lower ionosphere. But still, it has been unclear what boundary condition should be applied within the quasi-electrostatic approach in the magnetosphere and how to transfer the boundary condition to the lower ionosphere and whether such a transfer is valid. This basic problem is due to the high anisotropy of the conductivity of the upper ionosphere and the magnetosphere and mathematically is expressed through the product of a small value of the electric field along the geomagnetic field line and a very big value of the conductivity along with it, in order to calculate the density of the electric current. Thus, within the framework of the electrostatics, the correct formulation of the upper boundary condition for the electric potential is doubtful.
Nevertheless, the processes in the ionosphere and the magnetosphere are principally dynamic, so the case of the quasi-electrostatics may be justified only at small spatial scales. To perform this approximation, the condition (25) must be satisfied. The maximal period T can be estimated as T ≤ 2·10 3 s, i.e., the variation time of the current sources is smaller than 1 hour. Therefore, for the validity of the quasi-electrostatics, the distance along the geomagnetic field line should be L z < 2·10 4 km. From this simple estimation, it is seen that in the case of the open geomagnetic lines the quasi-electrostatic approximation is doubtful, whereas for the closed geomagnetic field lines it could be valid.
Generally, the dynamic approach should be used to consider the problem of penetration of the ULF electric field to the ionosphere. The simulations within the quasielectrostatic approximation should correspond to the results of the dynamic simulations. The dynamic method is based on the Maxwell equations and proper boundary conditions for the electromagnetic field. Because the physical processes possess various temporal scales, to make the qualitative analysis as clear as possible, it is rather better to apply the spectral approach, i.e., within the frequency domain. Our dynamic simulations have been proven that in the case of open geomagnetic field lines the quasi-electrostatic approximation does not provide adequate results. Moreover, the results of simulations obtained within the dynamic method do not depend on the position of the upper boundary within the magnetosphere when the heights are ≥1000 km, which proves the validity of the dynamic approach. The boundary conditions for the EM field are the absence of the ingoing magnetohydrodynamic waves. Contrary to this, in the case of open geomagnetic lines, the results of the quasi-electrostatic simulations depend on the position of the upper boundary condition and the penetrated electric field in the ionosphere F-layer is several orders smaller than one obtained in the dynamic simulations.
In the case of the open geomagnetic field lines, the dynamic method for the penetration of the ULF electric fields from near-Earth sources to the upper ionosphere and magnetosphere is similar to the consideration of penetration of the magnetosphere perturbations due to the solar wind to the lower ionosphere and the atmosphere. Namely, those perturbations propagate along the magnetic field line as the magnetohydrodynamic waves, and only near the Earth's surface the quasi-static approximation is valid and used [69].
The physical reason for the significant difference between the results obtained on the basis of the dynamic and quasi-static approaches (Figures 5 and 7) in the case of open magnetic field lines is as follows. In the dynamic approach, we use the radiation condition (Equations (11) or (12)) as the "upper boundary condition". This condition is set at the heights of the upper ionosphere or magnetosphere, where, for the frequency range under consideration, inhomogeneities in the medium are so smooth on the wavelength scale that the reflection for a wave emitted from the studied area can be neglected with sufficient accuracy [70]. This condition corresponds to the principle of causality [71], in particular, in the presence of (external current) sources only in the lower atmosphere. It turns out that the field excited by a given source in the lower atmosphere and penetrating into the ionosphere and magnetosphere is negligibly dependent on the location of the boundary at which the "upper boundary condition" (of radiation) is set in the upper ionosphere or magnetosphere (see Figure 5 and its description in Section 3). As for the quasi-static approach, it is applicable, if (25) is satisfied. Obviously, in the case of open geomagnetic field lines, a similar condition is violated at a sufficiently large value of L z , for example, at L z = 25,000 km. This is confirmed by the corresponding modeling illustrated in Figure 7.
In addition to this, we can also offer the following visual explanation of the significant difference between the results obtained on the basis of the dynamic and quasi-static approaches. If the radiation condition under the dynamic approach is formulated at a certain value of the coordinate z = L z1 , then we can assume that at z > L z1 there are, respectively, only waves radiated from the considered region 0 ≤ z ≤ L z . If we now set the radiative boundary condition at z = L z2 > L z1 , then in both regions L z1 < z ≤ L z2 , z > L z2 , in fact, only outgoing waves will be present, since the reflection in the region L z1 < z ≤ L z2 is negligible. Therefore, the field in the region z < L z1 , including the lower and middle ionosphere, will be practically insensitive to the location of the "upper boundary" selected in the upper ionosphere or magnetosphere (at z = L z1 or z = L z2 ). The situation is significantly different when the quasi-static approximation is used. In this approximation, there is neither delay nor radiation. As a result, when using the quasi-static approximation, the entire region 0 ≤ z ≤ L z automatically belongs to the "region of influence" in the solution, even at very large values of L z . Therefore, the value L z , which determines the place where the "upper boundary condition" is specified, significantly affects the solution when applying the quasi-electrostatic approximation. This indicates the inapplicability of the quasi-electrostatic approximation in the case of open geomagnetic field lines.
In the case of the closed geomagnetic field lines both the dynamic approach and quasi-electrostatic one yield practically similar results for the penetration of the electric field into the ionosphere F-layer. Therefore, it is natural to set the appropriate boundary conditions on the surface of the Earth, in the regions from which the geomagnetic field lines emerge and where they return. Boundary conditions are set in a similar way when using the quasi-static approximation. Both in the quasi-electrostatics and in the dynamic approach the boundary conditions are zero values of the tangential electric field on the Earth surface; in the quasi-electrostatics, they are equivalent to the zero-electric potential.
It was mentioned in [46] the difference between the pointed above configurations of the geomagnetic field lines for the electrostatic simulations and the necessity to use different boundary conditions in the lower ionosphere. The upper boundary conditions used in [46,47] correspond qualitatively and by mathematical structure to the upper boundary condition derived in the present paper, see Appendix C, Equation (A53).
In Section 3 it has been mentioned the mathematical limitation for the lowest frequency that can be used in our model of the dynamic simulations. This mathematical limitation seems not occasional and has a physical sense with respect to an attempt to "move" into the region of lower frequencies on the basis of the dynamic approach. Namely, with decreasing frequency, one should expect that the solution for the electromagnetic field can be represented in terms of electromagnetic potentials [66,71,72]. Alternatively, such a solution can be presented in terms of magnetic and electrical Hertz potentials, which in the limit of low frequencies are transformed into electrostatic and magnetostatic potentials [71]. With this approach, it is quite obvious that the solution of the dynamic problem, in principle, cannot reach the purely electrostatic limit, since in this case "the magnetostatics is lost". In any case, this must necessarily happen in the presence of low-frequency currents in the system, as in our case. In a geophysical experiment, this corresponds to the presence of a magnetic ULF field, excited by the corresponding currents, which can be measured at ground-based and satellite observatories.
Despite the fact that in this work a different approach is used with the reduction of Maxwell equations in the frequency domain to a system of two equations of the second order (8) for the Fourier amplitudes of two transverse components of the electric field, it seems obvious that also, in this case, two corresponding functions, associated with the transverse components of the electric field, cannot in the limit tend to a single function, namely, the electrostatic potential. The aforementioned mathematical difficulty with an increase in the required amount of RAM with a decrease in the source frequency is probably a mathematical reflection of the physical problem of the transition to quasi-electrostatics and quasi-magnetostatics. In the future, it will be of interest to consider the question of the "magneto-electrostatic" limit for the electrodynamic solution in the LEAIM system.
Our dynamic simulations, supplemented by quasi-electrostatic simulations (for the case of closed geomagnetic field lines) have demonstrated that the values of the electric field are, by the order of value, of about several, namely (1-5) mV/m in the lower ULF range of frequencies (10 −3 s −1 ≤ ω ≤ 5 · 10 −2 s −1 ) in the ionosphere F-layer at the heights z = 200-600 km and in the lower magnetosphere. Moreover, the corresponding simulations may be done even at such low frequencies, as 10 −4 s −1 ≤ ω < 10 −3 s −1 , but only formally, if one would ignore the fact, that the parameters of the atmosphere and ionosphere change essentially during the period (2π/ω ∼ 20hours, ω ∼ 10 −4 s −1 ) of corresponding current source oscillations.
The results for the penetration of electric field into the ionosphere are obtained using the models of external current sources localized in the lower atmosphere with a maximum value of 0.1 µA/m 2 , generated by tropical cyclones [54,56,57] or seismogenic (earthquake preparation) processes [1,6,24]. Such electric fields, by the order of value, are corresponding to the ionospheric observations above tropical hurricanes [61,73] and the region of the preparation of strong earthquakes [40]. Such values of the electric field correspond also to the observations of seismogenic variations in electron concentration and TEC [6,41]. Therefore, the simulated values of the electric field are of a qualitative agreement with the results of observations. Note also that the value of the extraneous current maximum of 0.1 µA/m 2 accepted in the calculations in this work is far from the maximum possible corresponding value known from the literature. Namely, there are theoretical estimates that substantiate the possibility of generating external currents in the lower atmosphere with values exceeding those mentioned above by one or even two orders of magnitude. Such estimates are given, in particular, in works devoted to modeling and observations of ionospheric effects caused by electric fields from surface current sources of a seismogenic nature (during the preparation of powerful earthquakes) [6,67,74,75] and sources associated with tropical cyclones [54,56].
The external current only mimics the real current generated in the lower atmosphere by all physical processes taking place in a self-consistent manner, but at the same time, it allows one to determine the field penetrating into the upper ionosphere/magnetosphere. To adequately determine the electric field in the lower atmosphere, a synergetic approach is required with consideration of all relevant physical processes, instabilities in the LEAIM system, their nonlinear limitation [9,15,34], as well as the corresponding boundary conditions, including ones in the lower atmosphere, for various current components associated with different types of charged particles, including aerosols and ions, and taking into account photochemistry, convection, hydration and other dynamic processes [1,4,6,24,32]. The synergetic approach would include, naturally, the self-consistent investigations with accounting for positive feedback between the penetrated EM field and changes in the ionosphere parameters. This may be a subject of future works.

Conclusions
The penetration of an electric field through the LEAIM system is of a dynamic nature, in particular, in the presence of external current sources exciting the electromagnetic fields in the lower atmosphere or lithosphere. This dynamic nature is due to the variation in time in both properties of the atmosphere, ionosphere, and the magnetosphere and of the external sources. Also, it should be taken into account the possible propagation of electromagnetic disturbances over large distances along the geomagnetic field lines. This propagation results in retardation of the electromagnetic field. Thus, the correct formulation of the problem should be done, generally speaking, on the basis of a dynamic approach to solving the Maxwell equations.
In this case, the corresponding dynamic boundary conditions for the electromagnetic field must also be formulated correctly. The ULF waves considered in this paper propagate to a large distance along the lines of the geomagnetic field. Two essentially different situations arise: (1) the geomagnetic field lines emerging from the Earth-Atmosphere boundary in the region of the location of an external electric current source are open, passing from the upper ionosphere into the deep magnetosphere, and do not return into the ionosphere (for sufficiently high latitudes); (2) after penetration from the ionosphere into the lower magnetosphere, the geomagnetic field lines return to the magnetically conjugated ionosphere and then again fall from the atmosphere into the Earth (lithosphere). Naturally, for the electromagnetic field at the "second end" of the line of force, where the numerical calculation completes, there are significantly different boundary conditions. In the case of open field lines, this is the condition of radiation, and in the case of closed field lines, the tangential component of the electric field is equal to zero, if the conductivity of the Earth is considered equal to infinity. It turns out that in the case of open geomagnetic field lines only the dynamic approach is applicable, but the quasi-static approach is not valid, and in the case of closed geomagnetic field lines the quasi-static approximation can also be used, since it gives results that qualitatively and quantitatively, by the order of value, coincide with the results obtained on the base of the dynamic approach. In the case of closed filed lines, the electric potential calculated in the upper part of the first ionosphere (located above the region where external sources are considered) can be transferred along the geomagnetic field lines to the upper part of the magnetically conjugated ionosphere, and then it will continue the calculation, in particular, in the atmosphere magnetically conjugated region.
The values of the electric field obtained in the dynamic approximation are of about several (1-5) mV/m in the ionosphere F-layer and in the lower magnetosphere. In this case, the results are similar to ones obtained within the quasi-electrostatic approximation simulations when the geomagnetic field lines return to the Earth's surface. In the region near the boundary between the atmosphere and the ionosphere z = 60 km the electric fields are several (1-5) V/m. The external electric current sources have a maximum value of 0.1 µA/m 2 and are localized in the lower atmosphere. The simulated values of the electric field are of a qualitative and even quantitative, by the order of value, agreements with the results of observations of electric fields and perturbations of the electron concentration/TEC, caused by the current developed in connection with tropical cyclones or the seismogenic processes (the preparation of the most powerful earthquakes).
Note that to the parameters used for the modeling presented in Figures 5 and 6, the estimation (24) of the penetration depth appeared to be suitable and corresponds qualitatively to the results presented in the Figures mentioned above.
Funding: Financial support for this research was partially provided by CONACYT-Fondos Sectoriales grant CONACYT project A1-S-29604 to Vsevolod Yutsis. Yuriy Rapoport, Asen Grytsai, Vasyl Ivchenko, Sergei Petrishchevskii and Valery Krivodubskij are thankful to the Ministry of Education and Science of Ukraine for the partial support of the present work, grant 20BF051-02 [0120U102178] "Wave processes and effects in active resonant layered plasma media and metamaterials". Alla Fedorenko is thankful to the National Research Foundation of Ukraine for the support of the work by grant 2020.02/0015 "Theoretical and experimental studies global perturbations of natural and artificial origin in the Earth-atmosphere system".
Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: The corresponding ionospheric data have been taken from the wellknown published handbooks, referred in the paper [18,51,52,[63][64][65]. There is implementation of the standard and approbated computed code without using the commercial software. The computer codes are the finite difference schemes, the Fast Fourier transform, the matrix factorization method, the correct approximation of the boundary conditions for different cases. The description of the computer codes is presented in the text of the paper in Appendices B and C, as well as in proper references [48][49][50].

Conflicts of Interest:
The authors declare no conflict of interest. The relations (A12)-(A17) completely determine the coefficients of the matrix equation (8). Assume that (a) the medium above the upper boundary z = L z is a homogeneous anisotropic but not gyrotropic layer F [52] with parameters that do not depend on z; this is justified by the fact, that the parameters in the upper ionosphere and magnetosphere change much slower than in E region of the ionosphere [15]; (b) there are no sources of ULF waves in the region z ≥ L z . Therefore the plane (outgoing) waves propagate in the region z ≥ L z in the positive direction of axis Z, and it is possible to find the connection between ∂E x,y /∂z and E x,y at z = L z + δ, where δ is infinitesimally small positive value. Therefore accounting for that E x,y are continuous at z = L z , we need to find the relationship between the derivatives ∂E x,y /∂z at z = L z − δ and z = L z + δ.
Tensorε in the F layer of the ionosphere in the coordinate frame x y z has the following form [51,52]:ε =   ε 1 0 0 0 ε 1 0 0 0 ε 3   . (A18) Let us use for the ULF waves in F region the approximation and account for that the following relationships for the vector components in the coordinate frames denoted with and without the primes (Figure 1): k x = k x cos θ + k z sin θ k z = −k x sin θ + k z cos θ E z = E z cos θ − E x sin θ ≈ 0 E x = E x cos θ + E z sin θ ≈ E x / cos θ.
In the matrix factorization method the following expressions are used: At the upper boundary in the magnetosphere the boundary conditions are, Equations (11) and (12)  At the first Earth's surface the boundary condition is And, once all the factorization coefficientsβ j , → r j ; j = N, . . . , 1 are already known, one can compute the transverse components of the electric field → F j , j = 1, . . . , N sequentially from (A41), as well as all another electromagnetic components.
The results of simulations do not depend on the step along OZ axis when this step is quite small, as usually h z < 1 km.