Whistler waves produced by monochromatic currents in the low nighttime ionosphere

We use a full-wave approach to find the field of monochromatic whistler waves which are excited and propagating in the low nighttime ionosphere. The source current is located in the horizontal plane and can have arbitrary distribution over horizontal coordinates. The ground-based horizontal magnetic field and electric field at 125 km are calculated. The character of wave polarization on the ground surface is investigated. The percentages of source energy supplied to the Earth-ionosphere waveguide and carried upward ionosphere are estimated. Received results are important for the analysis of ELF/VLF emission 5 phenomena observed both on the satellites and on the ground.


Introduction
ELF/VLF waves, which propagate in the ionosphere in whistler mode, are an important part of the ionosphere dynamics.
Such waves can be emitted by various natural phenomena such as atmospheric lightning discharges and volcanic eruptions, magnetospheric chorus and hiss. Artificial ELF/VLF waves have been produced by ground based transmitters and by modulated 10 HF heating of the ionosphere current system responsible for S q variations or auroral electrojet, which is by a now well known technique.
Several numerical methods have been developed for calculating of whistler wave fields in the Earth's ionosphere (Pitteway, 1965;Wait, 1970;Bossy, 1979;Nygre'n, 1982;Budden, 1985;Nagano et al., 1994, Yagitani et al., 1994Shalashov and Gospodchikov, 2011). One of the main difficulties is numerical instabilities caused by a large imaginary part of the vertical wave 15 number. General full-wave analysis, including the problem of numerical 'swamping' of the evanescent wave solutions, was made, for example, by Nygre'n (1982), Nagano et al. (1994), Budden (1985). Full wave calculation of ELF/VLF propagation from a dipole source located in the lower ionosphere has been made by Yagitani et al. (1994). The idea of recursive calculation of mode amplitudes was developed and used for an arbitrary configuration of the radiating sources by Lehtinen and Inan (2008).
Nevertheless, finding fields created by both natural and artificial ELF/VLF radiating sources is still very actual. 20 In this paper, we use numerical methods to find the field of a whistler wave generated and propagating in low night ionosphere. We use a technique known as the two-point boundary-value problem for ordinary differential equations (Kierzenka and Shampine, 2001). Using this technique in early work (Bespalov and Mizonova, 2017;Bespalov et al., 2018) has provided numerically stable solutions of a complete system of wave equations for arbitrary altitude profiles of plasma parameters and for arbitrary angles of wave incidence. Here, we find a wave field created by a monochromatic source current located in the low night ionosphere. As an example of calculations we use current distributions similar to those simulated by HF heating of the auroral electrojet (Payne et al., 2007). The obtained results are important for analysis of the ELF/VLF emission phenomena observed both in the ground-based observatories and on board of satellites.

2 Basic equations
We consider a whistler wave which is excited and propagating in the layer 0 ≤ z ≤ z max of the non-homogeneous stratified ionosphere. We choose a coordinate system with vertical upward z axis and x, y axes in horizontal plane, suppose that plasma parameters depend on coordinate z, plane z = 0 corresponds to the ground surface, above the boundary z = z max ionosphere plasma is close to homogeneous, the ambient magnetic field B 0 belongs to the y, z plane and is inclined at an angle ϑ to the z 10 axis. We assume that external currents have monochromatic dependence on time and flow in the source plane z = z s (1) At first, we use the Fourier composition of the source current density over the horizontal coordinates and wave electric and magnetic fields at each altitude and find field amplitudes E(n ⊥ , z), H(n ⊥ , z) corresponding to the horizontal wave vector component k ⊥ = k 0 n ⊥ , k 0 = ω/c.
Here we use SI units for E and modified units for H = Z 0 H SI (Budden, 1985), where Z 0 = µ 0 /ε 0 is the impedance of free 15 space. Then we write the Maxwell equations where c is the speed of light in free space andε is permittivity tensor, which yield a set of four equations for the horizontal components E ⊥ (n ⊥ , z), H ⊥ (n ⊥ , z) (in a case of source-free medium see, e.g., Budden, 1985;Bespalov et al., 2018;Mizonova, 2019) dF/dz =MF + Z 0 Jδ(z − z s ) .
Here we have taken into attention that the horizontal refractive index of the wave propagating through the stratified medium is conserved due to Snell's law. In Eq. (5) F(n ⊥ , z) and J(n ⊥ , z) are four-component column vectors M is a matrix of which the elements m ij are expressed in terms of components of the transverse wave vector k ⊥ = k 0 n ⊥ , ε, η, g are elements of the permittivity tensor which depends on the z coordinate Mizonova, 2017, Bespalov et al., 2018), ε zz = εsin 2 ϑ + ηcos 2 ϑ.

