Electrostatic electron cyclotron instabilities near the upper hybrid layer due to electron ring distributions

A theoretical study is presented of the electrostatic electron cyclotron instability involving Bernstein modes in a magnetized plasma. The presence of a tenuous thermal ring distribution in a Maxwellian plasma decreases the frequency of the upper hybrid branch of the electron Bernstein mode until it merges with the nearest lower branch with a resulting instability. The instability occurs when the upper hybrid frequency is somewhat above the third, fourth, and higher electron cyclotron harmonics, and gives rise to a narrow spectrum of waves around the electron cyclotron harmonic nearest to the upper hybrid frequency. For a tenuous cold ring distribution together with a Maxwellian distribution an instability can take place also near the second electron cyclotron harmonic. Noise-free Vlasov simulations are used to assess the theoretical linear growth-rates and frequency spectra, and to study the nonlinear evolution of the instability. The relevance of the results to laboratory and ionospheric heating experiments is discussed.


Introduction
The interest in cyclotron harmonic phenomena in magnetized plasma [1][2][3][4][5][6][7][8][9] were triggered by the observation of emissions, absorption, and resonances at harmonics of the electron cyclotron frequency in laboratory experiments and satellite sounding experiments [10,11]. Solar radio emissions [12,13], and some spectral features of radiation near electron cyclotron harmonics in ionospheric heating experiments [14] have also been attributed to electron cyclotron instabilities. Ringlike ion distributions may also be formed due to neutral beam injection [15,16], charge-exchange collisions in drifting plasma [17,18], solar wind interactions with comets [19,20], etc, and result in the excitation of electrostatic low-frequency waves. In general, maser instabilities due to loss-cone or ring-like electron distributions give rise to excitations in the Z (slow X) mode or upper hybrid mode branch for pe ce ω ω [13], while radiation in the fast X mode branch near ce ω dominates in the opposite limit. The latter has been studied extensively in the framework of auroral kilometric radiation [21], solar and stellar radio bursts [22], etc, and via simulations and laboratory experiments [23][24][25]. There are also anisotropydriven electromagnetic and electrostatic waves for different sets of parameters [26]. Recent laboratory experiments [27,28] using a mirror confined plasma have shown a number of different instabilities and emissions, including oscillations near the second electron cyclotron harmonic attributed to the presence of thermal ring distributions [12]. Theoretical and numerical studies have been carried out of thermal ring electron distributions [30,31] and of a maser-beam instability of Bernstein waves [29]. Particle-in-cell simulations have been used to study Bernstein waves [32] and instabilities for combinations of ring and core populations [33][34][35]. It has been recognized that the inclusion of a cold core distribution lowers the threshold for instability compared to a pure ring distribution [33,34,36,37].
The aim of this paper is to carry out a theoretical and numerical study of the electrostatic electron cyclotron instability and to discuss its relevance to spectral features of stimulated electro magnetic emissions escaping the plasma in laboratory [27,28] and ionospheric heating [14] experiments. A parallelized Vlasov code [38,39] in two spatial and two velocity dimensions (x, y, v x , v y ), plus time, is used to carry out noise-free simulations to assess the theoretical results and to study the nonlinear saturation of the instabilities.

