Two‐Dimensional Vlasov Simulations of Fast Stochastic Electron Heating in Ionospheric Modification Experiments

Ionospheric heating experiments using high‐frequency ordinary (O)‐mode electromagnetic waves have shown the induced formation of magnetic field‐aligned density striations in the ionospheric F region, in association with lower hybrid (LH) and upper hybrid (UH) turbulence. In recent experiments using high‐power transmitters, the creation of new plasma regions and the formation of descending artificial ionospheric layers (DAILs) have been observed. These are attributed to suprathermal electrons ionizing the neutral gas, so that the O‐mode reflection point and associated turbulence is moving to a progressively lower altitude. We present the results of two‐dimensional (2‐D) Vlasov simulations used to study the mode conversion of an O‐mode pump wave to trapped UH waves in a small‐scale density striation of circular cross section. Subsequent multiwave parametric decays lead to UH and LH turbulence and to the excitation of electron Bernstein (EB) waves. Large‐amplitude EB waves result in rapid stochastic electron heating when the wave amplitude exceeds a threshold value. For typical experimental parameters, the electron temperature is observed to rise from 1,500 K to about 8,000 K in a fraction of a millisecond, much faster than Ohmic heating due to collisions which occurs on a timescale of an order of a second. This initial heating could then lead to further acceleration due to Langmuir turbulence near the critical layer. Stochastic electron heating therefore represents an important potential mechanism for the formation of DAILs.


Introduction
Ionospheric heating facilities such as the High-frequency Active Auroral Research Program (HAARP) (Milikh et al., 1994) array in Alaska, USA, the Sura (Boiko et al., 1985) installation in Vasilsursk, Russia, and the European Incoherent Scatter superheater facility (Barr & Stubbe, 1991) in Tromsø, Norway, have been used for a rich array of experiments to study the interaction between large-amplitude high-frequency (HF) electromagnetic waves and ionospheric plasma. Large-amplitude HF ordinary (O)-mode pump waves injected into the overhead ionosphere can excite Langmuir and upper hybrid (UH) turbulence via parametric instabilities (Thidé et al., 2005;Eliasson, 2013;Najmi et al., 2015). The stimulated electromagnetic emissions escaping the plasma contains a rich spectrum of O-mode polarized sidebands (Leyser, 2001;Frolov et al., 2001;Grach, 1999. Recently, the formation of descending artificial ionospheric layers (DAILs) (Mishin & Pedersen, 2011;Mishin et al., 2016) have been observed. These are attributed to the ionization of neutral gas by suprathermal electrons (Mishin & Pedersen, 2011) accelerated by the ionospheric turbulence Yoon et al., 2012), leading to the reflection of the O-mode pump wave at a progressively lower altitude (Eliasson et al., 2012). This process is more efficient if the electrons are preheated so that a sufficient portion of the electrons in the tail of the distribution function can be resonantly accelerated by the Langmuir turbulence. A possible mechanism is that the electrons are bulk heated at the UH layer, then diffuse toward the critical layer where electron acceleration via Langmuir turbulence takes place. As an alternative acceleration mechanism, when the pump frequency is near an electron cyclotron harmonic nf ce (where n is an integer), it is predicted by quasi-linear theory that the excitation of UH and electron Bernstein (EB) waves can result in electron acceleration and the formation of high-energy tails in velocity perpendicular to the ambient magnetic field (Grach, 1999). Collisional Ohmic heating by mode-converted UH waves is conventionally considered for electron bulk heating near the UH layer (Gurevich et al., 1995). However, in the weakly collisional ionospheric F region, Ohmic heating occurs on a comparatively long timescale (of the order of seconds,) and stochastic Table 1 Ambient Plasma Parameters in the 2-D Vlasov Simulations: Run I Uses the Pump Frequency 0 = 21.72 × 10 6 s −1 (f 0 = 3.46 MHz) Corresponding to the Resonance Frequency of the N = −1, R = 1 Mode (cf. Figure 2), and Run II Uses 0 = 21.96 × 10 6 s −1 (f 0 = 3.49 MHz) Corresponding to that of the N  heating could potentially be a more efficient mechanism (Najmi et al., 2016). Stochastic heating of charged particles can occur transverse to a magnetic field in the presence of large-amplitude electric field gradients (Balikhin et al., 1993;McChesney et al., 1987;Stasiewicz et al., 2000) exceeding a threshold amplitude above which the individual particle orbits become unstable and diverge with time. A one-dimensional (1-D) Vlasov simulation study (Najmi et al., 2016) has demonstrated the viability of such a heating scheme at the UH layer, whereby electrons are efficiently energized via stochastic heating by short-wavelength, large-amplitude EB waves generated through the parametric decay of trapped UH waves in a magnetic field-aligned density striation (FAS). Submeter, supersmall FAS are believed to be produced by electrostatic processes (Gurevich & Zybin, 2006;Milikh et al., 2008;Najmi et al., 2014) within a few seconds of the O-mode pump wave being initiated and later evolve to be a few meters across but tens of kilometers in length along the magnetic field lines (Kelley et al., 1995;Franz et al., 1999).
In this paper we present the results of 2-D Vlasov simulations conducted to study the collisionless heating of electrons by the electrostatic turbulence resulting from trapped UH eigenmodes in a submeter scale FAS having circular cross section. The trapped UH eigenmodes are driven resonantly by an oscillatory left-hand circularly polarized dipole field representing the O-mode wave. The coupling between UH and LH waves leads to parametric instabilities and turbulence, and to the excitation of short-wavelength EB waves. Fast stochastic heating of the electrons occurs when the EB waves reach a threshold amplitude.

