Linear and Nonlinear Plasma Processes in Ionospheric HF Heating

: Featured observations of high frequency (HF) heating experiments are ﬁrst introduced; the uniqueness of each observation is presented; the likely cause and physical process of each observed phenomenon instigated by the HF heating are discussed. A special point in the observations, revealed through the ionograms, is the competition between the Langmuir parametric instability and upper hybrid parametric instability excited in the heating experiments and the impact of the natural cusp at foE (the peak plasma frequency of the ionospheric E region) on the competition. The ionograms also infer the generation of Langmuir and upper hybrid cavitons. Ray tracing theory is formulated. With and without the appearance of large-scale ﬁeld-aligned density irregularities in the background ionosphere, ray trajectories of the ordinary mode (O-mode) and extraordinary mode (X-mode) sounding pulses are calculated numerically. The results explain the artiﬁcial Spread-F recorded by the digisondes in the heating experiments. Parametric instabilities, which are the directly relevant processes to achieve effective heating of the ionospheric F region, are formulated and analyzed. The threshold ﬁelds and growth rates of Langmuir and upper hybrid parametric instabilities are derived as the theoretical basis of many radar observations and electron-plasma wave interactions. Harmonic cyclotron resonance interaction processes between electrons and upper hybrid waves are introduced. Formulation and analysis are presented. The numerical results show that ultra-energetic electrons are generated. These electrons enhance airglow at 777.4 nm as well as cause ionization. Physical processes leading to the generation of artiﬁcial ionization layers are discussed. The nonlinear Schrodinger equation governing the nonlinear evolution of Langmuir waves and upper hybrid waves are derived and solved. The nonlinear periodic and solitary solutions of the equations are obtained. The localized Langmuir and upper hybrid waves generated by the HF heater form cavitons near the HF reﬂection layer and near the upper hybrid resonance layer, which induce bumps in the virtual height spread of the ionogram trace similar to that induced by the density cusp at E-F1 transition layer; the down-going Langmuir waves and upper hybrid waves evolve into nonlinear periodic waves propagating along the magnetic ﬁeld, which backscatter incoherently the sounding pulses to cause downward virtual height spread.


Introduction
Using ionosphere as an unbounded plasma laboratory to study linear and nonlinear plasma processes, instigated by high frequency (HF) pump waves, has been a highly active research area over the past four decades [1][2][3]. Experiments were conducted by transmitting powerful HF waves from the ground to observe Ionospheric modification via remotesensing instruments. A major facility has been built in Gakona, Alaska for conducting ionospheric heating experiments, as part of the High-Frequency Active Auroral Research Program (HAARP) [4]. The HAARP HF transmitting system has a rectangular planar array of 180 elements, which consist of a low band (2.8 to 7.6 MHz) and a high band (7.6 to 10 MHz) crossed dipole antenna in each element. Each crossed dipole radiates circularly polarized wave up to 20 kW, so that the HAARP HF transmitter radiates circularly polarized waves, in the frequency band from 2.8 MHz to 10 MHz, up to 3.6 MW. The antenna gain increases from 15 dB to 30 dB with an increase of the radiating frequency from 2.8 MHz to 10 MHz. An effective radiated power (ERP) up to 93 dBw (at 10 MHz) is available in heating experiments, which explore modification effects on the bottom-side of the ionosphere as illustrated in Figure 1. Through the observations of various heating induced phenomena presented in Section 2, considerable advances are toward the theory of wave-wave and wave-particle interactions. Electromagnetic waves interact with charged particles; only electrons can effectively respond to the fast oscillation of the HF wave electric field. Through elastic and inelastic collisions with charged particles (mainly electrons), neutral particles can also be indirectly affected by the presence of HF waves. In D and E regions, electron-neutral collision frequency is higher than electron-electron and electron-ion collision frequencies. Although neutral particles can share considerable wave energy with charged particles and rapidly thermalize the indirectly absorbed wave energy, the background temperature elevation is small due to high neutral-particle density. Nevertheless, a long heating period in the daytime could still make an impact on the ionization balance, which involves the processes of (1) photoionization by the solar illumination, (2) electron-heating-caused reduction of the recombination coefficient, and (3) enhancement of the electron attachment to the oxygen molecules. Moreover, electrojet appeared in these regions can be modulated by intensity modulated HF heating waves to become a virtual antenna, which radiates ELF/VLF waves [5]. Applying an HF heating facility to establish an ionospheric virtual ELF/VLF transmitter for undersea communications is another focus of research.
However, most of the HF heating experiments were focused on the F region modification on the bottom side of the ionosphere, where the HF heating wave induces significant electron quiver motion in the region below and near the HF reflection height; parametric instabilities are excited to convert EM wave to plasma (ES) waves [6]. This is done by employing O-mode heating wave with frequency less than the maximum cutoff frequency (FoF2) of the ionosphere to confine the HF wave in the bottom side of the ionosphere. The excitation of parametric instabilities is essential because the collision (electron-electron and electron-ion) processes cannot efficiently absorb the electromagnetic (EM) wave energy delivered to the F region of the ionosphere, and a fast conversion of EM wave into electrostatic (ES) waves (via parametric instabilities) minimizes the loss of the heating power due to reflection back to the ground. The excited plasma waves interact with electrons as well as are interacting among themselves; as a result, new phenomena have been observed in the experiments, which are called for theoretical interpretations.
The present work is aimed at providing theoretical foundation for the understanding of experimental observations in active wave-ionosphere interaction and the underlying In Arecibo HF heating experiments, overshoot [7,8] of the intensity and down shifting [9,10] of the originating height in time of the HFPLs were observed. It is realized that the relevant parametric instabilities prefer to excite Langmuir waves along the geomagnetic field. On the other hand, the Langmuir waves ascribed to the HFPLs, which must have an oblique propagation angle of 40 • at the Arecibo site, are not the preferred ones. Thus, these Langmuir waves are suppressed in time by those preferred ones having much larger growth rates (mode competition nonlinear damping). A spectral distribution of Langmuir waves, excited by the HF heater via parametric instabilities, interact with background electrons. The non-resonance interaction with bulk electrons (wave phase velocity is much larger than the electron velocity) in the velocity distribution broadens the distribution to have an equivalent temperature T eff = T e0 + T w , where T e0 is the electron temperature of unperturbed background plasma and T w = 2 0 ∑ k ω 2 pe |E 1k | 2 /n 0 ω 2 kr is the time average wave energy per electron. The matching height of the Langmuir wave decreases with the increase of the effective electron temperature, which is proportional to the total spectral intensity of the Langmuir waves.
Measurements also show asymmetry [11] of the HFPL distribution on the negative side (corresponding to up-going plasma lines) and on the positive side (down-going plasma lines). It is realized that the spectral intensity on the positive side of the HFPLs is ascribed to the excited down-going plasma waves as well as the excited up-going plasma waves after being reflected. Both OTSI and PDI were excited very close to the HF reflection height, and the reflected back plasma waves did not attenuate significantly when combining with the down-going plasma waves.

