Excitation of planetary electromagnetic waves in the inhomogeneous ionosphere

In this paper we develop a new method for the analysis of excitation and propagation of planetary electromagnetic waves (PEMW) in the ionosphere of the Earth. The nonlinear system of equations for PEMW, valid for any height, from D to F regions, including intermediate altitudes between D and E and between E and F regions, is derived. In particular, we have found the system of nonlinear one-fluid MHD equations in theβ-plane approximation valid for the ionospheric F region ( Aburjania et al. , 2003a, 2005). The series expansion in a “small” (relative to the local geomagnetic field) non-stationary magnetic field has been applied only at the last step of the derivation of the equations. The small mechanical vertical displacement of the media is taken into account. We have shown that obtained equations can be reduced to the well-known system with Larichev–Reznik vortex solution in the equatorial region (see e.g. Aburjania et al. , 2002). The excitation of planetary electromagnetic waves by different initial perturbations has been investigated numerically. Some means for the PEMW detection and data processing are discussed.


Introduction and formulation of the problem
In recent years a considerable effort has been made for treating waves propagation problems in the inhomogeneous Earth's ionosphere (see e.g.Clark et al., 1971;Rapoport et al., 2004Rapoport et al., , 2009;;Sorokin and Fedorovich, 1982;Alperovich and Fedorov, 2007).In particular, the theory of "neutralloaded" MHD waves in the ionospheric plasma (including motion of neutrals (see e.g.Sorokin and Fedorovich, 1982)), as well as the phenomena of "mass loading" in space plasmas (Szegö et al., 2000) has been developed.Vortex structures in the Earth atmosphere, ionosphere and magnetosphere have been studied by Onishchenko et al. (2008Onishchenko et al. ( , 2013)); Onishchenko and Pokhotelov (2012); Saliuk and Agapitov (2013); Saliuk et al. (2012); Verkhoglyadova et al. (2001).Earlier studies have shown that MHD waves, propagating along a magnetic field line through the magnetosphere between the magneto-conjugated ionospheres could transform into horizontally propagating hydromagnetic waves (see e.g.Tolstoy, 1967).Herron (1966) has found that these waves can be captured into the ionospheric waveguide at the Alfvén speed minimum in the F region.Further development of these ideas brought the concept of vortex planetary electromagnetic waves (Kaladze et al., 2003;Khantadze, 1973;Aburjania et al., 2002Aburjania et al., , 2003aAburjania et al., , 2005;;Aburjania and Chargazia, 2011;Rapoport et al., 2011Rapoport et al., , 2012c)).The developed theory of PEMW is supported by a number of observational results of planetary-scale electromagnetic perturbations (see e.g.Aburjania and Chargazia, 2011;Burmaka et al., 2006).However, further theoretical and experimental investigations were necessary to reveal both "slow" and "fast" PEMWs in the ionosphere region (Khantadze, 1973;Aburjania et al., 2003a;Kaladze et al., 2003;Aburjania and Chargazia, 2011).Southwood and Kivelson (1993) have found that travelling ionospheric vortices (TIVs) can be divH = 0; divU = 0; (3) Here c is the speed of light; U , ρ 0 are the velocity and density of media at height z of the chosen β-plane, respectively; 0 is the rotation parameter; γ is the effective mechanical damping parameter; E, H are the electric and magnetic fields; Ẽ = E + c −1 [U × H ] is the electric field in the moving media; j is the electric current density; F A = c −1 j × H is the Ampere force.
Here, we apply approximations of an incompressible fluid and absence of stratification due to gravity, respectively to perturbations.Nevertheless, hydrostatic equilibrium (in other words stratification) is accounted for with steady-state parameters such as density of the neutral and ion components.We use, in general, the approximation used before in the set of papers devoted to PEMW (Aburjania et al., 2003a, b;Kaladze et al., 2003).Practically, stratification is accounted for, but only respectively to stationary distribution of the density, as a result of hydrostatic steady state (Gill, 1982).We account for this, in the frames of the method of frozen coefficients, by comparing the structures at 300 km and 600 km with different parameters (see Sect. 6).At the same time, because the mechanical motion of the perturbations is slow comparatively to the sound velocity, approximation of incompressibility is used respectively to perturbations.As a result, gravity force is included only into the hydrostatic steady state and not into the dynamics of perturbations.This is why the term proportional to (free-fall acceleration) is absent in the RHS of Eq. ( 15) for velocity perturbation U z .We included into the consideration the vertical velocity perturbation U z , but vertical displacement depends only on X, Y and is still small enough, providing that change of the ionospheric parameters due to vertical displacement is negligibly small, not to violate the condition of incompressibility for the perturbations.Respectively, the relation of "zero divergence" is valid and the corresponding presentation of the velocity X, Y components through ONE function ξ (see Eq. 12) are applicable for the horizontal components of velocity, even if a small vertical component corresponding to relatively small vertical displacement is also included in our model.Therefore, our model is internally consistent.The errors made by assuming that the magnetic and rotation axes coincide is of the same order as the errors inherent in the β-plane approximation (see e.g.Aburjania et al., 2005;Kaladze et al., 2003).For the above scales of disturbances, the consequences of the misalignment of both axes, which is about 0.2 rad, are beyond the accuracy of the β-plane approximation and, therefore, cannot change the results qualitatively.Note that the variation of ion mass, reflecting the changing in composition, is included into the model through the elements of conductivity tensor in the Eq. ( 4).
Before introducing the terms of conductivity tensor in the Eq. ( 4), let us specify the magnetic field.We define the velocity U and total magnetic field H by the relations: Here U 0 and H 0 are the stationary velocity and geomagnetic field; U and H are the small perturbations of the velocity and magnetic field; h, h 0 are the unit vectors that determine the directions of the total and stationary geomagnetic field; H 0 and H = (H 0 + H ) 2 are the values of stationary geomagnetic and total magnetic field; h determines the direction of the non-stationary part of the magnetic field.We also assume that √ h 2 1.By taking into account the geometry of the model (see Fig. 1), dipole approximation of the geomagnetic field, and Eq.(3), the stationary (wind) velocity and magnetic field can be written as follows: Here x, y and z are directed along the parallel, meridian and vertical, respectively; λ = 0 and λ = π/2 are the latitudes which correspond to the equator and pole, respectively (see Fig. 1).In approximation z/R 1 the components of the stationary geomagnetic field H 0y and H 0z , and the rotation parameter 0 can be represented as (see e.g.Kaladze et al., 2003) Note that the chosen coordinate dependence of the stationary velocity in Eqs. ( 6) also satisfy the second equation in Eqs.(3).Therefore, as it follows from the relations for the stationary geomagnetic field Eqs. ( 7), with accuracy to z/R 1, we obtain (Kaladze et al., 2003) The elements of conductivity tensor are defined by the wellknown formulas (see e.g.Guglielmi and Pokhotelov, 1996;Sorokin and Fedorovich, 1982;Kaladze et al., 2003): Here σ || , σ P , σ H are the parallel, Pedersen, and Hall conductivities, respectively; ω H e,i are the electron and ion gyrofrequencies (defined by the value of the total magnetic field H ); ν e is the sum of collision frequencies of electrons with neutrals and ions; ν i is the collision frequency of ions with neutrals.Considering Eq. ( 4) as algebraic equation for electric field E and by taking into account the second equation in Eqs.(2), the electric field in the moving media Ẽ can be expressed in terms of the magnetic field H as and, therefore, we can exclude the electric field from the system of MHD equations for PEMW.Here The coefficients q 1 and q 2 correspond to the "transverse" electric field (defined by the Hall current) and to the "longitudinal" electric field, respectively.The dependence of the elements of conductivity tensor σ eff , coefficients q 1,2 on height are shown in the Fig. 2. The limiting case of isotropic D region of the ionosphere corresponds to the low altitudes, where q 1 = 0, q 2 = 0.In the E region (where q 1 q 2 ), we arrive at the gyrotropic medium (Guglielmi and Pokhotelov, 1996)  The values of coefficients q 1,2 vs. height, z, are plotted at (b). q 1 and q 2 are coefficients in the expressions, corresponding to the "transverse" electric field, determined by Hall current and "longitudinal" electric field (see Eq. 10).
anisotropic and non-gyrotropic (Fig. 2a).The present model is also applicable for intermediate heights (i.e. for the heights between D and E, and between E and F regions), where the relationships mentioned above between q 1 and q 2 are not satisfied (see also Fig. 2b).The elements of conductivity tensor (see Eq. 9) depend on total magnetic field H (see Eq. 8) which is the function of the non-stationary magnetic field H (see Eq. 5).Therefore, the electric current Eq. ( 4) and electric field Eq. ( 10) become nonlinear.Due to such nonlinearity, hereafter we will call Eq. ( 4) as "nonlinear Ohm's law".In the linear limit (i.e. in Eq. ( 4)) h is replaced by h 0 and the elements of the conductivity tensor ( σ ) are replaced by corresponding linear values, the Eq.(4) become linear.In this case, hereafter, we will call Eq. ( 4) as "linear Ohm's law".Practically, it corresponds to replacement of H by H 0 in Eqs. ( 9) in all dependent on magnetic field terms.The β-plane approximation require implementation of "the method of frozen coefficients" (see e.g.Aburjania et al., 2003a;Kaladze et al., 2003) and satisfying of the following assumptions Second and third relations in Eq. ( 11) correspond to the case of neglecting gravity stratification.However, by taking into account equations for the stationary geomagnetic field, that is, Eqs. ( 7) and ( 8), we should require ∂H 0y,z /∂z = 0 and, therefore, ∂h 0 /∂z = 0, ∂h /∂z = 0 (see also Eqs. 5 and 7).By using the conditions of incompressibility, that is, Eq. ( 3) and first two relations in Eq. ( 11), the horizontal components of the non-stationary velocity and magnetic field can be presented as where ⊥ = ∂ 2 ∂x 2 + ∂ 2 ∂y 2 .Note that in contrast to the previous PEMW models (see Aburjania et al., 2002Aburjania et al., , 2003a;;Aburjania and Chargazia, 2011;Aburjania et al., 2005;Kaladze et al., 2003) in this case the perturbed velocity U z is not equal to zero.By taking the z component of the first equation in Eqs.(2) and then taking the z component after applying operation curl to the same equation, and taking into account the second relation in Eqs. ( 12) yields equations of the form: The electric field in the moving media Ẽ in Eqs. ( 13) and ( 14) is defined by the Eq. ( 10).Similarly to Eqs. ( 13) and ( 14), one can obtain the corresponding equations from Eq. (1), by taking into account Eq. ( 11), the first relation in Eq. ( 12) and the approximation of "frozen coefficients": Note that with such approximations (see also the last equation in Eq. 11) the pressure vanishes from Eqs. ( 15) and ( 16).Also, we assume that non-stationary parts are revealed from the LHS of Eqs. ( 13) and (14).The system of Eqs. ( 13)-( 16), along with the first and second equations in (12), complete the close nonlinear system for PEMW.This system is valid for any heights from the D to F regions of the ionosphere.Due to ionospheric motions, additional terms associated with an electric field and Ampere force, which is physically associated with the Lorentz force, come up in the nonlinear equations for the magnetic field and velocity (see Eqs. 13 and 15).
Another nonlinear term is associated with convective nonlinearity.In the β-plane and frozen coefficients approximations, the vertical curl operations are applied to Eqs. ( 13) and ( 14) yielding Eqs. ( 15) and ( 16), respectively.Therefore, the inhomogeneity contributes, along the meridional and vertical directions, in the linear and nonlinear parts of Eqs. ( 14) and ( 15) simultaneously.Herewith, in the approximation of an incompressible fluid, the entire nonlinearity has the vector form and described by the Jacobian of the corresponding variables (see also the last part of Sect.3).It can be explained as follows.An incompressible medium is free of shear perturbation, and therefore, some possible distortion is described by the effective "torsion".Qualitatively such a "mechanical" interpretation remains valid, despite the presence of electromagnetic forces, since the Ampere force in Eqs. ( 15) and ( 16) is similar to the Coriolis force, while the magnetic field is similar to the vector rotation.
In the next sections we will analyse the consequences following from the system of Eqs. ( 13)-( 16).Firstly, we will use an approximation of the "linear Ohm's law".The nonlinearities in the second terms in the RHSs of Eqs. ( 13) and ( 14) and in the first and second terms in the RHSs of Eqs. ( 15) and ( 16) are taken into account.Next the procedure of derivation of the total nonlinearity, including one in the "Ohm's law" and terms with the electric field in Eqs. ( 13) and ( 14) will be outlined in general.Finally, the developed numerical method and computational results for the near-polar region in the "linear Ohm's law" approximation will be presented.

