On the role of ion-scale whistler waves in space and astrophysical plasma turbulence

Competition of linear mode waves is studied numerically to understand the energy cascade mechanism in plasma turbulence on ion-kinetic scales. Hybrid plasma simulations are performed in a 3-D simulation box by pumping large-scale Alfvén waves on the fluid scale. The result is compared with that from our earlier 2-D simulations. We find that the whistler mode is persistently present both in the 2-D and 3-D simulations irrespective of the initial setup, e.g., the amplitude of the initial pumping waves, while all the other modes are excited and damped such that the energy is efficiently transported to thermal energy over non-whistler mode. The simulation results suggest that the whistler mode could transfer the fluctuation energy smoothly from the fluid scale down to the electron-kinetic scale, and justifies the notion of whistler turbulence.


Introduction
Turbulence in space plasmas is fundamentally different from that in ordinary, neutral fluids in that linear mode waves or electromagnetic waves can potentially be a carrier of the fluctuation energy in the spectral domain toward higher wavenumbers.One may speak of weak turbulence if a perturbative approach is valid such that the linear modes play an important role in the energy cascade, and strong turbulence if eddies play an important role in the energy cascade.Another difference of space plasma turbulence from ordinary fluid turbulence is that the inertial range or the energy cascade mechanism can be subdivided into magnetohydrodynamic (MHD) or fluid scale, ion-kinetic scale at wavelengths around the ion inertial length or gyro-radius, and electron-kinetic scale at wavelengths around the electron inertial length or gyro-radius.For example, the transition from the MHD scale to the ion-kinetic scale is about 100 to 1000 km in the solar wind, and that from the ion to electron-kinetic scale is of the order of 10 km.
Naively speaking, there are four fluctuation types that may exist as a linear-mode wave in the ion-kinetic range for a Maxwellian velocity distribution function: whistler, ion Bernstein, kinetic Alfvén, and kinetic slow modes.Figure 1 shows the dispersion relations of these modes at a propagation angle of θ kB = 85 • in a low-beta plasma (β = 0.05) using the numerical algorithm developed by Gary (1993).These four modes have the following relevance or applications in space plasmas.
-Whistler mode exists not only as a right-hand mode in the parallel direction (to the mean magnetic field) in the cold plasma treatment, but also in the quasiperpendicular limit.The whistler mode may be regarded as a kinetic extension of the MHD fast mode ( Gary, 1986).In the dispersion relation diagram, the frequencies increase monotonously as a function of the wavenumbers and extend across the ion gyro-frequency, its harmonics, and the lower-hybrid frequency (at about 43 times higher than the ion gyro-frequency but well below the electron gyro-frequency).There are numerous observations of the whistler mode in various regions in near-Earth space: upstream waves of the Earth bow shock (Narita and Glassmeier, 2005), magnetospheric chorus (Katoh, 2014), magnetotail right-hand waves (Tsurutani and Smith, 1984;Tsurutani et al., 1985), lion roar waves in the magnetosheath and the dayside magnetosphere (Baumjohann et al., 1999(Baumjohann et al., , 2000)), as well as waves at the dayside magnetopause (Vaivads et al., 2007) and in the outflow region from magnetic reconnection (Eastwood et al., 2009).
-Ion Bernstein mode appears as a breakup of the whistler mode at larger propagation angles from the mean magnetic field and are essentially a resonance mode at the ion gyro-frequency (fundamental mode) and harmonics.There are only few observation cases in space plasmas within a frequency range between the fundamental and several harmonics of the ion gyro-frequency such as fluctuations in the solar wind (Perschke et al., 2013(Perschke et al., , 2014) ) and in the magnetic reconnection outflow region (Narita et al., 2016a).The ion Bernstein modes are clearly visible in direct numerical simulations, and offer various kinds of wave-wave interactions for the parametric instability (Jenkins et al., 2013).
-Kinetic Alfvén mode is a kinetic extension of the MHD Alfvén mode, and can also be obtained as a quasiperpendicular limit of the ion cyclotron mode.Various spacecraft measurements in the solar wind indicate the kinetic Alfvén mode at intermediate frequencies (between 0.1 and 100 Hz in the spacecraft frame) of so-lar wind turbulence (Bale et al., 2005;Sahraoui et al., 2010).
-Kinetic slow mode is a kinetic extension of the MHD slow mode (Zhao et al., 2014;Narita and Marsch, 2015) and its existence is suggested by a pressure-balance structure in the solar wind (Yao et al., 2013).Kinetic slow mode is obtained as a quasi-perpendicular limit of the ion acoustic waves and becomes less damped at highly oblique angles at around 85 • and larger.Both kinetic extensions of the Alfvén and slow modes have the lowest frequencies (and therefore represent nearly standing structures) and stay well below the ion gyrofrequency.
We address the importance of the wave dispersion relations in plasma turbulence as they can be the primary channel of the energy cascade mechanism in the inertial range.In this paper, we report that the whistler mode is the most persistent mode that survives both in 2-D and 3-D magnetized plasmas as fluctuations evolve into turbulence.Namely, whistler turbulence is the most appropriate picture to describe plasma turbulence on the ion-kinetic scale.We present a study on the competition of linear wave modes toward turbulence on the ion-kinetic scales.Ion-scale turbulence represents a transition from MHD turbulence to electron-scale turbulence, and may exhibit both wave-wave interactions and wave-particle interactions.The ion-scale spectral domain may therefore be regarded as the dispersive-dissipative range.