Competition between Langmuir PDI and Upper Hybrid PDI
Parametric instabilities excited in the HF heating experiments have been monitored by UHF/VHF radars [1,3], which receive return signals backscattered by the plasma waves. Because of the imposed Bragg backscattering conditions, the exploration of parametric instabilities by UHF/VHF backscatter radars in HF heating experiments is limited to HFPLs and HFILs, which do not represent the spectra of the plasma waves excited by parametric instabilities. Moreover, UHF/VHF radars do not detect magnetic field-aligned waves (i.e., wavevectors are closely perpendicular to the magnetic field), such as the upper and lower hybrid waves, because those waves do not backscatter radar pulses. On the other hand, the O-mode heating wave propagates through the upper hybrid resonance layer before reaching the reflection height in an over-dense ionosphere; thus, Langmuir and upper hybrid parametric instabilities can be excited simultaneously by a HF heating wave [12,13]. Langmuir and upper hybrid waves exert ponderomotive and thermal pressure forces on electrons [14,15], which induce nonlinear feedback on the waves. Thus, intense Langmuir and upper hybrid waves evolve to nonlinear waves, which have periodic and solitary envelopes; solitary waves press out local plasma in a self-consistent way to form cavitons [16,17].
A digisonde [18] was used to explore Langmuir and upper hybrid parametric instabilities excited in the HF heating experiments. The experiment was conducted with 2 min on and 2 min off on 20 November 2009 from 21:00 to 23:04 UT. In the on period, the polarization of the heating wave was switched alternately with O mode and X mode. Echo spread with bump(s) occurs only after O mode heating. The spread was fading away in the subsequent off period and eliminated further by the X-mode heater. In essence, it was 2 min on and 6 min off. The Sun was above the HAARP horizon for the entire experiment period. Therefore, there was no precondition on the background plasma for each O mode heating period. The experimental observations of exciting Langmuir and upper hybrid waves were manifested by bumps in the virtual spread of the ionogram trace, which are located closely below the HF heater frequency and the upper hybrid resonance frequency [12,13]. The time development of bumps in the ionogram trace, which signifies the competition between Langmuir PDI and upper hybrid PDI, was observed [12,13]. This is demonstrated in a sequence of 15 ionograms, presented in Figure 2a  The results show the time change of the virtual height spread as well as the development of bumps next to the HF reflection layer and upper hybrid resonance layer. In Figure 2a-j, bumps, indicated by arrows, appear simultaneously in the regions below 3.2 MHz (i.e., HF reflection height) and 2.88 MHz (upper hybrid resonance height), respectively. Based on the locations of these speculative bumps, one may infer that Langmuir PDI and upper hybrid PDI were excited simultaneously by the HF heater. In the early experimental period (from 21:10 to 21:52 UT), when the natural virtual height bump at E-F1 transition layer was still relatively small, the ionograms in Figure 2a-f show that the Langmuir bump was rising while the upper hybrid bump was dropping. As the experiment proceeded, the natural virtual height bump at E-F1 transition layer was rising; it seems to favor upper hybrid PDI in competing with Langmuir PDI. As shown in Figure 2g-j, the Langmuir bump was weakening and shifting down slightly in frequency while the upper hybrid bump was strengthening. As the experiment proceeded further, the natural virtual height bump at the E-F1 transition layer rose to become very high as shown in Figure 2k-o, and the Langmuir bump was suppressed completely, leaving the strengthened upper hybrid bump to stand out.
The development infers that the Langmuir PDI was suppressed by the upper hybrid PDI. The likely process is the pump depletion. First, the density ledge (cusp) at foE (the peak plasma frequency of the ionospheric E region) increased the loss of the HF heater in the E-F1 transition region; next, the heating power was further drained in the upper hybrid resonance region before the heating wave reaches the HF reflection height.
The impact of the natural cusp at foE on the competition between Langmuir PDI and upper hybrid PDI is further evidenced by the ionograms presented in Figure 3. The ionograms in Figure 3a,b were acquired on 16 November 2009 at 30 s before the heating experiment started and at the beginning of the moment when the O mode heater turned off after on for 2 min, respectively. As shown, the heating induced considerable virtual height spread as well as a bump located slightly below the HF reflection height. The ionograms in Figure 3c,d were acquired, on 20 November 2009, at 90 s after the heating experiment ended and at the beginning of the moment when the O mode heater turned off after on for 2 min, respectively. As shown, the heating also induced considerable virtual height spread as well as a bump but located slightly below the upper hybrid resonance layer. Ionograms in Figure 3a,c represent the background conditions in the two heating periods  The ionograms in Figure 3 indicate that the natural cusp was small on 16 November, thus, Langmuir PDI prevailed. On the other hand, on 20 November when the natural cusp was large, Langmuir PDI was suppressed and upper hybrid PDI prevailed. Because the density cusp introduces anomalous loss of the HF heater, the upper hybrid PDI was not as strongly excited as the Langmuir PDI excited on 16 November. Consequently, the upper hybrid bump in Figure 3d is much smaller than the Langmuir bump in Figure 3b.
Monostatic UHF radar only detect backscattered lines, such as HFPLs and HFILs; on the other hand, parametric instabilities excite plasma waves with broad spectra, and their evolution to the nonlinear states cannot be revealed by the UHF radar. Moreover, due to the field-aligned nature of the upper hybrid waves, these waves also cannot be detected directly by the UHF radar. The distinctive virtual height bumps in Figures 2 and 3 suggest that the digisonde can be a key diagnostic instrument to explore the nonlinear evolution of Langmuir waves and upper hybrid waves excited in the HF heating experiment.

Airglow Enhancement
Superthermal electron flux excites atomic oxygen via collision to produce airglow, which is emissions mainly at 630 nm and 557.7 nm. The minimum electron energies to stimulate airglow at 630 nm and 557.7 nm are 3.1 eV and 5.4 eV, respectively. Airglow enhancements were observed in the O-mode HF heating experiments. Electron plasma waves excited by parametric instabilities accelerate electrons through Doppler-shifted resonance interaction (quasi-linear diffusion) to reach the energy thresholds [6] for stimulating those airglows.
The enhancement of airglow at 777.4 nm was also observed as the O-mode HF heating wave was transmitted near the second and third harmonic of the electron gyro frequency [19,20]. The minimum electron energy to excite 777.4 nm emissions is 10.7 eV [19], which well exceeds the super thermal energy range and the capacity of quasi-linear diffusion.
Harmonic cyclotron resonances with electrons require the aid of finite Larmour radius effect [21], which is proportional to kv e /Ω e , where k and v e , and Ω e are the transverse wavenumber and electron speed, and the electron gyrofrequency. The finite Larmour radius effect works to shift down the wave frequency to the fundamental cyclotron resonance frequency for Doppler-shifted cyclotron resonance interaction, as well as to provide a positive feedback to the interaction.
Thus, heating at the harmonic cyclotron resonances directly by the HF heating wave is not effective because the wavenumber of the HF heating wave is small. On the other hand, an O-mode HF heating wave can excite parametric instabilities (upper-hybrid OTSI and PDI) in and below the upper hybrid resonance layer to produce short scale upper hybrid waves, which can interact effectively with electrons at cyclotron harmonic resonance. In Section 5, the theory of harmonic cyclotron resonance interaction procedure will be elaborated, and numerical results will be presented to show the generation of very-highenergy electrons by the upper hybrid waves [22,23].

Energetic Electron Flux
In the O-mode HF heating experiments, electron fluxes in the energy range from 10 to 25 eV were detected in situ by a probe in rocket [24] as well as on the ground by UHF radar inferred by an ultra-upshifted frequency band [25] and by the CCD camera inferred by the airglow enhancement at 777.4 nm [19].
As described in Section 2.1, the spectral lines in the radar returns have a fixed wavelength equal to half of the radar wavelength; thus, an ultra-upshifted frequency band signifies that the radar scatterers are moving at ultra-high phase speeds. These high-speed electron waves, which are not plasma modes, are virtually formed by the wakefields of ultra-high energy electron bunches.
Plasma waves excited by the parametric instabilities are likely to be responsible for the electron acceleration to such a high energy level. Langmuir waves implement quasi-linear diffusion through a Doppler-shifted resonance interaction to produce superthermal electrons. Upper hybrid waves further accelerate super-thermal electrons to ultra-energy level via Doppler-shifted harmonic cyclotron resonance interaction, which will be demonstrated by the numerical simulations of Section 5.