Theory
Electrostatic waves in a magnetized electron plasma with stationary ions are governed by the Vlasov-Poisson system and where e is the magnitude of the electron charge, m e is the electron mass, B 0 is the external magnetic field, 0 ε is the vacuum electric permittivity, and n 0 is the total equilibrium electron number density. Linear electron Bernstein waves are found by linearizing the unknown variables as f | | , and Fourier analyzing the Vlasov-Poisson system of equations by assuming that f 1 and 1 φ are where k is the wave vector and ω is the frequency.
We will concentrate on electron Bernstein waves propagating perpendicularly to the magnetic field B B z 0 0 =ˆ so that x y r x y = + and , where x , y , and ẑ are the units vectors along the x, y and z axes. Hence the wave vector component parallel to B 0 is set to zero, k z = 0. To obtain the dispersion relation for a distribution of electrons, we follow the formalism of [5]. The starting point is the deltafunction electron velocity ring distribution where δ is Dirac's delta function, v z0 is the particle velocity along the z-axis, parallel to the magnetic field, and all electrons have the same speed v 0 ⊥ transversely orbiting the magnetic field. We have denoted v v v x y For perpendicular propagation (k z = 0), the equilibrium distribution function (3) then yields the electron susceptibility on integral form [5] k v sin sin J 2 cos 2 d sin , is the electron plasma frequency, eB m ce e 0 / ω = is the electron gyrofrequency, J 0 is the Bessel function of the first kind of order zero, and k k k To obtain the susceptibility for a different velocity distribution, it now suffices to multiply equation (4) by a weighting function depending on the shape of the velocity distribution (normalized to unity) and to integrate over velocity space.
For a Maxwellian equilibrium distribution where the integration over v z0 has been carried out. By using the Bessel function identity [40] the integration over v 0 ⊥ can be carried out, with the result [5] exp sin sin sin exp cos d , . The integral form of the susceptibility is more amenable to numerical evaluation than forms containing infinite sums of Bessel functions.
For a thermal ring distribution of suprathermal electrons on the form [4,12,14] F we have the susceptibility where the integration over v z0 has been carried out. By using the identity which is obtained by differentiating both sides of equation (7) with respect to 2 γ − , the integration over v 0 ⊥ in equation (10) 2 . For a plasma consisting of a sum of a Maxwellian core distribution and thermal and delta-function ring distributions, the dispersion relation governing the linear growth-rate of the instability can now be written is the relative number density fraction of the thermal and delta-function ring electrons, respectively. To obtain the complex-valued wave frequency i I R ω ω ω = + from the dispersion relation (13), where R ω is the real frequency and I ω the growth-rate, we carry out the integrals over ψ numerically by means of a sum representation and solve the dispersion relation iteratively.