Persistence of whistler mode
We use the method of hybrid plasma simulations in a 3-D setup for the following reasons.MHD or Hall-MHD simulations cannot properly resolve wave-particle interactions for ions.Also, particle-in-cell simulations are numerically too demanding to perform 3-D turbulence simulations with a high mass ratio from electrons to ions.Hybrid simulations treat ions as charged particles (strictly speaking, superparticles using the particle-in-cell algorithm) and electrons as a massless, finite-pressure fluid as a charge-neutralizing background.
Our numerical studies follow those presented in Verscharen et al. ( 2012) and Comişel et al. (2013Comişel et al. ( , 2015) ) for a 2-D setup.Using the hybrid plasma simulation code AIKEF (Müller et al., 2011), we solve a set of equations of motion for ions (treated as superparticles in the particle-in-cell algorithm) and the Maxwell equations in a self-consistent way.The equations are time-integrated for 3-D vectorial components (e.g., for particle velocity, electric field, and magnetic field) using a finite-element method in the coordinate space.Electrons are treated as a finite-pressure massless fluid and serve as a charge-neutralizing background.
The simulation box size used in the 2-D run has L ⊥1 = 250d i and L = 250d i in the perpendicular and parallel direc- tions to the mean magnetic field, respectively (Comişel et al., 2015).Here, d i = V A / p denotes the ion inertial lengths for protons.The Alfvén speed V A and the proton gyro-frequency p are used to estimate the inertial length.The 3-D run is different from the 2-D run in that the third direction (the ⊥ 2 axis) is introduced and the box size is extended.The box size for the 3-D run is L ⊥1 = L ⊥2 = L = 128d i in the two perpendicular directions and the parallel direction, respectively.The simulation box is spanned by a 3-D mesh with a size of d i /4 to resolve the ion gyro-motion.The boundary condition is periodic in all directions.
We set in our numerical study, an ion beta of β i = 0.1.This value of beta is a moderate reachable beta in the computational setup that we use.Simulation studies using the particle-in-cell algorithm (for ions in our case) become increasingly more expensive in a numerical sense at higher values of beta because one has to put an even larger number of superparticles to minimize the numerical noise.Furthermore, in the 3-D spatial setup, the number of superparticles required for the simulation is more demanding as the number scales to a cubic law of the mesh point per spatial axis.All the variables relevant in the simulations are normalized using the beta, the magnetic field magnitude, and the Alfvén speed.
As a pump for turbulence excitation, low-frequency electromagnetic waves are isotropically set in the simulation box as an initial condition.In the 2-D setup, the value of ion beta was β i = 0.05 and a number of 20 wave modes have been used, spanning 50 equal propagation directions.The wave amplitudes are derived with the Kolmogorov power-law scaling (i.e., with the spectral index −5/3) up to a cutoff at 20 % of the inertial length wavenumber, kd i ≤ 0.2; see e.g., Verscharen et al. (2012).A number of 18 wave modes are set for the 3-D setup without imposing any power law scaling.Their corresponding wavenumbers are k ⊥1 = 0, ±2π/L ⊥1 , k ⊥2 = 0, ±2π/L ⊥2 , and k = ±2π/L .The initial wave phases are chosen as random.The wave frequencies are derived from the MHD Alfvén waves, ω = ±k V A .No additional pump waves are added during the simulation run, nor is the value of beta reset.The equation of motion and the Maxwell equations are time-integrated in an alternate fashion.The time step of the fast Fourier transformation applied in the time domain is 1 −1 p gyro-period while the time range is about 100 −1 p gyro-periods.The initial spectrum is established such that the magnetic field fluctuations in the coordinate space is 1 % of the mean magnetic field and this value is implicitly assumed in our simulation results unless a higher amplitude of 10 % is specified.Verscharen et al. (2012) discussed the regime of turbulent cascade and the role of the fluctuation amplitude of the initial pumping waves in their study at ion beta 0.05.The authors noticed that the turbulent cascade is more pronounced at larger pumping-wave amplitudes (10 % of the mean magnetic field).The magnetic energy spectrum in the wavenumber-frequency domain shows significant amount of energy only along the dispersion relation for the whistler mode.
The spectral decay in the wavenumber domain for the 3-D simulation run is given in Fig. 2 at three different times: t p = 100, t p = 200, and t p = 300.While the spectrum is rather isotropic in the perpendicular wavenumber domain (Fig. 2, upper row), the spectral decay is flatter in the perpendicular direction and steeper in the parallel direction (Fig. 2, bottom row).
We are aware that at this small 1 % initial amplitude, it is difficult to clearly distinguish the flow of the energy from MHD scales to the kinetic scales.Figure 3a  wavenumber spectrum of the magnetic field fluctuations obtained at time T = 500 −1 p (solid black line).The gray line indicates the spectral curve from a simulation run without any initial pump.The black dashed line represents the slope of the Kolmogorov spectrum (−5/3).At larger wavenumbers, the overlap of the spectra for the initial MHD pumping waves and the thermal noise can have a physical interpretation like a mixture between fluctuations originating in the Alfvén wave excitation and the thermal noise manifested by the solar wind plasma.As an additional comment, in order to minimize the consequences of the numerical noise, we carried out highly demanding computational runs by using more than 10 10 superparticles in the simulation box.
In the 2-D simulation, the energy spectrum develops by showing various dispersion relations.At a time of 600 ion gyro-periods, whistler mode, ion Bernstein modes, and ion cyclotron mode are clearly visible in the wavenumberfrequency spectrum (as a slice of along the k ⊥1 axis) in Fig. 4, obtained by using data from Comişel et al. (2013).This earlier study has found that the frequencies follow first for the linear modes, and then deviate or become broadened from the linear modes.
In the 3-D simulation, in contrast to the 2-D case, the energy spectrum develops primarily by showing the whistler mode branch.Strong ion Bernstein modes no longer appear.The low-frequency modes such as the ion cyclotron mode and the kinetic slow mode cannot be clearly identified because of the broad frequency distribution in a wider range of the wavenumbers.Figure 5a displays a slice of the magnetic energy spectrum along the k ⊥1 axis.The fluctuation energy is axi-symmetrically distributed around the mean magnetic field direction.Only the whistler mode can be mainly identified in the sliced spectrum along the k ⊥2 axis during the simulation run.