Artificial Ionization Layers (AILs)
Optical emissions observed by Pedersen et al. [26,27] indicated that new ionization layers were produced by the O-mode HF heating wave at frequency (~2.85 MHz) set near double the electron gyro frequency. Digisonde [18] echoes, shown in Figure 4a, provided evidence to support the generation of artificial ionization layers, indicated by an arrow. The generation and development of AILs during the heating experiment are seen by a sequence of seven ionograms recorded from 01:59:40 to 02:03:00, as presented in Figure 4b. Ionization energy of the atomic oxygen "O" is 13.6 eV. Experiments [28][29][30] were further conducted during twilight and early evening hours in Alaska local time, when the photoionization was weak. The wave-ionosphere interaction occurs in the region around 230 to 250 km, where the O + ions are dominant. Experimental results show that the enhanced optical emissions (inferring ionization enhancement) descend in the background F-region ionosphere [31] and relatively thin artificial ionization structures, seen directly in the ionograms of the digisonde, are emerging from the lower F-region ionosphere and descending ( Figure 4b) to settle at the base of the F-region (Figure 4a).
In later experiments with the heating wave frequencies set at 4.34 MHz and at 5.8 MHz, which are around the third and fourth harmonic of the electron gyro frequency, digisonde ionograms show that AILs also emerge from the base of the ambient F region as relatively thin layers, like those formed with the heating wave frequency set near the second harmonic of the electron gyro frequency. Experiments found that AILs were preferred to be generated with the heating wave frequency tuned either slightly above a harmonic of the electron gyro frequency or far below the second harmonic gyro frequency.
The power transfer from wave to the electron moving along the magnetic field, i.e., P = −e E z v ez , depends on the phase of E z which varies with time due to the mismatch frequency ∆ω = 0. Although the phase is distributed randomly from 0 to 2π, there are electrons of appropriate phases that will gain energy from the wave, i.e., P > 0 and v ez increases. In the case of ∆ω > 0, ∆ω decreases to give a positive feedback to enhance waveelectron interaction. When ∆ω drops to a negative value, the feedback of the interaction becomes negative; thus, a sufficient large ∆ω 0 is necessary to give adequate positive feedback interaction period for the generation of energetic electrons. On the other hand, if ∆ω 0 is too large, the initial interaction will be too weak, and the available interaction spatial length will be too short to generate energetic electrons. When the heating wave frequency is far below the second harmonic of the electron gyro frequency, the excited Langmuir waves can resonantly interact with electrons through Doppler-shifted cyclotron resonance; at fundamental cyclotron resonance, the interaction does not need the aid of the finite Larmour radius effect, and positive feedback to the interaction can still be sufficient at large frequency mismatch.
The numerical results to be presented in Section 5 show that the parametrically excited upper hybrid waves can energize electrons to exceed 13.6 eV, the ionization energy of the atomic oxygen "O", through the harmonic cyclotron resonance interaction [22,23]. Moreover, the parametrically excited upper hybrid waves have a frequency band, which it makes accessible for the match of cyclotron harmonic resonance over an altitude region. Thus, electrons, while moving downward, can be accelerated continuously along a slightly increasing geomagnetic field, through the Doppler-shifted harmonic cyclotron resonance interaction with a band of spatially distributed upper hybrid waves. It results in a major ionization occurrence at the bottom of the F region. This explains the descending feature of the enhanced optical emissions in the evolution of AILs and the emergence of AILs as relatively thin artificial ionization structures at the bottom of the F region.

Artificial Spread-F
Digisonde is HF radar, which probes the electron density distribution in the bottom side of the ionosphere. O-mode and X-mode sounding pulses with carrier frequency swept from 1 to 10 MHz are transmitted for sounding echoes, which are then recorded in an ionogram. An O-mode or a X-mode sounding echo represents the backscatter of a corresponding sounding pulse from a layer of the ionosphere, where the electron plasma density N(h) of the layer matches the O-mode cutoff density N cO = (f/9000) 2 cm −3 or X-mode cutoff density N cX = f(f − f e )/(9000) 2 cm −3 set by the carrier frequency f of the respective sounding pulse, where h and f e are the virtual height of the layer and the electron gyrofrequency. The virtual height h is determined by the time delay τ of the echo to be h = cτ/2, where c is the speed of light in free space. At HAARP site, N(h) can be converted to a true height profile N(h, ) by a profile conversion (NHPC) algorithm [32], which is available in the software program SAO Explorer [33].
The digisonde radiates at a large cone angle, each sounding pulse can be decomposed into many rays, which have different ray trajectories, and only backscattered rays can return to the digisonde and are recorded as the ionogram echoes.
In the unperturbed ionosphere, only a few rays, which are close to the vertical transmission, are backscattered. The virtual height traces of the sounding echoes in the ionogram have narrow virtual height spreads. When the ionosphere is perturbed, it also affects ray trajectories. In the presence of large-scale field-aligned density irregularities (FAIs with scale lengths of a few hundreds of meters to kilometers), ray trajectories [34] can be significantly modified as elaborated in Section 3.
In general, multiple incident rays (at different oblique angles) are backscattered to the digisonde receiver to produce multiple sounding echoes at the same radar frequency but at different return times [35], resulting in the spread of the virtual height traces. In the naturally perturbed ionosphere, the spread usually appears in the frequency band corresponding to the F region of the ionosphere and is thus termed "Spread-F".
In HF heating experiments conducted on November 13, 2009 from 21:00 to 23:00 UTC (12:00 to 14:00 local time), spread-F induced by the O-mode HF heater at 5.75 MHz (slightly above the fourth harmonic of the electron gyrofrequency) has also been observed and is termed "artificial spread-F". This is demonstrated in Figure 5 comparing a pair of ionograms without ( Figure 5a) and with heating effect (Figure 5b) [34]. As shown, the artificial spread-F extends from 2.4 MHz (foF1~2.38 MHz) to foF2~6.4 MHz; on the other hand, the parametric instabilities excited by the O-mode heating wave of 5.75 MHz occur only locally near the HF reflection height (at 5.75 MHz) and around upper hybrid resonance region (at 5.57 MHz). Although the large-scale FAIs could be generated through the filamentation of the HF heating wave [36], the FAIs responsible for the virtual height spread appearing in Figure 5b are likely generated by the thermal instability [37] driven by the downward and upward heat flow from the heat sources (converted from the plasma waves excited by the parametric instabilities) located near the HF reflection region and near the upper hybrid resonance region.

Ionization Enhancement
In the 20 November 2009 experiment, although the O-mode heating wave frequency of 3.2 MHz was not close to an electron harmonic cyclotron resonance frequency, artificial ionization enhancement was observed together with artificial spread-F [34]. Experiments were conducted around local solar noon when the photoionization was strong and the wave-electron interaction occurred in the lower F region (<180 km) of the ionosphere, where the electron-ion effective recombination coefficient depends strongly on the electron temperature T e [38]. Anomalous electron heating through parametric instabilities and thermal diffusion it reduces the recombination coefficient to change the balance between the photoionization and the recombination loss over a large region. As a result, there is an electron density enhancement in the heated region below~180 km. This is demonstrated in Figure 6, in which electron density distributions at times (a) 21:29:30UT, after the O-mode heater off more than 90 s and (b) 21:28:00UT, off within 10 s, are presented for comparison. As shown, the distribution at time (b) (i.e., blue curve) has a higher density in the entire modified region than that at time (a) (i.e., red curve). The percentage of the electron density increase exceeds 10% over a region extended in height (> 30 Km) from below to above the HF reflection height.