Numerical setup and simulation results
We carry out a set of Vlasov simulations using combinations of Maxwellian and ring distributions, and compare the results with solutions of the linear dispersion relation. The plasma parameters used for the 9 simulation runs are listed in table 1, where it is indicated in which figures the respective numerical result are presented. An equilibrium distribution function consistent with the dispersion relation (13) is The numerical simulation method is based on the Vlasov equation Fourier transformed in velocity space [38] where the solver is parallelized [39] by (i) integrating the Vlasov equation for different particle species in parallel on separate processors, and (ii) domain decomposition of each particle species in the Fourier transformed velocity space and configuration space. The communication between processors is done by means of message passing interface (MPI). The initially Maxwellian and ring distributions are treated as different electron species in the simulations. Expressions for the initial conditions in the Fourier transformed velocity space are given in the appendix. As described in [38], a 4th-order compact difference scheme is used to calculate derivatives in the Fourier transformed velocity space, a pseudo-spectral method is used in configuration space with periodic boundary conditions, and the standard 4th-order Runge-Kutta scheme is used to advance the solution in time. The code is run in electrostatic mode to lessen the stability constraint (Courant condition) on the time step.
The hot thermal ring distribution is typically one order of magnitude wider in velocity than the Maxwellian distribution, while the opposite holds in the Fourier transformed velocity space. To resolve both the Maxwellian and ring populations accurately in the Fourier transformed velocity space, we have used the following numerical parameters: the simulation box is along the x-direction is  (d)), we use a density fraction of 10% for the ring distributions and 90% for the Maxwellian distribution. Figures 1-3 show the real frequencies and growth rates obtained by solving the dispersion relation (14) for a sum of a Maxwellian core distribution and a thermal ring distribution, as well as the frequency spectra and estimated growth-rates from the corresponding Vlasov simulations (see runs 1-4 in table 1). The spectra and growth-rates are estimated for times (indicated in the figure captions) before nonlinear saturation of any instability. Figures 1-3 use an upper hybrid frequency somewhat above the 4th, 3rd and 2nd electron cyclotron harmonic, respectively, with 4   ) ω space are concentrated to the linear Bernstein modes, and that the growth-rate of the instability is consistent with linear theory. The growth-rate is estimated in the simulations by assuming that the amplitude of each wavenumber component of the electric field is growing as t exp I ( ) ω for the growth-rate I ω , and hence its logarithm grows linearly with time as A t where A is a constant; a linear regression scheme 4 is used to estimate I ω for each wavenumber to produce the plots in figures 1(d) and 2(d). Compared to the purely Maxwellian case (indicated with dashed lines in figures 1(a) and 2(a)), the upper hybrid branch the decreases its frequency until it touches the nearest lower electron Bernstein branch. An important result is here that while the upper hybrid frequency has been chosen slightly above an electron cyclotron harmonic harmonic, the frequency of the unstable waves is close to the harmonic of the electron cyclotron frequency, which is significantly below the upper hybrid frequency. Theoretical and simulation results using the delta-function ring distribution in the initial conditions (see runs 5-9 in table 1) are shown in figures 4-7. The ring distribution leads to a very rapid growth of the waves, with a growth-rate about one order of magnitude larger than for the corresponding case using a thermal ring distribution (see figures 4(d) and 5(d)). The large growth-rate leads to an increase of the amplitude a factor two in each wave period, leading to a somewhat blurred picture of the frequency spectra in figures 4(c) and 5(c).        / ω ω = ) shows unstable wave modes for three groups of wavenumbers with frequencies near the 4th electron cyclotron harmonic, and in addition an instability where the 2nd electron Bernstein branch merges with the 3rd branch. Also for 3.1 uh ce / ω ω = in figure 5 there are two groups of unstable waves near the 3rd cyclotron harmonic. Figure 6 shows that in contrast to the case of a thermal ring distribution there is also an instability for 2.1 uh ce / ω ω = using the deltafunction ring distribution where the the 1st and 2nd electron Bernstein branches merge. When the upper hybrid frequency is decreased to 2.0 uh ce / ω ω = (see figure 7(a)), the two modes no longer merge, and the instability disappears. The case of a pure ring distribution ( 1 η = δ , 0 th η = ) in figure 7(b) is also stable with respect to the perpendicular electron cyclotron instability.
We have here mostly concentrated on cases where the thermal ring distribution is a relatively small fraction of the total electron distribution. The opposite case, when the ring distribution has a larger density [33] or when the core distributions is absent [5], leads in general to instabilities between electron cyclotron harmonics (half-harmonic radiation) at points in the spectrum where the electron Bernstein branches merge, with maximum growth-rates for waves propagating at oblique angles to the magnetic field. It is interesting that the inclusion of a cold core distribution lowers the threshold for instability compared to a pure ring distribution [33,34,36,37]. As an example, the pure delta-function ring distribution becomes unstable only for 6.62 2.57 pe ce ce ω ω ω ≈ [5], while we see in figure 6 that the combination of a tenuous delta-function ring distribution and a Maxwellian core distribution is unstable for 2.1 uh ce On the other hand, when the delta-function ring distribution is replaced by the thermal ring distribution, in figure 7(a), the system is again stabilized for 2.1 uh ce . This is consistent with the study of electron cyclotron harmonic instabilities [42] where a large enough cold component stabilizes the low electron cyclotron harmonic bands. We have not been able to find an instability near the second electron cyclotron harmonic using the thermal ring distribution. Recent laboratory experiments [27,28], on the other hand, have observed particularly strong bursts of electromagnetic emissions when the upper hybrid frequency is in the vicinity of the second electron cyclotron harmonic. It could indicate the presence of a cold electron population and a relatively cold ring distribution. It seems from figures 1 and 2 that the electron cyclotron instability is easier to excite when the upper hybrid frequency is somewhat above the third and higher electron cyclotron harmonics where the threshold for instability is lower. Figure 8 shows the exponential growth in time and nonlinear saturation of the wave electric field for different combinations of thermal and delta-function ring distributions with Maxwellian distributions corresponding to the parameters in figures 1, 2, 4 and 5. For the thermal ring distribution in figures 8(a) and (b), the electric wave field reaches a saturation level that stays throughout the simulation, while for the thermal ring-distribution in figures 8(c) and (d), the rapidly growing electric field peaks at saturation after which the electric field amplitude is weakly decreasing with time. There is no significant difference in electric field amplitudes comparing 3 / ω ω = (bottom panels). As noted in previous studies, the nonlinear saturation of the instability leads to a loss of energy of the ring distributions and an energy gain (heating) of the core distribution [33]. This is consistent with the volume-rendered 5 phase-space plots in figures 9 and 10 for the thermal ring and delta-function ring distributions, respectively, for the case The instability gives rise to large amplitude waves, in particular visible in the initially cold Maxwellian distribution of electrons (right-hand columns of figures 9 and 10). The large amplitude waves modify the distribution functions via wave-particle interactions. Most notable are the changes in the tenuous ring distributions, seen in the right-hand panels of figures 9 and 10, since these are the minority species in our simulations: in general, the instability leads to that a portion of the ring-distributed electrons move towards lower speeds to fill in the minimum in the total distribution function. The wave turbulence also accelerates a small portion of the initially cold Maxwellian distribution, as is visible in the right-hand columns in figures 9 and 10. Once marginal stability has been obtained, the wave amplitude stops growing in time and either maintains the saturation level or slowly decreases. The saturation process may involve resonance broadening of the cold component [43] as well as nonlinear Landau damping [45] and quasilinear diffusion of the hot component.

