Numerical Simulation of Very-Low-Frequency Waves Passing through the Magnetoactive Plane-Layered Plasma of Earth’s Lower Ionosphere

Abstract—The article describes an algorithm for calculating the characteristics of the passage of low-frequency electromagnetic waves through the magnetoactive plane-stratified lower ionosphere. The algorithm is based on numerical integration of wave equations for the case of a plane wave incident on the ionosphere from below, taking into account the plasma’s ionic composition. As an illustration, the basic characteristics of passage of waves in the 0.05–1 kHz frequency range are calculated. It is shown that the power transmission coefficient of these waves is quite large, ranging from 0.2–0.4 for most of the frequencies used. It is shown that in the daytime near a frequency of ~300 Hz, the transmission coefficient in the daytime has a maximum and is approximately double the nighttime values.

1 INTRODUCTION

Numerical simulation of remote recording of low-frequency waves in the magnetosphere and on Earth’s surface from artificial disturbances in the ionosphere is of great importance (see, e.g., [1–3]). This is due to the commissioning of a number of powerful ground-based high-frequency (HF) transmitters, the radiation of which is modulated at a low frequency (LF). The radiation of such transmitters (“heating stands”) is intensively absorbed in the lower ionosphere, creating a nonlinear demodulation effect and the emission of LF waves into the ionosphere and the Earth–ionosphere waveguide [4]. More than ten heating stands are functioning around the world; the EISCAT (Norway) and HAARP (Alaska) stands are the most powerful and are equipped with scientific equipment [5, 6].

To date, it is possible to consider that the theory of very- and extremely-low-frequency (VLF and ELF) wave propagation in the near-Earth space, in particular, in the Earth–ionosphere waveguide (resonator), has basically been constructed [7–10]. The corresponding algorithms are often insufficiently illustrative physically and require considerable effort to interpret the results. For this reason, the problem can be divided into related parts, for each of which relatively simple numerical algorithms are constructed. This primarily relates to the algorithms for constructing ray trajectories in the ionosphere and magnetosphere [11–14]. Even in the Earth–ionosphere waveguide, where at low frequencies the geometric optics approximation is formally not feasible near the upper “wall,” it is possible to construct effective numerical algorithms using the “jump method” [15–17]. In this case, the coefficients of reflection from the lower inhomogeneous anisotropic ionosphere must be calculated independently.

The calculation of such coefficients is part of the general solution to the problem of electromagnetic waves passing through a locally plane-layered anisotropic plasma of the lower ionosphere.

The objective of this work is to describe one of the possible algorithms for solving such a problem. By way of illustration, the calculations were performed for VLF (50–1000 Hz) case, which has recently experienced increased attention in diagnostics of the state of Earth’s crust using remote sensing methods, as well as in communication problems.

2 BASIC RELATIONS

Mathematically, the problem of electromagnetic (EM) waves passing through a locally plane-layered anisotropic plasma (the ionosphere) was formulated a relatively long time ago [18]. However, in attempts to numerically calculate the EM field strength along the path and the corresponding reflection and transmission coefficients, some fundamental difficulties arose associated with integrating a system of homogeneous differential equations with exponential numerical instability of the solution [18, 19]. Here, we briefly describe an efficient algorithm for solving this problem, taking into account the ionic composition of the plasma.

To calculate the field strength, we use the Maxwell equations together with a model of the medium in the form of an inhomogeneous cold multicomponent magnetoactive plasma. Such a medium is described by the dielectric constant tensor ${\hat {\varepsilon }}$ (the magnetic permeability of the ionospheric plasma is equal to the magnetic permeability of a vacuum μ₀) The tensor ${\hat {\varepsilon }}$ in the coordinate system in which the z axis is directed along the external geomagnetic field has the form [20]

$$\hat {\varepsilon } = \left[ {\begin{array}{*{20}{c}} S&{ - iD}&0 \\ {iD}&S&0 \\ 0&0&P \end{array}} \right],$$

(1а)

where

$$\begin{gathered} R = 1 + \sum {{{{{X}_{k}}} \mathord{\left/ {\vphantom {{{{X}_{k}}} {({{Y}_{k}} - {{U}_{k}})}}} \right. \kern-0em} {({{Y}_{k}} - {{U}_{k}})}}} ; \\ L = 1 - \sum {{{{{X}_{k}}} \mathord{\left/ {\vphantom {{{{X}_{k}}} {({{Y}_{k}} + {{U}_{k}})}}} \right. \kern-0em} {({{Y}_{k}} + {{U}_{k}})}}} ; \\ P = 1 - \sum {{{{{X}_{k}}} \mathord{\left/ {\vphantom {{{{X}_{k}}} {{{U}_{k}}}}} \right. \kern-0em} {{{U}_{k}}}}} ; \\ \end{gathered} $$

(1b)

$$\begin{gathered} S = {{\left( {R + L} \right)} \mathord{\left/ {\vphantom {{\left( {R + L} \right)} 2}} \right. \kern-0em} 2};\,\,\,\,D = {{\left( {R--L} \right)} \mathord{\left/ {\vphantom {{\left( {R--L} \right)} 2}} \right. \kern-0em} 2}; \\ {{U}_{k}} = 1--{{i{{\nu }_{k}}} \mathord{\left/ {\vphantom {{i{{\nu }_{k}}} \omega }} \right. \kern-0em} \omega }, \\ \end{gathered} $$

(1c)

ν_k is the effective collision frequency of k particles, f is the wave frequency, ω = 2πf, and i is the imaginary unit. Summation is carried out for k charged particles making up the plasma, taking into account the sign of the charge in terms of Y_k:

$${{Y}_{k}} \equiv {{{{f}_{{{\text{H}}k}}}} \mathord{\left/ {\vphantom {{{{f}_{{{\text{H}}k}}}} f}} \right. \kern-0em} f};\,\,\,\,{{X}_{k}} \equiv {{({{{{f}_{{{\text{p}}k}}}} \mathord{\left/ {\vphantom {{{{f}_{{{\text{p}}k}}}} f}} \right. \kern-0em} f})}^{2}},$$

(2)

where f_Hk and f_pk are the gyro and plasma frequency of a k particle: k = 1 (electrons), 2, … (ions). For electrons, we have