Artificial Cusp
A virtual height bump usually occurs at frequencies near foE in the presence of an F1 layer, as shown in Figure 5a; the true height profile (the solid line plot) indicates that there is a density ledge (cusp) at foE, which retards the propagation of sounding signals with frequencies f close to foE, causing spread of the sounding echoes in virtual height in the form of a distinct bump. Spread-F normally extends over a frequency band in the ionogram and is recognized to be caused by the large-scale FAIs in the plasma. However, localized anomalous echo spread appearing as a bump in the ionogram trace has also been observed in the O-mode HF heating experiments conducted on 16 November 2009. This is exemplified in Figure 7a Figure 7b); the similarity suggests that a heater-stimulated ionization ledge (cusp) appears near the upper hybrid resonance layer, which causes the virtual height spread in the form of a bump peaked slightly above the upper hybrid resonance layer [12]. This ledge can be created by the thermal pressure force arising from the Langmuir waves (i.e., the localized Langmuir waves form a caviton), which are excited parametrically by the HF heating wave. The matching height of the Langmuir wave drops due to the increase of the effective electron temperature, which is proportional to the total spectral intensity of the Langmuir waves.
The experiment was conducted around local solar noon when the photoionization was strong and the wave-electron interaction occurred in the lower F region (<180 km of the true height) of the ionosphere, where the electron-ion effective recombination coefficient depends strongly on the electron temperature T e [38]. The heat diffusion from the Langmuir PDI region downward and upward excites thermal instability to generate large-scale FAIs [37] as well as to cause ionization enhancement over a large region. An ionization enhancement, like that presented in Figure 6, was observed. On the other hand, the combined ionogram presented in Figure 7a shows that the heating wave induces a downward spread of the ionogram trace in the region from 2.1 to 2.8 MHz, different from that shown in Figure 5b in which the enhanced Spread-F appears with an upward-spread of the virtual height in the ionogram trace. It is likely the spread in Figure 7a is caused by a different mechanism; the down-going Langmuir waves evolve into nonlinear periodic plasma waves propagating along the magnetic field, which backscatter incoherently the sounding pulses to give rise to downward spread of the ionogram trace.
In sum, these observed phenomena are ascribed to various linear and nonlinear wave-plasma and wave-wave interactions instigated by the Langmuir waves, upper hybrid waves, and density irregularities, generated by the HF heating waves. Theoretical formulation and analysis of parametric instabilities and instigated interaction processes are presented in Section 4 as the theoretical basis of HF heating induced phenomena.

Ray Tracing
One approach to characterize wave propagation in the ionosphere is "ray tracing". It is applicable when the wave front extends uniformly over several wavelengths and the inhomogeneity scale lengths of the medium are large in comparison to the wavelength, and particularly, the wave has a dominant frequency. In this situation, the wave may be treated as a ray and its trajectory is tracked to explore wave propagation and the ionospheric plasma. The ionosphere is stratified into layers and the Snell's law of refraction is applied at the interface of two adjacent layers to setup the ray trajectory equation. This is illustrated by applying the arrangement shown in Figure 8, in which an inhomogeneous medium with refractive index n(z) is approximated by a series of plane slabs of thickness ∆z, where ∆z → 0 and each slab i has a uniform refractive index n i . As shown in Figure 8, a ray is incident from a uniform medium (e.g., below the ionospheric plasma) of refractive index n 0 at an angle, θ 0 , with respect to the vertical axis into this inhomogeneous medium; the Snell's law n j sin θ j = n j+1 sin θ j+1 is applied at the interface of two adjacent slabs j and j+1 to determine the moving direction of the ray; n 0 < n 1 < n 2 . . . is assumed in the plot (however, in the case of propagating upward into the ionospheric plasma, it is n 0 > n 1 > n 2 . . . .).
The path of the ray in each slab is a straight line, and one can identify the relations: ∆z/∆x i = cotθ i and n i sin θ i = n 0 sin θ 0 , where ∆x i is the horizontal displacement of the ray after transit through slab i. In the limit of ∆z → 0, these two relations become dz/dx = cotθ and n sin θ = n 0 sin θ 0 , and a trajectory equation is derived to be If n 2 can be modelled by a second-order polynomial, i.e., n 2 (z) = a + bz + cz 2 , (1) can be integrated analytically; otherwise, it is integrated numerically.

General Formulation of Ray Trajectory Equations
Consider a general case that wave propagation is governed by a dispersion equation with the generic form [39] G(k, ω; r, t) = 0.
We now introduce a generic variable "τ" and take a total τ derivative on (2), it yields The four terms in (3) are arranged into two groups to be ∂G ∂t where the first group of terms is related to the time variation of the media, while the second group related to the spatial variation of the media. Because (4) is obtained from a general approach and the spatial variation and temporal variation are separable, its general solution requires that the following two relations be satisfied simultaneously: and ∂G ∂t where dr dτ = dk dτ and dt dτ = dω dτ because r and k and t and ω are independent to each other; thus (5) and (6) deduce to the following relations In principle, (2) can be solved to obtain, for instance, ω = ω(k; r, t); in other words, there are only three independent variables in (2). Choose r, t, and k to be independent variables and ω a dependent variable, i.e., G(k, ω; r, t) = G(k, r, t; ω(k, r, t)) then take partial derivatives of (2) with respective to the three independent variables k, r, and t, respectively, it obtains the following three relationships ∂G ∂t With the aid of (8a) and (7d), (7a) becomes With the aid of (8b) and (7d), (7b) becomes With the aid of (8c) and (7d), (7c) becomes Equations (9) to (11), subject to a set of initial conditions, determine the ray trajectory in the phase space (i.e., r-k space). In essence, these are the Hamilton's equations of motion with ω and k to be the Hamiltonian and momentum of the ray.

Ray Trajectories of Sounding Pulses-Spread-F
This ray tracing technique is applied to explain the virtual height spread (e.g., Figure 5b) of the ionogram trace observed in the HF heating experiments. Digisonde (ionosonde), an HF radar, is a remote-sensing device for monitoring the plasma density profile of the bottom-side ionosphere. It transmits sounding pulses upward and records the return echoes in an ionogram as the virtual height trace. Digisonde radiates at a large cone angle, each sounding pulse can be decomposed into many rays, which have different ray trajectories, and only backscattered rays can return to the digisonde and are recorded as the ionogram echoes. In the unperturbed ionosphere, only a few rays, which are close to the vertical transmission, are backscattered. Thus, the virtual height traces of the sounding echoes in the ionogram have narrow virtual height spreads. However, when the ionosphere is perturbed in the presence of large-scale field-aligned density irregularities, ray trajectories can be significantly modified as elaborated in the following.
A carton of large-scale FAIs appearing in the F region of the ionosphere is presented in Figure 9a; x and z axes are in the horizontal and vertical directions. A sounding pulse of the digisonde is decomposed into a bunch of rays as also shown in Figure 9a. Only plasma affects ray trajectory, thus, from the ground to the reference layer, z = 0, rays are assumed to propagate in a free space, where the plasma density is still low and FAIs do not extend to. Ray trajectories from z = 0 up to their reflection layers and returning to z = 0 are determined via ray tracing technique. In this region, the background plasma has a density distribution N(z) = N 0 1 + z L , where N 0 and L are the plasma density at the reference layer, z = 0, and the inhomogeneity scale length of the background plasma, respectively. FAIs are represented by a single sinusoidal function with a spatial period d and an amplitude ∆N, i.e., N I (x, z) = ∆N sin κ(x sin θ d + z cos θ d ), where κ = 2π/d and θ d is the magnetic dip angle; θ d = 50 • , 75.6 • , and 78 • at Arecibo, HAARP, and EISCAT, respectively. In the model, plasma is uniform in the y direction, only two-dimensional trajectories on the x-z plane are considered in the numerical analysis; (9) and (10) become and Because the ionospheric plasma and FAIs can be stationary in the propagation period of a ray, ∂ω/∂t = 0, (11) leads to dω/dt = 0; thus ω = ω(k x , k z , x, z) is a constant of motion. With the aid of the dispersion relation ω = ω(k x , k z , x, z), right-hand side terms of the Equations in (12) and (13) can be expressed explicitly.
Digisonde transmits linearly polarized sounding pulses, which can be decomposed to RH-and LH-polarizations representing O-mode and X-mode sounding pulses. The dispersion equations for obliquely propagating O-mode and X-mode heating waves are given approximately by and where the parallel and perpendicular (to the magnetic field) components, k and k ⊥ of the wavevector are related to the x and z components, k x and k z of the wavevector by k = −k x cos θ d + k z sin θ d and k ⊥ = k x sin θ d + k z cos θ d ; ω 2 p = ω 2 p0 [1 + z L + ∆N N 0 sin κ(x sin θ d + z cos θ d )] and ω 2 p0 = N 0 e 2 m 0 . With the aid of (14) and (15), right hand side terms of the equations in (12) and (13) can be expressed explicitly for the O-mode and X-mode sounding rays, respectively.
In the numerical calculation, FAIs are represented by a single sinusoidal function with a spatial period of 200 m and an amplitude ∆N; the background plasma density has a linear scale length L = 30 km. Figure 9b demonstrates the change of the trajectory of a vertically incident O-mode ray of 2.8 MHz in the presence of FAIs with the amplitude ∆N increased from 0 to 0.1 N 0 , where N 0 is the plasma density at a reference layer, z = 0, having a plasma frequency ω p0 = 2.26 MHz (i.e., 3.2/ √ 2 MHz). In the absence of FAIs (∆N = 0), the vertically incident ray is backscattered. However, in the presence of FAIs (∆N/N 0 = 0.01 to 0.1), the vertically incident ray is not backscattered anymore.
A similar demonstration for the vertically incident x-mode ray of 2.8 MHz is presented in Figure 9c. As shown, the vertically incident ray is backscattered also only in the absence of FAIs (∆N = 0). In fact, the vertically incident ray is ducted by the FAIs (for ∆N/N 0 ≥ 0.03), the return signal propagates along the geomagnetic field, rather than propagating vertically downward. The x and z axes in Figure 9b,c are normalized to k 0 = ω 0 /c, where f 0 = ω 0 /2π is the frequency of the sounding pulse. For example, at 2.8 MHz, 1/k 0 = c/ω 0 = 17.05232 m; thus, the (x, z) axes for 2.8 MHz ray are normalized to 17.05232 m. The normalization is inversely proportional to the frequency, thereby it decreases to 13.64 m for 3.5 MHz ray.
Therefore, in the presence of FAIs, the traces in an ionogram are not contributed by the returns of the vertically incident rays. In order to achieve backscatter, the incident direction of the ray at the reference layer, z = 0, must have an off vertical angle. The trajectories of backscatter rays in the presence of 15% FAIs for both O-and X-mode sounding signals from 3.5 to 4 MHz are determined. 15% FAIs are relatively large for the general situation; it is adopted in the simulations for a better contrast. The results are shown in Figure 10, the backscatter trajectories in the absence of FAIs are also plotted for comparison. As shown in Figure 10a, the trajectories of the O-mode backscatter rays [3.5 (red), 3.6 (green), 3.7 (black), 3.8 (pink), 3.9 (blue), 4 (brown) MHz] are modified significantly by the FAI from those (light pink) in the absence of FAI, however the incident off-vertical angles are very small (< 1 • ); on the other hand, as shown in Figure 10b, the off-vertical angles of the x-mode backscatter rays in the presence of FAI are very large, and this off-vertical angle (> 15 • ) increases as the frequency of the sounding signal decreases. The off-vertical angle also increases with the increase of the FAI amplitude. In general, multiple incident rays (at different oblique angles) are backscattered to the digisonde receiver to produce multiple sounding echoes at the same radar frequency but at different return times [35], resulting in the spread of the virtual height traces. In the HF heating experiments, large-scale FAIs are generated by the O-mode heating waves, causing a spread of the virtual height traces like the natural Spread-F; the heating-wave-induced Spread-F is termed "Artificial Spread-F".