Discussion
In conclusion, we have carried out a theoretical study of the electron cyclotron instability due to the sum of a Maxwellian core distribution and a thermal ring or deltafunction ring electron distribution. The theoretical results are supported by Vlasov simulations. For a tenuous thermal ring distribution, the main result is that the electrostatic instability perpendicular to the magnetic field occurs when the upper hybrid frequency is slightly above one of the electron cyclotron harmonics, starting with the third electron cyclotron harmonic. For a sum of a Maxwellian distribution and a tenuous delta-function ring distribution, an instability occurs also near the second harmonic of the electron cyclotron frequency. Hence, the recent experimental observations [27,28] of emissions near the double resonance where the upper hybrid frequency is in the vicinity of the second electron cyclotron harmonic, could indicate the presence of a cold electron component and a cold ring distribution of electrons. It has also been suggested that the broad upshifted maximum feature of the stimulated electromagnetic emissions (SEE) in ionospheric heating experiments could be explained by an electron cyclotron instability when the upper hybrid frequency is slightly above one of the electron cyclotron harmonics [14], and where a small portion of the electrostatic wave energy is mode converted to ordinary mode waves escaping the plasma as SEE. However, since the instability takes place by the merging of the upper hybrid branch of the electron Bernstein wave and the next lower branch, giving rise to oscillations at or below the electron cyclotron harmonic, this would require a separate mechanism to increase the frequency of these oscillations to form the broad upshifted maximum. Experimental observations [46][47][48] show that the peak frequency BUM ω approximately obeys n BUM c e 0 0 ω ω ω ω − = − when the transmitted frequency 0 ω is above one of the electron cyclotron harmonics, n ce 0 ω ω > . An alternative model for the broad upshifted maximum is a 4-wave parametric decay scenario [49,50] involving the transmitted pump wave, the upper hybrid and nearest lower electron Bernstein branch, and the lower hybrid waves, which seems to account for the observed frequency upshift.

Acknowledgments
Discussions with Thomas Leyser at the Swedish Institute of Space Physics are gratefully acknowledged. This work was supported by the Engineering and Physical Sciences Research Council (EPSRC), U.K., Grant no. EP/M009386/1. Simulation data supporting the figures are available at http://dx.doi. org/10.15129/56448d9e-adb0-4d2b-afdb-029165a40f54.

Appendix. Initial conditions in Fourier transformed velocity space
The initial conditions for the Fourier transformed electron distribution functions to be used in the numerical code [39] are obtained by using the Fourier transform pair x y x y x x y y x y 2 (A.1b) which applied to equations (15a)-(15c) give the Fourier transformed Maxwellian, thermal ring, and delta-function ring distribution function, respectively, as