$${{f}_{{{\text{pe}}}}} = {{({{{{e}^{2}}{{N}_{{\text{e}}}}} \mathord{\left/ {\vphantom {{{{e}^{2}}{{N}_{{\text{e}}}}} {4{{\pi }^{2}}{{\varepsilon }_{0}}{{m}_{{\text{e}}}}}}} \right. \kern-0em} {4{{\pi }^{2}}{{\varepsilon }_{0}}{{m}_{{\text{e}}}}}})}^{{{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0em} 2}}}} \approx 8.97N_{{\text{e}}}^{{{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0em} 2}}}\,\,{\text{kHz}},$$

(3)

where N_e is measured in cm^–3; e, m_e are the electron charge and mass; ε₀ is the dielectric constant of the vacuum. At the F2 layer maximum, we have in the daytime f_pe ≈ 10 MHz; in the nighttime, f_pe ≈ 2–5 MHz, depending on the conditions (time of day, season, geographical location, and solar activity):

$${{f}_{{{\text{He}}}}} = {{e{{B}_{0}}} \mathord{\left/ {\vphantom {{e{{B}_{0}}} {2\pi {{m}_{{\text{e}}}}}}} \right. \kern-0em} {2\pi {{m}_{{\text{e}}}}}}.$$

(4)

Within the ionosphere, the spatial structure of the geomagnetic field vector $\overrightarrow {{{B}_{0}}}$ is well described by the model of a point magnetic dipole located near the center of the Earth with the axis inclined at an angle to Earth’s axis of rotation [21]. In analogy with geographical coordinates, the concepts of geomagnetic latitude Φ and longitude Λ are introduced. At height h and latitude Φ, the gyrofrequency can be estimated as follows:

$${{f}_{{{\text{He}}}}} \approx 876.0{{(1 + {h \mathord{\left/ {\vphantom {h {{{R}_{0}}}}} \right. \kern-0em} {{{R}_{0}}}})}^{{ - 3}}}{{(1 + 3{{\sin }^{2}}\Phi )}^{{{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0em} 2}}}}\,\,{\text{kHz}},$$

(5)

where R₀ is the average radius of the Earth (~6370 km). For k ions, in formulas (3), (4), the appropriate corresponding values of the concentration N_k and mass m_k should be used.

Let the origin of the right Cartesian rectangular coordinate system be located at O on Earth’s surface and the z axis be directed vertically upward. A plane EM wave with a wave vector $\vec {k}$ is incident from below onto the horizontal nonuniform (along the z axis) layer of the plasma at angle θ with the vertical. The upper and lower boundaries of the layer correspond to heights z₁ and z₂. We direct the x axis so that the vector $\vec {k}$ is in the plane XZ. The direction of the local geomagnetic field is determined by angle ψ of the vector ($- \overrightarrow {{{B}_{0}}}$) with the vertical and the given azimuth angle A, which is counted from the x axis in the direction of the projection $\overrightarrow {{{B}_{0}}}$ onto the Earth’s surface. Counting goes counterclockwise looking at the surface from top to bottom. The angle ψ is determined from the relation cot ψ = 2tan Φ [21]. To describe the direction $\overrightarrow {{{B}_{0}}}$, we use its directing cosines {l, m, n}, where

$$\begin{gathered} l = - \sin \psi \cos A,\,\,\,\,m = - \sin \psi \sin A, \\ n = - \cos \psi . \\ \end{gathered} $$

This choice of reference system makes it possible to record any variable wave quantity $\vec {V}\left( {x,y,z,t} \right)$ in the form

$$\vec {V}\left( {x,y,z,t} \right) = \vec {V}\left( z \right)\exp (i\omega t--i{{k}_{0}}\Gamma z),$$

where k₀ = ω/c is the wavenumber in vacuum; Γ = sinθ. Omitting everywhere this exponential factor, we can transform the system of partial differential equations into a system of ordinary first-order differential equations [18]:

$${{{\text{d}}\vec {e}(z)} \mathord{\left/ {\vphantom {{{\text{d}}\vec {e}(z)} {{\text{d}}z}}} \right. \kern-0em} {{\text{d}}z}} = - i{{k}_{0}}\hat {T}\vec {e}(z),$$

(6)

where $\vec {e}$ is a column vector of the form {E_x, –Ey, Z₀H_x, Z₀H_y} artificially created from the electrical field components E_x, E_y and the magnetic field components H_x, H_y of the wave; here, the magnetic field components for equalizing the dimensions are multiplied by the impedance of the vacuum Z₀ = (µ₀/ ε₀)^1/2 = 120π. The components of a 4 × 4 matrix ${\hat {T}}$ are expressed via the components of a 3 × 3 electrical susceptibility matrix ${\hat {M}} = ~\,\,{\hat {\varepsilon }} - ~\,\,\hat {1}$ ($\hat {1}~$ is an identity matrix):

$$\begin{gathered} {{T}_{{11}}} = {{ - \Gamma {{M}_{{31}}}} \mathord{\left/ {\vphantom {{ - \Gamma {{M}_{{31}}}} {\left( {1 + {{M}_{{33}}}} \right)}}} \right. \kern-0em} {\left( {1 + {{M}_{{33}}}} \right)}};\,\,\,\,{{T}_{{12}}} = {{\Gamma {{M}_{{32}}}} \mathord{\left/ {\vphantom {{\Gamma {{M}_{{32}}}} {\left( {1 + {{M}_{{33}}}} \right)}}} \right. \kern-0em} {\left( {1 + {{M}_{{33}}}} \right)}}; \\ {{T}_{{14}}} = {{\left( {{{C}^{2}} + {{M}_{{33}}}} \right)} \mathord{\left/ {\vphantom {{\left( {{{C}^{2}} + {{M}_{{33}}}} \right)} {\left( {1 + {{M}_{{33}}}} \right)}}} \right. \kern-0em} {\left( {1 + {{M}_{{33}}}} \right)}}; \\ {{T}_{{31}}} = --{{M}_{{21}}} + {{{{M}_{{23}}}{{M}_{{31}}}} \mathord{\left/ {\vphantom {{{{M}_{{23}}}{{M}_{{31}}}} {\left( {1 + {{M}_{{33}}}} \right)}}} \right. \kern-0em} {\left( {1 + {{M}_{{33}}}} \right)}}; \\ {{T}_{{32}}} = {{C}^{2}} + {{M}_{{22}}} - {{{{M}_{{23}}}{{M}_{{32}}}} \mathord{\left/ {\vphantom {{{{M}_{{23}}}{{M}_{{32}}}} {\left( {1 + {{M}_{{33}}}} \right)}}} \right. \kern-0em} {\left( {1 + {{M}_{{33}}}} \right)}}; \\ {{T}_{{34}}} = {{\Gamma {{M}_{{23}}}} \mathord{\left/ {\vphantom {{\Gamma {{M}_{{23}}}} {\left( {1 + {{M}_{{33}}}} \right)}}} \right. \kern-0em} {\left( {1 + {{M}_{{33}}}} \right)}}; \\ {{T}_{{44}}} = {{ - \Gamma {{M}_{{13}}}} \mathord{\left/ {\vphantom {{ - \Gamma {{M}_{{13}}}} {\left( {1 + {{M}_{{33}}}} \right)}}} \right. \kern-0em} {\left( {1 + {{M}_{{33}}}} \right)}}; \\ {{T}_{{41}}} = 1 + {{M}_{{11}}} - {{{{M}_{{13}}}{{M}_{{31}}}} \mathord{\left/ {\vphantom {{{{M}_{{13}}}{{M}_{{31}}}} {\left( {1 + {{M}_{{33}}}} \right)}}} \right. \kern-0em} {\left( {1 + {{M}_{{33}}}} \right)}}; \\ {{T}_{{42}}} = - {{M}_{{12}}} + {{{{M}_{{13}}}{{M}_{{32}}}} \mathord{\left/ {\vphantom {{{{M}_{{13}}}{{M}_{{32}}}} {\left( {1 + {{M}_{{33}}}} \right)}}} \right. \kern-0em} {\left( {1 + {{M}_{{33}}}} \right)}}; \\ {{T}_{{13}}} = {{T}_{{21}}} = {{T}_{{22}}} = {{T}_{{24}}} = {{T}_{{33}}} = {{T}_{{43}}} = {\text{ }}0;\,\,\,\,{{T}_{{23}}} = {\text{ }}1, \\ \end{gathered} $$