5
To solve the system (5), (6) we define four boundary conditions. We write two of them on the plane z = 0 assuming the ground surface to be perfect conductive We write two other conditions on the plane z = z max excluding wave energy coming from above. To clarify them we express the field vector column F above the boundary z = z max as sum of four wave modes Here A j = const, k zj are four roots of local dispersion relation and P j are four corresponding local polarization vectors. We 10 mention that values k zj and vectors P j are the solution of Eqs. (5), (8) for homogeneous plasma without sources. Assuming that indices 2 and 4 correspond to coming from above propagating and non-propagating wave modes (imaginary parts of k z2 and k z4 are negative) we write Solving the set of Eq. (5) with known source current density (6) and boundary conditions (7), (9) we can find the field of plane wave with horizontal wave vector k ⊥ = k 0 n ⊥ in the layer 0 < z < z max . Then, we use the inverse transform to calculate the total field.
We take into account that out of the plane z = z s the source current density (6) is equal to zero, so Eq. (5) becomes To solve Eq. (11) in layers 0 ≤ z < z s and z s < z ≤ z max we apply packaged solver Mathlab's bvp4c and use a method of known as the two-point boundary-value problem for ordinary differential equations (Kierzenka and Shampine, 2001). The solver solution starts with an initial guess supplied at an initial mesh points and changes step-size to get the specified accuracy.

5
At first we find two linearly independent solutions F 1 and F 2 of Eq. (11) in the layer 0 ≤ z < z s completing the boundary condition (7) on the plane z = 0 for arbitrary conditions on the plane z = z s − 0. For example, we use four conditions Then we write the general solution in the layer 0 ≤ z ≤ z s as sum Similarly, we find two linearly independent solutions F * 1 and F * 2 of Eq. (11) in the upper layer z s < z ≤ z max completing the boundary condition (9) on the plane z = z max . By arbitrary conditions on the plane z = z s + 0 we have A 2 = 0 , A 4 = 0, write the general solution in the layer z s < z ≤ z max as sum Integrating Eq. (5) over z coordinate in a narrow layer (z s − 0, z s + 0) we find a condition connecting solutions (12) and (13) 15 That condition yields four algebraic equations for the coefficients a, a * , b, b * . Thus, finding those coefficients we obtain the wave field F (6) in the layer 0 ≤ z ≤ z max . In particular, the fields on the ground surface can be expressed as where n 0z = 1 − n 2 ⊥ , E x ,y are electric strength components of incident wave in coordinate system with z -axis along ambient magnetic field. The wave polarization near the ground surface is determined by The electric field at the altitude z = z max can be expressed as The coordinate dependence of the wave field can be found from the inverse Fourier transform (10). The vertical energy flux (Poynting vector) is and the total power of source is Now we present the results of numerically computed solution of the set (5), (6) under model conditions for the nighttime ionosphere.

The ionosphere data and calculation results
In the calculations, we use the altitude profiles of the plasma density shown in Fig. 1a, and the collision frequencies be- x and y components j x = j y .
We assume that currents occupy a volume which has a shape of a horizontal pancake and use for calculations a Gaussian distributions over x and y coordinates J(r ⊥ ) = J max exp −x 2 /2L 2 x − y 2 /2L 2 y with characteristic horizontal sizes L x = 12 km, L y = 70 km. Corresponding current distribution (2) in Fourier-space is also Gaussian and has a form 20 J(n ⊥ ) = J 0 exp −k 2 0 L 2 x n 2 x /2 − k 2 0 L 2 y n 2 y /2 , where J 0 = 2πL x L y J max . We present the results of field calculation in Fourierspace (Figs. 2-3) and in coordinate space (Fig. 4). The dependences of amplitude of horizontal magnetic field H ⊥ (n x , n y )/E 0 on ground surface z = 0 and amplitude of electric field E(n x , n y )/E 0 at altitude z = z max are presented in Figs. 2a,b. The field  Figure 4a shows current density normalized by the value J max . Figure 4b shows electric field at the altitude z = z max . Figure 4c shows horizontal magnetic field on ground surface z = 0. Both electric and magnetic fields are normalized by the value Z 0 J max . Figure 4d shows polarization parameter φ. The arrow shows the current direction.

