Generation of auroral turbulence through the magnetosphere–ionosphere coupling

The shear Alfvén waves coupled with the ionospheric density fluctuations in auroral regions of a planetary magnetosphere are modeled by a set of the reduced magnetohydrodynamic and two-fluid equations. When the drift velocity of the magnetized plasma due to the background electric field exceeds a critical value, the magnetosphere–ionosphere (M–I) coupling system is unstable to the feedback instability which leads to formation of auroral arc structures with ionospheric density and current enhancements. As the feedback (primary) instability grows, a secondary mode appears and deforms the auroral structures. A perturbative (quasilinear) analysis clarifies the secondary growth of the Kelvin–Helmholtz type instability driven by the primary instability growth in the feedback M–I coupling. In the nonlinear stage of the feedback instability, furthermore, auroral turbulence is spontaneously generated, where the equipartition of kinetic and magnetic energy is confirmed in the quasi-steady turbulence.


Introduction
Plasma turbulence ubiquitously found in fusion, space, and astrophysical plasmas plays leading roles in mixing and transport of particles, momentum and energy. Solar wind is one of the most popular examples of space plasma turbulence (for example, see [1]), where the magnetohydrodynamic (MHD) turbulence [2] is considered to be dominant in a large scale. Because of the high-β value, the background (or mean) magnetic field intensity in the solar wind is the same order as the turbulence fluctuations. In contrast, turbulence intensity in low-β plasma is considered to be much weaker than the background field magnitude. It is typically found in fusion plasmas as well as in planetary magnetospheres. Our interest in the present study is generation of the low-β MHD turbulence in the auroral magnetosphere.
The MHD plasma turbulence is supposed to involve nonlinear interactions of shear Alfvén waves. In the Goldreich-Sridhar (G-S) theory [3], interaction of oppositely propagating waves with an equal amplitude constitutes an elementary process for the turbulence cascade, where the energy spectrum of the Alfvénic turbulence is predicted to be µk 5 3 . Here,k denotes the perpendicular component of the wavenumber vector with respect to the ambient magnetic field. The Alfvénic turbulence in the planetary magnetosphere would be a typical example of the low-β MHD turbulence, where a strong mean magnetic field makes the turbulence quasi two-dimensional as supposed in the G-S model. Indeed, some satellite observations in the Earth's magnetosphere report turbulent fluctuations of the electric and magnetic fields in the auroral region [4], and some types of aurorae are considered to be closely related to the observed Alfvén waves [5,6]. However, generation mechanism of the Alfvénic turbulence in the auroral magnetosphere remains to be answered. The auroral turbulence is, thus, an interesting application of the low-β MHD turbulence from both the viewpoints of plasma physics and space physics.
Aurorae are commonly observed in magnetized planets, such as Earth, Jupiter, and Saturn, etc. The magnetosphere-ionosphere (M-I) coupling system [7] is the platform of auroral phenomena. Auroral brightening is caused by precipitation of keV order electrons to the ionosphere where neutral atoms and molecules are excited. However, physics mechanisms of spontaneous growth of auroral arc structures and their dynamics are still open issues. One of the possible mechanisms of auroral growth is the feedback instability proposed by Sato in 1970s [8], where the shear Alfvén waves propagating along the planetary magnetic field resonantly interact with ionospheric density waves. The feedback instability growth is associated with enhancement ofÉ B shear flows and ionospheric density and field-aligned current fluctuations. Spatial structures and intensity of the auroral current system expected from the feedback instability are qualitatively consistent with auroral arc observations [8]. So far, a variety of linear analyses of the feedback instability have been carried out with the linear Alfvén wave response [9], the three-dimensional MHD model [10,11], the twofluid model [12], the gyrokinetic descriptions [13], the ionospheric density cavity [14], the dipole geometry [15], and so on. Nonlinear simulations of the feedback instability such as [16,17] are, however, still limited, but essential to investigate deformation of auroral arc structures and transition to the auroral turbulence.
In the present study, we investigate the nonlinear evolution of the feedback instability and transition to the Alfvénic turbulence in the M-I coupling. Deformation of the linear eigenmode structure of the feedback instability was demonstrated in our previous simulation, while its detailed mechanism remains to be investigated. Here, we apply a perturbative analysis to the M-I coupling, and have found a critical amplitude of the primary mode for growth of the secondary instability. When the primary instability is continuously driven by the external energy source, that is, a cross-field (É B) plasma flow, further nonlinear growth of instabilities leads to transition to turbulence. It also motivated us to develop a new reduced MHD simulation code with a spectral method. In the latter part of this article, we report the first simulation result of Alfvénic turbulence generation in the auroral M-I coupling.
Remainder of this paper is organized as follows. The M-I coupling model is elucidated in section 2. The perturbative analysis of the secondary instability is given in section 3. The nonlinear simulation results of the turbulence transition will be shown and discussed in section 4. The obtained results are summarized in section 5.