Nonlinear equations for PEMW at arbitrary altitudes from D to F regions in the "linear Ohm's law" approximation
Let us use the approximation of "linear Ohm's law".By introducing the scales L, t 0 , U 0 for the length, time and velocity, respectively, and the scales H sc , ζ sc , A zsc , zsc , ξ sc , 0sc for the normalisation of H x,y,z , ζ , A z , U z , z , ξ , 0y,z , respectively, the normalisation could be presented in the following form: Here Next, let us introduce dimensionless parameter α 0 = U 0 t 0 /L and the Alfvén speed V A0 = H sc / √ 4πρ 0 .Note, that α 0 = 1, when L = U 0 t 0 .By omitting "dash" above normalised variables we obtain the system of normalised equations in β-plane coordinates.This system includes four evolutional equations for ζ , H z , U z , z , and two Poisson equations for A z and ξ and has the following form: Here Equation ( 17) along with Eqs. ( 18)-( 24) form a complete set of nonlinear equations describing PEMW at any height between D and F .In particular, we can apply this theory to the heights below 150 km, where, as it follows from Fig. 2b, the difference of q 1 from 0 and q 2 from −1 is significant.Coefficients C ij are presented in Appendix A. F 1Qζ,2Qv describe the electromagnetic losses.The terms F NLζ,NLH,NLU,N L contain vector nonlinearities in the form of Jacobians.Note that in Eqs. ( 18)-( 23), some coefficients (see Eqs. A1-A5) contain derivatives of vector components h 0 , H 0 and 0 with respect to y and z.These coefficients describe the contribution of the curvature of the Earth in the electromagnetic and mechanical forces.
4 Nonlinear equations for PEMW in the "linear Ohm's law" approximation.Application to the "polar" and "equatorial" regions The application to the near-pole region (i.e. for the limit λ → π/2) is shown in Appendix B (see Eqs. B1-B5).The numerical results for this region are presented in Sect.6.In the limiting case λ → 0 (i.e. for the "equatorial" region) from the system of Eqs. ( 17), accounting for relations ( 18)-( 24) we obtain the system of normalised Eqs.(C1)-(C5) which is shown in Appendix C. Note, that in the approximation e. losses are absent), the Eqs.( C1)-(C5) reduces to the system for z , H z and ξ : This system is identical to the system obtained by Aburjania et al. (2002).Therefore, Eqs. ( 25)-( 27) yield vortices, described by the Larichev-Reznik solution, including the Bessel function in the core and the McDonald function at the periphery of the nonlinear structure (Aburjania et al., 2002).The numerical simulation of the excitation of PEMW at the middle longitudes is out of scope of this paper.Nevertheless, due to the presence of large stationary horizontal magnetic field, the corresponding estimates of the spatio/temporal scales can be found by applying Eqs.(C1)-(C5) for the near-equatorial region for both.For example, at Z = 300 km, the characteristic temporal scale is 10 3 s, the spatial scale is a few thousand km, and the velocity scale is a few km s −1 .These values are in qualitative agreement with observational results of nearly longitudinal propagation of slow MHD/PEMW disturbances (Burmaka et al., 2006;Aburjania and Chargazia, 2011).