(7)

where C = cos θ.

In the chosen coordinate system, it is easy to show that the elements of the matrix $\hat {M}$ have the form

$$\begin{gathered} {{M}_{{11}}} = \left( {1 - {{l}^{2}}} \right)S + {{l}^{2}}P - 1; \\ {{M}_{{12}}} = lm\left( {P - S} \right) + inD; \\ {{M}_{{13}}} = ln\left( {P - S} \right) - imD; \\ {{M}_{{21}}} = lm\left( {P - S} \right) - inD; \\ {{M}_{{22}}} = \left( {1 - {{m}^{2}}} \right)S + {{m}^{2}}P - 1; \\ {{M}_{{23}}} = mn\left( {P - S} \right) + ilD; \\ {{M}_{{31}}} = ln\left( {P - S} \right) + imD; \\ {{M}_{{32}}} = mn\left( {P - S} \right) - ilD; \\ {{M}_{{33}}} = \left( {1 - {{n}^{2}}} \right)S + {{n}^{2}}P - 1; \\ \end{gathered} $$

(8)

where the parameters S, P, D are determined in Eqs. (1) taking into account the ionic composition of the ionospheric layer and the effective frequency of particle collisions in the layer.

The missing components E_z and H_z in Eq. (6) are found explicitly from the Maxwell equations:

$${{E}_{z}} = - {{\left( {\Gamma {{Z}_{0}}{{H}_{y}} + {{M}_{{31}}}{{E}_{x}} + {{M}_{{32}}}{{E}_{y}}} \right)} \mathord{\left/ {\vphantom {{\left( {\Gamma {{Z}_{0}}{{H}_{y}} + {{M}_{{31}}}{{E}_{x}} + {{M}_{{32}}}{{E}_{y}}} \right)} {\left( {1 + {{M}_{{33}}}} \right)}}} \right. \kern-0em} {\left( {1 + {{M}_{{33}}}} \right)}},$$

(9a)

$${{H}_{z}} = {{\Gamma {{E}_{y}}} \mathord{\left/ {\vphantom {{\Gamma {{E}_{y}}} {{{Z}_{0}}}}} \right. \kern-0em} {{{Z}_{0}}}}.$$

(9b)

To isolate the incident and reflected waves in the waveguide at z < z₂, we write the relations between the components of the fields in this range:

$$\begin{gathered} {{H}_{x}} = \mp C{{E}_{y}};\,\,\,\,{{H}_{y}} = {{ \pm {{E}_{x}}} \mathord{\left/ {\vphantom {{ \pm {{E}_{x}}} {(C{{Z}_{0}})}}} \right. \kern-0em} {(C{{Z}_{0}})}}; \\ {{H}_{z}} = {{ \mp \Gamma {{H}_{x}}} \mathord{\left/ {\vphantom {{ \mp \Gamma {{H}_{x}}} C}} \right. \kern-0em} C} = {{\Gamma {{E}_{y}}} \mathord{\left/ {\vphantom {{\Gamma {{E}_{y}}} {{{Z}_{0}}}}} \right. \kern-0em} {{{Z}_{0}}}};\,\,\,\,{{E}_{z}} = \mp C{{E}_{y}}; \\ \end{gathered} $$

(10)

where the upper sign corresponds to the ascending wave, and the lower one, to the descending.

System (6) consists of four ordinary first-order linear differential equations. The system is homogeneous, i.e., having four linearly independent solutions ${{\vec {e}}_{j}}$ (j = 1–4); any other solution (6) is a linear combination of these solutions. It is easy to show that solutions ${{\vec {e}}_{j}}$ in a homogeneous medium correspond to the classical magnetoionic modes, conventionally called ordinary and extraordinary waves, which propagate upward and downward [18].

Let us assume that for z > z₁, the parameters of the ionosphere are changed little by the wavelength. In this case, when a wave is incident on the layer from below, only ordinary and extraordinary waves are directed upward at the upper boundary of the layer. Let these be waves j = 1, 2. For z = z₁, solution (6) can be represented as $\vec {e}$ = a₁${{\vec {e}}_{1}}$ + a₂${{\vec {e}}_{2}}$, where a_1,2 are complex coefficients. Due to the linearity and homogeneity of system (6), such a representation of the solution is valid for any height z < z₁; i.e., the coefficients a_1,2 are independent of z. Therefore, it is necessary to find only two linearly independent solutions (6) corresponding to the ascending ordinary and extraordinary waves. To do this, we calculate the initial values ${{\vec {e}}_{{1,2}}}$(z₁) in a homogeneous medium above the inhomogeneous layer. The existence of such solutions requires fulfillment of the well-known condition det(${\hat {T}}$ – q$\hat {1}$) = 0, where the coefficient q describes the change in a plane wave along the z axis as exp(–ik₀qz). Knowing the values q_1,2 at height z₁, we can integrate (6) downward for the corresponding eigenvectors ${{\vec {e}}_{{1,2}}}$(z) of the matrix ${\hat {T}}$.