Model
In the present study, we employ the reduced MHD model [18] for the magnetosphere with the uniform magnetic field B 0 and the homogeneous plasma density r 0 for simplicity. This is the simplest MHD model for the magnetosphere with a strong planetary magnetic field and low plasma pressure, where the shear Alfvén law plays a key role in the plasma dynamics. The ionosphere is filled with a partially ionized plasma with low temperature, and its vertical scale length is much shorter than the magnetospheric one. The height-integrated (two-dimensional) two-fluid equations are adopted to the ionosphere model. The M-I coupling model with the slab geometry is the same as used in the previous work [16].
The governing equations for the magnetosphere are given by the vorticity equation and the induction equation as follows: . Fluctuations of the vorticity and the magnetic flux function are denoted by ω and ψ, respectively, which are related to the stream function (electric potential) f and the parallel current j  , that is The unit vector parallel to B 0 is denoted by b 0 . The other notations are the uniform electric field E 0 perpendicular to B 0 , the viscosity ν, the resistivity η, and the permeability m 0 . It should be noted that the largescale convection in the magnetosphere is modeled by the equilibrium plasma flow ofÉ B B 0 2 which is a driving source of the feedback instability.
Model equations for the ionosphere are given by the continuity equations for ions and electrons with the quasi-neutrality condition, where ion dynamics is governed by the ion-neutral collisions and the perpendicular electric field while electrons are magnetized. In the followings, we assume that B 0 is in the vertical direction to the auroral ionosphere, and that the ionospheric electric field is given by the electric potential, f. Then, one finds the height-integrated continuity equations for the number density and the ionospheric current defined on the two-dimensional plane perpendicular to B 0 , that is where n and n 0 are perturbed and equilibrium parts of the ionospheric electron number density. The other symbols mean the effective ionospheric height h, the elementary charge e, the recombination rate α, the Pedersen mobility of ions m P , and the particle diffusion coefficient D. Here, we assume the linearized recombination term, a n n 2 0 . The magnetospheric and ionospheric equations are coupled through the continuity of the perpendicular electric field given by f -B 0 and the field-aligned current j  . In other words, a solution of the ionospheric equations provides the boundary condition for the magnetospheric MHD equations. Solutions of the above set of equations are simply described on the coordinates moving with the equilibrium flow, where the advection terms in equations (1), (2), and (5) are absorbed into the time derivative terms. Thus, the real frequency calculated on the moving frame is Doppler-shifted by^· k V E B , wherek is the perpendicular wavenumber.
In the followings, we employ the rectangular coordinates, ( ) x y z , , , with B 0 pointing the positive z direction. In the present model, x and y are assumed to be horizontal, and z is taken for the vertical coordinate, while ( ) x y z , , corresponds to the azimuthal, the latitudinal (at the foot point of field lines), and the field-aligned coordinates of the dipole geometry, respectively. The Poisson brackets, x . The background electric field E 0 is set in the y direction, and thus,V E B is parallel to the x axis. It should be noted that solutions of the above set of equations are classified into even and odd parity modes with respect to the magnetospheric equatorial plane, as the equilibrium is symmetric in the northern and southern hemispheres. In the present study, we consider the even parity for f (and ω) and the odd parity for ψ (and j  ), say,  (5) and (6) as theÉ B drift term due to E 0 . In the followings, we consider the same parameters as those used in [16], such that, the distance between the ionosphere and the magnetospheric equator = l L 1000 , the uniform external electric field = -E B V