An inclusion of nonlinearity of the ionospheric current: the case of F region
Let us outline the procedure of inclusion of the current nonlinearity into Eqs.( 3)-( 5).Estimations show that the condition of the weak nonlinearity (i.e.H /H 0 1) is satisfied with a good accuracy.Therefore, we can expand the secondorder nonlinearity in terms of the small parameter H /H 0 .
The unit vector of the whole magnetic field Eq. ( 5) and Pedersen conductivity Eq. ( 9) are modified as follows: where σ −1 P0 is the corresponding value of Pedersen conductivity for the linear case, µ P is a multiplier of order of unity.In the framework of β-plane approximation and "method of frozen coefficients" we are forced to put hereafter ∇σ P0 = 0 and ∇µ P = 0. To simplify, let us consider the F region, where σ eff → σ P , q 1 → 0, q 2 → −1 (see also Fig. 2).The nonlinear additional term to the electric field, δ ẼNL , can be revealed from the RHS of the Eq. ( 10): where By substituting the relations ( 28) and (29) into Eqs.( 13) and ( 14), we obtain the additional nonlinear terms in the RHS of the first two equations in the system (17), for example, Note the corresponding nonlinearities are derived, but due to bulkiness of expressions they are not presented here.

Numerical model and the results of modelling
To analyse generation and propagation of PEMW in the nearpole region numerically, we used the normalised set of equations:   heights within F region of the ionosphere.At higher height, i.e. at Z = 600 km the characteristic time for retain PEMW is smaller, i.e. < 10 2 s.This is due to the smaller value of the Pedersen conductivity σ P ∼ 10 3 s −1 and higher Alfvén speed.Excitations of PEMW at daytime at Z = 600 km are negligibly small and not shown here.
The results are presented for the initial normalised value of the magnetic field Hz = 0.01, what corresponds to a dimensional value H z ∼ 100 nT, that is possible for PEMW vortices (Aburjania et al., 2005).The chosen value of normalised vortex ζ = (curlH) z ∼ 0.01 means that the corresponding vertical current is, by order of magnitude, J z ∼ 10 −7 A/m 2 .It should be noted that the considered vertical currents are smaller than currents accompanying some typical events of the magnetosphere-ionosphere coupling.An example of such a current is one accompanying magnetic ULF disturbances in the night time auroral region, caused by magnetic impulse events (Pilipenko et al., 1999).Further, as it is seen from the comparison of Figs. 5 and 6, the values of H z and ζ at Z = 300 km, after passing corresponding characteristic time, are, depending on a type of excitation, comparable to or few times larger than at Z = 600 km (compare Figs. 3, 4, 5a-c and Figs.6a-d.The same concerns also the values of horizontal velocities (compare Fig. 5f with Fig. 6d).Note that the Alfvén velocity, used for normalisation of the velocity components, is about one order of value larger at altitude Z = 600 km than at Z = 300 km (see Table 1).Respectively, the time scale t 0 is much less at altitude Z = 600 km (t 0 ∼ 150s) than at Z = 300 km (t 0 ∼ 1400s).The shapes of the excited PEMW also depend on the type of initial This system consists of four evolution equations for ζ , H z , U z , and z , and two Poisson equations for A z and ξ .Note that the corresponding system of equations, which includes the coefficient C ij in the limit λ → π/2, is presented in the Appendix B. All functions are set to zero at the boundaries of numerical domain.The full numerical box is L x = 20R km wide and L y = 20R km high (where R is the Earth radius).
We have chosen a sufficiently big numerical domain to avoid an nonphysical numerical influence (i.e.possible reflection) of the boundaries.To solve the Poisson equations, we applied fast Fourier transform (FFT) (Press et al., 1997).The evolution equations have been solved by splitting with respect to the physical factors (Marchuk, 1994).This method can be called a spectral -grid Arakawa type method.Specifically, the first fractional step in splitting is related to the diffusion-like process (terms ⊥ ζ and ⊥ H z in the first and second Eq.31).To find the function ζ we solved: ∂ζ /∂t = C 1a ⊥ ζ + F 1 .Before applying FFT, the function F 1 has been calculated with finite differences.We kept the same numerical grid (x, y) for the finite differences and for spectral FFT.Roache (1998) and Arakawa (1997) have shown that simplest approximations of the spatial derivatives in the Jacobians are not conservative and may lead to essential errors at long times.Therefore, the Arakawa approximation for the Jacobians has been used only at the first fractional step.The second fractional step is related to the convection terms in all evolution equations.For example, for the function ζ it is ∂ζ /∂t = − U 0x ∂ζ /∂x + U 0y ∂ζ /∂y .Since the wind velocities U 0x,y are at least one order smaller than the Alfvén velocity, the simplest corner-like upwind scheme can be applied there (see e.g.Roache, 1998).
This method demonstrated a good stability and flexibility.This is further development of our techniques proposed before for modelling of solitons and bullets in nonlinear gyrotropic structures (see e.g.Rapoport et al., 2002;Slavin et al., 2003;Zaspel et al., 2001), and solitons and strongly nonlinear waves in metamaterial structures (see e.g.Boardman et al., 2010Boardman et al., , 2011a, b;, b;Rapoport et al., 2012a, for details).The conductivities dependence on height and time of the day has been calculated based on data from Alperovich and Fedorov (2007) (see Fig. 2a).Based on that, we found coefficients q 1 , q 2 which are included in the RHS of the first four parts of Eqs.(17) (see also relations ( 18)-( 24) and Appendix A).The height dependence of the q 1 , q 2 is shown in Fig. 2b.
During each simulation PEMW was excited by initial perturbation of the vertical components of magnetic vortex ζ and magnetic field H z .The spatial distributions of these functions in the form of Gaussian are as follows:  bations.In particular, an initial perturbation in the form of a bell-like vertical component of the magnetic vortex (Fig. 4a) leads to the formation of the same shape of vertical magnetic field (Fig. 4c).An excitation by means of the initial bell-like vertical component of the magnetic field (Fig. 3a) results in a quadrupole magnetic vortex (Fig. 3c).The combined excitation by means of initial vertical components of both magnetic field and vortex (Figs.5a-c Fig. 5g demonstrates a new and interesting behaviour of velocity U z at Z = 300 km, during night time conditions, with a remarkable nonlinearity and due to a presence of combined source; that is, if both H z and ζ are present at t = 0. We use the initial perturbations of normalised H z and ζ with maximum values equal to 0.05, which corresponds to a magnetic field of order of 1000 nT.However this rather extreme value is still in a range of possible amplitudes for PEMW vortex structures and in agreement with evaluations by Aburjania et al. (2005).The corresponding value of the vertical current is of order of J ∼ 3 × 10 −7 A/m 2 .Based on results which are shown on Fig. 5g, we can conclude, that under such con-ditions, a vertical displacement of a medium during a characteristic time of order of 10 3 s reaches the value of order of 100 m.However, the vertical displacement corresponding to the combined sources at Z = 600 km is negligibly small (not shown here).Note that in case of combined source (i.e. both H z and ζ), the nonlinear 'effective source' of excitation of U z in the near-pole region is provided, which is evident from Eq. (B3).It is further recalled that, at the middle latitudes, both linear and nonlinear excitations of the vertical velocity component are possible (see the third Equation in (17) and relation ( 22)).A rough estimation shows that the vertical displacement at height Z = 300 km reaches a value of order of 1 km.We have found that, within the ionosphere F region (i.e. at Z = 300 km) the pulse-like PEMW can retain their structure up to 10 3 s.At the higher altitudes (Z = 600 km) PEMW could only be observed during the night time and they preserve their shapes up to 10 2 s due to the much smaller Pedersen conductivity in the magnetosphere than in the ionosphere.The magnetic field of PEMW has the characteristic values: H z of order one to few nT, H x,y of order of few to 10 nT.These values are well above the sensitivity of the modern magnetometers, and therefore, they can be detected (Prattes et al., 2011, see also Section 7).
Here ζ (t = 0) max , H z (t = 0) max and x ζ , y ζ are the maximum values and widths of spatial distributions in X, Y directions for ζ and H z , respectively.x 0H , y 0H are the shift values of the maximum of H z distribution from the coordinate system origin (which coincides with the pole) and θ is the rotation angle in respect to the X axis.If θ is equal to zero, the corresponding axis coincides with X axis.
The results of numerical simulations for parameters of ionosphere at different heights (see Table 1) are shown in Figs.3-6.The initial perturbations with bell-like spatial distribution have been taken for ζ , H z , or both variables simultaneously.The spatial scale is the Earth's radius, that is, R = 6.4 × 10 3 km and the temporal scale is t 0 = R/V A0 ≈ 1.5 × 10 3 s for Z = 300 km.At the night-time and Z = 600 km the temporal scale is t 0 = R/V A0 ≈ 1.5 × 10 2 s.The wind velocity U 0x,y = 0.01V A0 , and their influence is negligible at the considered time intervals.The propagation of the excited PEMWs depend significantly on daytime, height and excitation sources.The calculation for the night-time is shown on Figs.3-6 (the daytime case is not shown); dependences on height for Z = 300 km and Z = 600 km are shown on Figs.3-5 and Fig. 6, respectively.The excitations of PEMW separately by vertical magnetic field H z and magnetic vortex ζ are presented in Figs.3-4.Figures 5 and 6 illustrate the results of calculations for the combined source.In these cases, initial perturbations of both vertical magnetic field and magnetic vortex component remain bell shapes, but the corresponding spatial distributions are shifted and rotated relatively to each other.Note that at Z = 300 km at  day conditions, the values of magnetic field perturbations are about half of order less than at night, because the Pedersen conductivity in F region at night is larger and the effective diffusion, respectively, smaller, than at daytime (Alperovich and Fedorov, 2007).The PEMW excited by the pulse-like initial perturbation can exist for about ∼ 500 s (see Figs. 3-5) and even up to 10 3 s, for the heights within F region of the ionosphere.At higher height (i.e. at Z = 600 km) the characteristic time for retain PEMW is smaller (i.e.< 10 2 s).This is due to the smaller value of the Pedersen conductivity σ P ∼ 10 3 s −1 and higher Alfvén speed.Excitations of PEMW at daytime at Z = 600 km are negligibly small and not shown here.
The results are presented for the initial normalised value of the magnetic field Hz = 0.01, what corresponds to a dimensional value H z ∼ 100 nT, that is possible for PEMW vortices (Aburjania et al., 2005).The chosen value of normalised vortex ζ = (curlH) z ∼ 0.01 means that the corresponding vertical current is, by order of magnitude, J z ∼ 10 −7 A m −2 .It should be noted that the considered vertical currents are smaller than currents accompanying some typical events of the magnetosphere-ionosphere coupling.An example of such a current is one accompanying magnetic ULF disturbances in the night-time auroral region, caused by magnetic impulse events (Pilipenko et al., 1999).Further, as it is seen from the comparison of Figs. 5 and 6, the values of H z and ζ at Z = 300 km, after passing corresponding characteristic time -depending on a type of excitation -are comparable to or few times larger than at Z = 600 km (compare Figs. 3,4,.The same also concerns the values of horizontal velocities (compare Fig. 5f with Fig. 6d).Note that the Alfvén velocity, used for normalisation of the velocity components, is about one order of value larger at altitude Z = 600 km than at Z = 300 km (see Table 1).Respectively, the timescale t 0 is much less at altitude Z = 600 km (t 0 ∼ 150 s) than at Z = 300 km (t 0 ∼ 1400 s).The shapes of Secondly we obgeneral set of equations for PEMW in a weakly onosphere, based on 'nonlinear Ohm's law'.Also, riable (vertical component of velocity) is included sponding nonlinear system.Within this approach, tral -grid Arakawa type method has been devele numerical realisation of the PEMW model.This lises splitting of the calculations into two steps l factors.The pure transfer is realised by a finitescheme at the second half-step, while all other efesented by means of two dimensional purely real t the first half-step.This method provides high staccuracy, in spite of the simultaneous presence of linearity and a number of derivatives in different ns.The analytical and numerical simulation yield and temporal characteristic scales of PEMW exhich, by the orders of values, agree with of ob- (Burmaka et al., 2006;Aburjania & Chargazia, 2011), and the connection between the spatial shapes of excited structures and the initial perturbations of PEMW.The developed model can be useful in the preparation of an experimental methodology for registration and identification of PEMWs and their sources by ground and space-based instruments.In the latter case, there are well-established methods.With regard to measurements of mechanical movements, they require additional effort, since speed is low, that is combined with the relatively small magnitudes of these displacements.Here it is appropriate to provide selected data on the possibilities of relevant methods for ionospheric disturbances measurements.The DOPE (Doppler Pulsation Experiment) HF Doppler sounder located near Tromso, Norway has a spatial resolution of the order of 4 km, for F-region reflection height of 250 km and a sounder frequency of 4.45 MHz (see e.g.Wright et al., 1997).The method of GPS differential TEC (see e.g.Borries et al., 2007) has a spatial resolution of about 50-100 m.At last, the CHInese MAGnetometer (CHIMAG) fluxgate magnetometer chain has an accuracy of 8 pT at a temporal resolution of 1 Hz (see e.g.Prattes et al., 2011).The spatial dimensions of the objects under consideration and vertical displacements are of the order of 1000 km and 1 km, respectively, and the magnetic fields are of the order of 10 nT at a frequency of about 0.01 Hz.In view of this it can be concluded that based on the multipoint synchro- the excited PEMW also depend on the type of initial perturbations.In particular, an initial perturbation in the form of a bell-like vertical component of the magnetic vortex (Fig. 4a) leads to the formation of the same shape of vertical magnetic field (Fig. 4c).An excitation by means of the initial bell-like vertical component of the magnetic field (Fig. 3a) results in a quadrupole magnetic vortex (Fig. 3c).The combined excitation by means of initial vertical components of both magnetic field and vortex (Fig. 5a-c) cases, besides a simple deformation of the horizontal components of the magnetic field and velocities (Fig. 5c, e, f) causes, also a qualitative change in the shape of ζ .Precisely, the spatial distribution of ζ possesses two "dips", around the central maximum (Fig. 5c). Figure 5g demonstrates a new and interesting behaviour of velocity U z at Z = 300 km, during night-time conditions, with a remarkable nonlinearity and due to a presence of combined source; that is, if both H z and ζ are present at t = 0. We use the initial perturbations of normalised H z and ζ with maximum values equal to 0.05, which corresponds to a magnetic field of order of 1000 nT.However this rather extreme value is still in a range of possible amplitudes for PEMW vortex structures and in agreement with evaluations by Aburjania et al. (2005).The corresponding value of the vertical current is of order of J ∼ 3 × 10 −7 A m −2 .Based on results which are shown on Fig. 5g, we can conclude, that under such conditions, a vertical displacement of a medium during a characteristic time of order of 10 3 s reaches the value of order of 100 m.However, the vertical displacement corresponding to the combined sources at Z = 600 km is negligibly small (not shown here).Note that in case of combined source (i.e. both H z and ζ ), the nonlinear "effective source" of excitation of U z in the near-pole region is provided, which is evident from Eq. (B3).It is further recalled that, at the middle latitudes, both linear and nonlinear excitations of the vertical velocity component are possible (see the third Eq. 17 and relation 22).A rough estimation shows that the vertical displacement at height Z = 300 km reaches a value of the order of 1 km.We have found that, within the ionosphere F region (i.e. at Z = 300 km) the pulse-like PEMW can retain their structure up to 10 3 s.At the higher altitudes (Z = 600 km) PEMW could only be observed during the night-time and they preserve their shapes up to 10 2 s due to the much smaller Pedersen conductivity in the magnetosphere than in the ionosphere.The magnetic field of PEMW has the characteristic values: H z of the order of one to a few nT, H x,y of the order of a few to 10 nT.These values are well above the sensitivity of the modern magnetometers, and therefore, they can be detected (Prattes et al., 2011, see also Sect.7).

Discussion and conclusions
The new nonlinear model of vortex PEMWs for arbitrary height, which does not necessarily belong to only one of the "typical" D, E or F regions, has been developed.In contrast to previous models, the present model is suitable for either of these regions, or to "intermediate" regions between D and E and E and F layers within the ionosphere.Secondly, we obtained the general set of equations for PEMW in a weakly nonlinear ionosphere, using "linear" and then "nonlinear" Ohm's law.Also, an extra variable (vertical component of velocity) is included in the corresponding nonlinear system.Within this approach, a new spectral -grid Arakawa type method has been developed for the numerical realisation of the PEMW model.This method utilises splitting of the calculations into two steps by physical factors.The pure transfer is realised by a finite-difference scheme at the second half-step, while all other effects are presented by means of two-dimensional purely real sine FFT at the first half-step.This method provides high stability and accuracy, in spite of the simultaneous presence of vector nonlinearity and a number of derivatives in different combinations.The analytical and numerical simulation yield the spatial and temporal characteristic scales of PEMW excitations, which, by the orders of values, agree with observations (Burmaka et al., 2006;Aburjania and Chargazia, 2011), and the connection between the spatial shapes of excited structures and the initial perturbations of PEMW.The developed model can be useful in the preparation of an experimental methodology for registration and identification of PEMWs and their sources by ground-and space-based instruments.In the latter case, there are well-established methods.With regard to measurements of mechanical movements, they require additional effort, since speed is low, that is combined with the relatively small magnitudes of these displacements.Here it is appropriate to provide selected data on the possibilities of relevant methods for ionospheric disturbances measurements.The DOPE (Doppler Pulsation Experiment) HF Doppler sounder located near Tromsø, Norway has a spatial resolution of the order of 4 km, for F region reflection height of 250 km and a sounder frequency of 4.45 MHz (see e.g.Wright et al., 1997).The method of GPS differential TEC (see e.g.Borries et al., 2007) has a spatial resolution of about 50-100 m.At last, the CHInese MAGnetometer (CHIMAG) fluxgate magnetometer chain has an accuracy of 8 pT at a temporal resolution of 1 Hz (see e.g.Prattes et al., 2011).The spatial dimensions of the objects under consideration and vertical displacements are of the order of 1000 km and 1 km, respectively, and the magnetic fields are of the order of 10 nT at a frequency of about 0.01 Hz.In view of this it can be concluded that based on the multipoint synchronised measurements using differential GPS TEC and chains of magnetometers and advanced data processing methods, it is possible to identify PEMWs in the ionosphere.
It should be noted that large-scale experiments to study the phenomenon of PEMW have not been carried out yet.And this, despite the fact that, according to Aburjania et al. (2005), PEMWs are natural ionospheric oscillations, and therefore, they are a primary link in the chain of the ionosphere's response to disturbances both from below and from above (earthquakes, hurricanes, high-frequency ionosphere heating, magnetosphere's and solar wind impacts).
The theoretical background of relevant experiments, therefore, includes (a) classification of various modes of hydromagnetic waves in the ionosphere (Aburjania et al., 2002(Aburjania et al., , 2003b;;Aburjania and Chargazia, 2011); (b) estimations of wave parameters in linear theory (Aburjania et al., 2002(Aburjania et al., , 2003b;;Rapoport et al., 2007;Aburjania and Chargazia, 2011); and (c) systematisation of estimates to prepare a measurement scheme (Aburjania et al., 2002(Aburjania et al., , 2003b;;Aburjania and Chargazia, 2011).In particular, it has been found that the phase speed is of the order of 1 km s −1 to 10 3 km s −1 and even higher, depending on height, while the wavelength is of order 10 3 to 10 4 km.The most promising estimate relates to the magnitude of magnetic disturbances caused by PEMW: it is of the order of up to 10 2 nT and in certain cases -up to 10 3 nT in situ for vortex structures propagating along the F layer (Aburjania et al., 2005).
There are a number of experiments devoted to measurements of the excitation of ULF hydromagnetic waves by different sources in the ionosphere (see e.g.Burmaka et al., 2006).Wave disturbances of the geomagnetic field associated with distant rocket launches were measured using the modified ionosonde "Basis", the fluxgate magnetometer, and the radar of incoherent scattering at the Ionosphere observatory, near Kharkov, Ukraine.In total, measurements of ionospheric wave disturbances from about 140 distant missile launches were carried out for about 10 yr.
The observations were made at various distances from the launching sites (Baikonur, Plesetsk, et al.): from 700 km to about 10 4 km.It was found that (a) perturbations with a speed of 10-20 km s −1 are quite often observed at distances of up to 2300 km from the rocket; (b) there were 3 groups of speeds of disturbances: 0.5-0.7 km s −1 , 2-3 km s −1 , and 10-25 km s −1 ; (c) there were 3 groups of periods of dominant geomagnetic micropulsations: about 6, 10, and 20 min, their amplitudes on the ground attained a value of 3-5 nT; (d) the launch-resulted waves were most pronounced at altitudes 150-350 km; (e) relative amplitudes of disturbances of the electron concentration reached 5-7 %.The analysis of such observations shows that the difficulty in identifying waves is that the frequencies of typical PEMW packets belong to the ULF band and may be buried in the noise that exists there.To be able to detect PEMW packets in this background, it seems possible to use sophisticated data analysis techniques that have worked well in the study of seismic and oceanic wave processes: directional spectra estimation (see e.g.Sun et al., 2005); study of a packet modification due to the coaction of nonlinearity, dispersion, and dissipation; using the concept of so called "coda waves" (see e.g.Sèbe et al., 2005).These techniques become more efficient when we take into account the internal multi-scale nature of the processes under study.For example, wavelet packet transform where appropriate statistical decision rule is applied to subset of wavelet coefficients (Lyubushin, 2007) gives a tool to detect complex changes, say, in magnetometer recordings.
Finally, it should be emphasised that the present state of theoretical and modelling studies in the area of PEMW puts on the agenda the question of systematic experiments on the identification of PEMWs under various geomagnetic conditions, on revealing their sources, correlations with space weather events and catastrophic events throughout the atmosphere-ionosphere system.This requires the organisation of measurements accounting for the large-scale (10 3 -10 4 km) character of the PEMWs on the base of a planetary size network of ionospheric/magnetospheric observatories.
The proposed theory can be extended by including the connection between particular coupling processes in the "magnetosphere-ionosphere-atmosphere-lithosphere (MIAL)" system and corresponding sources of full-spectrum MHD waves/PEMW in the form of the effective "external currents/fields".
Here we sum up the results of the present work.
1.A new nonlinear analytical model of PEMW, valid in the local β-plane approximation for arbitrary altitudes in D, E, and F ionospheric regions, as well as for intermediate altitudes is developed.We propose for the first time for the PEMW the method of consistently utilising the "nonlinear Ohm's law" and the series expansion in the relatively small non-stationary part of the magnetic field.
Yu. Rapoport et al.: Excitation of planetary electromagnetic waves in the inhomogeneous ionosphere 2. A new numerical method for the simulations of the evolution of PEMW is developed.This is a highly stable and efficient hybrid finite difference-spectral method, based on splitting by physical factors.
3. The shapes of the excited PEMW depend on the types of initial conditions or sources.In particular, the excitation of PEMW by an initial bell-like vertical component of the magnetic vortex leads to the formation of a bell-like vertical magnetic field.The excitations by means of an initial bell-like vertical component of the magnetic field results in a quadrupole magnetic vortex.
4. We have found that diffusion processes are rather important in the spreading of PEMW excitation in the F region of the ionosphere, while in most previous papers, such as Aburjania et al. (2005), this factor was underestimated.The vertical component of the velocity of PEMW is not essential for an altitude of Z = 600 km.At altitude Z = 300 km this component provides a vertical displacement of ∼ 100 m on a timescale ∼ 10 3 s at the high latitudes, while the estimate of vertical displacement is ∼ 1 km at the middle latitudes, where both linear and nonlinear wave excitation of PEMW is possible.
5. The components of the magnetic field of PEMW are of the order of one to a few nT for the vertical component and from a few to tens of nT for the horizontal component; therefore, it can be measured by modern instruments.

Appendix C Equations for the near-equator region
In the near-equator region, we obtain the following set of normalised equations from the system of Eqs. ( 17) accounting for the relations ( 18)-( 24) in the limiting case λ → 0:

Fig. 1 .
Fig. 1.Coordinate system of the model.X axis pierces the plane of figure.R is the radius of the Earth, 0 is the angular velocity of the Earth rotation; λ = 0 and λ = π/2 correspond to equator and pole, respectively.

Fig. 2 .
Fig.2.The dependence of ionosphere conductivities σ 0 , σ H , σ P and σ eff as a function of height, z, are shown in (a).The values of coefficients q 1,2 vs. height, z, are plotted at (b). q 1 and q 2 are coefficients in the expressions, corresponding to the "transverse" electric field, determined by Hall current and "longitudinal" electric field (see Eq. 10).
Fig. 3: Excitation of PEMW by means of vertical component of magnetic field.Spatial distributions of normalised components of magnetic field (H x , H y , H z ) and ζ at different times are shown at hight Z = 300 km (night time).X and Y are normalised by R = 6.4 × 10 3 km.Components of magnetic field H are normalised by H 0 = 3 × 10 4 nT.Dimensionless size of the source are ∆x = 0.25 and ∆y = 0.17; ∆X 0H = ∆Y 0H = 0 and θ = 0. PEMW was excited by initial perturbation of H z (H z (t = 0) max = 0.01).

Fig. 3 .
Fig. 3. Excitation of PEMW by means of vertical component of magnetic field.Spatial distributions of normalised components of magnetic field (H x , H y , H z ) and ζ at different times are shown at hight Z = 300 km (night-time).X and Y are normalised by R = 6.4 × 10 3 km.Components of magnetic field H are normalised by H 0 = 3 × 10 4 nT.Dimensionless size of the source are x = 0.25 and y = 0.17; X 0H = Y 0H = 0 and θ = 0. PEMW was excited by initial perturbation of H z (H z (t = 0) max = 0.01).
Fig. 4: Excitation of PEMW by means of vertical component of magnetic vortex.The spatial distributions of normalised ζ and components of magnetic field (H x , H y , H z ) are taken at hight Z = 300 km (night time).Amplitudes of initial excitations are equal to ζ(t = 0) max = 0.01 and H z (t = 0) max = 0. Dimensionless size of the source are ∆x ζ = ∆y ζ = 0.17.The scaling of X, Y and normalised values of H are the same as in Fig. 3.
) cases, besides a simple deformation of the horizontal components of the magnetic field and velocities (Figs.5c, e, f) causes, also a qualitative change in the shape of ζ.Precisely, the spatial distribution of ζ possesses two 'dips', around the central maximum (Figs.5c).

Fig. 4 .
Fig. 4. Excitation of PEMW by means of vertical component of magnetic vortex.The spatial distributions of normalised ζ and components of magnetic field (H x , H y , H z ) are taken at hight Z = 300 km (night-time).Amplitudes of initial excitations are equal to ζ (t = 0) max = 0.01 and H z (t = 0) max = 0. Dimensionless size of the source are x ζ = y ζ = 0.17.The scaling of X, Y and normalised values of H are the same as in Fig. 3.
Yu. Rapoport et al.: Excitation of planetary electromagnetic waves in the inhomogeneous ionosphere ort: Excitation Of Planetary Electromagnetic Waves In The Inhomogeneous Ionosphere 11 (a) ζ(t = 0) (b) ζ(t = 0.06) (c) Hx(t = 0.06) (d) Ux(t = 0.06) spatial distributions of normalised ζ, H x and U x components of magnetic field and velocity are taken at hight m (night time).PEMW was excited by vertical component of magnetic vortex ζ z (ζ(t = 0) max = 0.01).The scaling d normalised value of H x is the same as in Figs. 3. The velocity component U x is normalised by V A = 46 km/s.sion and Conclusions onlinear model of vortex PEMWs for arbitrary ich does not necessarily belong to only one of the , E or F regions has been developed.In contrast us models, the present model is suitable for either ions, or to 'intermediate' regions between D and d F layers within the ionosphere.

Fig. 6 .
Fig. 6.The spatial distributions of normalised ζ , H x and U x components of magnetic field and velocity are taken at hight Z = 600 km (night-time).PEMW was excited by vertical component of magnetic vortex ζ z (ζ (t = 0) max = 0.01).The scaling of X, Y and normalised value of H x is the same as in Fig. 3.The velocity component U x is normalised by V A = 46 km s −1 .

Table 1 :
The ionosphere night time parameters.

Table 1 .
The ionosphere night-time parameters.