Discussion and outlook
Why does only the whistler mode survive during the plasma evolution into turbulence in the 3-D coordinate space and the other modes do not?In fact, one may expect that the ions should be heated by the dissipation of the left-hand modes.In the 3-D simulation, two different scenarios could explain the missing ion Bernstein (IB) modes: (1) the IB modes are excited but the resonance with the ions is so efficient that the energy of the IB modes goes immediately into thermal energy and (2) the IB modes are not excited in the 3-D setup and the ion heating is not significantly occurring.Consequently, we first search for evidence of increasing ion temperature.The total temperature increase in the simulation box at the latest time of the 3-D run (t ∼ 700 −1 p ) is about 8 % of the initial temperature.This temperature increase along the parallel direction cannot account for the particle heating provided by the dissipation of the IB modes.For example, in the 2-D scenario, the perpendicular temperature increase is about 70 % of the initial temperature; see e.g., Verscharen et al. (2012).
For a better clarification of the above proposed scenarios, the evidence for wave damping is investigated by checking the occurrence of the left-hand mode.We apply a method of decomposition of the right-hand (R) and left-hand (L) modes based on the Stokes parameters {I, Q, U, V } from the magnetic field.The decomposition into R and L modes is meaningful for field-aligned and oblique wave propagation.The procedure is described in Appendix A. Each component of the fluctuating field is completed from real into complex values by shifting a phase of 90 • using the Hilbert transformation.The phase information in each component is used to determine the rotation sense of fluctuation.The method is calibrated by exciting the ion cyclotron instability.A 1-D hybrid simulation is performed starting with a high initial temperature anisotropy; see e.g., Gary and Saito (2003).The left-hand mode is assigned to the strongest branch resulting in the decomposition method.Figure 6a shows the energy of the right-hand mode E R (solid black line) and the left-hand mode E L (solid gray line) with respect to the elapsed time in the simulation.Both curves evolve smoothly at the same level until time t p ∼ 150, then they quickly begin growing in magnitude.At about time t p ∼ 300, the left-hand mode becomes clearly stronger than the right-hand mode.The difference between the two modes initially increases but at later times (not shown in Fig. 6.), the left-hand modes slowly achieve a decaying phase.We interpret this result as a coupling between the MHD Alfvén waves and the fast mag-netosonic waves at time t p ∼ 150.The energy is pumped from MHD scale in the system and both the R and L modes rapidly start to grow until time t p ∼ 1000.The oscillations seen at the latest time in both R and L modes are provided by the vibration of the mean magnetic field induced by the MHD alfvénic wave.The oscillations have amplitudes at 1 • from the mean magnetic field direction while the mean magnetic field is now tilted at an angle of 3.5 • .In this regime, the wavenumber-frequency spectrum (not shown here) is closer to that one obtained by Verscharen et al. (2012) by using 10 % amplitude of the pumping waves; the persistent linear wave modes are the whistlers.The other left-hand modes have been decayed.
The right-and left-hand modes obtained by decomposing the magnetic field fluctuations from the 3-D setup are also plotted in Fig. 6a by using the same color convention (black line for R mode, gray line for L mode).E R and E L oscillate due to the initial MHD excitation and smoothly increase at later times due to a weak inclination of the mean magnetic field.Their profiles do not encounter the exaltation of the 2-D analogous modes.The circular polarization of the waves is weaker than in the 2-D setup and the waves are much more linearly polarized.
An additional 3-D run is carried out for a further investigation on the role played by the pumping Alfvén waves for the persistency of the whistler modes.In order to attain a clear energy cascade from the MHD scale to the ion-kinetic scale, an asymmetric box is set with a diminished length on the direction parallel with respect to the mean magnetic field (L = 64d i ).The current number of superparticles in the computational cell is doubled, from 200 superparticles to 400 superparticles, the ion beta parameter is decreased from β i = 0.1 to β i = 0.05, and the amplitude of the initial pumping waves is raised at a value of 10 % from the mean magnetic field.Figure 3b shows, by using the incremental grayscale color nuances, the 1-D reduced power spectra along the perpendicular wavenumber axis at times 0, 100, 200, 300, 400, and 500 ion gyro-periods.At the initial time, the power spectrum is given at larger scales by the solid, light gray line.At a time of 100 gyro-periods, the energy is transported at lower scales by quickly decreasing until a minimum value is reached at the wavenumber kV A / p ≈ 1.At later times (t > 100 −1 p ), the fluctuation level is higher than the noise level while the spectral slope is attending the Kolmogorov value of −5/3.The turbulent fluctuations reach a quasi-stationary state (in terms of the energy spectra) by 500 ion gyro-periods.
The wavenumber-frequency spectrum is shown at a time of 500 gyro-periods in Fig. 5b.The whistler mode is the wave mode solely excited at the perpendicular direction with respect to the mean magnetic field.The other weak modes observed in Fig. 5a (in particular the second harmonic of Bernstein mode) are suppressed.The result of the decomposition method is shown in Fig. 6c.Left-and right-hand modes are clearly excited in the sense that their profiles are uncorrelated during the time evolution of plasma turbulence.Figure 6b shows, for comparison purposes, the result obtained for the 2-D setup using equivalent driving amplitudes.At the time when the energy is about to cascade at smaller scales (t > 100 −1 p ), the strength of both left-and right-hand modes starts to grow.In contrast with the former 3-D setup, Fig. 7 shows a temperature increasing in the perpendicular direction of about 15 % from the initial value.The coupling between oblique whistler modes and oblique Bernstein modes has been proved by Markovskii et al., (2010) to manifest as a proper mechanism for heating the protons in a 2-D hybrid simulation by using initial fast magnetosonic waves.A similar process, seen as an indirect heating by means of the whistler waves, may occur in our second 3-D setup, and then the first assumed scenario at the beginning of this section is available.Nevertheless, whistler mode is persistent irrespective of dimensionality or initial setup because all other modes are excited and immediately damp such that the energy is efficiently transported to thermal energy over non-whistler mode.
We summarize the 3-D simulation results as follows.By applying lower driving amplitude waves, the IB modes or other modes are weakly excited -the opposite of the 2-D result.The decomposition method shows weak evidence of circular polarized waves at parallel and oblique propagation angles with respect to the mean magnetic field.When higher driving amplitude waves are used, the suppression of the linear modes except for the whistler modes is observed, as shown in the 2-D simulation of Verscharen et al. (2012).Actually, the detailed process of the suppression and the appearance of the outstanding whistler mode (fast mode) are still unclear even in the present study using the 3-D simulation results.But it is worth noting that the outstanding whistler waves, which is similar to those in the 2-D simulation of Verscharen et al. ( 2012), also appears in the 3-D simulation, even though the only smallest wavenumber is given as the initial condition of the 3-D setup.This suggests that the appearance of the whistler mode itself is relatively robust for the difference of the initial shape of wavenumber spectra.In contrast, the difference between the 2-D and 3-D setup in the case of -Figure 7. Time evolution of the perpendicular (black) and parallel (gray) temperature of the protons for the 3-D setup using higher driving amplitude waves.lower driving amplitude wave possibly occurs due to the difference of the initial conditions.From this point of view, the difference of temperature in the case of higher driving amplitude waves can also come from the difference of the initial conditions.To clarify these points and scenarios, more parameter studies are necessary.
As an alternative scenario, the azimuthal degree of freedom blocks the development of the electrostatic component such that the IB modes (which have a larger electrostatic component) are suppressed and the whistler mode (which is an electromagnetic mode) can evolve.Ganguli et al. ( 2010) discussed on the 3-D character of the electromagnetic whistler turbulence in the intermediate frequency range above the proton gyro-frequency and below the electron gyro-frequency at low beta plasmas.The authors' opinion is that the wave-particle scattering can convert electrostatic waves with low group velocity into electromagnetic waves with large group velocity.This process can convect energy away from the region with the result of a wave energy loss without significant local particle heating.
The lessons from the 3-D hybrid simulation of plasma turbulence can be resumed as follows.The whistler mode is persistent in hybrid simulation, and could fill the gap of the energy cascade process between MHD turbulence and electron-scale turbulence.The existence of whistlers as the only persistent mode justifies the previously proposed scenario that whistlers are the pump or a major energy carrier for electron-kinetic scale plasma turbulence (Saito et al., 2008(Saito et al., , 2010;;Gary et al., 2012;Chang et al., 2013).In fact, Narita et al. (2016b) has recently reported new findings of the whistler turbulence in the solar wind at electron scales for oblique propagation angles.The other candidate wave modes with quasi-perpendicular wavevectors such as kinetic Alfvén mode, kinetic slow mode, and IB modes are either not excited by wave-wave couplings or strongly damped.The IB mode exists only in a 2-D domain.
We note possible applications of our findings.The persistency of whistler mode motivates us to construct a phenomenological model for whistler turbulence as formulated earlier by Narita and Gary (2010) and Saito et al. (2010).For example, it would be interesting to extend the critical balance hypothesis to the ion-kinetic scale by balancing the eddy turnover time and the whistler wave scatter time.On the other hand, search for the ion-scale whistler and the IB modes is a suitable task to understand the turbulent heating process in the inner heliosphere in the upcoming Solar Orbiter and Solar Probe Plus missions.

Data availability
Data from our 2-D and 3-D hybrid simulations supporting the results presented in this paper are stored at the Institut für Theoretische Physik -Technische Universität Braunschweig.Data can be obtained by writing to the following email addresses: h.comisel@tu-braunschweig.de or comisel@spacescience.ro.

Figure 1 .
Figure 1.Dispersion relations calculated for a propagation angle of 85 • for low-beta plasma (β = 0.05), showing whistler modes (WH), ion Bernstein fundamental mode (IB1) and harmonics (IB2), kinetic Alfvén mode (KA), and kinetic slow mode (KS) under the condition of the damping rate magnitude less than the wave frequency.Dashed line shows ion cyclotron mode (IC) for a propagation angle of 75 • .

Figure 2 .
Figure 2. Wavenumber spectrum of the magnetic field fluctuations for the 3-D plasma simulation during the time evolution into turbulence.
Figure3.(a) One-dimensional wavenumber spectra of magnetic field fluctuations for a cut along the perpendicular-1 axis at time 500 −1 p for the actual run (black) and for a run without initial pump (gray).(b) Reduced magnetic energy spectrum (normalized to their initial value) obtained for higher amplitudes of the pumping waves at a time step of 100 −1 p plotted by an incremental level of gray, starting with time 0 (light gray) and ending at time 500 −1 p (black).As reference, the dashed line is the Kolmogorov value of the spectral slope −5/3.

Figure 4 .
Figure 4. Wavenumber-frequency spectrum for magnetic field fluctuations in a 2-D plasma simulation at a time of 600 ion gyroperiods.

Figure 5 .
Figure 5. Wavenumber-frequency spectrum for magnetic field fluctuations in 3-D plasma simulations at a time of 500 ion gyroperiods.(a) and (b) display slices along the k ⊥1 axis obtained from using lower and higher amplitudes of the initial pumping waves.

Figure 6 .
Figure 6.(a) Energy (normalized to the initial magnetic field energy) of the right-hand mode (black) and the left-hand mode (gray) of the magnetic field fluctuations calculated for the 2-D and 3-D setups using lower amplitudes of the initial pumping waves.(b) and (c) show results obtained from the 2-D and the 3-D setups using higher amplitudes of the initial pumping waves, respectively.