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 finite 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 proportion in which source energy supplies the Earth–ionosphere waveguide or flows upward can be adjusted by distribution of the source current. Received results are important for the analysis of ELF/VLF emission phenomena observed both on the satellites and on the ground.

1 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 high-frequency (HF) heating of the ionosphere current system responsible for S_q variations or the auroral electrojet, which is by now a well-known technique. ELF/VLF waves modulated by the HAARP heating facility can be injected in the Earth–ionosphere waveguide as far as 4400 km (Moore et al., 2007) and also into space (Inan et al., 2004).

Several numerical methods have been developed for calculating 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., 1994; Shalashov and Gospodchikov, 2011). One of the main difficulties is numerical instabilities caused by a large imaginary part of the vertical wave number. General full-wave analysis, including the problem of numerical “swamping” of the evanescent wave solutions, was performed, for example, by Nygre'n (1982), Nagano et al. (1994), and Budden (1985). A traditional approach in full-wave analysis is to divide a stratified ionosphere into a number of thin horizontal and homogeneous slabs and then connect the solutions in each slab by applying the boundary conditions. Such a technique has been used by Yagitani et al. (1994) to study ELF/VLF propagation from an infinitesimal dipole source located in the lower ionosphere. 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 relevant.

In this paper, we use numerical methods to find the field of ELF/VLF waves, which have been produced in a low nighttime ionosphere. On the one hand, significant inhomogeneity of plasma parameters, strong wave mode attenuation and effect of wave mode transformation (for example, whistler to vacuum electromagnetic) in the low-altitude nighttime ionosphere make the problem considered difficult enough and fundamentally important. On the other hand, it has practical significance, as an example, for interpretation of numerous experimental results on HF heating which modulate natural ionospheric currents at altitudes of 60–100 km.

In calculations, 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 in the stratified ionosphere for arbitrary angles of wave incidence from above. Here, we find a wave field created by a monochromatic source current located in the low night ionosphere. We examine the influence of peculiarities of current distribution on the proportion in which source energy supplies the Earth–ionosphere waveguide or flows upward. 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 onboard satellites.

2 Basic equations

We consider a whistler wave which is excited and propagating in the layer $\mathrm{0}\le z\le z_{max}$ of the non-homogeneous stratified ionosphere. We choose a coordinate system with a vertical upward z axis and x and y axes in the 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, and the ambient magnetic field B₀ belongs to the y–z plane and is inclined at an angle ϑ to the z axis. We assume that external currents have monochromatic dependence on time and flow in the source plane z=z_s:

$$\begin{array}{}\text{(1)}& \mathit{j}\left({\mathit{r}}_{\perp },z,t\right)=\mathit{J}\left({\mathit{r}}_{\perp }\right)\mathit{\delta }\left(z-z_s\right)e^{-i\mathit{\omega }t}.\end{array}$$

At first, we use the Fourier composition of the source current density over the horizontal coordinates,

$$\begin{array}{}\text{(2)}& \mathit{J}\left({\mathit{n}}_{\perp },z\right)=\int \mathit{J}\left({\mathit{r}}_{\perp },z\right)e^{-i{k}_{\mathrm{0}}{\mathit{n}}_{\perp }{\mathit{r}}_{\perp }}{k}_{\mathrm{0}}^{\mathrm{2}}\mathrm{d}{\mathit{r}}_{\perp },\end{array}$$

and wave electric and magnetic fields at each altitude

$$\begin{array}{}\text{(3)}& {\displaystyle }\begin{array}{rl}& \mathit{E}\left({\mathit{n}}_{\perp },z\right)=\int \mathit{E}\left({\mathit{r}}_{\perp },z\right)e^{-i{k}_{\mathrm{0}}{\mathit{n}}_{\perp }{\mathit{r}}_{\perp }}{k}_{\mathrm{0}}^{\mathrm{2}}\mathrm{d}{\mathit{r}}_{\perp },\\ & \mathit{H}\left({\mathit{n}}_{\perp },z\right)=\int \mathit{H}\left({\mathit{r}}_{\perp },z\right)e^{-i{k}_{\mathrm{0}}{\mathit{n}}_{\perp }{\mathit{r}}_{\perp }}{k}_{\mathrm{0}}^{\mathrm{2}}\mathrm{d}{\mathit{r}}_{\perp },\end{array}\end{array}$$