Numerical Setup
Vlasov simulations are carried out in two spatial (x, y) and two velocity (v x , v y ) dimensions with the simulation region defined as a transverse plane (perpendicular to the ambient magnetic field) having a small FAS of circular cross section. The simulations solve the 2-D Vlasov equation governing the evolution of the distribution function f for a charged particle species (= e for electrons and i for electrons) in position and velocity space, plus time, is the force and q and m are the particle charge and mass, respectively, with q i = e and q e = −e where e is the unit charge. The ambient magnetic field B 0 =ẑB 0 is directed along the z axis, perpendicularly to the simulation plane, and E ext = E 0 [x cos( 0 t) +ŷ sin( 0 t)] is an external left-hand polarized oscillatory dipole field representing the O-mode (L-O mode) pump wave with frequency 0 and amplitude E 0 . Herex,ŷ, andẑ are unit vectors in the x, y, and z directions. We use a pump wave amplitude E 0 = 2.0 V/m, with a hyperbolic tan rise time of 30 μs. For numerical efficiency, we adopt the electrostatic approximation E = −∇ for the self-consistent electric field, where the electrostatic potential is obtained from Poisson's equation The Vlasov equation is solved numerically using a Fourier transform method in velocity space (Eliasson, 2002(Eliasson, , 2010. The 2-D Vlasov simulation code is parallelized by means of domain decomposition (Daldorff & Eliasson, 2009), using the Message Passing Interface (MPI) for communication and is run using 72 cores on the Atoms, Beams and Plasmas Group's high performance computing (HPC) cluster as well as at the Archie-WeSt supercomputer facility, based at the University of Strathclyde. The simulations have a 2-D spatial domain of 3 × 3 m with a uniform grid of 200 × 200 intervals, yielding the spatial grid size Δx = Δy = 1.5 cm. Periodic conditions are used in space, with pseudospectral methods employed to approximate spatial derivatives. In velocity space, a grid of 80 × 80 mesh cells is used to resolve the velocity (v x , v y ) domain for both the ions and electrons. The electron distribution function is represented in the range ±17.9v Te0 for v x and v y where the electron thermal speed is v Te0 = √ k B T e0 ∕m e . The corresponding range for the ion distribution function is ±13.2v Ti0 for v x and v y , where the ion thermal speed is v Ti0 = √ k B T i0 ∕m i . The simulation duration is 260 μs over 200,000 time steps, using a fixed time step Δt = 1.3 ns.
A summary of the plasma parameters used in the simulations are provided in Table 1. To increase the numerical efficiency while giving sufficient separation in frequency between the high-frequency UH and EB waves and the low-frequency LH waves, we use the proton-to-electron mass ratio m i ∕m e = 1, 836, which is a factor of 16 smaller than the ion-to-electron mass ratio of the atomic oxygen majority species at the F 2 layer. The initial electron and ion temperatures are set to T e = 1, 500 K and T i = 1, 000 K, respectively, corresponding to typical values in the ionospheric F region (Kelley, 1989). We assume an ambient ion and electron number density of n 0 = 1.27 × 10 11 m −3 , giving the electron and ion plasma frequencies pe = √ n 0 e 2 ∕( 0 m e ) = 20.09 × 10 6 s −1 and pi = √ n 0 e 2 ∕ 0 m i = 4.691 × 10 5 s −1 . Using an applied geomagnetic field B 0 = 5.17 × 10 −5 T gives the electron and ion cyclotron frequencies ce = eB 0 ∕m e = 9.09 × 10 6 s −1 and ci = eB 0 ∕m i = 4.952 × 10 3 s −1 . The corresponding lower hybrid (LH) resonance frequency is LH = [ −2 pi +   Centered at x = y = 0 is a small 2-D striation (cf. Figure 1) represented by the initial ion and electron number densities where the normalized striation amplitude is = 0.1 and the transverse size is D str = 0.5 m. Trapped UH waves have frequencies below the ambient UH frequency but above the local UH frequency 21.11×10 6 s −1 at the bottom of the striation. For above the ambient UH frequency there is also a continuum of freely propagating UH waves. In the inhomogeneous plasma of the ionosphere, trapped UH waves are driven resonantly at quantized altitudes where the pump frequency matches the resonance frequency of the trapped eigenmodes (Dysthe et al., 1982;Mjølhus, 1998;Eliasson & Papadopoulos, 2015;Eliasson & Leyser, 2015). A field-structure analysis, based on the fluid model and method of Eliasson and Leyser (2015), is conducted of UH modes trapped in the FAS. In the analysis, the electrostatic potential is assumed to be proportional to (r) exp(iN ) exp(−i t) where r = √ x 2 + y 2 is the radial coordinate, = arctan(y∕x) the azimuthal coordinate, and N = 0, ±1, ±2, … the azimuthal mode indices. For trapped waves with below the ambient UH frequency, the eigenvalue problem with appropriate boundary conditions then gives discrete radial modes with mode indices R = 1, 2, … corresponding to the number of extrema of (r). The radial profiles and resonance frequencies of the first four trapped eigenmodes are shown in Figure 2.
In a similar manner as for 1-D striations (Mjølhus, 1998), the number of discrete trapped eigenmodes in a 2-D FAS can be estimated as (Eliasson & Leyser, 2015) where = 2 ∕( 2 − 4 2 ce ) ≈ 3.3 is a correction to the electron thermal pressure derived from kinetic theory, and De = v Te ∕ pe = 7.2×10 −3 m is the electron Debye length, giving M ≈ 12 trapped eigenmodes for = 0.1 and D str = 0.5 m. Of these, purely radial modes (N = 0) are not excited by a dipole field, and the L-polarized dipole field primarily excites modes with N < 0. The number of eigenmodes grows rapidly as the square of the striation size D str , leading to a continuum of trapped modes for large-scale striations. For example, a striation having a depth of = 0.05 and a width at half maximum of 15 m (Kelley et al., 1995;Franz et al., 1999), corresponding to D str ≈ 8 m would support M ≈ 1, 200 eigenmodes, approaching a continuum. Furthermore, ion density fluctuations due to turbulence would further blur the distinction between trapped eigenmodes into a continuous spectrum of trapped waves. It should be noted that taking into account electromagnetic effects (not covered in our electrostatic model), there would also be a leakage of Z mode (slow X mode) waves (Eliasson & Leyser, 2015;Hall & Leyser, 2003;Mjølhus, 1983), which, however, is important only for trapped UH waves in very small-amplitude striations and not for the large-amplitude striation studied here.

Simulation Results
We present the results of two simulation Runs (cf. Table I) using pump frequencies corresponding to the resonance frequencies of two different trapped UH modes. In Run I, we choose a pump wave frequency 0 = 2.389 ce = 21.72 × 10 6 s −1 which is at the resonance frequency of the N = −1, R = 1 mode (cf. Figure 2), while in Run II, we use 0 = 2.415 ce = 21.96 × 10 6 s −1 , corresponding to the resonance frequency of the N = −1, R = 2 mode. The pump frequencies in Runs I and II are somewhat above the local UH frequency UH = 21.59 × 10 6 s −1 at the locations where n e = 0.95n 0 (indicated in Figure 1).  Figure 3a, resonant excitation of the N = −1, R = 1 mode is clearly visible, resulting in that the x component of the electric field at the center of the density striation (x = y = 0) increasing in amplitude to about 30 V∕m at t = 0.07 ms. This is followed by a sharp decrease in amplitude until t = 0.12 ms, after which the amplitude again rises until t = 0.18 ms and again decreases, in a quasi-periodic manner. Figures 3b and 3c show the spatial profile of E x at t = 0.06 ms and t = 0.08 ms, respectively, which is near the time of maximum amplitude of E x in Figure 3a. At t = 0.06 ms in Figure 3b, the characteristic structure of the N = −1, R = 1 trapped UH wave is evident. Notably at t = 0.08 ms in Figure 3c, radially propagating small-scale EB waves escaping the FAS are also visible. Referring to Figure 4, it is clear that the sudden decrease in the amplitude of the trapped UH waves near t = 0.08 ms in Figure 3a coincides with the appearance of radial banding in the electron and ion spatial density profiles (Figures 4a  and 4b, respectively), indicative of the generation of radially propagating short-wavelength EB and LH waves. The approximate wavelength of the EB waves is 8.5 cm, corresponding to a wave number k x = ±75 m −1 . The coupling between the pump wave and the trapped UH wave decreases due to the formation of a cavity in the density profile at the center of the striation, seen in Figure 4b, resulting in a detuning between the pump frequency and the resonance frequency of the trapped N = −1, R = 1 mode, and an abrupt decrease of the wave electric field. The cavity formation is the result of the ponderomotive pressure of the high-amplitude UH waves acting on the electrons, which, in turn, pull the ions through the low-frequency space charge electric field. Figure 4d shows the frequency-wave number power spectrum, where E x has been Fourier transformed over the two spatial coordinates x, y, and time t, and the result plotted in logarithmic scale (dB), 20 log 10 |Ê x (f , k x , k y )|. Due to the rotational symmetry of the system, the turbulence is almost gyrotropic, and we have chosen to plot the data in Figure 4d as a slice at k y = 0. The obtained spectrum is consistent with the ones produced in the earlier 1-D simulations (Najmi et al., 2016), where only one spatial coordinate was used. It is seen in Figure 4d that strongly defined spectral components are concentrated at the second EB mode around f = 3.4 MHz, k x = ±75 m −1 , correlated with LH waves at f ≈ 100 kHz, k x = ±75 m −1 . The coupling between the pump wave and the trapped UH wave is clearly evident at f = 3.44 MHz at small wave numbers |k x | ≲ 10 m −1 , and the excitation of LH waves at small wave numbers |k x | < 10 m −1 at frequency f ∼ 35 kHz indicates that a parametric coupling between UH and LH waves takes place. There is also evidence of excitations of first EB mode waves at around f = 1.95 MHz, k x = ±65 m −1 . The coupling between UH, LH, and the second branch EB waves can be explained via three-wave parametric decay scenarios  (Najmi et al., 2016), while the coupling to the first EB branch could potentially be due to four-wave couplings involving the third EB branch, two first EB waves, and the trapped UH wave.

Run I:
We next turn to the collisionless heating of electrons by the large-amplitude electrostatic waves. The kinetic electron temperature perpendicular to the magnetic field is calculated from the 2-D electron distribution function as where v ex,ey = (1∕n e ) ∫ ∫ v x,y f e dv x dv y are the mean velocity components and n e = ∫ ∫ f e dv x dv y the electron number density. The kinetic temperature coincides with the conventional temperature for a Maxwellian electron distribution and approximately represents the "temperature" for non-Maxwellian distributions as long as they do not deviate too much from a Maxwellian distribution (Sckopke et al., 1983). Figure 5a shows the spatiotemporal evolution of T e plotted through the center of the striation at y = 0, while Figure 5b shows the maximum of T e over the whole simulation domain as a function of time. The trapped UH waves excite LH waves resulting in localized adiabatic compression of the electron density and low-amplitude variations of T e in the interior of the striation before t = 0.08 ms (cf. Figure 5a). These fluctuations correspond to LH waves with a frequency of ∼ 31 kHz. Following the parametric decay of the UH wave to EB waves at t ≈ 0.07 ms, there is a spike in the maximum T e to ∼2,150 K, correlated with the appearance of EB waves in the corresponding spatiotemporal evolution of E x in Figure 4c. This is followed by a decrease of the temperature fluctuations when the amplitude of the trapped UH wave decreases. Figure 5c shows a snapshot in time of adiabatic temperature fluctuations in the x-y plane due to the long wavelength LH waves, and Figure 5c shows smaller scale temperature fluctuations due to short-wavelength LH waves. The short-wavelength oscillations seen in Figure 5a for t > 0.17 ms and in Figure 5d are consistent with LH waves having the frequency ∼100 kHz and a wavelength of ∼ 0.1 m. In general, there is a relatively slow increase in the mean electron temperature in Run I.

Run II: N = −1, R = 2
Results from simulation Run II (cf. Table I) are shown in Figures 6-8. In this case, we tuned the pump frequency to the resonance frequency of the N = −1, R = 2 trapped UH mode (cf. Figure 2). In general, trapped UH modes with higher R have their maximum amplitude moving away from the center to the edge of the density striation (Eliasson & Leyser, 2015). Figure 6a shows the temporal evolution of E x at x = 0.25 m, y = 0, corresponding to the position of peak electric field (gradient of the electrostatic potential) for the N = −1, R = 2 trapped UH mode (see Figure 2). It is seen in Figure 6a that the amplitude of E x rises to a peak of about 25 V∕m at t ≈ 0.065 ms followed by a sharp decrease at t ≈ 0.075 ms, after which there is a rapid variation of the amplitude in a random fashion. Figures 6b and 6c show the spatial profile of E x at t = 0.05 ms and t = 0.075 ms, respectively, corresponding to when the amplitude of E x is near its maximum (cf. Figure 6a). At t = 0.05 ms (Figure 6b), the profile of E x is consistent with the N = −1, R = 2 trapped UH wave, while at t = 0.075 ms (Figure 6c) there are also short-wavelength EB waves propagating radially away from the density striation. The spatiotemporal evolution of the electron and ion number densities seen in Figures 7a and 7b, respectively, show that there is significant UH and LH wave turbulence at both the edges and interior of the density striation for t > 0.075 ms. This is consistent with the nonlinear saturation of a parametric instability involving trapped UH waves and LH waves. Visible in Figures 7a and 7b, there are also short-wavelength EB and LH waves radiating outward from the density striation at t ≈ 0.075 ms. The EB waves are also visible in the electric field in Figure 7c. Significant breakup of the density striation boundary occurs at t ≈ 0.08 ms with a reduction in both the trapped UH and EB wave amplitudes. There is also a visible increase in the spatial scale of the trapped UH wave E x profile in Figure 7c due to a collective increase in radial dimensions of the density striation. This results in an effective reduction in the N = −1, R = 2 eigenmode resonance frequency, and a detuning relative to the pump frequency. From t = 0.125 ms onward, there is also evidence of larger scale modulations in the ion density ( Figure 7b) that resemble wavefronts propagating away from the striation. These oscillations are consistent with long wavelength LH waves having the frequency f LH ≃ 31 kHz. The wave turbulence leads to a broad and diffuse frequency-wave number power spectrum of E x , shown in Figure 7d, with a rich array of multiwave couplings between the trapped UH wave, the first, second, and third EB modes, and the LH wave mode. The strongest coupling with the first EB mode occurs at k x ≈ ±60 m −1 , and there is a peak in the second EB mode coupling at k x ≈ ±75 m −1 .
Interestingly, there is a significant increase of the electron temperature in Run II, as illustrated in Figure 8. In Figure 8a it is seen that there is a sudden rise in electron temperature within the density striation at t ≈ 0.08 ms, with small-scale localized peaks in temperature. As the edge of the density striation is modulated and collectively expands in radius, the region of electron heating also expands. There are localized regions within the density striation where T e > 6,000 K. Figure 8b shows that the maximum electron temperature over the simulation domain rises abruptly by ∼3,000 K at t ≈ 0.08 ms, after which it fluctuates around T e ≈ 5, 000 K. At t ≈ 0.20 ms, there is a second rapid temperature rise to T e = 8,000 K followed by fluctuations around this temperature. As seen in the 2-D temperature plots in Figures 8c and 8d, the temperature is distributed as "hot spots" within the striation, with several thousands of degrees kelvin between the hotter and cooler regions.
In order to diagnose the variation in parametric decay/coupling to the first and second EB modes as a function of time, we performed a spectral analysis applying Gaussian temporal windows to the electric field data and performing fast Fourier transforms (FFTs) on the resultant data subsets. We applied four separate temporal windows numbered 1-4, and centered respectively at t = 0.035 ms, 0.075 ms, 0.145 ms, and 0.195 ms, where all windows had a full width at half maximum of 0.028 ms. The corresponding wave spectra are presented in Figure 9. Gaussian windows 2 and 4 centered at t = 0.075 ms and t = 0.195 ms, respectively, are of particular interest, as these correspond to the periods of sharp temperature rise in Figure 8b. Window 1 at t = 0.035 ms shows the initial wave spectrum prior to significant temperature rise, and for window 3 at t = 0.145 ms we see the spectral content during the intermediate saturation period prior to secondary temperature rise. With reference to Figure 9, the wave spectrum for the first Gaussian window centered around t = 0.035 ms (Figure 9a) has a region of intense trapped UH wave excitation at around 3.5 MHz (close to the pump frequency) and ranging from k = 0 → ±20 m −1 . There is also some evidence of weak coupling/decay to the second EB mode at higher wave numbers (k > 50 m −1 ) and LH content mirroring the wave number range of the trapped UH waves and second EB coupling. This is indicative of three-wave parametric decay of the pump-coupled UH wave to frequency-downshifted UH waves mediated by LH waves, and to second EB waves and LH waves. Using the Gaussian window centered at t = 0.075 ms, Figure 9b shows a large increase in the coupling to the second EB wave at k = ±75 m −1 . This coincides with the steep rise in electron temperature observed in Figure 8b, suggesting the large increase in the second EB wave excitation is a potential driver of the electron heating. Although coupling to the first EB branch is present, it is very weak in comparison to that of the second EB wave. A very different picture is seen in Figure 9c using the third Gaussian window centered at t = 0.145 ms. During this saturation period in electron temperature, there is an increased coupling to the 1st EB mode for k = ±(50 → 90) m −1 , and some coupling to the third EB wave. There is also a reduction in the coupling to the second EB branch. When we reach the region of secondary temperature rise at around t = 0.18 ms covered by the fourth Gaussian window, we see in Figure 9d that the spectral/wave number range of coupling to the first EB mode has increased further, while there is no evidence of increased or restored coupling to the second EB mode at k = ±75 m −1 . This suggests that the electron temperature rise at t ≈ 0.18 ms is not associated with the second EB mode but with the first EB mode. The temporal window analysis therefore suggests that both the first and second EB modes are directly associated with electron heating, albeit at different times in Run II.

Analysis of Results: Stochastic Electron Heating
In order to determine whether the observed heating of electrons in simulation Runs I and II is associated with stochastic heating by large-amplitude EB waves, we consider the phenomenology of stochastic particle Figure 10. The normalized wave amplitude A (defined in equation (7)) as a function of x (at y = 0) and t for (a) Run I and (b) Run II.
heating and the associated threshold amplitude in the electrostatic wave field. Some of the earliest observations of stochastic electron and ion heating were made in theta pinch plasma devices (Demchenko et al., 1971) and magnetic mirror traps (Namba & Kawamura, 1970). These initial considerations which primarily considered electron cyclotron coupling (Namba & Kawamura, 1970;Kawamura et al., 1971) and ion/electron heating in turbulent plasma with cross-field current flow (Demchenko et al., 1971) were followed by studies of stochastic ion heating due to electrostatic LH waves propagating both oblique (Smith & Kaufman, 1975) and perpendicular (Fukuyama et al., 1977;Karney, 1978) to a magnetostatic field. Stochastic heating of ions and electrons by perpendicular collisionless shocks (Balikhin et al., 1993;Stasiewicz et al., 2000) and ion stochastic heating by drift Alfvén waves (McChesney et al., 1987) have also been studied. In general, any scenario whereby a charged particle is accelerated by a sharp electric field gradient (perpendicular to a magnetostatic field) can lead to stochastic motion and associated rapid heating of the plasma species.
The equations of motion for a charged particle subject to an electrostatic plane wave propagating perpendicularly to a static magnetic field can be reduced to the dimensionless model equation (McChesney et al., 1987) where X represents the particle's position normalized by k −1 , is time normalized by −1 ci , A = (mk 2 )∕(qB 2 0 ) is the normalized wave amplitude of the electrostatic wave, Ω = ∕ c is the wave frequency normalized to the cyclotron frequency, and is the electrostatic potential. For Ω < 1 and small v it was noted in reference (McChesney et al., 1987) that particle orbital motion becomes chaotic and divergent for |A| ≳ 1; therefore, this represents the threshold of stochasticity (other regimes of stochasticity exist for Ω ≫ 1 (Karney, 1978) and near cyclotron harmonics where Ω takes integer values (Fukuyama et al., 1977)). At the threshold |A| ≳ 1, the particle displacement induced by the oscillatory electrostatic field of the wave becomes comparable to the wavelength of the electrostatic wave, resulting in breakdown of the polarization drift approximation McChesney et al. (1987). An analogous expression, was used by Balikhin et al. (1993) and Stasiewicz et al. (2000) (restricted to 1-D) as a criterion for the onset of stochastic motion of charged particles in large-amplitude electric field gradients in collisionless shocks and nonlinear structures. As for the previous case for monochromatic electrostatic waves, when |A| ≳ 1, the orbits of neighboring particles become exponentially divergent in time, leading to rapid stochastic heating of the charged particles.
To compare the conditions for onset of stochastic heating in simulation Runs I and II, we plot the normalized electric field gradient A in equation (7) (using q = −e and m = m e ) as a function of x (at y = 0) and t in Figure 10. The color map is limited to |A| ≤ 1 for clarity. For the N = −1, R = 1 case in Run I, stochasticity could potentially set in (A ≥ 1) during two short time intervals t = 0.06 ms to 0.08 ms and t = 0.17 ms to 0.19 ms (cf. Figure 10a). However, as seen in Figure 5, there was no significant electron heating for this case. Looking at the spatiotemporal evolution of A for the N = −1, R = 2 case (Run II), it is seen in Figure 10b that there is a steady increase in the amplitude of A with time within the density striation until t = 0.065 ms when high-amplitude, small-scale radial bands appear in A within the striation and radiating outward from the striation. The amplitude increase of A and wave breakup at t ≈ 0.085 ms seen in Figure 10b coincides approximately with the sharp rise of T e to about 5, 000 K seen in Figure 8 at t ≈ 0.8 ms. The turbulent structure in A remains present until t ≈ 0.17 ms, where an amplitude increase again appears at x = ±0.25 m. This coincides with a second rise of T e to 8000 K at t = 0.18 ms seen in Figure 8. Hence, it is likely that the sharp temperature rises observed in the N = −1, R = 2 simulation (Run II) are due to stochastic heating by large-amplitude short-wavelength EB waves excited at the first and second EB branches.

Conclusions
In summary, the mode conversion of a left-hand polarized O-mode pump wave to upper hybrid (UH) waves trapped in a small-scale density striation, and subsequent electrostatic wave turbulence and stochastic electron heating have been investigated by 2-D Vlasov simulations. The trapped UH wave undergoes parametric instabilities giving rise to short-wavelength electron Bernstein waves and lower hybrid waves. When exceeding a threshold value, electric field gradients lead to stochastic motion of the electrons and a rapid increase in temperature. Hence, in addition to Ohmic heating due to electron collisions with neutrals and ions, stochastic heating may lead to efficient conversion of wave energy to kinetic energy of electrons. With application to ionospheric heating experiments, stochastic electron heating could therefore be an important mechanism in the observed creation of new plasma regions and the formation of ionization fronts/descending artificial ionospheric layers (DAILS) when the energy of the fastest electrons exceeds the ionization potential of the neutral gas. Future works will explore alternative acceleration scenarios occurring when the pump frequency is near a cyclotron harmonic of the charged particles/electrons. Quasi-linear theory (Grach, 1999) and studies of stochastic heating (Fukuyama et al., 1977) predict that the resonant wave-particle interactions can then result in the formation of high-energy tails in velocity space perpendicular to the ambient magnetic field.