Parametric Instabilities Excited in HF Heating Experiments
Plasma supports high-frequency and low-frequency modes; in the absence of external sources, these modes present in plasma as thermal fluctuations. A large-amplitude, highfrequency wave, E p ω 0 , k p , (either electromagnetic (EM) or electrostatic (ES)), can be a pump wave, which excites high-and low-frequency plasma waves concurrently [39]. The electric field of the high frequency pump wave sets up a quiver motion in the electron plasma at velocity v eq (t). On the other hand, electrons and ions oscillate together to maintain quasi-neutrality, i.e., n es (ω s , k s ) ∼ = n is (ω s , k s ) = n s (ω s , k s ), in the low frequency plasma wave field; it facilitates the low-frequency plasma wave to buildup plasma density perturbation n s (ω s , k s ). A coupling of the pump wave and the lowfrequency plasma wave induces a high-frequency space charge current density, −en s v eq , in the electron plasma. Such a current drives beat waves, E(ω, k) and E ω , k , with the wavevectors and frequencies satisfying the matching conditions: − ω s and k p = k + k s = k − k s These beat waves, in turn, also couple with the pump wave to induce in the electron plasma a low-frequency nonlinear force, which has the same frequency and wavevector as the density perturbation, n s (ω s , k s ), to drive this density perturbation. When the couplings generate large enough positive feedback to overcome linear losses of the coupled waves, instability is excited to grow the coupled plasma waves exponentially in the expense of the pump wave energy. This is called "parametric instability"; a pump wave, E p ω 0 , k p , decays to two sidebands, E(ω, k) and E ω , k , through coupling to a low-frequency decay mode, n s (ω s , k s ). Parametric coupling is imposed by the frequency and wavevector matching conditions; further, the instability requires the pump electric field intensity to exceed a threshold. When the decay mode, n s (ω s , k s ), has a finite oscillation frequency, two sidebands cannot satisfy the same dispersion relation concurrently. The frequencyupshifted sideband, E ω , k , has large frequency discrepancy from the plasma mode and can be disregarded to simplify the parametric coupling to be three-wave interaction.
Parametric excitation of Langmuir/upper hybrid waves, φ(ω, k), and low-frequency plasma waves, n s (ω s , k s ), by the electromagnetic or Langmuir/upper hybrid pump waves, E p ω 0 , k p , in the bottom-side of the ionosphere, are explored in the following, where E p, φ, and n s denote electric field of a pump wave, electrostatic potential of the Langmuir/upper hybrid sideband, and density perturbation of a low frequency decay mode, respectively. Langmuir waves can have large oblique propagation angles (with respect to the geomagnetic field B 0 = −ẑB 0 ), upper hybrid waves have near 90 • propagation angles, and low-frequency plasma waves include ion acoustic/lower hybrid waves as well as purely growing modes.
The most likely parametric instabilities excited directly by a HF heating wave in the bottom-side of the ionosphere include 1) oscillating two-stream instability (OTSI), and 2) parametric decay instability (PDI), in both mid-latitude and high-latitude ionosphere [4,40,41]. The sidebands are Langmuir waves and upper hybrid waves. The instabilities with Langmuir waves as sidebands must compete with those excited in the upper hybrid resonance region, where the upper hybrid waves are the sidebands of the instabilities. The wavenumber k 0 of the HF heating wave is much smaller than the wavenumbers of the electrostatic sidebands and decay modes, thus a dipole pump, i.e., k 0 = 0, is generally assumed. As the excited Langmuir/upper hybrid sidebands grow to large amplitudes, they also become pump waves to excite parametric instabilities [42,43].
When the theory is used to explain the observations, it is noted that HFPLs and HFILs are monitored by UHF and VHF backscatter radars; those lines are associated with plasma waves propagating oblique to the geomagnetic field at an angle conjugate to the magnetic dip angle. Moreover, upper hybrid waves, lower hybrid waves, and field-aligned purely growing modes cannot be detected.
The coupled mode equation for the Langmuir/upper hybrid sideband φ(ω, k) is derived from the electron continuity and momentum equations, together with the Poisson's equation, to be [44] where represents a filter, which only keeps terms with the same phase function as that of the function φ on the left-hand side; Ω e = eB 0 /m, ω p = n 0 e 2 /m 0 1/2 , v te = (T e /m) 1/2 , and m, are the electron cyclotron frequency, plasma frequency, thermal speed, and mass, respectively. (16) is derived from the fluid equations, however the kinetic effect of electron Landau damping is included phenomenologically in the collision damping rate, i.e., ν eh = ν e + ν eL , which is the effective electron collision frequency, where ν e = ν en + ν ei , a sum of the electron-neutral elastic collision frequency, ν en , and the electron-ion and ν eL = (π/2) 1/2 ω 2 ω 2 p /k z k 2 v 3 te exp −ω 2 /2k 2 z v 2 te , a phenomenological collision frequency to represent electron Landau damping.
The low-frequency wave fields move both electrons and ions, hence the formulation of the coupled mode equation for the low-frequency decay mode, n s (ω s , k s ), includes electron and ion fluid equations. Because electrons and ions move together in the low frequency field, the formulation is simplified by introducing quasi-neutral condition: n si ∼ = n se = n s . The ion fluid equations are similar to the electrons except that the electron mass m is changed to the ion mass M, the charge −e changed to e, and the collision terms in the electron and ion fluid equations are modeled differently by ν ei (v e − v i ) + (ν en + ν eL ) v e and ν ie (v i − v e ) + ν i v i , respectively, where ν i = ν in + ν iL , ν in is the ion-neutral collision frequency, ν iL ∼ = (π/2) 1/2 ω 2 s /k z V s (T e /T i ) 3/2 exp −ω 2 s /2k 2 z v 2 ti represents ion Landau damping on the low frequency decay mode, and V s = (T e /M) 1/2 .
Neglecting the electron inertial term and the ion convective term in the momentum equations, these two momentum equations are combined into a one fluid equation, where the relations ν ie = m M ν ei and mΩ e = MΩ i are applied. This equation is then combined with the ion continuity equation to derive a coupled mode equation, for the low-frequency (ion acoustic/lower hybrid) decay mode [45] where Ω i is the ion cyclotron frequency, C s = [(T e + 3T i )/M] 1/2 the ion acoustic speed, and M the ion (O + ) mass; (m/M)ν e << ν i is applied; ν es = ν e + ν eLs , and ν eLs = (π/2) 1/2 ω 2 s ω 2 p /k z k 2 v 3 te exp −ω 2 s /2k 2 z v 2 te . The coupling terms a p = v e ·∇v e and J B = n e v e arise from plasma nonlinearities. The linear responses of the electron density and velocity to the total high frequency wave fields are used to present these coupling terms.
In the following, the coupled mode algebraic Equations (21) and (22) are employed to analyze parametric instabilities excited in the ionospheric HF heating experiments.