Opening up the determinant, taking into account the values of the elements of the matrix ${\hat {T}}$ (7), we obtain the well-known complex fourth-degree algebraic equation with respect to q (the “Booker quartic”) [18]. Figure 1 schematically shows the location of the roots of the quartic q_j (j = 1–4) on the complex plane in the LF case, when the difference of the roots for the four magnetoionic modes is the most illustrative. For ascending waves, Re(q)> 0. It can be seen from the figure that for the ordinary wave, |Im(q₁)| $\gg$ Re(q₁), and this wave attenuates relatively quickly in an inhomogeneous medium. For the extraordinary wave, conversely, |Im(q₂)| $\ll$ Re(q₂), attenuation is relatively small. For this reason, whistlers extend over great distances from hemisphere to hemisphere in the ionosphere and magnetosphere. Therefore, this LF mode is often called whistling.

Fig. 1.

Diagram of roots q_j ( j = 1–4) on complex plane.

Determining the roots of the quartic seems a simple task, since it is well known that a fourth-degree algebraic equation with one unknown has an analytical solution. However, the corresponding formulas contain embedded square and cubic roots, as well as differences of quantities that are close in value. In addition, the imaginary and real parts of the roots can vary greatly in absolute value, while the roots themselves can approach the coordinate axes at distances comparable to the calculation accuracy, which in practice can lead to gross computational errors. We propose an efficient numerical algorithm for calculating the roots of the quartic, based on a simplified version of the simplex method well known in linear programming.

(1) On the surface q in the complex space, we construct a regular triangle near the origin with a characteristic size much smaller than the refractive index of the wave ${{\vec {e}}_{2}}$ at height z₁. As the characteristic size, we choose, e.g., the radius of a circumscribed circle.

(2) We check the q values at the vertices of the triangle and find the smallest.

(3) We move the center of the triangle to this vertex and reduce its characteristic size.

With cyclic repetition of procedure 2, 3, the triangle inevitably “hits” one of the roots, stopping when the required accuracy is achieved.

(4) The degree of the polynomial decreases by one, and procedures 1–4 are repeated until all four roots are calculated.

We have long used this simple algorithm in this problem and have never produced an erroneous result.

For numerical integration (6), we use the stable numerical Runge–Kutta–Merson integration method with a variable step and control of the integration error [22], although other methods [23] can be effectively applied.

Then we find the initial values ${{\vec {e}}_{{1,2}}}$(z₁) for integrating system of equations (6) as the eigenvectors of the matrix ${\hat {T}}$ at height z₁ with the known eigenvalues q_j:

$$(\hat {T} - {{q}_{j}}\hat {1}){{\vec {e}}_{j}} = 0\,\,{\text{at}}\,\,j = 1,2.$$

(11)

Linear system of equations (11) is easily solved up to a constant factor, e.g., E_y:

$${{E}_{x}} = - \alpha {{E}_{y}};\,\,\,\,{{Z}_{0}}{{H}_{x}} = - q{{E}_{y}};\,\,\,\,~{{Z}_{0}}{{H}_{y}} = - \beta {{E}_{y}},$$

(12)

where the index j is not given,

$$\begin{gathered} \alpha = {{[{{T}_{{34}}}{{T}_{{42}}} + ({{q}^{2}}--{{T}_{{32}}})({{T}_{{44}}}--q)]} \mathord{\left/ {\vphantom {{[{{T}_{{34}}}{{T}_{{42}}} + ({{q}^{2}}--{{T}_{{32}}})({{T}_{{44}}}--q)]} {[{{T}_{{31}}}({{T}_{{44}}}--q){\text{ }}--{{T}_{{34}}}{{T}_{{41}}}]}}} \right. \kern-0em} {[{{T}_{{31}}}({{T}_{{44}}}--q){\text{ }}--{{T}_{{34}}}{{T}_{{41}}}]}}, \\ \beta = {{ - [{{T}_{{12}}} + \alpha ({{T}_{{11}}}--q)]} \mathord{\left/ {\vphantom {{ - [{{T}_{{12}}} + \alpha ({{T}_{{11}}}--q)]} {{{T}_{{14}}}}}} \right. \kern-0em} {{{T}_{{14}}}}}. \\ \end{gathered} $$

At low frequencies, an ordinary wave strongly attenuates along the z axis; therefore, when integrating from top to bottom, an ordinary wave ${{\vec {e}}_{1}}$(z) rapidly increases in comparison to an extraordinary wave ${{\vec {e}}_{2}}$(z) (Fig. 1). Due to the finite number of places in the mantissa of real numbers in a computer, this effect of a strong increase in one solution can lead to loss of linear independence; therefore, it must be periodically restored during integration of Eqs. (6) taking advantage of the fact that orthogonal vectors are always linearly independent [24]. Upon reaching the amplitude of the solution ${{\vec {e}}_{1}}$(z) of a given threshold, we form a vector ${{\vec {e}}_{{20}}}$ = ${{\vec {e}}_{2}}$ + $\xi {{\vec {e}}_{1}}$ and require that the scalar product ($\vec {e}_{1}^{*}$, ${{\vec {e}}_{{20}}}$) be equal to 0. Hence, we find the factor

$$\xi = {{ - \left( {\vec {e}_{1}^{*},{{{\vec {e}}}_{2}}} \right)} \mathord{\left/ {\vphantom {{ - \left( {\vec {e}_{1}^{*},{{{\vec {e}}}_{2}}} \right)} {\left( {\vec {e}_{1}^{*},{{{\vec {e}}}_{1}}} \right)}}} \right. \kern-0em} {\left( {\vec {e}_{1}^{*},{{{\vec {e}}}_{1}}} \right)}},$$

where the asterisk indicates complex conjugation. Since the vector ${{\vec {e}}_{{20}}}$ is a linear combination of solutions (6), it is also a solution that is linearly independent with ${{\vec {e}}_{1}}$. Note that in a strongly inhomogeneous medium, the concept of magnetoionic mode loses its original physical sense. As a result of the above numerical integration procedure, we obtain for z = z₂ two formally correct linearly independent solutions (6). However, as a result of applying the orthogonalization procedure, these solutions together have no physical sense. Using them, it is possible to construct all the characteristics of the wave field in empty space for z ≤ z₂ [23].

First, it is possible to separate the waves below the layer into ascending and descending parts. From expressions (7) and (8), it is easy to obtain q = ±C; the upper sign corresponds to ascending waves ($\vec {U}$), and the lower sign, to descending ($\vec {D}$), i.e.,