and find field amplitudes $\mathit{E}({\mathit{n}}_{\perp },z)$, $\mathit{H}({\mathit{n}}_{\perp },z)$ corresponding to the horizontal wave vector component ${\mathit{k}}_{\perp }=k_{\mathrm{0}}{\mathit{n}}_{\perp }$, $k_{\mathrm{0}}=\mathit{\omega }/c$. Here we use SI units for E and modified units for H=Z₀H_SI (Budden, 1985), where $Z_{\mathrm{0}}=\sqrt{{\mathit{\mu }}_{\mathrm{0}}/{\mathit{\epsilon }}_{\mathrm{0}}}$ is the impedance of free space. Then we write the Maxwell equations

$$\begin{array}{}\text{(4)}& \begin{array}{rl}& \mathrm{\nabla }\times \mathit{H}=Z_{\mathrm{0}}\mathit{j}-i{k}_{\mathrm{0}}\widehat{\mathit{\epsilon }}\mathit{E},\\ & \mathrm{\nabla }\times \mathit{E}=i{k}_{\mathrm{0}}\mathit{H},\end{array}\end{array}$$

where c is the speed of light in free space and $\widehat{\mathit{\epsilon }}$ is the permittivity tensor, which yield a set of four equations for the horizontal components ${\mathit{E}}_{\perp }({\mathit{n}}_{\perp },z)$, ${\mathit{H}}_{\perp }({\mathit{n}}_{\perp },z)$ (in a case of a source-free medium, see, e.g., Budden, 1985, Bespalov et al., 2018, and Mizonova, 2019)

$$\begin{array}{}\text{(5)}& \mathrm{d}\mathit{F}/\mathrm{d}z=\widehat{\mathbf{M}}\mathit{F}+Z_{\mathrm{0}}\mathit{I}\mathit{\delta }(z-z_s).\end{array}$$

Here we have taken into account that the horizontal refractive index of the wave propagating through the stratified medium is conserved due to Snell's law. In Eq. (5) $\mathit{F}({\mathit{n}}_{\perp },z)$ and $\mathit{I}({\mathit{n}}_{\perp },z)$ are four-component column vectors

$$\begin{array}{}\text{(6)}& \mathit{F}=\left(\begin{array}{l}{E}_x\\ E_y\\ H_x\\ H_y\end{array}\right),\mathit{I}=\left(\begin{array}{l}{n}_x{J}_z/{\mathit{\epsilon }}_{zz}\\ n_y{J}_z/{\mathit{\epsilon }}_{zz}\\ J_y-J_z(\mathit{\eta }-\mathit{\epsilon })\sin \mathrm{2}\mathit{\vartheta }/\mathrm{2}{\mathit{\epsilon }}_{zz}\\ -J_x+J_{z}ig\sin \mathit{\vartheta }/{\mathit{\epsilon }}_{zz}\end{array}\right).\end{array}$$

$\widehat{\mathbf{M}}$ is a matrix of which the elements m_ij are expressed in terms of components of the transverse wave vector ${\mathit{k}}_{\perp }=k_{\mathrm{0}}{\mathit{n}}_{\perp }$, and ε, η, and g are elements of the permittivity tensor which depends on the z coordinate (Bespalov and Mizonova, 2017; Bespalov et al., 2018), ${\mathit{\epsilon }}_{zz}=\mathit{\epsilon }{\sin }^{\mathrm{2}}\mathit{\vartheta }+\mathit{\eta }{\cos }^{\mathrm{2}}\mathit{\vartheta }$.

To solve the system in Eq. (5), in Eq. (6) we define four boundary conditions. We write two of them on the plane z=0 assuming the ground surface to be perfectly conductive:

$$\begin{array}{}\text{(7)}& E_x\left(z=\mathrm{0}\right)=\mathrm{0},E_y\left(z=\mathrm{0}\right)=\mathrm{0}.\end{array}$$

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 the sum of four wave modes

$$\begin{array}{}\text{(8)}& \mathit{F}\left(z\right)=\sum _{j=\mathrm{1}}^{\mathrm{4}}{A}_j{\mathit{P}}_j\exp \left(i{k}_{zj}\left(z-z_{max}\right)\right).\end{array}$$

Here A_j=const, k_zj are four roots of the local dispersion relation and P_j are four corresponding local polarization vectors. We mention that values k_zj and vectors P_j are the solutions of Eqs. (5) and (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

$$\begin{array}{}\text{(9)}& A_{\mathrm{2}}=\mathrm{0},A_{\mathrm{4}}=\mathrm{0}.\end{array}$$

Solving the set of Eq. (5) with known source current density Eq. (6) and boundary conditions Eqs. (7) and (9), we can find the field of the plane wave with horizontal wave vector ${\mathit{k}}_{\perp }=k_{\mathrm{0}}{\mathit{n}}_{\perp }$ in the layer $\mathrm{0}<z<z_{max}$. Then, we use the inverse transform

$$\begin{array}{}\text{(10)}& {\displaystyle }\mathit{E}\left({\mathit{r}}_{\perp },z\right)=\int \mathit{E}\left({\mathit{n}}_{\perp },z\right)e^{i{k}_{\mathrm{0}}{\mathit{n}}_{\perp }{\mathit{r}}_{\perp }}{\displaystyle \frac{\mathrm{d}{\mathit{n}}_{\perp }}{{\left(\mathrm{2}\mathit{\pi }\right)}^{\mathrm{2}}}},\text{(11)}& {\displaystyle }\mathit{H}\left({\mathit{r}}_{\perp },z\right)=\int \mathit{H}\left({\mathit{n}}_{\perp },z\right)e^{i{k}_{\mathrm{0}}{\mathit{n}}_{\perp }{\mathit{r}}_{\perp }}{\displaystyle \frac{\mathrm{d}{\mathit{n}}_{\perp }}{{\left(\mathrm{2}\mathit{\pi }\right)}^{\mathrm{2}}}},\end{array}$$

to calculate the total field.

3 Description of the solution algorithm

We take into account that out of the plane z=z_s, the source current density Eq. (6) is equal to zero, so Eq. (5) becomes

$$\begin{array}{}\text{(12)}& \mathrm{d}\mathit{F}/\mathrm{d}z=\widehat{\mathbf{M}}\mathit{F}.\end{array}$$

To solve Eq. (11) in layers $\mathrm{0}\le z<z_s$ and $z_s<z\le z_{max}$, we apply packaged solver MATLAB bvp4c and use a method 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 point and changes step size to get the specified accuracy.

At first we find two linearly independent solutions F₁ and F₂ of Eq. (11) in the layer $\mathrm{0}\le z<z_s$ completing the boundary condition Eq. (7) on the plane z=0 for arbitrary conditions on the plane $z=z_s-\mathrm{0}$. For example, we use four conditions $E_x\left(z=\mathrm{0}\right)=\mathrm{0}$, $E_y\left(z=\mathrm{0}\right)=\mathrm{0}$, $E_x\left(z=z_s-\mathrm{0}\right)=E$, and $E_y\left(z=z_s-\mathrm{0}\right)=\mathrm{0}$ for the solution F₁ and four conditions $E_x\left(z=\mathrm{0}\right)=\mathrm{0}$, $E_y\left(z=\mathrm{0}\right)=\mathrm{0}$, $E_x\left(z=z_s-\mathrm{0}\right)=\mathrm{0}$, and $E_y\left(z=z_s-\mathrm{0}\right)=E$ for the solution F₂, where E is constant. Then we write the general solution in the layer $\mathrm{0}\le z\le z_s$ as the sum

$$\begin{array}{}\text{(13)}& \mathit{F}=a{\mathit{F}}_{\mathrm{1}}+b{\mathit{F}}_{\mathrm{2}}.\end{array}$$

Similarly, we find two linearly independent solutions ${\mathit{F}}_{\mathrm{1}}^{*}$ and ${\mathit{F}}_{\mathrm{2}}^{*}$ of Eq. (11) in the upper layer $z_s<z\le z_{max}$ completing the boundary condition Eq. (9) on the plane $z=z_{max}$. By arbitrary conditions on the plane $z=z_s+\mathrm{0}$, we have A₂=0, A₄=0, $E_x\left(z=z_s+\mathrm{0}\right)=E$, and $E_y\left(z=z_s+\mathrm{0}\right)=\mathrm{0}$ and A₂=0, A₄=0, $E_x\left(z=z_s+\mathrm{0}\right)=\mathrm{0}$, and $E_y\left(z=z_s+\mathrm{0}\right)=E$, respectively. We write the general solution in the layer $z_s<z\le z_{max}$ as the sum

$$\begin{array}{}\text{(14)}& \mathit{F}=a^{*}{\mathit{F}}_{\mathrm{1}}^{*}+b^{*}{\mathit{F}}_{\mathrm{2}}^{*}.\end{array}$$

Integrating Eq. (5) over the z coordinate in a narrow layer $\left(z_s-\mathrm{0},z_s+\mathrm{0}\right)$, we find a condition connecting solutions Eqs. (12) and (13):

$$\begin{array}{}\text{(15)}& \mathit{F}\left(z=z_s+\mathrm{0}\right)-\mathit{F}\left(z=z_s-\mathrm{0}\right)=Z_{\mathrm{0}}\mathit{I}.\end{array}$$

That condition yields four algebraic equations for the coefficients a, a^*, b, and b^*. Thus, finding those coefficients, we obtain the wave field F Eq. (6) in the layer $\mathrm{0}\le z\le z_{max}$. In particular, the fields on the ground surface can be expressed as

$$\begin{array}{}\text{(16)}& \begin{array}{rl}& H_{\perp }\left(z=\mathrm{0}\right)=\sqrt{{\left|H_x\left(z=\mathrm{0}\right)\right|}^{\mathrm{2}}+{\left|H_y\left(z=\mathrm{0}\right)\right|}^{\mathrm{2}}},\\ & E_{y^{\prime }}\left(z=\mathrm{0}\right)=\left(\sin \mathit{\vartheta }-{\displaystyle \frac{n_y}{n_{\mathrm{0}z}}}\cos \mathit{\vartheta }\right)n_x{H}_y\left(z=\mathrm{0}\right)-\\ & -\left(\left(\mathrm{1}-n_y^{\mathrm{2}}\right){\displaystyle \frac{\cos \mathit{\vartheta }}{n_{\mathrm{0}z}}}+n_y\sin \mathit{\vartheta }\right)H_x\left(z=\mathrm{0}\right),\\ & E_{x^{\prime }}\left(z=\mathrm{0}\right)={\displaystyle \frac{\mathrm{1}-n_x^{\mathrm{2}}}{n_{\mathrm{0}z}}}{H}_y\left(z=\mathrm{0}\right)+{\displaystyle \frac{n_x{n}_y}{n_{\mathrm{0}z}}}{H}_x\left(z=\mathrm{0}\right),\end{array}\end{array}$$

where $n_{\mathrm{0}z}=\sqrt{\mathrm{1}-n_{\perp }^{\mathrm{2}}}$, $E_{x^{\prime },y^{\prime }}$ are electric strength components of the incident wave in the coordinate system with the z^′ axis along the ambient magnetic field. The wave polarization near the ground surface is determined by

$$\begin{array}{}\text{(17)}& \mathrm{\Pi }=E_{y^{\prime }}\left(z=\mathrm{0}\right)/E_{x^{\prime }}\left(z=\mathrm{0}\right).\end{array}$$

The electric field at the altitude $z=z_{max}$ can be expressed as

$$\begin{array}{}\text{(18)}& \begin{array}{rl}& E\left(z=z_{max}\right)=\\ & =\sqrt{{\left|E_x\left(z=z_{max}\right)\right|}^{\mathrm{2}}+{\left|E_y\left(z=z_{max}\right)\right|}^{\mathrm{2}}+{\left|E_z\left(z=z_{max}\right)\right|}^{\mathrm{2}}}.\end{array}\end{array}$$

The coordinate dependence of the wave field can be found from the inverse Fourier transform Eq. (10). The vertical energy flux (Poynting vector) is

$$\begin{array}{}\text{(19)}& S_z=(\mathrm{2}{Z}_{\mathrm{0}}{)}^{-\mathrm{1}}\mathit{\text{Re}}{\left[{\mathit{E}}_{\perp }^{*},{\mathit{H}}_{\perp }\right]}_z,\end{array}$$

and the total power of the source is

$$\begin{array}{}\text{(20)}& P={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\text{Re}}\int \left(j_x{E}_x^{*}\left(z=z_s\right)+j_y{E}_y^{*}\left(z=z_s\right)\right)\mathrm{d}x\mathrm{d}y.\end{array}$$

Now we present the results of the numerically computed solution of the set Eqs. (5) and (6) under model conditions for the nighttime ionosphere.

4 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 between charged and neutral particles shown in Fig. 1b. The plasma density data are taken from International Reference Ionosphere (Bilitza and Reinisch, 2007), available from CCMC at https://ccmc.gsfc.nasa.gov/modelweb/models/iri2016_vitmo.php (last access: 28 April 2021), and correspond to 68⁰ N, 25⁰ E; 4 September 2019, 00:30 LT. The collision frequency data are taken from the book of Gurevich and Shvarcburg (1973). The angle of magnetic field inclination with respect to the axis z is equal to ϑ=168⁰. We use the value $z_{max}=\mathrm{125}\mathrm{km}$ as the upper boundary of the solution. At this altitude a typical spatial scale of plasma inhomogeneity exceeds 70 km, and it is much more than the wavelength, which in the considered case is on order 20 to 25 km. As an example, we calculate the fields created by varying at frequency 3 kHz and flowing at the altitude z_s=80 km horizontal current.

Figure 1 Nighttime ionosphere model: (a) the electron plasma density and (b) the collision frequency between electron and neutral particles (solid line) and collision frequency between ion and neutral particles (dashed line).

We assume that currents occupy a volume which has a shape of a horizontal pancake and use for calculations a Gaussian distribution over x and y coordinates

$$\begin{array}{}\text{(21)}& J\left({\mathit{r}}_{\perp }\right)=J_{max}\exp \left(-x^{\mathrm{2}}/\mathrm{2}{L}_x^{\mathrm{2}}-y^{\mathrm{2}}/\mathrm{2}{L}_y^{\mathrm{2}}\right).\end{array}$$

The corresponding current distribution Eq. (2) in Fourier space is also Gaussian and has a form J(n_⊥) = $J_{\mathrm{0}}\exp \left(-k_{\mathrm{0}}^{\mathrm{2}}{L}_x^{\mathrm{2}}{n}_x^{\mathrm{2}}/\mathrm{2}-k_{\mathrm{0}}^{\mathrm{2}}{L}_y^{\mathrm{2}}{n}_y^{\mathrm{2}}/\mathrm{2}\right)$, where $J_{\mathrm{0}}=\mathrm{2}\mathit{\pi }{k}_{\mathrm{0}}^{\mathrm{2}}{L}_x{L}_y{J}_{max}$. We calculate the wave field in N=1000 points with steps and then use inverse fast Fourier transform (Cooley and Tukey, 1965) to find its coordinate dependence. We present the results of field calculation in Fourier space (Figs. 2–3) and in coordinate space (Fig. 4). The dependences of amplitude of horizontal magnetic field $H_{\perp }(n_x,n_y)/E_{\mathrm{0}}$ on ground surface z=0 and amplitude of electric field $E(n_x,n_y)/E_{\mathrm{0}}$ at altitude $z=z_{max}$ are presented in Fig. 2a, b. The field values are normalized by the value E₀=Z₀J₀. Polarization ellipse parameters ϕ, $E_{y^{\prime }}=E_{x^{\prime }}{e}^{-i\mathit{\varphi }}$, and $\log \left|E_{y^{\prime }}/E_{x^{\prime }}\right|$ are presented in Fig. 2c, d. Positive values of phase ϕ correspond to right-hand polarization, typical for whistler waves, and negative values of phase ϕ correspond to left-hand polarization. Positive values of parameter $\log \left|E_{y^{\prime }}/E_{x^{\prime }}\right|$ mean that the polarization ellipse is elongated predominantly along the y axis, and negative values of parameter $\log \left|E_{y^{\prime }}/E_{x^{\prime }}\right|$ mean that the polarization ellipse is elongated predominantly along the x axis. Examples of altitude dependences of normalized electric and magnetic fields corresponding to different horizontal refractive indices are presented in Fig. 3. The level z=z_s of the source action is marked by a dotted line. Figure 4 shows contour maps of fields created by horizontal currents (20) with characteristic horizontal sizes L_x≃12 km, L_y≃70 km and equal x and y components J_x=J_y current density (Fig. 4a), electric field at altitude $z=z_{max}$ (Fig. 4b), horizontal magnetic field on ground surface z=0 (Fig. 4c) and polarization parameter ϕ (Fig. 4d). Current density is normalized by the value J_max, both electric and magnetic fields are normalized by the value Z₀J_max, and the arrow shows the current direction. Examples of ground-based horizontal magnetic field and electric field at altitude $z=z_{max}$ corresponding to the coordinate x=0, characteristic horizontal sizes of source L_y≃70 km, L_y≃12 km and different distributions of J_x,y are presented in Fig. 5.

Figure 2 Fields in the z=0 plane: (a) characteristic horizontal sizes of the horizontal magnetic field, (b) the electric field at the height z=125 km, (c) the polarization ellipse parameters $\mathit{\varphi }\left(E_{y^{\prime }}=E_{x^{\prime }}{e}^{i\mathit{\varphi }}\right)$, and (d) $\log \left|E_{y^{\prime }}/E_{x^{\prime }}\right|$.

Figure 3 Altitude dependences of the amplitudes of horizontal electric and magnetic fields.

5 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 the horizontal plane and to have in this plane an arbitrary finite-space distribution. At first, we consider a plane wave with the horizontal components of the refractive index n_⊥ generated by the current $J\left({\mathit{r}}_{\perp }\right)\sim e^{i{k}_{\mathrm{0}}{\mathit{n}}_{\perp }{\mathit{r}}_{\perp }}$. The set of wave equations in the layer $\mathrm{0}<z<z_{max}$ for each n_⊥ component is completed by boundary conditions assuming a perfect conductivity of the ground surface and excluding wave energy coming on the upper boundary $z=z_{max}$ from above. The mathematical 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. When the dependencies of source current and wave field in n_⟂ space are finite functions with discretized values, the fast Fourier transform algorithm can be used for numerical calculations. Inverse fast Fourier transform yields space dependence of the wave field.

Figure 4 Source currents and field space distributions: (a) space distributions of source currents, (b) the electric field at z=125 km, (c) the horizontal magnetic field at z=0, and (d) the polarization ellipse parameters ϕ, $E_{y^{\prime }}=E_{x^{\prime }}{e}^{-i\mathit{\varphi }}$ at z=0.

Figure 5 Current density (red line) and field space distributions at x=0 for two different ratios $J_y/J_x$, outlined by green ($J_y/J_x=\mathrm{1}$) and blue ($J_y/J_x\simeq -i$) colors. Left plot column (a) corresponds to the $L_y=\mathrm{70}\mathrm{km}\ge \mathrm{1}/k_{\mathrm{0}}$ source scale and right plot column (b) corresponds to the $L_y=\mathrm{12}\mathrm{km}\le \mathrm{1}/k_{\mathrm{0}}$ source scale.

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 a Gaussian distribution of source current density over x and y coordinates. We mention that the model of a plane source current can also be effective in a more general case of a current layer with small thickness $\mathrm{\Delta }z\ll {\mathit{\lambda }}_z\sim \mathrm{60}\mathrm{km}$. The ground-based horizontal magnetic field and the electric field at 125 km are calculated in both Fourier (n_x, n_y) and coordinate (x, y) spaces, since a wave can achieve the ground surface in penetrating mode in case $n_{\perp }\le \mathrm{1}$. The magnetic field $H_{\perp }({\mathit{n}}_{\perp },z=\mathrm{0})$ is noticeably non-zero for Fourier components with $n_{\perp }<\mathrm{1}$ and is practically equal to zero for Fourier components with ${\mathit{n}}_{\perp }\gg \mathrm{1}$. Waves with $n_{\perp }<\mathrm{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_{\perp }({\mathit{n}}_{\perp },z=\mathrm{0})$ can increase by 2 or 3 times (see Fig. 2a). The polarization parameter ϕ of such components can be negative, similarly to the left-hand-polarized waves (see Fig. 2c). A typical size of the spot and a character of polarization of the ground-based magnetic field $H_{\perp }\left(r_{\perp },z=\mathrm{0}\right)$ depend on the source distribution. If the horizontal size of the radiating currents exceeds a value $L\ge \mathrm{1}/k_{\mathrm{0}}\sim \mathrm{15}\mathrm{km}$ (for the wave of frequency 3 kHz), then the Fourier components with $n_{\perp }\le \mathrm{1}$ dominate in the total field. In that case, the magnetic field $H_{\perp }\left(r_{\perp },z=\mathrm{0}\right)$ is predominantly localized under the source. If the horizontal size of radiating currents is small enough $L\le \mathrm{1}/k_{\mathrm{0}}$, then Fourier components $n_{\perp }\sim \mathrm{1}$ can make a notable contribution to the total field. The magnetic field $H_{\perp }\left(r_{\perp },z=\mathrm{0}\right)$ occupies a spot $\sim \mathrm{1}/k_{\mathrm{0}}$ but can sometimes be registered far from the source (see Figs. 4c, 5b). The field above the source is determined by both direct radiation from the source and the reflected radiation. Because of the effective night reflection of ELF/VLF waves, a spot of the electric field at 125 km can exceed the size of the source by several times (Fig. 4b).

At altitudes of source and above, the two ELF/VLF wave modes Eq. (8) are weakly damped and have right-hand polarization, and another two modes are evanescent and have left-hand polarization. At altitudes ∼60–70 km and below all four modes transform into vacuum electromagnetic ones. Source currents can excite both right-hand- and left-hand-polarized modes. Hence, current distribution properties specify field distribution, field polarization, and proportion in which source energy supplies the Earth–ionosphere waveguide or flows upward. As an example, Figs. 5a, b show coordinate dependencies of source density, fields $E\left(x=\mathrm{0},y,z=\mathrm{0}\right)$, $H_{\perp }\left(x=\mathrm{0},y,z=\mathrm{0}\right)$ and polarization parameter ϕ in cases of large $L_y=\mathrm{70}\mathrm{km}\ge \mathrm{1}/k_{\mathrm{0}}$ (a) and small $L_y=\mathrm{12}\mathrm{km}\le \mathrm{1}/k_{\mathrm{0}}$ (b) sources for two different ratios $J_y/J_x$ in distribution (20), outlined by green and blue colors. In the first case of ratio $J_y/J_x=\mathrm{1}$ (blue lines on the plots), source currents excite predominantly right-hand-polarized modes. So the polarization parameter at z=0 is mainly positive under the source. Approximately 50 %–60 % of the source energy is carried upward, 30 %–40 % is adsorbed, and about 10 % is supplied to the Earth–ionosphere waveguide. This can be useful for modification processes in the plasma magnetic trap (Bespalov and Trakhtengerts, 1986). Of course, the electromagnetic ELF/VLF radiation of ionospheric currents themselves is not sufficient for a noticeable modification of the Earth's electron radiation belts. However, under quiet conditions in the nighttime magnetosphere, these emissions can ensure the radiation belt transition through the threshold of the cyclotron instability. This process can be accompanied by a significant precipitation of energetic electrons into the ionosphere and other geophysical manifestations.

In the second case of ratio $J_y/J_x\simeq -i$ ratio (green lines on the plots), source currents excite predominantly left-hand-polarized modes. Respectively, the polarization parameter at z=0 becomes predominantly negative. In that case, regardless of source size, approximately 10 %–12 % of the energy is supplied to the Earth–ionosphere waveguide, and about 90 % is absorbed and not carrying upward to the magnetosphere. The proportion in which source energy supplies the Earth–ionosphere waveguide or flows upward mainly depends on polarization of the source current.

If the horizontal size of radiating currents is small enough, $L\le \mathrm{1}/k_{\mathrm{0}}$, a ground-based magnetic field can be noticeably non-zero far (a few thousand kilometers) from the source (see Fig. 5b). Propagation of modulated ELF/VLF signals in the Earth–ionosphere waveguide far from the source and also into space (Inan et al., 2004) is observed experimentally. The HAARP heating facility in Gakona, Alaska, injects ELF/VLF waves in the Earth–ionosphere waveguide as far as 4400 km (Moore et al., 2007).

The current distributions used as an example in our calculation and presented in Fig. 4 can be similar to electrojet currents modulated in the D region by the HAARP HF heating facility (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 the ambient electrojet field is directed along the x axis, then j_x≈j_y. The maximum surface density of modulated currents Eq. (1) has an order $J_{max}\sim {\mathrm{10}}^{-\mathrm{6}}$–${\mathrm{10}}^{-\mathrm{5}}\mathrm{A}{\mathrm{m}}^{-\mathrm{1}}$. Using in our calculations the magnitude of current density $J_{max}\sim \mathrm{5}\cdot {\mathrm{10}}^{-\mathrm{6}}\mathrm{A}{\mathrm{m}}^{-\mathrm{1}}$ yields the total power of the source ∼36 W, the ground-based horizontal magnetic field under the source $B_{\perp }\sim \mathrm{1}\mathrm{pT}$, and the electric field at the altitude of 125 km above the source $E\sim \mathrm{400}\mathrm{\mu }\mathrm{V}{\mathrm{m}}^{-\mathrm{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 of 125 km is $\sim \mathrm{3.2}\mathrm{nW}{\mathrm{m}}^{-\mathrm{2}}$, and total power is ∼17 W. About half of the source energy is carried upward, approximately 20 % of the energy is supplied to the Earth–ionosphere waveguide, and approximately 30 % of the energy is absorbed.

6 Conclusions

We find a field of monochromatic whistler waves which are excited and propagating in the low nighttime ionosphere. Using a MATLAB boundary-value problem solver enables us to find numerically stable solutions of a full set of the wave equations applying to conditions of an inhomogeneous ionosphere at altitudes below 125 km. Above this altitude the ionosphere plasma is slightly inhomogeneous, and hence approximate methods are suitable. As an example, this calculation technique is applied to the problem of ELF/VLF wave radiation from modulated HF-heated electrojet currents. At first we consider a plane wave with a 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 an inverse fast Fourier transform to find a total field, get the dependencies of a wave field at 0 and 125 km, and analyze the type of wave polarization on the ground surface. The proportion in which source energy supplies the Earth–ionosphere waveguide, absorbs or flows upward depends on the altitude profile of the ionosphere plasma. Besides, the spatial distribution of radiating energy can be regulated by modulated currents. Depending on their properties, radiating energy can predominantly flow into space or inject into the Earth–ionosphere waveguide far from the source. The obtained values are in good agreement with ground and satellite observations and known calculation results. Using a model of the plane horizontal source, currents can be generalized for the arbitrary altitude source distribution.