OTSI and PDI near the HF Reflection Height
The RH circularly polarized HF heating wave propagates to the region near the reflection height, and it converts to the O-mode with the electric field E p ∼ =ẑE 0z cos ω 0 t, where E 0z ∼ = √ 2E 0 .

OTSI-Excitation of Langmuir
Sidebands φ ± (ω, ±k) Together with Purely Growing Density Striations n s (ω s = iγ s , k s ) by the HF Heating Wave The parallel (to the magnetic field) component k sz of the wavevector of the small-scale purely growing mode n s (ω s = iγ s , k s ) is not negligibly small, and the density striations are mainly driven by the parallel component of the ponderomotive force induced by the high frequency wave fields. Apply the condition Ω e Ω i ξ~Ω 2 e k 2 z /k 2 >> k 2 C 2 s to simplify (22) and set ω s = iγ s and ω = ω 0 + iγ s , (21) and (22) are combined to be (23), the threshold field is obtained to be It shows that the threshold field of OTSI varies with the propagation angle θ and wavelength λ of the Langmuir sidebands as well as the location of excitation (i.e., ∆ω 2 ). The threshold field (24) has the minimum when the (k,θ) lines are excited in a preferential layer at altitude h(k,θ), where ω 2 When the instability is excited at the matching height h of its Langmuir sideband (k, θ), i.e., ω r = ω kθ = ω 2 p + 3k 2 v 2 te + Ω 2 e sin 2 θ 1/2 and ω sr = kC s , the minimum threshold field is obtained to be Equations (25) and (27) show that the layers of exciting OTSI and PDI move downward and the threshold fields increase as the oblique propagation angles θ of the OTSI and PDI lines increase. When the heating wave field E 0 is large, the Langmuir sidebands excited by the OTSI and PDI will have angular (θ) and spectral (k) spread lines, distributed in altitude layers, where are higher and narrower for the OTSI excitation layer than that for the PDI layer.

OTSI-HF Heating Wave Decaying to Upper Hybrid Sidebands and Field-Aligned Density Irregularities
The O-mode HF heating wave can access the upper hybrid resonance layer, which is located below the O-mode HF reflection height. In high-latitude regions, such as at Tromso, Norway and Gakona, Alaska, RH circularly polarized heating wave can be transmitted along the geomagnetic field. In the region near the upper hybrid resonance layer, the heating wave field still remain at RH circular polarization, i.e., E p = 1 2 x + iŷ E 0 exp(−iω 0 t) + c.c.; the heating wave decays to two upper hybrid sidebands φ + ω ± , k ± = ±x k and a field-aligned purely growing mode n s ω s = iγ s , k s = −xk , where ω ± = ω = ω 0 + iγ s .

PDI-Decay of HF Heating Wave to an Upper Hybrid Sideband φ(ω, k) and a Lower
Hybrid Decay Wave n s (ω s , k s ) A dispersion equation is derived, by combining (21) and (22) where ν i << ν es is assumed. The minimum threshold field E updi m of the instability excited at the matching height of its upper hybrid sideband φ(k, θ) is obtained to be where ν e = ν ei + ν en , ν es = ν e + ν eLs , and ν eLs = (π/2) ; in other words, the preference region of exciting upper hybrid PDI is below the upper hybrid resonance layer.

Height Separations of the Instability Layers in the Bottom-Side of the Ionosphere
Let h r /ω pr , h p /ω pp , and h u /ω pu be the HF reflection height, Langmuir PDI height, and upper hybrid resonance layer height, and the corresponding electron plasma frequency at each height, respectively. In the small region from the upper-hybrid resonance layer to the HF reflection layer, a linear increasing plasma density profile n = n 0 (1 − ∆h/L) is assumed, where n 0 is the electron density at h r (i.e., ω 2 pr = ω 2 0 ), L is the linear scale length, and ∆h is the distance below h r (i.e., ∆h = h r − h). With the aids of the HFPL dispersion relation and the definition of the upper-hybrid resonance frequency, this leads to ω 2 pr − ω 2 pp ∼ 12k 2 R v 2 te = ω 2 0 ∆h p /L and ω 2 pr − ω 2 pu ∼ Ω 2 e = ω 2 0 (∆h u /L); then, ∆h p = h r − h p and ∆h u = h r − h u are derived to be and where k R is the wavenumber of the UHF radar.

Impact of Double Resonances on Parametric Excitation of Upper Hybrid PDI
Double resonances is the situation where the upper-hybrid sideband of the upperhybrid PDI excited in the upper hybrid resonance region is also in resonance with electrons at the Doppler-shifted nth harmonic cyclotron resonance, i.e., ω ∼ ω uk = nΩ e + k z v te , where the upper-hybrid frequency ω uk = ω 2 0 + Ω 2 e + 3k 2 v 2 te 1/2 . The frequency matching condition ω 0 = ω + ω s of the upper-hybrid PDI, can be satisfied if ω 0 > nΩ e but close to nΩ e . The lower-hybrid wave frequency at HAARP is given to be f Lk A frequency sweep method can be applied to achieve a double resonance situation. It slowly sweeps the heating wave frequency from nΩ e to a frequency covering the lowerhybrid wave frequency f Lk and the Doppler-shifted frequency ∆f = k z v te /2π. The linear velocity response of the electrons in the upper-hybrid wave field E u = −ikφ + c.c. is given to be v u = − ke/m ω 2 − Ω 2 e ωx + iΩ eŷ φ + c.c.. At the Doppler-shifted cyclotron harmonic resonance, ω = nΩ e + k z v te , there is an additional resonance response, ascribed to the finite Larmour radius effect; it is derived to be v ur = i(ke/mωv e ) n 2 Λ n (β e )/β e ωx + iΩ eŷ φ + c.c., where β e = k 2 v 2 te /Ω 2 e , Λ n (β e ) = I n (β e ) exp(−β e ), and I n (β e ) is the modified Bessel function of the first kind. The total electron velocity response to the upper-hybrid wave becomes v ut = − ke/m n 2 − 1 Ω 2 e ] 1 − i n n 2 − 1 Ω e Λ n (β e )/β e ν e ωx + iΩ eŷ φ + c.c.; the electron ponderomotive force, driving the lower-hybrid waves, is then enhanced by g 2 = |v ut | 2 /|v u | 2 . It lowers the threshold field, given in (31), by a factor g = 1 + n n 2 − 1 Ω e Λ n (β e )/2β e ν e 2 1/2 and elevates the growth rate by g 2 ; thus, the upper-hybrid and lower-hybrid waves will be intensified as the finite Larmour radius factor β e = k 2 v 2 te /Ω 2 e is not too small. For β e ≥ 0.1, it requires that the wavelength λ u of the upper-hybrid wave is less than 0.3 m, i.e., λ u ≤ 0.3 m.

