HF wave propagation and induced ionospheric turbulence in the magnetic equatorial region

The propagation and excitation of artificial ionospheric turbulence in the magnetic equatorial region by high‐frequency electromagnetic (EM) waves injected into the overhead ionospheric layer is examined. EM waves with ordinary (O) mode polarization reach the critical layer only if their incidence angle is within the Spitze cone. Near the critical layer the wave electric field is linearly polarized and directed parallel to the magnetic field lines. For large enough amplitudes, the O mode becomes unstable to the four‐wave oscillating two‐stream instability and the three‐wave parametric decay instability driving large‐amplitude Langmuir and ion acoustic waves. The interaction between the induced Langmuir turbulence and electrons located within the 50–100 km wide transmitter heating cone at an altitude of 230 km can potentially accelerate the electrons along the magnetic field to several tens to a few hundreds of eV, far beyond the thresholds for optical emissions and ionization of the neutral gas. It could furthermore result in generation of shear Alfvén waves such as those recently observed in laboratory experiments at the University of California, Los Angeles Large Plasma Device.


Introduction
The interaction of large-amplitude high-frequency (HF) radio waves with the ionospheric plasma in the vicinity of the dip equator is significantly different than the well-studied interactions in the polar and middle latitude regions. Previous theoretical study [Erukhimov et al., 1997] proposed a variety of effects including the creation of quasi-periodic structures with vertical periodicity of approximately half the wavelength of the electromagnetic (EM) wave and extending horizontally along the magnetic field lines due to Ohmic heating of the plasma, as well as creation of a virtual antenna at ULF/ELF/VLF frequencies by modulation of the D/E conductivity of the equatorial electrojet similar to the modulation of the polar electrojet in the auroral region [Rietveld et al., 1984[Rietveld et al., , 1987[Rietveld et al., , 1989Papadopoulos et al., 1990Papadopoulos et al., , 2005Moore, 2007;Payne et al., 2007]. Advantages of the equatorial over the polar electrojet modulation were noted by Papadopoulos et al. [2005], while Erukhimov et al. [1997] emphasized the excitation of low-frequency Alfvén waves in the 1 Hz frequency range. EM waves exceeding a threshold amplitude could also drive artificial ionospheric turbulence (AIT) by the nonlinear interaction of ordinary (O) mode wave with the plasma electrons. In this case, since the O mode has its electric field along the ambient magnetic field lines when it reaches the critical layer, it could excite AIT associated with large amplitude Langmuir waves and ion acoustic waves propagating along the magnetic field. Experiments at the Large Plasma Device (LAPD) at University of California, Los Angeles (UCLA) [Van Compernolle et al., 2006;Wang et al., 2016] have also shown the formation of suprathermal electrons accompanied by the excitation of Alfvén waves when large-amplitude microwaves are injected into the plasma perpendicular to the magnetic field. Particle-in-cell simulations [Tsung et al., 2007] and subsequent laboratory experiments indicate that accelerated electrons due to resonant absorption of the O mode waves, and nonlinear processes involving wave collapse may explain the observation. Ionospheric high-latitude heating experiments at the High-Frequency Active Auroral Program (HAARP) have also demonstrated that AIT correlated with optical emissions and the formation of artificial plasma layers [Pedersen et al., 2009[Pedersen et al., , 2010 can be attributed to high-energy suprathermal electrons accelerated by the turbulence near the O mode critical layer to energies large enough to ionize the neutral gas [Mishin and Pedersen, 2011;Eliasson et al., 2012].
instabilities driven by the O mode to Langmuir and ion acoustic waves at different altitudes is calculated. The acceleration of electrons by the turbulence is calculated in section 5, and its implications for both ionospheric and laboratory conditions are discussed. Finally, conclusions are drawn in section 6.

Ionospheric Model and Ray Tracing
The results of this study are based on the use of an ionospheric model whose parameters are listed in Table 1. For simplicity we use a Gaussian vertical profile of the ionospheric layer of the form n 0 z ð Þ ¼ n 0; max exp h i , and we adopt a coordinate system such that the horizontal magnetic field B 0 is directed along the x axis.
To lowest order, the propagation of the HF O and extraordinary (X) modes can be estimated using a simple ray-tracing model [e.g., Whitham, 1974] dk where the wave frequency ω(k, r) is governed by the Appleton-Hartree dispersion relation based on a cold fluid plasma model [e.g., Stix, 1992] k 2 c 2 where Δ ¼ ω 2 ce sin 4 θ þ 4ω À2 ω 2 À ω 2 pe 2 cos 2 θ ! 1=2 and θ is the angle between k and B 0 , obtained by the cosine formula cos θ = k Á B 0 /(|k||B 0 |). Here c is the speed of light in vacuum, ω pe ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi n 0 e 2 =ϵ 0 m e p is the electron plasma frequency, and ω ce = eB 0 /m e is the electron cyclotron frequency, where e is the unit charge, m e is the electron mass and ϵ 0 is the electric vacuum permittivity. The upper sign (+) is for the O mode and the lower sign (À) is for the X mode. For numerical convenience, equation (3) can be written as Figures 1 and 2 show a comparison between the propagation of O and X mode waves injected at different angles of incidence for ω = ω 0 . The O mode is converted to Z mode (slow X mode) waves at two angles of incidence, χ = ± χ S , where the Spitze angle for a horizontal magnetic field (α = 90°) is given by [Mjølhus, 1990]  and where Y = ω ce /ω 0 . Within the Spitze cone, |χ| < χ S , the O mode wave reaches the critical altitude where ω = ω pe , where there is a strong interaction between the EM wave and electrostatic Langmuir waves propagating along the magnetic field. Since Y 2 ≈ 0.011 ≪ 1, the upper hybrid (UH) layer z UH is only ≈ L n Y 2 ≈ 500 m below the critical layer z O (cf. Figure 2 and Table 1). However, since the wave electric field of the O mode is almost parallel to the ambient magnetic field for incidence angles within the Spitze cone, the coupling to UH and lower hybrid waves is expected to be less important. The fast X mode wave is reflected at an altitude about 5 km below the critical layer, and therefore, its coupling to electrostatic waves is expected to be weak.   Figure 3 shows the absorption coefficient (absorbed intensity divided by incident intensity) given in equation (18) of Mjølhus [1990]. The normalized (by ω 0 /c) components of the wave vector of the O mode are N x = sin χ cos ϕ and N y = sin χ sin ϕ. The O mode is efficiently absorbed only in a narrow region (the radio window) within 1 or 3 degrees around the Spitze angle χ = χ S = 17.95°and ϕ = 0, where the O mode is converted to Z mode waves. As seen in the close-up in Figure 2, the mode converted Z mode wave propagates up to an altitude about 1 km above the critical layer of the O mode wave, after which it turns and propagates down to the electrostatic resonance at the UH layer. Outside the Spitze region (e.g., for 21°in Figure 2), the O mode is reflected below the critical layer.

Standing Wave Pattern and Swelling Near the O Mode Turning Point
To investigate the standing wave pattern and swelling of an O mode wave near its turning point, we carry out full-wave simulations for different angles of incidence using the model of Eliasson et al. [2015], with a linearly polarized (in the x-z plane) EM wave injected from the bottom side at z = 200 km with a reference amplitude of 1 V/m. A finite electron temperature of 3000 K is used, which has significance only for short wavelength electrostatic waves. As the EM wave reaches the turning point within the Spitze region, there is a swelling of the x component of the electric field, and a steady state Airy pattern of the standing wave is formed at time t = 1 ms, as seen in Figure 4. The wave amplitude at the first Airy maximum is about 5-8 times larger than the injected amplitude. On the other hand, both the y and z components (not shown here) are only weakly excited within the Spitze cone. Hence, it is expected that the main nonlinear interaction in this region is Langmuir turbulence due to the coupling of the O mode wave to Langmuir and ion acoustic waves propagating along the magnetic field lines, while the coupling to UH and lower hybrid waves propagating across the magnetic field lines is likely to be less significant.
At the Spitze angle χ = χ S = 17.95°( Figure 4) there are large amplitude mode-converted Z mode waves reaching an altitude of 240 km, about 1 km above the critical layer of the O mode wave. This is consistent with the ray-tracing results in Figure 2. As the Z mode wave is reflected and propagates downwards, it turns into a short wavelength electrostatic UH wave as it passes the UH resonance layer z UH where ω 0 = ω UH . Large amplitude vertical (z) component of the electric field (not shown) for χ = 17.95°could potentially lead to nonlinear excitations of UH and lower hybrid turbulence in a small region (cf. Figure 3) near the Spitze angle.

Threshold for Artificial Ionospheric Turbulence
A vertically (along z) injected O mode wave has the wave electric field parallel to the ambient magnetic field in the x direction. Above a threshold amplitude, the O mode will be unstable and decay to high-frequency Langmuir waves and low-frequency ion density fluctuations. The swelling of the EM wave seen in Figure 4 near the O mode turning point significantly boosts the nonlinearity. In general, the four-wave oscillating two-stream instability (OTSI) dominates at higher altitudes near the critical layer, while at lower altitudes the three-wave parametric decay instability (PDI) dominates [e.g. Mjølhus et al., 1995]. Figure 5 shows a schematic of the two processes. The OTSI excites nonresonant, purely growing ion density fluctuations and saturates nonlinearly by strong turbulence associated with wave localization and collapse (nucleation and burnout cycles) [e.g., Russell et al., 1988], while the PDI excites resonant Langmuir and ion acoustic waves associated with weak turbulence where the involved waves approximately obey the respective linear dispersion relations.  The growth rates of the PDI and OTSI can be obtained from a model involving large amplitude electric field oscillations coupled with low-frequency ion fluctuations via the ponderomotive force acting on the electrons [e.g., Papadopoulos et al., 1974;Freund and Papadopoulos, 1980]. Parallel to the magnetic field, the dispersion relation for the growth rate of the instability is (see details in Appendix A) Te k 2 1 þ ω 2 pe represent the sidebands of the Langmuir wave and E 0 is the local amplitude of the O mode pump electric field. The perturbation wave number and frequency is denoted k 1 and ω 1 = ω R + iγ, respectively, where ω R is the real frequency and γ is the growth rate. Here is the ion acoustic speed, where m i is the ion mass (we use atomic oxygen ions) and k B is Boltzmann's constant. Near the critical layer, where ω 0 ≈ ω pe , the OTSI dominates, while the PDI dominates at lower altitudes where the frequency mismatch Δω = ω 0 À ω pe is larger.
In equation (6), ν e = ν ei + ν eL and ν i = ν iL are electron and ion damping rates due to collisions and Landau damping. The electron-ion collision frequency due to Coulomb collisions may be estimated by De is the plasma parameter. In the numerical work below, we use ν ei = 10 3 s À 1 , which is consistent with the ionospheric parameters in Table 1. The relative importance of collisions increases with density, which has to be taken into account when comparing ionospheric and laboratory experiments. The electron and ion Landau damping rates are approximated by [e.g., Krall and Trivelpiece, 1986] where R = T e /T i is the electron-to-ion temperature ratio and λ De = v Te /ω pe is the electron Debye length. The used T e = 3000 K and T i = 1000 K gives significant ion Landau damping. Electron Landau damping of Langmuir waves becomes significant when k 1 λ De ≳ 0.2, which for our parameters (with λ De ≈ 4 mm) is for wave numbers k 1 ≳ 50 m À 1 corresponding to wavelengths smaller than about 13 cm.
On longer length scale and time scale, the coupling to collisional modes and the plasma drift due to gradients [Stenflo, 1985] may become important for the formation of quasi-periodic structures in the vertical direction [e.g., Erukhimov et al., 1997] and other temperature-driven effects. However, while thermal instabilities in ionospheric heating experiments may develop on time scales of the order of a second [e.g., Gurevich, 1978], Langmuir turbulence usually develops on millisecond time scales. We therefore neglect the effect of plasma temperature changes in the present treatment. Figures 6-9 show numerical solutions of equation (6) as a function of perturbation wave number at different altitudes, where estimates of the O mode pump electric field E 0 = |E x | at different altitudes are taken from the standing wave pattern for vertical incidence (i.e., for 0°in Figure 4). Since the length scale of the O mode wave is about 2-3 orders of magnitude larger than that of the Langmuir waves, we assume that E 0 (as well as n 0 ) is locally independent of space in the numerical solution of equation (6). As seen in Figures 6 and 8, the instability takes place in a quasi-periodic pattern with maximum growth rate at the Airy maxima of the O mode wave. Figures 6 and 7 show the growth rate of the instabilities for an injected O mode amplitude of E EM = 1 V/m. It is seen in Figures 6b and 6c that the instability takes place in a region covering about 4-5 km below the critical layer. This is in contrast to the high-latitude cases [e.g., Eliasson and Stenflo, 2008;Eliasson et al., 2012Eliasson et al., , 2015, where the region of Langmuir turbulence covers only a fraction of a kilometer.  With a growth rate of the order 10 4 s À 1 (cf. Figure 6b) the instability will typically saturate nonlinearly within a few milliseconds after switch on of the heating. Near the critical layer, the instability gives rise to Langmuir waves and ion density fluctuations at wave number of the order k 1 ∼ 5-15m À 1 (cf. Figure 6c) corresponding to wavelengths of the order 0.5-1 m, while at lower altitudes the instability takes place for larger wave numbers (shorter wavelengths) of the generated waves. The electron Landau damping becomes strong enough to quench the instability about 4.5 km below the critical layer, where the wave number of the Langmuir waves is sufficiently large, k 1 ≈ 45 m À 1 (cf. Figure 6c) corresponding to a wavelength of about 15 cm.
The line plots in Figure 7 show that the OTSI has a higher the growth rate than the PDI at the first Airy maximum (z = 238.85 km) where the frequency mismatch Δω = ω 0 À ω pe = 7.7 × 10 3 s À 1 (1.2 kHz), while the PDI has a larger growth rate at the lower altitude 2 km below the critical layer (z = 236.98 km) where Δω = 1.2 × 10 6 s À 1 (190 kHz). While the OTSI gives rise to purely growing ion fluctuations (cf. Figure 7b), the PDI is an oscillating instability giving rise to ion acoustic waves with frequencies near the unperturbed ion acoustic frequency.
While the relatively large injected O mode amplitude of E EM = 1 V/m is believed to be typical for large facilities such as HAARP with access to large transmitted powers and antenna gains, such high amplitudes may be out of reach for proposed  portable heating facilities in the equatorial region. It is therefore of interest to investigate a case of significantly lower amplitude. As an example, we take an injected amplitude of 0.2 V/m, corresponding to a factor 25 smaller power than for 1 V/m. Most importantly, it is seen in Figure 8 that the excitation of AIT is still possible at this relatively low O mode amplitude. This is due to the swelling of the injected O mode near the critical layer with a local increase of the amplitude by a factor 5-8 compared to the free space amplitude, above the threshold for instability. The instability takes place in a somewhat smaller, 3.5 km wide region, and with an order of magnitude smaller growth rate. For this case, the PDI has in general a larger growth rate than the OTSI (Figure 9), and the OTSI quickly vanishes at lower altitudes (cf. Figure 9c).

Nonlinear Evolution and Formation of Suprathermal Electrons
To study the turbulence induced by the large amplitude EM wave and the resulting electron acceleration, we carry out local simulations using a two-dimensional (2-D) generalized Zakharov model in the x-z plane, in which the envelope of the HF electrostatic waves is nonlinearly coupled to the slowly varying ion density fluctuations via the ponderomotive force acting on the electrons. Such local turbulence simulations have previously been used to study turbulence and stimulated EM emissions from unmagnetized and magnetized plasma [Mjølhus et al., 1995[Mjølhus et al., , 2003]. The electron continuity and momentum equations of motion for the complex valued envelopes of the high-frequency electron density and velocity are ∂n e ∂t ¼ iω 0 n e À ∇Á n s v e ð Þ and ∂v e ∂t x. The envelope of the HF electrostatic potential, ϕ, is obtained from Poisson's equation The external electric x is a dipole field (constant in space) along the x axis, representing the O mode wave, where we will use the values of E 0 discussed in Figures 7 and 9. The influence of the low-frequency fluctuations on the high-frequency dynamics is through the slow time scale electron and ion density n es = n is = n s in equation (9). The low-frequency ion motion is governed by the ion continuity and momentum equations and ∂v s ∂t where the second term on the right-hand side is due to the ponderomotive force acting on the electrons [see, e.g., Thornhill and Ter Haar, 1978]. At equilibrium, n s = n 0 , where n 0 is the local density given by the altitudedependent ionospheric profile. Figures 10 and 11 show numerical solutions of equations (9)-(13) for the large and small amplitude of the pump electric field used in Figures 7 and 9, respectively. The spatial derivatives as well as the Landau damping effects are accurately approximated in the simulations using pseudo spectral methods, in which the formulas (7) and (8) are used to calculate electron and ion Landau damping in Fourier space [cf. Mjølhus et al., 2003]. The quadratic nonlinearities are dealiased using a two-third rule based on zero padding in Fourier space before multiplication [e.g., Gumerov et al., 2011]. A two-dimensional domain is used, resolved by 250 grid points with periodic boundary conditions in both the x and z direction. The domain size is taken to be 10 × 10 m in Figures 10a-10c and 11a-11c and 3 × 3 m in Figures 10d-10f and 11d-11f to accurately resolve the turbulence. A standard fourth-order Runge-Kutta scheme is used to advance the solution in time, with a time step Δt = 10 À 8 s À 1 . Small-amplitude density perturbations (random numbers) of the order 10 7 m À 3 are added to n e to seed the instability. The instability in general first develops with electrostatic waves along the magnetic field lines, consistent with equation (6). After the initial nonlinear one-dimensional saturation, the turbulence in some cases enters a second phase of OTSI-generated 2-D nucleation and burnout cycles similar to that discussed by Russell et al. [1988] and others. The simulations initially show exponential growth of the PDI and OTSI instabilities in time, after which a transition to steady state turbulence takes place. The results in Figures 10 and  11 are recorded when the turbulence has reached steady state.
For the large amplitude injected O mode, the OTSI at the higher altitude at first Airy maximum of the O mode wave (cf. Figure 10a-10c) leads to turbulence characterized by nucleation and burnout cycles [cf. Russell et al., 1988], in which large amplitude localized electric fields self-trapped in ion density cavities collapse and dissipate. At the lower altitude (cf. Figure 10d-10f), the PDI leads to the excitation of propagating ion acoustic waves, in addition to the OTSI-generated nucleation and burnout cycles.
For the smaller injected O mode amplitude E EM = 0.2 V/m is shown in Figure 11. The turbulence near the first Airy maximum, shown in Figures 11a-11c, is characterized by PDI-generated wave turbulence together with only a few nucleation and burnout events due to the OTSI, giving rise to transient large-amplitude electric fields correlated with ion cavities. In this process, ion acoustic waves are excited and propagate away radially giving rise to a few circular patterns visible in Figure 11c. At the lower altitude, displayed in Figures 11d-11f, the instability is dominated by the PDI, while the OTSI is absent (cf. Figure 9d), and the nonlinear saturation is characterized by weak turbulence involving propagating Langmuir and ion acoustic waves, but no nucleation and burnout events are visible.
The x components of the HF electric fields in Figures 10 and 11 are used to estimate the acceleration and transport of electrons along the magnetic field lines by the turbulent electric field. When the electrons have sufficiently high velocities, they feel an almost constant electric field when they pass over localized electric field envelopes. This leads to a random walk process and a diffusion of the electron distribution in velocity space. As a crude approximation, this process is modeled by a one-dimensional Fokker-Planck equation for the averaged electron distribution [e.g., Sagdeev and Galeev, 1969] with the diffusion coefficient where W ω 0 ; k x ð Þ¼ ΔE 2 x =Δk x is the spectral energy density of the electric field per wave number Δk and ΔE 2 is the differential squared electric field. Here W is in V 2 /m and is normalized such that the total spectral energy equals the mean squared electric field where L x is the width of the turbulent region. Extensions of the model to multiple dimensions leads to a set of integrodifferential equations [e.g., Vedenov et al., 1961;Ivanov et al., 1975]. A convective term has been added to the left-hand side of equation (14) to account for the transport of electrons through a turbulent layer of finite width. To numerically construct the diffusion coefficient (15), the energy spectrum of E x is calculated along the x direction and is averaged in the y direction over the simulation domain. An average is also taken over 10 different times covering 10 À 4 s. The resulting diffusion coefficients are shown in Figure 12 [2006] where the width of the interaction region of 10 cm is approximately three wavelengths of the 9 GHz microwave. In Figure 13a we see that the electrons passing through this region form significant suprathermal tails with the fastest electrons having a velocity of approximately 0.5 × 10 7 m/s, corresponding to about 70 eV energy. This is in line with the experiment [Van Compernolle et al., 2006] where the highest energy of the electrons were about 75 eV. For a 10 times wider region shown in Figure 13b, which may be more representative of ionospheric experiments, the highest energy suprathermal electrons reach about 10 7 m/s, corresponding to about 280 eV. The line plots in Figure 13c show the distribution function of the electrons streaming out of the turbulent region for the respective cases. It should be noted that electrons above about 2 eV energy lead to optical emissions by exciting the 1-D state of atomic oxygen, while energies above 12 eV leads to ionization of molecular oxygen [e.g., Pedersen et al., 2009;Rees, 1989]. The suprathermal electrons far exceed these thresholds, and hence, both optical emissions and the formation of artificial plasma near the Airy maxima of the O mode wave are expected for this case. For the case of a weaker injected O mode amplitude of 0.2 V/m, strong turbulence is observed only at the higher altitude near the first Airy maximum (cf. Figures 11a-11c), and the corresponding diffusion coefficient derived from the turbulent field in Figure 12c is relatively small amplitude. The Fokker-Planck simulation in Figure 15a shows the excitation of relatively weak suprathermal tails. For the wider (1 km) region in Figure 15b, a small fraction of electrons reaches a velocity of about 2 × 10 6 m/s, corresponding to about 11 eV energy. The turbulence would likely give rise to optical emissions but not to significant ionization of the neutral gas. At the lower altitude (cf. Figures 11d-11f) the weak wave turbulence gives rise to a narrow wave spectrum and to diffusion coefficient which only has two narrow peaks (cf. Figure 12d), and for this case no significant suprathermal tails are created in the electron distribution function (not shown). The injected amplitude of E EM = 0.2 V/m may be regarded as a lower bound for inducing artificial ionospheric turbulence strong enough to efficiently accelerate electrons. Thermal expansion near the O mode Airy maxima [Erukhimov et al., 1997] could here lead to periodic patterns extending along the magnetic field lines.

Conclusions
We have carried out an extensive investigation of the propagation and nonlinear interaction between large-amplitude EM waves and the ionospheric plasma in the magnetic equatorial region. The equatorial region has unique features that distinguish it from the midlatitude and auroral regions. Since the magnetic field is horizontal, O mode waves injected vertically and within the Spitze region will be linearly polarized near the critical layer with the wave electric field parallel to the magnetic field. This leads to the excitation of Langmuir turbulence in a several kilometers wide region below the critical layer and efficient electron acceleration along the magnetic field lines due to resonant interactions with the turbulent Langmuir waves.
Simulations of the nonlinear evolution of the instability and of the interaction between the electrons and the turbulent electric field show that suprathermal tails of electrons are created by the turbulence and can reach energies of several tens to a few hundreds of eV, which is far above the thresholds for optical emissions and ionization of the neutral gas. The suprathermal electrons could also excite shear Alfvén waves such as those observed in recent experiments at the UCLA LAPD device [Van Compernolle et al., 2006;Wang et al., 2016].