$$\vec {U}\sim \exp ( - i{{k}_{0}}Cz)\,\,{\text{and}}\,\,\vec {D}\sim \exp ( + i{{k}_{0}}Cz).$$

We take an electric field $\vec {E}$, then $\vec {E}$ = $\vec {U}$ + $\vec {D}$ and ${\text{d}}\vec {E}$/dz = –ik₀C$\vec {U}$ + ik₀C$\vec {D}$. The expression for the derivative can also be obtained from Eq. (6). Taking into account relations (10), we obtain

$$\begin{gathered} {{U}_{x}} = {{({{E}_{x}} + C{{Z}_{0}}{{H}_{y}})} \mathord{\left/ {\vphantom {{({{E}_{x}} + C{{Z}_{0}}{{H}_{y}})} 2}} \right. \kern-0em} 2}; \\ {{U}_{y}} = {{({{E}_{y}}--{{{{Z}_{0}}{{H}_{x}}} \mathord{\left/ {\vphantom {{{{Z}_{0}}{{H}_{x}}} C}} \right. \kern-0em} C})} \mathord{\left/ {\vphantom {{({{E}_{y}}--{{{{Z}_{0}}{{H}_{x}}} \mathord{\left/ {\vphantom {{{{Z}_{0}}{{H}_{x}}} C}} \right. \kern-0em} C})} 2}} \right. \kern-0em} 2}; \\ {{U}_{z}} = {{ - \Gamma {{U}_{x}}} \mathord{\left/ {\vphantom {{ - \Gamma {{U}_{x}}} C}} \right. \kern-0em} C}; \\ \end{gathered} $$

$$\begin{gathered} {{D}_{x}} = {{({{E}_{x}}--C{{Z}_{0}}{{H}_{y}})} \mathord{\left/ {\vphantom {{({{E}_{x}}--C{{Z}_{0}}{{H}_{y}})} 2}} \right. \kern-0em} 2}; \\ {{D}_{y}} = {\text{ }}{{({{E}_{y}} + {{{{Z}_{0}}{{H}_{x}}} \mathord{\left/ {\vphantom {{{{Z}_{0}}{{H}_{x}}} C}} \right. \kern-0em} C})} \mathord{\left/ {\vphantom {{({{E}_{y}} + {{{{Z}_{0}}{{H}_{x}}} \mathord{\left/ {\vphantom {{{{Z}_{0}}{{H}_{x}}} C}} \right. \kern-0em} C})} 2}} \right. \kern-0em} 2}; \\ {{D}_{z}} = {{\Gamma {{D}_{x}}} \mathord{\left/ {\vphantom {{\Gamma {{D}_{x}}} C}} \right. \kern-0em} C}, \\ \end{gathered} $$

where C = cos θ.

Second, it can be shown that the linear combination ${{\vec {U}}_{{\text{p}}}}$ = ${{\vec {U}}_{2}}$ + $\zeta {{\vec {U}}_{1}}$, where ζ = –($\vec {U}_{1}^{*}$${{\vec {U}}_{2}}$)/($\vec {U}_{1}^{*}$${{\vec {U}}_{1}}$), is a “penetrating” wave, which generates the least attenuating extraordinary wave above the level z₁. Linearly independent solutions ${{\vec {e}}_{1}}$ and ${{\vec {e}}_{{\text{p}}}}$ to Eq. (6) lower than height z₂ have a simple physical sense: they have, respectively, the smallest and largest coefficients of power transmission through an inhomogeneous plasma layer z₂ – z₁, and ${{\vec {e}}_{{\text{p}}}}$ = ${{\vec {e}}_{2}}$ + $\zeta {{\vec {e}}_{1}}$.

Third, the vectors ${{\vec {e}}_{1}}$ and ${{\vec {e}}_{{\text{p}}}}$ can be used as a basis for obtaining the characteristics of a known plane wave incident on the layer from below (reflection and transmission coefficients, as well as power). Traditionally, incident plane waves with two independent polarizations are considered. Vertical polarization means that the electric field strength vector $\vec {E}$ of the wave lies in the plane (x, z), in which lies the vector $\vec {k}$. In the case of horizontal polarization, the vector $\vec {E}$ is perpendicular to this plane, i.e., parallel to the y axis. We denote these vectors as ${{\vec {E}}_{\parallel }}$ and ${{\vec {E}}_{ \bot }}$, respectively. Since the waves of these conditionally selected polarizations should be linear combinations of the basis solutions (${{\vec {e}}_{1}}$, ${{\vec {e}}_{{\text{p}}}}$), from the simplest geometric considerations, it is easy to obtain the expansion coefficients ${{\vec {E}}_{\parallel }}$ and ${{\vec {E}}_{ \bot }}$ for waves going up (↑) and down (↓), which will obviously be the same for these ranges: $\vec {E}_{\parallel }^{ \uparrow }$ = A_∥${{\vec {U}}_{1}}$ + B_∥${{\vec {U}}_{{\text{p}}}}$,

(13a)

$$\begin{gathered} {{A}_{\parallel }} = {{ - {{E}_{\parallel }}C{{U}_{{{\text{p}}y}}}} \mathord{\left/ {\vphantom {{ - {{E}_{\parallel }}C{{U}_{{{\text{p}}y}}}} {({{U}_{{1y}}}{{U}_{{{\text{p}}x}}} - {{U}_{{1x}}}{{U}_{{{\text{p}}y}}})}}} \right. \kern-0em} {({{U}_{{1y}}}{{U}_{{{\text{p}}x}}} - {{U}_{{1x}}}{{U}_{{{\text{p}}y}}})}}; \\ {{B}_{\parallel }} = {{{{E}_{\parallel }}C{{U}_{{1y}}}} \mathord{\left/ {\vphantom {{{{E}_{\parallel }}C{{U}_{{1y}}}} {({{U}_{{1y}}}{{U}_{{{\text{p}}x}}}--{{U}_{{1x}}}{{U}_{{{\text{p}}y}}})}}} \right. \kern-0em} {({{U}_{{1y}}}{{U}_{{{\text{p}}x}}}--{{U}_{{1x}}}{{U}_{{{\text{p}}y}}})}}. \\ \end{gathered} $$

(13b)

If the amplitude E_∥ of an incident wave with vertical polarization is given, then the structure of the field below the layer can be determined from relations (13). Similarly, if the amplitude of the field E_⊥ of an incident wave with horizontal polarization is given, then