Instabilities under Double Resonance Situation
It was shown by Huang and Kuo [46] that when the heating wave frequency ω 0 is near the third harmonic electron cyclotron resonance frequency, i.e., ω 0 ∼ 3Ω e , upper-hybrid OTSI is suppressed by the parametric excitation of the electron Bernstein sidebands and small-scale FAIs. As a result, the source upper hybrid waves of the downshifted maximum (DM) lines of the stimulated electromagnetic emissions (SEEs) [47] is effectively suppressed; it explains the quenching of the DM feature in the SEE spectrum observed in the ionospheric heating experiments. However, the upper-hybrid PDI, which decays the HF pump wave into upper-hybrid sidebands and lower hybrid decay modes, can still be excited in the region below the upper-hybrid resonance region. As ω 0 < 3Ω e , upper-hybrid (OTSI and PDI) instabilities are excited in the upper-hybrid resonance region as well as in a large region underneath. In this case, upper-hybrid instabilities deplete the HF heating power considerably before the HF heating wave reaches the PDI layer, thus the Langmuir PDI is suppressed.
On the other hand, as the heating wave frequency ω 0 is ramping-up swept to exceed nΩ e (i.e., ω 0 > nΩ e ), the frequency matching condition, ω 0 = ω uk + ω Lk , of the upperhybrid PDI and the gyro-harmonic resonance condition ω uk = nΩ e + k z v te , can be satisfied simultaneously in the same altitude region, where ω uk = ω 2 uh + 3k 2 v 2 te 1/2 and ω 2 uh = ω 2 pu + Ω 2 e . Although double resonance suppresses upper-hybrid OTSI, it enhances the excitation of upper-hybrid PDI as discussed in Section 4.4; however, the upper-hybrid PDI zone is narrow; the instability will not be able to suppress the Langmuir PDI/OTSI. Langmuir PDI generated both up-going and down-going HFILs; but only the downgoing HFILs will continuously propagate downward. On the other hand, the strongly excited upper-hybrid sidebands propagate downward and upward as electrostatic pump waves [48] that can cascade to down-going Langmuir sidebands and up-going ion-acoustic waves as well as up-going Langmuir sidebands and down-going ion-acoustic waves over an additional altitude range below h = h r . In terms of the wavevector matching conditions for the generation of up-going/down-going HFILs, the cascade processes are illustrated in Figure 11a,b, respectively [49]. Upward propagating upper-hybrid waves stop at a height near the HF reflection layer (at h = h r , as shown in Figure 11b).
The Langmuir sidebands in the two cascade instabilities have different wavenumbers and the oblique angles to the magnetic field. The cascade process involving up-going Langmuir sidebands (Figure 11b) has a higher threshold field (due to larger oblique angle) and lower growth rate (due to smaller wavenumber) than the other one involving downgoing Langmuir sidebands (Figure 11a). Because the mode characteristic of the HFILs is not sensitive to the background plasma density variation, down-going HFILs mainly attribute to those generated via Langmuir PDI.

Harmonic Electron Cyclotron Resonance Interactions
The motion of an electron, starting at origin and interacting with the excited upper hybrid wave field, E(r, t) = 1 2π x k x +ẑk z /k E k cos (k x x + k z z − ωt) dk, in a uniform magnetic fieldẑB 0 , is considered, where |k z | < |k x | and E k is an even function of k; thus, the initial wave field amplitude at origin is given by E u0 = E(0, 0) = E k dk/2π. The trajectory of the electron is governed by the equations: dr(t)/dt = v(t) and dv(t)/dt = −eE(r, t)/m − Ωv(t) ×ẑ, where Ω = eB 0 /mc is the electron cyclotron frequency. The trajectory equations are integrated to be where At Doppler-shifted cyclotron resonance, ω = nΩ + k z v z1 , the trajectory would be that determined by the resonant (secular) terms together with small perturbations oscillating at high frequencies due to the non-resonant interactions. v z1 = (ω − nΩ)/k z is the harmonic cyclotron shifted wave phase velocity along the magnetic field and the resonant electrons propagate along the magnetic field in the same direction as that of the wave when ω > nΩ.

Second Harmonic Cyclotron Resonance Case (n = 2)
We keep only the J 1 term in (38) to (42) and J 1 (w)≈ w/2 for |w| < 1; (38) and (39) are converted to the differential equations where a 2 = (k x eE kx /4mΩ) cos[Φ + k z (v z0 − v z1 )t] dk/2π and b 2 = (k x eE kx /4mΩ) After the spectral integration, we obtain a 2 = γ 2 sin[k z (v z0 − v z1 )t + 2θ 0 ] and b 2 = −γ 2 cos[k z (v z0 − v z1 )t + 2θ 0 ], where γ 2 = k 2 m eE u0 v 0 /4mΩ 2 sin θ 0 ; the even function property of E k is applied and k m is close to the wavenumber of the spatial spectral peak of the upper hybrid wave. The value of γ 2 depends on the initial velocity (speed and phase angle) of the electron and the field intensity. The difference between the parallel velocities v z0 and v z1 causes detuning on the resonance, which results in the reduction of the acceleration. We consider the optimal case that v z0 = v z1 = (ω − nΩ)/k z ; thus a 2 and b 2 are constant parameters, and (43) and (44) are solved analytically. Equations (43) and (44) can be decoupled to become two second-order differential equations: Equations in (45) are solved with the initial conditions: The resonance electrons have the energy where 0 = mmv 2 ⊥0 /2 is the initial electron transverse kinetic energy. The acceleration time depends on the available interaction time τ, which is determined by the shortest Coulomb collision time τ e, and the transit time, τ T , governed by the size of the interaction region. Thus, τ = τ e τ T /(τ e + τ T ). Since v z is assumed to be small, τ T >> τ e , andτ ≈ τ e .
In the numerical evaluations, the parameter values used in the (1) second harmonic case are E u0 = 0.015 V/m, λ m = 0.5 m, and the interaction time τ = τ ei τ en /(τ ei + τ en ), where τ ei = 0.  It is realized that the geomagnetic field varies slightly with the altitude and it takes time and interaction length to energize electrons in the resonance interaction. Theoretical results [23] indicate that the excited upper hybrid waves have a frequency bandwidth. In the parametric coupling, the upper hybrid wave frequency is f uk = f 0 − f Lks , where f 0 is the heating wave frequency and the lower hybrid frequency f Lks = f LH ξ 1/2 . In HAARP, the lower hybrid resonance frequency is f LH ∼ (Ω e Ω i ) 1/2 /2π ∼ 8.3 kHz. For ξ = 1 to 10, f Lks ∼ 8.3 to 26.3 kHz. Hence, f uk has a bandwidth of 18 kHz, which covers the change of the ionospheric electron gyro frequency over an altitude range of over 10 km. This interaction length is crucial for producing energetic electron flux for ionizations.
The plasma frequency increases faster than the decrease of the gyro frequency in the upward direction. Thus, the resonance interaction region is extended only downward from the upper hybrid resonance layer and with the presence of the downward propagating upper hybrid waves. This explains why the artificial ionization layers descend during their development as observed in electron cyclotron harmonic resonance heating experiments.