Perturbative analysis
The eigenmode analyses for the linearized set of equations have been carried out for the slab [16] and the dipole [15] configurations. Linearizing equations (1)- (6), and using the periodic boundary condition for x and y, one finds the dispersion relation of the shear Alfvén waves coupled with the ionospheric density wave. For the slab geometry with the uniform background (B 0 , E 0 and r 0 ), the linear eigenfunction for n h = = 0 is given by where A L is an arbitrary constant and k  is a complex valued wavenumber along B 0 . The parallel wavenumber is related to the eigenfrequency, such that, W = k V L A  . The perpendicular wavenumber is denoted by = ( ) k k k , x y . In the above, the two-dimensional ionosphere and the magnetospheric equatorial plane are set at z=0 and l, respectively.
If the convection electric field E 0 exceeds a critical value, the feedback instability is switched on, and amplitude of the unstable eigenmode grows exponentially in time with a sinusoidal mode pattern in the x-y plane (for example, see figure 4 (a)). When the mode amplitude exceeds a critical level, nonlinear coupling with other modes influences the linear instability growth. Actually, our previous simulation of the feedback instability demonstrates deformation of the mode structure in a nonlinear phase of the instability growth [16], where rolling-up of vortex sheets driven by the Kelvin-Helmholtz (K-H) type instability is observed on the magnetospheric equatorial plane. In this section, we investigate the secondary growth of the K-H type instability driven by the primary unstable mode.
Let us suppose that the perturbed quantities, such as ω, ψ, and n, are represented by a sum of a linear eigenmode with a finite amplitude and deviations from it. For example w w w y y y , a n d , 9 L L L where the subscript L denotes the linear eigenfunction with the eigenfrequency W L , and · means a deviation from it. Substituting equation (9) into equations (1)- (6), and subtracting the linearized equations for w L , y L , and n L on the moving frame withV E B , one finds Second, the eigenmode solution is substituted into (10)- (13), where the mode amplitude A L is fixed as a given parameter, while the wave propagation is introduced in terms of the real part of the eigenfrequency W ( ) . This perturbative (or quasilinear) approach provides us an instantaneous estimate of the secondary mode growth for a given value of A L . Starting with small amplitude random perturbations of w y (˜˜˜) n , , , we observed the exponential growth of the secondary mode, if A L exceeds a critical value. The mode amplitude is normalized so that the maximum of w |˜| on the equatorial plane at z=l is equal to unity. For clarity of the visualization, we artificially reduced the plotting scale in the z direction. In this computation we introduced a finite value of the viscosity and the resistivity, that is, n h m = =´-LV 4 10 0 7 A as same as those used in the nonlinear simulation [16]. On the equatorial plane shown in figure 1(a), vortex streets with positive and negative signs accompanied with the secondary instability are formed along the sinusoidal wave pattern of the primary mode. The secondary mode structure with the vortex streets is similar to that typically seen in the K-H type instability, and can cause deformation of the vortex sheet of the primary mode. The vorticity distribution on the ionospheric plane shown in figure 1(b) differs from that found on the equatorial plane, and is not a simple projection from the magnetosphere because of the nonlinear mode coupling. The V-shaped vortex pattern on the ionosphere was also identified in the nonlinear simulation of the feedback instability [16]. Superposition of w to the linear eigenfunction of the primary mode w L is also shown in figures 1(c) and (d) where the vorticity distribution defined by w w +0.25 L is plotted. The secondary mode structure propagates with the primary one in the perpendicular direction to the ambient magnetic field. Therefore, growth of the secondary mode causes deformation of the primary mode structures , and leads to rolling-up of the vortex sheet.
A parameter survey for the primary mode amplitude, A L , shows existence of a critical amplitude. Growth rates of the secondary instability are plotted in figure 2, where the horizontal line shows the primary mode growth rate of As shown in figure 2, the secondary instability grows for >´-A B V L 2.5 10 L 5 0 A , and its growth rate monotonically increases with A L beyond the growth rate of the primary one. It suggests the following scenario for saturation of the feedback instability growth in the nonlinear simulation of the M-I coupling. When the primary mode amplitude exceeds the critical value (that is, 0 A in the present case), the secondary instability is switched on. As A L increases with the  feedback instability growth, the linear growth rate of the secondary instability is enhanced, and can be even larger than that of the primary one. It enables a rapid growth of the secondary instability, and its amplitude becomes comparable to the primary mode. Finally, the secondary mode interacts with the primary one, leading to the saturation of the feedback instability growth. Indeed, a nonlinear growth of ω is found when the averaged amplitude was w á ñ » B V L 0.01 0 A in the nonlinear simulation (see

Nonlinear simulation
The secondary instability investigated in section 3 leads to nonlinear saturation of the feedback instability growth and turbulence transition. The previous simulation of the feedback M-I coupling, however, could not be continued long after the saturation of the primary instability growth due to a numerical difficulty. We have newly developed a spectral code which numerically solves equations (1)-(6) by means of the Fourier spectral method for the perpendicular coordinates while the finite difference is applied to the field-aligned coordinate. Here, nonlinear terms in equation (6), such as f  · n and f { } n , are neglected, assuming small amplitudes of density perturbations | | n n 0  , which greatly simplifies the spectral code. Other nonlinear terms are computed in the real space by means of the fast Fourier transform with the '3/2-rule' for de-aliasing. In the followings, we discuss the nonlinear simulation result for the same physical parameters shown above, where  ( ) 84, 84 Fourier modes are employed for the perpendicular wavenumber space = ( ) k k k , and small random perturbations.
In figure 3 shown are time histories of the total kinetic (electric) and magnetic energy integrated over the simulation domain, that is, ò 2 , respectively. The linear instability growth in an early phase of the simulation is followed by peaks of the kinetic and magnetic energy at » t L V 8000 A , which is consistent with the previous result [16]. The magnetic energy dominates the kinetic one at the peak level. After the saturation of the initial instability growth, the fluctuation energy drops, and then, keeps a nearly constant level, where the instability drive due to the fixed external electric field E 0 balances with the finite dissipation due to the viscosity and the resistivity. In the quasi-steady turbulence state, the magnetic and kinetic energy is almost the same as seen in figure 3. The energy equipartition is a typical feature of the Alfvénic turbulence.
Spatial distributions of vorticity at t=5000, 7000, 9000, and 20 000 L V A are plotted in figures 4(a)-(d), respectively. During the linear growth phase of the feedback (primary) instability, one finds a sinusoidal mode pattern in the x-y plane (see figure 4(a)). When the secondary instability grows, the vortex sheets are rolled-up due to the K-H type mode as shown in figure 4(b). Further deformation of vortex sheets leads to spiral and folding motions, and turbulent fluctuations develop as clearly recognized in figure 4(c). In the quasi-steady state (see figure 4(d)), the Alfvénic turbulence is continuously driven through the M-I coupling with the constant externalÉ B flow.
In this problem, the ionospheric equation with E 0 as well as the initial condition introduces anisotropy on the outer scale with respect to the direction ofk . We confirmed, however, turbulent fluctuation spectra in  where Dk denotes the grid spacing of the discrete wavenumber space. In figure 5, one finds that, even from the initial monotonic spectrum of the linear eigenmode (but with very weak perturbations), the fluctuation spectrum broadens into a wide wavenumber space, and is followed by an exponentially decaying tail on the high  , which are integrated for z and shell-averaged on thek space. k side. The fluctuation spectrum on the lowk side (roughly^ k L 50) seems to be~k 2 or slower, but further studies for cases with much weaker dissipation are necessary for evaluation of the power law spectrum. Nevertheless, the present simulation result clearly demonstrates that the Alfvénic turbulence can be generated through the feedback instability in the M-I coupling system, which is a possible cause of the broad band Alfvén wave fluctuations observed in the auroral magnetosphere.

Summary
The secondary instability driven by the linear feedback (primary) unstable modes and the transition to turbulence have been examined by means of the M-I coupling model. The perturbative (quasilinear) analysis has shown that there is a critical amplitude of the primary mode for the secondary instability growth, and that the secondary mode structure has the vortex streets on the magnetospheric equatorial plane and the V-shaped pattern on the ionosphere. The mode structures and the critical amplitude obtained by the perturbative analysis are consistent with the nonlinear simulation results of the feedback instability in the M-I coupling system.
A nonlinear simulation of the feedback instability has also been carried out by means of the reduced MHD code with the spectral method. Transition to the Alfvénic turbulence through the secondary instability growth is successfully simulated, where the quasi-steady turbulent state is realized through the M-I coupling with a continuous instability drive due to theÉ B convection. The equipartition of the kinetic and magnetic energy is found in the Alfvénic turbulence. While broadening of the energy spectrum of turbulence fluctuations into the higher wavenumber space is observed, more elaborate simulation studies for cases with higher resolution and weaker dissipations are necessary to discuss the inertial subrange and the power law spectrum.
In this paper, we have demonstrated the first simulation result of the auroral turbulence generation through the M-I feedback interaction. It is shown that the Alfvénic turbulence can be generated spontaneously in the M-I coupling system in case with theÉ B plasma convection continuously driven in the magnetosphere, which may provide a plausible explanation to the origin of Alfvénic turbulence observed in auroral regions [4]. A more detailed nonlinear simulation study on the Alfvénic turbulence is currently in progress with higher numerical resolution, and the results will be reported elsewhere.