(14а)

$$\begin{gathered} {{A}_{ \bot }} = {{ - {{E}_{ \bot }}{{U}_{{{\text{p}}x}}}} \mathord{\left/ {\vphantom {{ - {{E}_{ \bot }}{{U}_{{{\text{p}}x}}}} {({{U}_{{1y}}}{{U}_{{{\text{p}}x}}}--{{U}_{{1x}}}{{U}_{{{\text{p}}y}}})}}} \right. \kern-0em} {({{U}_{{1y}}}{{U}_{{{\text{p}}x}}}--{{U}_{{1x}}}{{U}_{{{\text{p}}y}}})}}; \\ {{B}_{ \bot }} = {{ - {{E}_{ \bot }}{{U}_{{1x}}}} \mathord{\left/ {\vphantom {{ - {{E}_{ \bot }}{{U}_{{1x}}}} {({{U}_{{1y}}}{{U}_{{{\text{p}}x}}}--{{U}_{{1x}}}{{U}_{{{\text{p}}y}}})}}} \right. \kern-0em} {({{U}_{{1y}}}{{U}_{{{\text{p}}x}}}--{{U}_{{1x}}}{{U}_{{{\text{p}}y}}})}}. \\ \end{gathered} $$

(14b)

If a wave, e.g., with horizontal polarization is incident on the anisotropic layer, then waves with both horizontal and vertical polarizations will be reflected. Therefore, it is customary to represent the reflection coefficient as a 2 × 2 matrix $\hat {R}$ [18]. Here, we define of the elements of the matrix:

$$\begin{gathered} {{R}_{{11}}} = {{H_{y}^{ \downarrow }} \mathord{\left/ {\vphantom {{H_{y}^{ \downarrow }} {H_{y}^{ \uparrow }}}} \right. \kern-0em} {H_{y}^{ \uparrow }}}\,\,{\text{for}}\,\,{{E}_{ \bot }} = 0; \\ {{R}_{{12}}} = {{{{Z}_{0}}H_{y}^{ \downarrow }} \mathord{\left/ {\vphantom {{{{Z}_{0}}H_{y}^{ \downarrow }} {E_{y}^{ \uparrow }}}} \right. \kern-0em} {E_{y}^{ \uparrow }}}\,\,{\text{for}}\,\,{{E}_{\parallel }} = 0; \\ \end{gathered} $$

(15а)

$$\begin{gathered} {{R}_{{21}}} = {{E_{y}^{ \downarrow }} \mathord{\left/ {\vphantom {{E_{y}^{ \downarrow }} {H_{y}^{ \uparrow }}}} \right. \kern-0em} {H_{y}^{ \uparrow }}}\,\,{\text{for}}\,\,{{E}_{ \bot }} = 0; \\ {{R}_{{22}}} = {{E_{y}^{ \downarrow }} \mathord{\left/ {\vphantom {{E_{y}^{ \downarrow }} {E_{y}^{ \uparrow }}}} \right. \kern-0em} {E_{y}^{ \uparrow }}}\,\,{\text{for}}\,\,{{E}_{\parallel }} = 0. \\ \end{gathered} $$

(15b)

From relations (10), (13), (14) and (15), we obtain

$${{R}_{{11}}} = {{({{U}_{{1y}}}{{D}_{{{\text{p}}x}}} - {{U}_{{{\text{p}}y}}}{{D}_{{1x}}})} \mathord{\left/ {\vphantom {{({{U}_{{1y}}}{{D}_{{{\text{p}}x}}} - {{U}_{{{\text{p}}y}}}{{D}_{{1x}}})} {({{U}_{{1x}}}{{U}_{{{\text{p}}y}}} - {{U}_{{1y}}}{{U}_{{{\text{p}}x}}})}}} \right. \kern-0em} {({{U}_{{1x}}}{{U}_{{{\text{p}}y}}} - {{U}_{{1y}}}{{U}_{{{\text{p}}x}}})}},$$

(16а)

$${{R}_{{12}}} = {{({{D}_{{1x}}}{{U}_{{{\text{p}}x}}} - {{U}_{{1x}}}{{D}_{{{\text{p}}x}}})} \mathord{\left/ {\vphantom {{({{D}_{{1x}}}{{U}_{{{\text{p}}x}}} - {{U}_{{1x}}}{{D}_{{{\text{p}}x}}})} {({{U}_{{1x}}}{{U}_{{{\text{p}}y}}} - {{U}_{{1y}}}{{U}_{{{\text{p}}x}}})C}}} \right. \kern-0em} {({{U}_{{1x}}}{{U}_{{{\text{p}}y}}} - {{U}_{{1y}}}{{U}_{{{\text{p}}x}}})C}},$$

(16b)

$${{R}_{{21}}} = {{C({{U}_{{{\text{p}}y}}}{{D}_{{1y}}} - {{U}_{{1y}}}{{D}_{{{\text{p}}y}}})} \mathord{\left/ {\vphantom {{C({{U}_{{{\text{p}}y}}}{{D}_{{1y}}} - {{U}_{{1y}}}{{D}_{{{\text{p}}y}}})} {({{U}_{{1x}}}{{U}_{{{\text{p}}y}}} - {{U}_{{1y}}}{{U}_{{{\text{p}}x}}})}}} \right. \kern-0em} {({{U}_{{1x}}}{{U}_{{{\text{p}}y}}} - {{U}_{{1y}}}{{U}_{{{\text{p}}x}}})}},$$

(16c)

$${{R}_{{22}}} = {{({{U}_{{1y}}}{{D}_{{{\text{p}}x}}} - {{U}_{{{\text{p}}y}}}{{D}_{{1x}}})} \mathord{\left/ {\vphantom {{({{U}_{{1y}}}{{D}_{{{\text{p}}x}}} - {{U}_{{{\text{p}}y}}}{{D}_{{1x}}})} {({{U}_{{1x}}}{{U}_{{{\text{p}}y}}} - {{U}_{{1y}}}{{U}_{{{\text{p}}x}}})}}} \right. \kern-0em} {({{U}_{{1x}}}{{U}_{{{\text{p}}y}}} - {{U}_{{1y}}}{{U}_{{{\text{p}}x}}})}}.$$

(16d)

Let us determine the polarization $\rho$ of the basic solutions:

$${{\rho }_{{\text{p}}}} = {{{{E}_{{{\text{p}} \bot }}}} \mathord{\left/ {\vphantom {{{{E}_{{{\text{p}} \bot }}}} {{{E}_{{{\text{p}}\parallel }}}}}} \right. \kern-0em} {{{E}_{{{\text{p}}\parallel }}}}},\,\,\,\,~{{\rho }_{1}} = {{{{E}_{{1 \bot }}}} \mathord{\left/ {\vphantom {{{{E}_{{1 \bot }}}} {}}} \right. \kern-0em} {}}{{E}_{{1\parallel }}},$$

(17)

because ${{\vec {e}}_{1}}$, ${{\vec {e}}_{{\text{p}}}}$ are orthogonal, |ρ_p||ρ₁| = 1 and argρ₁ = arg ρ_p + π. From expansions (13) and (14), we obtain

$${{\rho }_{{\text{p}}}} = {{C{{U}_{{{\text{p}}y}}}} \mathord{\left/ {\vphantom {{C{{U}_{{{\text{p}}y}}}} {{{U}_{{{\text{p}}x}}}}}} \right. \kern-0em} {{{U}_{{{\text{p}}x}}}}},\,\,\,\,~{{\rho }_{1}} = {{C{{U}_{{1y}}}} \mathord{\left/ {\vphantom {{C{{U}_{{1y}}}} {{{U}_{{1x}}}}}} \right. \kern-0em} {{{U}_{{1x}}}}}.$$

(18)

The transmission coefficients of the electric field of the wave are easily obtained, taking into account that only one penetrating wave passes to the level z₁. We obtain two coefficients:

$$\begin{gathered} {{t}_{\parallel }} = {{{{E}_{{{\text{p}}x}}}({{z}_{1}})} \mathord{\left/ {\vphantom {{{{E}_{{{\text{p}}x}}}({{z}_{1}})} {{{E}_{\parallel }}}}} \right. \kern-0em} {{{E}_{\parallel }}}}\,\,{\text{for}}\,\,{{E}_{ \bot }} = 0; \\ ~{{t}_{ \bot }} = {{{{E}_{{{\text{p}}x}}}({{z}_{1}})} \mathord{\left/ {\vphantom {{{{E}_{{{\text{p}}x}}}({{z}_{1}})} {{{E}_{ \bot }}}}} \right. \kern-0em} {{{E}_{ \bot }}}}\,\,{\text{for}}\,\,{{E}_{\parallel }} = 0. \\ \end{gathered} $$

(19)

These expressions can be rewritten in simpler form:

$${{t}_{ \bot }} = {{C{{E}_{{{\text{p}}x}}}({{z}_{1}})} \mathord{\left/ {\vphantom {{C{{E}_{{{\text{p}}x}}}({{z}_{1}})} {{{U}_{{{\text{p}}x}}}({{\rho }_{p}} - {{\rho }_{1}})}}} \right. \kern-0em} {{{U}_{{{\text{p}}x}}}({{\rho }_{p}} - {{\rho }_{1}})}} = {{ - {{t}_{\parallel }}} \mathord{\left/ {\vphantom {{ - {{t}_{\parallel }}} {{{\rho }_{1}}}}} \right. \kern-0em} {{{\rho }_{1}}}}.$$

(20)

In conclusion, we present the expressions for the power transmission coefficients. The power of the penetrating wave mode incident on the layer along the z axis is

$${{P}_{{z2}}} = C{{({{{\left| {{{U}_{{{\text{p}}x}}}} \right|}}^{2}} + {{{\left| {{{U}_{{{\text{p}}y}}}} \right|}}^{2}} + {{{\left| {{{U}_{{{\text{p}}z}}}} \right|}}^{2}})} \mathord{\left/ {\vphantom {{({{{\left| {{{U}_{{{\text{p}}x}}}} \right|}}^{2}} + {{{\left| {{{U}_{{{\text{p}}y}}}} \right|}}^{2}} + {{{\left| {{{U}_{{{\text{p}}z}}}} \right|}}^{2}})} {2{{Z}_{0}}}}} \right. \kern-0em} {2{{Z}_{0}}}}.$$

Above the layer for z = z₁, the vertical part of the power is formed by the initial field of the whistling mode ${{\vec {e}}_{{\text{p}}}}$:

$${{P}_{{z1}}} = {{{\text{Re}}\left( {{{E}_{x}}H_{y}^{*} - {{E}_{y}}H_{x}^{*}} \right)} \mathord{\left/ {\vphantom {{{\text{Re}}\left( {{{E}_{x}}H_{y}^{*} - {{E}_{y}}H_{x}^{*}} \right)} 2}} \right. \kern-0em} 2}.$$

The vertical power transmission coefficient D_pz is determined as

$$\begin{gathered} {{D}_{{{\text{p}}z}}} = {{{{P}_{{z1}}}} \mathord{\left/ {\vphantom {{{{P}_{{z1}}}} {{{P}_{{z2}}}}}} \right. \kern-0em} {{{P}_{{z2}}}}} \\ = {{{{Z}_{0}}\operatorname{Re} \left( {{{E}_{x}}H_{y}^{*} - {{E}_{e}}H_{x}^{*}} \right)} \mathord{\left/ {\vphantom {{{{Z}_{0}}\operatorname{Re} \left( {{{E}_{x}}H_{y}^{*} - {{E}_{e}}H_{x}^{*}} \right)} {C\left( {{{{\left| {{{U}_{{{\text{p}}x}}}} \right|}}^{2}} + {{{\left| {{{U}_{{{\text{p}}y}}}} \right|}}^{2}} + {{{\left| {{{U}_{{{\text{p}}z}}}} \right|}}^{2}}} \right)}}} \right. \kern-0em} {C\left( {{{{\left| {{{U}_{{{\text{p}}x}}}} \right|}}^{2}} + {{{\left| {{{U}_{{{\text{p}}y}}}} \right|}}^{2}} + {{{\left| {{{U}_{{{\text{p}}z}}}} \right|}}^{2}}} \right)}}. \\ \end{gathered} $$

(21)

Let us determine the total power transmission coefficient D_p:

$${{D}_{{\text{p}}}} = C{{D}_{{{\text{p}}z}}}\sqrt {1 + s_{x}^{2} + s_{y}^{2}} ,$$

(22)

where the relative “side” power fluxes s_{x, y} at height z₁ are calculated as the ratios of the corresponding components of the Umov–Poynting vector to its z-component:

$${{s}_{x}} = {{\operatorname{Re} \left( {{{E}_{y}}H_{z}^{*} - {{E}_{z}}H_{y}^{*}} \right)} \mathord{\left/ {\vphantom {{\operatorname{Re} \left( {{{E}_{y}}H_{z}^{*} - {{E}_{z}}H_{y}^{*}} \right)} {\operatorname{Re} \left( {{{E}_{x}}H_{y}^{*} - {{E}_{y}}H_{x}^{*}} \right)}}} \right. \kern-0em} {\operatorname{Re} \left( {{{E}_{x}}H_{y}^{*} - {{E}_{y}}H_{x}^{*}} \right)}},$$

$${{s}_{y}} = {{\operatorname{Re} \left( {{{E}_{z}}H_{x}^{*} - {{E}_{x}}H_{z}^{*}} \right)} \mathord{\left/ {\vphantom {{\operatorname{Re} \left( {{{E}_{z}}H_{x}^{*} - {{E}_{x}}H_{z}^{*}} \right)} {\operatorname{Re} \left( {{{E}_{x}}H_{y}^{*} - {{E}_{y}}H_{x}^{*}} \right)}}} \right. \kern-0em} {\operatorname{Re} \left( {{{E}_{x}}H_{y}^{*} - {{E}_{y}}H_{x}^{*}} \right)}},$$

the components of the field strength of the wave E_z and H_z are determined from relations (9a) and (9b).

3 CALCULATION RESULTS

Below, we everywhere use the numerical values of plasma parameters corresponding to the subpolar ionosphere at latitude 60° during the equinox with moderate solar activity [25]. Table 1 gives the percentage of the ionic components in the layer [z₁, z₂].

Table 1.   Parameters of daytime and nighttime ionosphere

The N_e(z) and ν_e(z) profiles were approximated by cubic splines. This method ensures the absence of “spurious” partial reflections in a strongly inhomogeneous layer of the lower ionosphere. For the daytime ionosphere, the 50–100 km layer was chosen; for the nighttime ionosphere, 75–150 km. At the upper boundary of the layer, the plasma frequency for electrons is 1268.2 kHz during the day and 1750.6 kHz at night; the electron gyrofrequencies are 1507.1 kHz during the day and 1472.7 kHz at night. For the ions, the corresponding frequencies are given in Table 1. The lower hybrid resonance frequencies f_LHR are 4.95 and 5.15 kHz, respectively.

The calculations were performed for θ = 10°, A = 180° (direction to magnetic south) and frequencies in the 50–1000 Hz range. The choice of these values was determined by the operating conditions of the HAARP station. However, it should be noted that the wavelength at a frequency of 50 Hz is 6000 km in the atmosphere. With an equator circumference of 40 000 km, six integer wavelengths fit on it. Such a wave process can only conditionally be considered “local,” and a further decrease in frequency requires consideration of the spherical Earth–ionosphere cavity in the form of a resonator, not a waveguide. Given, however, that in the lower ionosphere at height z₁ the refractive index of waves is relatively large, we obtain a value of the wavelength in the medium of ~30–40 km. This value can be considered “local.”

Figure 2 shows the calculated frequency dependence of the refractive index n of the whistling mode at z = z₁ for the daytime and nighttime ionosphere models. The n value for night is higher than the daytime values because the nighttime height z₁ is significantly greater than the daytime one.

Fig. 2.

Frequency dependences of refractive index n of whistling mode at z = z₁ for daytime (curve 1) and nighttime (curve 2) ionosphere models.

Figure 3 shows the frequency dependence of the total power transmission coefficient D (which in fact corresponds to D_p in the text) for the daytime and nighttime ionosphere models. For the ELF range, the daytime values D, as a rule, are much less than the nighttime ones due to collisional wave absorption at heights of 50–70 km in the daytime. Figure 3 shows that such a pattern begins to develop only at frequencies of f > 500 Hz. At lower frequencies, ion resonances in the daytime form a characteristic transmission maximum at frequencies of f ~ 30 Hz, which corresponds to approximately half the gyrofrequency of the lightest ion N⁺. Note that the D value at frequencies of 0.2–1 kHz is not small either during the day or night and almost reaches a value of 0.4. Excluding the influence of the geomagnetic field, the D value would be zero.

Fig. 3.

Frequency f dependence of power transmission coefficient D for daytime (curve 1) and nighttime (curve 2) ionosphere models.

Figure 4a shows the frequency dependence of the moduli of the matrix components ${\hat {R}}$ of the reflection coefficient (see (15), (16)) for day (solid curves) and night (dashed line). We present the f dependences on the quantities R₁₁ (curves 1, 4), R₁₂ (curves 2, 5), and R₂₂ for the day (curve 3). For the selected azimuth value, the figure shows (curves 1, 3) that the R₁₁ and R₂₂ values are very close, or almost coincide (as are the values of nondiagonal elements R₁₂ and R₂₁, of which the curve is given only for the quantity R₁₂). For this reason, Fig. 4b shows only the dependences of the arguments R₁₁(f), curves 1, 4; and arguments R₁₂(f), curves 2, 5.

Fig. 4.

Frequency f dependences of moduli (a) and arguments (b) of matrix components ${\hat {R}}$ of reflection coefficient for day (solid curves) and night (dashed line) for quantities R₁₁ (1, 4), R₁₂ (2, 5), R₂₂ (3).

Figure 4 shows that diagonal elements dominate in the chosen frequency range ${\hat {R}}$, and at night, the argument of the matrix elements is relatively stable. In the daytime, with increasing frequency, when moving from the region of ionic vibrations to the region of electronic waves, the argument of diagonal elements changes significantly (Fig. 4b, curve 1). Curve 3 for argument R₂₂ in Fig. 4b is absent, since it virtually coincides with curve 1. We do not give similar dependences of the polarization of the reflected wave as a whole, since the ionosphere plays the role of an almost ideal polarization filter in the used frequency range, so the polarization of the reflected wave is circular.

4 CONCLUSIONS

We have described an algorithm for calculating the characteristics of the passage of LF electromagnetic waves through the magnetoactive plane-layered lower ionosphere with parameters similar to those of the ionosphere above the HAARP station, Alaska. The algorithm is based on numerical integration of the wave equations for the case of a plane wave incident on the ionosphere from below, taking into account the ionic composition of the plasma. As an illustration of operation of the algorithm, the main characteristics of passage of the wave through the lower ionosphere in the 0.05–1 kHz frequency range are calculated.

The power transmission coefficient of these LF waves is quite large, in the range 0.2–0.4 for most of the frequencies used. In the daytime near a frequency of ~300 Hz, the transmission coefficient in the daytime has a maximum and is approximately double the nightly values. In the entire studied frequency range, the polarization of the reflected wave is circular with a good degree of accuracy.