Down-Going Upper Hybrid Waves and Artificial Ionization Layer (AIL) Location
Upward-propagating upper-hybrid waves reflect at a height near the HF reflection layer (h = h r ) to become down-propagating upper hybrid waves. These waves (ω uk , k =xk +ẑk z ) continuously propagate down from the upper-hybrid resonance layer, h = h u , along the geomagnetic field, to a height h, and it evolves to (ω uk , k 1 =xk 1 +ẑk 1z ); its wavevector changes the magnitude and direction (i.e., the oblique angle, θ 1 = tan −1 (k 1 /k 1z ), to the geomagnetic field) due to the dropping of the plasma density.
In the propagation, the wave frequency and the horizontal component of the wavevector do not change, i.e., ω uk (h u ) = ω uk (h) and kcosθ = k 1 sin(θ 1 − θ), where θ is the conjugate of the geomagnetic dip angle α~76 • , which set up two equations for k 1 andθ 1 to be These two equations are solved to obtain k 1z (h) = k 1 cos θ 1 = k 2 1 − k 2 cos 2 θ 1/2 cos θ − (k/2) sin 2θ, which increases with an increase of k 1 and with the decrease of θ as waves propagate down; upper-hybrid waves then convert linearly into oblique propagating Langmuir waves [49].
The process is presented in Figure 13, in which the blue curve, ω/ω p = 1 + Ω 2 e + 3k 2 v 2 te /ω 2 p 1/2 , and red curve, ω/ω p = ω 2 p (h a ) + Ω 2 e sin 2 θ a + 3k 2 v 2 te 1/2 /ω p are the dispersion curves of the upper-hybrid wave (k = k x , i.e., |k z /k x | << 1) and the Langmuir wave (k = k 1 ) at the parametric excitation layer (at h = h p ) and at the AIL (at h = h a ), respectively, where ω p = ω p h p and θ a is the oblique angle (to the geomagnetic field) of the Langmuir wave at AIL location (h = h a ).
As the upper-hybrid wave propagates downward, its frequency does not change, but the wavenumber increases and the wavevector inclines toward the magnetic field. As shown in Figure 13, the cross shifts horizontally to the right to the circle (o) on the red curve to become an obliquely propagating Langmuir wave. The downward change of the wave location can also be represented by rescaling the vertical axis by the ratio . The dot curve intersects with the blue curve at a point also marked with a cross (x), which then represents the new mode appearing in the AIL. It shows that the upper-hybrid wave converts to an obliquely propagating Langmuir wave with k L = k L x sin θ a +ẑ cos θ a . Because the horizontal component of the wavevector does not change, i.e., k L sin(θ a − θ) = k m cos θ, thus k L = k m cos θ/ sin(θ a − θ), where k m is the initial wavenumber of the upper-hybrid wave in a range (k 1 , k 2 ). The converted Langmuir wave further cascade if it still has a large amplitude.

Generation of Energetic Electrons
The excited upper-hybrid and lower-hybrid waves accelerate electrons to higher energy. As the heating wave frequency ramps up, the harmonic cyclotron resonance region moves downward because the electron cyclotron frequency increases with lowering the altitude; thus, the double resonance layer shifts downward for shorter wavelength upper hybrid waves; it extends the generation region of the upper-hybrid waves, and also enhances the harmonic cyclotron resonance interaction between the excited waves and the electrons because shorter wavelength upper-hybrid waves provide stronger finite Larmour radius effect, as discussed in Section 5.
Those electrons at Doppler shifted cyclotron harmonic resonance interaction with the excited waves (i.e., ω uk = nΩ e + k 1z v ez ), can be energized to high energy while moving down along the geomagnetic field if the interaction persists long enough. Constant v ez is assumed because the cyclotron harmonic resonance interaction mainly accelerates the gyration speed of the electron. However, if resonance interaction occurs at h = h u , i.e., ω uk = nΩ e + k z0 v ez , where k z0 = k z (h u ), the resonance condition will be detuned as the electron moves downward, and the magnitude |∆ω| of the detuning frequency ∆ω = ω uk − [nΩ e + k 1z (h)v ez ] = − [ k 1z (h) − k z0 ]v ez increases with the downward distance. It means that the interaction becomes weaker in time. On the other hand, if ω uk > nΩ e + k z v ez at h = h u , i.e., ω 0 > nΩ e by, e.g., 100 kHz, then ∆ω = ω uk − [nΩ e + k 1z (h)v ez ] = ∆ω 0 − [ k 1z (h) − k z0 ]v ez , where ∆ω 0 = ω uk − (nΩ e + k z0 v ez ) > 0, which is a gain factor; ∆ω 0 decreases as the electron moves downward with the wave. Hence, the interaction is continuously enhanced, and reaches cyclotron harmonic resonance at h = h i , where nΩ e + k 1z (h i )v ez = ω uk , and resonant electrons can gain considerable energy from the upper hybrid waves.
The power transfer from wave to the electron moving along the magnetic field, i.e., P || = −e E z v ez , depends on the phase of E z which varies with time due to ∆ω = 0.
Although the phase is distributed randomly from 0 to 2π, there are electrons of appropriate phases that will gain energy from the wave, i.e., P || > 0 and v ez increases; consequently, ∆ω decreases to give a positive feedback to enhance wave-electron interaction. When ∆ω drops to a negative value, the feedback of the interaction becomes negative; thus, a sufficient large ∆ω 0 is necessary to give adequate positive feedback interaction period for the generation of energetic electrons. On the other hand, if ∆ω 0 is too large, the initial electron-wave interaction will be too weak to lock in the positive feedback in the interaction, or the available interaction spatial distance will be too short to generate energetic electrons. A frequency ramp-up sweep serves the purpose of setting the proper gain factor for achieving continuous positive feedback during the wave-electron interaction in the available altitude range.

Spatial Bunching to Buildup Energetic Electrons Density in a Thin Layer for AIL Production
Only a small fraction of electrons in a height layer can achieve cyclotron harmonic resonance with the upper-hybrid waves at h = h i ; as the upper-hybrid waves propagate down, more electrons along the path can meet the condition to be in cyclotron harmonic resonance with the upper-hybrid waves. These electrons can bunch together in a layer, where significant ionizations occur to form an AIL.
Because the cyclotron harmonic resonance interaction mainly accelerates the electron gyration speed, constant v ez is adopted to illustrate the bunching phenomenon in Figure 14. As shown, down-going electrons initially located higher with slightly larger initial v ez (or at the same location but starting later), and those located lower with smaller initial v ez (or at the same location but starting earlier), can bunch together at h = h a . Although the three bunches between two dash lines in Figure 14 show the accumulation of energetic electrons at three different times, the time delay (difference) of the bunch converts to the spatial difference of the starting location. In essence, Figure 14 also illustrates that groups of energizing electrons starting at different altitudes can bunch together at the same location and time. Such a bunching process builds up the density of energetic electrons at h = h a , where an AIL is produced. 1+2bH/a 2 < 1 2 , where φ 2 10 = a b √ 1 + 2bH/a 2 − 1 , (71) is normalized to be where y 1 = d 2 dη 2 1 y 1 . The solution of (73a) is a Jacobi elliptic (cosine amplitude) function cn (η 1 , k 1 ); thus, φ p1 (ξ) = φ 10 cn (η 1 , k 1 ), which is a symmetric alternate function. √ 1 + 2bH/a 2 + 1 , (71) is normalized to be y 2 + 1 − 2k 2 2 y 2 + 2k 2 2 y 3 2 = 0 , where y 2 = d 2 dη 2 2 y 2 . (73b) has the same form as (73a), its solution is also a Jacobi elliptic function cn (η 2 , k 2 ). Thus, φ p2 (ξ) = φ 20 cn (η 2 , k 2 ); likewise, it is a symmetric alternate function. , where φ 2 30 = − a b √ 1 + 2bH/a 2 + 1 , (71) is normalized to be where y 3 = d 2 dη 2 3 y 3 . The solution of (73c) is a Jacobi elliptic (delta amplitude) function dn (η 3 , k 3 ); thus, φ p3 (ξ) = ξ 30 dn (η 3 , k 3 ), which is a non-alternate periodic function.