Discussion
We use a full-wave approach to find the field of monochromatic whistler waves which are excited and propagating at night in the strongly inhomogeneous low ionosphere. A source current is assumed to be located in horizontal plane and to have in this plane an arbitrary space distribution. At first, we consider a plane wave with the horizontal components of the refractive index n ⊥ generated by the current j(r ⊥ ) ∼ e ik0n ⊥ r ⊥ . The set of wave equation in the layer 0 < z < z max for each n ⊥ -component is 15 completed by boundary conditions assuming a perfect conductivity of ground surface and excluding wave energy coming on the upper boundary z = z max from above. The method known as the two-point boundary-value problem for ordinary differential equations (Kierzenka and Shampine, 2001) is applied to find the solutions of wave equations above and below the source plane. Then we connect these solutions using source current distribution. Inverse Fourier transform yields space dependence of the wave field.
As an example, we calculate the fields created by varying at a frequency of 3 kHz and flowing at z s = 80 km horizontal current, with coincided x and y components of current density j x = j y . We use a Gaussian distributions of source current density over x and y coordinates with characteristic horizontal sizes L x = 12 km, L y = 70 km. We mention that the model of a plane source current can also be effective in a more general case of current layer with small thickness ∆z λ z ∼ 60 km.
The ground-based horizontal magnetic field and the electric field at 125 km are calculated both in Fourier (n x , n y ) and coordinate (x, y) spaces. Since a wave can achieve the ground surface in penetrating mode in case n ⊥ ≤ 1. The magnetic field H ⊥ (n ⊥ , z = 0) is noticeably non-zero for Fourier components with n ⊥ < 1 and is practically equal to zero for Fourier 5 components with n ⊥ 1. Waves with n ⊥ < 1 are right hand polarized (see Fig. 2c), which is typical for whistlers. However, if the horizontal component of the refractive index has an order of unit, then the magnetic field H ⊥ (n ⊥ , z = 0) can increase two or three times (see Fig. 2a). The polarization parameter φ of such components can be negative similar to the left hand polarized waves (see Fig. 2c). We mention that if the horizontal size of radiating currents is small enough L ≤ 1/k 0 ∼ 15 km (for wave with frequency 3 kHz) then Fourier components with n ⊥ ∼ 1 can make a noticeable contribution into the whole field 10 and change the character of polarization.
The magnetic field H ⊥ (r ⊥ , z = 0) is predominantly localized under the source (Fig. 4c). The electric field at 125 km occupies a spot of several times larger size (Fig. 4b). The polarization parameter φ at z = 0 is positive almost everywhere, but becomes negative at large distances from the source across the current direction (Fig. 4d). We mention that the used in our calculation current distributions can be similar to electrojet currents modulated in D-region by the HAARP HF heating facility 15 (Keskinen and Rowland, 2006;Payne et al., 2007). For example, according to data collected during an experimental campaign run in April 2003 and results of numerical simulations (Payne et al., 2007;Lehtinen and Inan, 2008), the maximum change in modulated conductivity occupies approximately 10 km over the height and occurs at altitude ∼ 80 km. Pedersen and Hall conductivities approximately coincide so if ambient electrojet field is directed along x axis, then j x ≈ j y . The maximum surface density of modulated currents has an order J max ∼ 10 −6 ÷ 10 −5 Am −1 . Using in our calculations the magnitude of current density J max ∼ 5 · 10 −6 Am −1 yields the total power of the source ∼ 36W, the ground-based horizontal magnetic field under the source B ⊥ ∼ 1pT and the electric field at the altitude of 125 km above the source E ∼ 400µVm −1 . The magnitude of the magnetic field is similar to the field measured at VLF sites in the immediate vicinity of the HAARP heating facility (Payne et al., 2007) and calculated by Lehtinen and Inan (2008). The maximum vertical energy flux (Poynting vector) at the altitude 5 of 125km is ∼ 3, 2nWm −1 and total power is ∼ 17W. The percentages of source energy supplied by the Earth-ionosphere waveguide and carried upward ionosphere depends on altitude profile of ionosphere plasma. If low boundary of ionosphere is sharp enough then sufficient part of the source energy is supplied to the Earth-ionosphere waveguide. In a case of smooth ionosphere low boundary sufficient part of the source energy is carried upward. By used for calculation altitude profile of plasma density about half of the source energy is carried upward, approximately twenty percent of the energy is supplied to the 10 Earth-ionosphere waveguide and approximately thirty percent of the energy is absorbed.

Conclusions
We find a field of monochromatic whistler waves which are excited and propagating in the low nighttime ionosphere. Using a technique known as the two-point boundary-value problem for ordinary differential equations enables to find numerically stable solutions of full set of the wave equations applying to conditions of inhomogeneous ionosphere at altitudes below 125 km. 15 Above this altitude the ionosphere plasma is slightly inhomogeneous, hence approximate methods are suitable. As example, this calculation technique is applied to the problem of ELF/VLF waves radiation from modulated HF-heated electrojet currents.
At first we consider a plane wave with known horizontal component of the refractive index, find a wave field and analyze a character of wave polarization on the ground surface. Then we use Fourier transform to find a total field. The ground-based magnetic field is predominantly located under the source. The wave field at 125 km occupies a spot of several times larger 20 size. The polarization almost everywhere is typical for whistlers (right hand), but at large distances from the source across the current direction can become left hand. About half of the source energy is carried upward, approximately twenty percent of the energy is supplied to the Earth-ionosphere waveguide and approximately thirty percent of the energy is absorbed. The obtained values are in a good agreement with ground and satellite observations and known calculation results. Using model of plane horizontal source currents can be generalized for the arbitrary altitude source distribution.

25
Data availability. The paper is theoretical and no new experimental data are used. The data are taken from International Reference Ionosphere model (Bilitza and Reinisch, 2007) (https://ccmc.gsfc.nasa.gov/modelweb/models/iri2016_vitmo.php).
All figures are obtained from numerical calculation in MATLAB codes.
Author contributions. VM produced the calculations, analyzed results, and wrote the paper. PB proposed the problem, discussed results, and wrote the paper.