Dispersion relation as a channel of plasma turbulence evolution

Abstract

We present evidence from ion-scale plasma turbulence simulations that normal-mode waves and sideband waves co-exist when a plasma evolves from wave fields into turbulence. Ion-scale fluctuations represent normal-mode fluctuations such as Bernstein and cyclotron modes, accompanied by sideband waves having moderate frequency mismatch with the normal modes. The evolution process is studied using the method of energy partition and frequency broadening. Wave evolution into turbulence exhibits distinct stages: growth phase of normal modes, saturation, and decay. The saturation time is of the order of 1000 ion gyroperiods, which is delayed at higher values of the plasma parameter beta. We conclude that the dispersion relations play the central role in plasma turbulence evolution, serving as the energy flux channel in the spectral domain.

Findings

Space and astrophysical plasmas are collisionless. Their turbulent state, particularly in the ion kinetic regime, is believed to play a fundamental role in transferring energy between an electromagnetic field and a plasma, serving as an effective dissipation mechanism. A useful way to describe plasma turbulence is to view it as a set of normal- and non-normal-mode waves. While normal-mode waves are characterized by dispersion relations associating frequencies with wavevectors and may be sustained for a long time, non-normal-mode waves are not strictly constrained to dispersion relations and may have a relatively shorter lifetime. The wave picture of turbulence should be, naively speaking, valid under a well-established large-scale magnetic field such that perturbative treatment is applicable. In fact, multi-spacecraft observations indicate that the magnetic field fluctuations represent normal modes or linear modes that have been predicted by the linear Vlasov theory (Dudok de Wit et al. 1995; Narita and Glassmeier 2005). Foreshock waves represent the low-frequency whistler mode excited by right-hand ion beam instability (Gary 1991). Magnetosheath fluctuations represent the mirror mode (Gary 1993). Furthermore, fluctuations near the location of magnetic reconnection show characteristics of the whistler mode (Eastwood et al. 2009); the whistler-mode chorus resonates with energetic electrons in the inner magnetosphere (Katoh 2014); and fluctuations in the solar wind show the dispersion relations for kinetic Alfvén waves and ion Bernstein modes (Perschke et al. 2013; Sahraoui 2010).

Of course, not all the fluctuating fields represent normal modes characterized by a dispersion relation, and wave-wave coupling leads to the excitation of sideband waves. A frequency mismatch with the dispersion relation naturally explains the mechanism of wavevector anisotropy in plasma turbulence (Gary 2013). Yet, the spacecraft observations cannot reveal the detailed mechanism of turbulence evolution or information on the time elapsed after pump wave excitation.

Here, we address a question that is essential for understanding the evolution of plasma turbulence - how does a plasma evolve into turbulence? Since the dispersion relations are known to depend on the value of beta, it is natural to anticipate that turbulence evolution might also be influenced by beta. We adopt a numerical approach to answer this question, and use the hybrid plasma code Adaptive Ion Kinetic Electron Fluid (AIKEF) (Müller et al. 2011). This direct numerical simulation extends the earlier works, ref. (Comişel et al. 2013; Verscharen et al. 2012), in two ways: (1) by achieving the longest possible simulation run to resolve turbulence evolution processes and (2) by extending the condition of beta to the highest possible values (up to beta 0.2). Both improvements require substantial computational resources, enabling one to survey systematically the influence of beta on turbulence evolution for, to our knowledge, the first time in the ion kinetic regime. The method using the hybrid code does not assume the existence of normal mode waves in the kinetic regime a priori. Wave fields are produced in the simulation box, and the energy spectra are determined and analyzed in the wavenumber-frequency domain at different times. The concept of energy partition and strength of broadening are used as the analysis tool for the normal-mode waves and sideband waves to track the time evolution of wave fields into turbulence.

Methods

Hybrid simulation

Direct numerical simulation is the best method for studying the role of dispersion relations in turbulence evolution. The AIKEF code has been developed for hybrid plasma simulation, and it solves the set of equations of motion for ions (as particles) and the Maxwell equations self-consistently using the finite-element method. Electrons are treated as a charge-neutralizing massless fluid. No gyration average (or gyro-kinetic treatment) is needed in this simulation. Hence, the hybrid code is advantageous in that one need not assume any statistical properties of the fluctuating field (e.g., Gaussian statistics), wave modes, or propagation directions. The code is disadvantageous with a weakness in that strong electrostatic effects due to charge localization would be immediately canceled out by the mobile electron fluid, but this effect can be safely neglected for our purpose.

Energy conservation under the AIKEF code simulation was verified for freely streaming particles in a periodic simulation box (or domain) in three dimensions with a total length of 32 ion inertial lengths and the ion bulk flow at 8 Alfvén speeds (Müller et al. 2011). Four hundred iterations (or four ion gyroperiods) are needed for one domain crossing. Particle mass (associated with the adaptive mesh refinement, not used in our simulation), particle energy, particle momentum, and electromagnetic energy were recorded at each iteration. These quantities are normalized to their initial values. Energy conservation has been confirmed with errors of less than 1% per gyration. Strictly speaking, the total particle energy decreases at a constant rate of approximately 0.015% per ion gyration, while the electromagnetic energy moderately increases at variable rates of approximately 0.005% per gyration.

We use the ion beta β _i defined as follows (in the SI units):

$$ \beta_{\mathrm{i}} = \frac{2 \mu_{0} n_{\mathrm{i}} k_{\mathrm{B}} T_{\mathrm{i}}}{B_{0}^{\ 2}}, $$

((1))

where μ ₀ is the permeability of free space, n _i is the ion number density, k _B is the Boltzmann constant, T _i is the ion temperature, and B ₀ is the magnitude of the mean magnetic field. The value of beta is set before the simulation run. The magnetic field magnitude is scaled to unity, and all the quantities relevant to the magnetic field magnitude are scaled according to the dimensionless unit. In our simulations, the ion beta is an input parameter. By combining the ion beta with the number density (which is also an input parameter), the temperature is determined from a given set of beta and density values. This temperature value is further used to generate the thermal population of ions for the simulation runs.

Ion-kinetic waves are generated numerically using the hybrid code in the two-dimensional spatial setup spanning the parallel and the perpendicular directions of the large-scale magnetic field. The simulation in this configuration allows one to track the wave evolution into turbulence directly in the spectral domain (spanning the wavenumbers and frequencies) without any additional complication caused by eddies originating from the advective nature of fluid turbulence. Vectorial quantities such as ion velocity, magnetic field, and electric field are treated as three-dimensional. The wave field is generated on the large possible spatial scale in the simulation box with the smallest possible mesh size to achieve the highest spectral resolution and range. The code is applied to a simulation box with an integration size of 250×250 ion inertial lengths (for protons), a mesh size of a quarter of the inertial length, spatially periodic boundaries, and large-scale constant magnetic field in the z-direction (hereafter the parallel direction). Fluctuating magnetic fields are normalized to the large-scale field magnitude. Initially, a superposition of 1000 Alfvén waves are excited on the large scales where the magnetohydrodynamic picture is valid. The initial or pump waves have random phase distribution with the cutoff wavenumber set to 20% of the inertial length wavenumber ($\frac {|\vec {k}| V_{\mathrm {A}}}{\Omega _{\mathrm {p}}} < 0.2$, where Ω _p and V _A denote the ion gyrofrequency and Alfvén speed, respectively). The amplitudes for the pump waves are determined from the model energy spectrum, which is isotropic, and obeys the Kolmogorov power-law scaling with the spectral index −5/3. The total magnetic field fluctuation amplitude (or the standard deviation) is set to 1% of the large-scale magnetic field. During the numerical run, no additional energy is given to the simulation box, and the value of beta is not reset. The equation of motion and the Maxwell equation are solved alternately at a time step of 0.01% of ion gyroperiods for a time length of up to 2,000 gyroperiods. Wave evolution is studied using three simulation runs, each with a different value of beta: β _i=0.05, β _i=0.1, and β _i=0.2. Four-hundred super-particles are set in the computation cell, which corresponds to 8.4×10⁸ particles in the entire simulation box.

Wave analysis

Spectral estimator

The numerical simulation generates fluctuations in plasma and magnetic field at wavelengths smaller than the ion inertial length. Figure 1 displays an example of the magnetic field structure at beta 0.05 at 700 ion gyroperiods. As is seen in the figure, the wave dynamics produces a filamentary structure along the large-scale magnetic field, confirming the previous results of the wavevector anisotropy (Comişel et al. 2013; Verscharen et al. 2012). The fluctuation data are obtained first in the space-time coordinates and are then transformed into the spectral domain spanning the wavevectors and frequencies by using the fast Fourier transform algorithm. The analysis below is limited to the perpendicular direction to the large-scale field.

Figure 1

Snapshot of magnetic field distribution. Snapshot of magnetic field distribution (normalized to the large-scale field magnitude B ₀) in the spatial coordinates spanning the perpendicular distance (the y component) and the parall el dist ance (the z component) in units of the ion inertial length. Data represent ion beta β _i=0.05 and at 700 ion gyroperiods (adopted from ref. (Comişel et al. 2014)).

The fluctuation data are obtained at each evolution time step for the three field components (x and y are in the direction perpendicular to the mean magnetic field. x is in the direction out of the simulation box or plane, and z is in the direction parallel to the mean field) at all the mesh grid points in the two-dimensional space. The simulation box therefore spans by the y and z directions. The obtained data are then Fourier-transformed into the three-dimensional energy spectra as a function of the frequency ω and the two components of the wavevectors k _y and k _z. We are interested in turbulence evolution in the direction perpendicular to the mean magnetic field, as the spectral energy is transported primarily in that direction. The parallel component of the wavevector is therefore set to zero in the data analysis, and we obtain the two-dimensional energy spectrum in the domain spanning the perpendicular components of the wavevectors and frequencies. This implies that one averages the three-dimensional spectra over the parallel spatial coordinate z to obtain results with statistical significance, rather than use the data at any specific value of z.

We limit our wave analysis to the propagation direction perpendicular to the mean magnetic field because most of the magnetic fluctuation energy is stored along the axis of the perpendicular wavevectors, as shown in the previous studies (Comişel et al. 2013; Verscharen et al. 2012). Wave analysis in other directions, parallel or oblique to the mean magnetic field, is also possible and worth pursuing. The present analysis provides information on the higher-order picture of turbulence evolution as most of fluctuation energy is in the perpendicular direction.

Figure 2 shows the two-dimensional spectra for the magnetic field fluctuations E(k _⊥,ω) obtained at an early evolution stage (200 ion gyroperiods in all the three beta cases, top three panels) and a late evolution stage (800, 1200, and 1600 gyroperiods at beta 0.05, 0.1, and 0.2, respectively). The fluctuation energy is investigated by slicing the energy spectra at $\frac {k_\| V_{\mathrm {A}}}{\Omega _{\mathrm {p}}} = 0$ with the accuracy $\left | \frac {\Delta k_\| V_{\mathrm {A}}}{\Omega _{\mathrm {p}}} \right | \le 4.0 \times 10^{-3}$; this accuracy limit stems from the finite simulation box size. The spectra exhibit an organization, that is, the spectral density has peaks along the dispersion relations for the ion Bernstein modes and the oblique ion cyclotron modes, including both the fundamental mode and harmonics in frequencies.

Figure 2

Snapshots of magnetic energy spectra in the spectral domain spanning the perpendicular components of the wavevectors. Snapshots of magnetic energy spectra in the spectral domain spanning the perpendicular components of the wavevectors k _⊥ to the large-scale magnetic field (normalized to the ion inertial scale using the Alfvén speed V _A and the ion gyrofrequency Ω _p) and the frequencies (normalized to Ω _p). Left column shows the snapshots for ion beta 0.05 at 64 ion gyroperiods (top, the earliest possible time for recording the wavenumber-frequency spectrum), 200 gyroperiods (middle, the intermediate state before the transition), and 800 gyroperiods (bottom, after the transition). Middle column shows the snapshots for ion beta 0.1 at time 64 gyroperiods (top), 200 gyroperiods (middle), and 1200 gyroperiods (bottom). Right column shows the snapshots for ion beta 0.2 at time 64 gyroperiods (top), 200 gyroperiods (middle), and 1600 gyroperiods (bottom).

Partition and broadening

For a quantitative analysis of the wave-field evolution, we extract two essential features as scalar quantities out of the spectra. The first is the energy partition of the normal mode m:

$$\begin{array}{@{}rcl@{}} E_{m} = \int_{0}^{\omega_{\text{max}}} \mathrm{d}\omega \int_{k_{\text{min}}}^{k_{\text{max}}} \mathrm{d}k_{\perp} E(k_{\perp},\omega) \delta \left(\omega - f_{m}(k_{\perp}) \right). \end{array} $$

((2))

Here, ω _max denotes the maximum frequency attained in the Fourier analysis $\left (\frac {\omega _{\text {max}}}{\Omega _{\mathrm {p}}} = 2\pi \right)$, and ω=f _m (k _⊥) is the dispersion relation of the m-th mode ω=f _m (k _⊥) evaluated in the kinetic range between $\frac {k_{\text {min}} V_{\mathrm {A}}}{\Omega _{\mathrm {p}}}=1$ and $\frac {k_{\text {max}} V_{\mathrm {A}}}{\Omega _{\mathrm {p}}}=6$. By filtering using the Dirac delta function, the partition quantity indicates the amount of wave energy stored in the normal modes. The second quantity is the frequency broadening around the dispersion relation for the mode m:

$$\begin{array}{@{}rcl@{}} \Delta_{m} = \frac{1}{\Omega_{\mathrm{p}}} \left\langle \int_{\omega_{-}}^{\omega_{+}} \mathrm{d}\omega \left(\omega - f_{m}(k_{\perp}) \right)^{2} E(k_{\perp},\omega) \right\rangle_{k_{\perp}}, \end{array} $$

((3))

where ω ₋ and ω ₊ are the lower and upper limits in the computation of frequency broadening. The moment calculation is conducted for the lower and upper limits such that the slice of the energy spectrum (or frequency spectrum) at a given wavenumber can be regarded as a locally single-peaked distribution around the frequency calculated for the theoretical dispersion relations (so as not to mix different modes in the moment calculation). The exact values are summarized in Table 1. We measure the spectral broadening in the frequency direction so that the dispersion relations are calculated as a function of wavevectors in the linear Vlasov theory. It is, however, worth mentioning that the spectral broadening can be measured in the wavenumber direction as an additional wave analysis tool. The angular bracket 〈⋯ 〉 denotes the operation of averaging over the wavenumbers. The broadening estimator is used as a measure of the strength of the sideband wave associated with the normal mode m. The method with the partition and broadening estimators can be applied to any other relevant mode, but for realizing the objective of this work, we select three representative modes that contain a major fraction of the fluctuation energy: fundamental ion Bernstein mode (IB1), second-harmonics of the ion Bernstein mode (IB2), and ion cyclotron mode (IC). We study the wave-field evolution using the two estimators for these dominant wave modes.

Table 1 Lower and upper limits of the frequency range used in the moment calculation

Results

Wavenumber-frequency spectra

The evolution of the space-time structure of fluctuations is studied by evaluating the two-dimensional energy spectra by spanning the perpendicular components of the wavevectors and the frequencies. Figure 2 shows the spectra at two different evolution stages at plasma beta 0.05 (left column), 0.1 (middle column), and 0.2 (right column). In the early evolution stage (at approximately 200 ion gyroperiods), the fluctuation energy is mostly constrained to a limited frequency range, forming line-shaped structures in the spectral domain that are identified as the realization of dispersion relations for the ion Bernstein mode and the oblique cyclotron mode. Figure 3 shows the dispersion relations derived from the linear Vlasov theory for reference. As wave dynamics proceeds, both fluctuation types (normal-mode waves and sideband waves) attain increasingly larger energy until the fluctuation growth saturates. After saturation, the fluctuation energy is stored more as sideband waves, indicated by the frequency broadening of the spectra around the dispersion relations. The bottom panels of Figure 2 display snapshots in the late evolution stage at approximately 800 gyroperiods at beta 0.05, 1200 gyroperiods at beta 0.1, and 1600 gyroperiods at beta 0.2. The growth of the sideband wave components are particularly visible for ion Bernstein modes (fundamental and second harmonic) at beta 0.1 and 0.2.

Figure 3

Model dispersion relations used for wave mode identification in the wavenumber-frequency spectra in Figure 2 . Dispersion relations are obtained by solving the linearized Vlasov equation for beta 0.1 and quasi-perpendicular wavevectors to the large-scale magnetic field. IB1 denotes the fundamental ion Bernstein mode, IB2 its second harmonic, and IC the ion cyclotron mode.

Normal modes and sideband waves

In order to resolve the transition to sideband wave growth and the saturation effect, a quantitative analysis is applied to the energy spectra by using the method of energy partition and frequency broadening. Three wave modes are investigated in detail: fundamental the ion Bernstein mode (IB1), second harmonic of ion Bernstein mode (IB2), and oblique ion cyclotron mode (IC). These modes contain the largest fraction of fluctuation energy in the spectra.

The left column in Figure 4 displays the time evolution of energy partition for the three modes obtained under different beta conditions. The evolution profile of the energy partition shows the following tendency in all the three modes: the partition increases in the early evolution phase, reaches a maximum value (or peak of the energy partition), and becomes diminished in the late evolution phase. The estimated saturation time is found to be dependent on both the modes and the value of beta.

Figure 4

Time evolution of energy partition for the three modes obtained under different beta conditions. Energy partition E _m (left column) and frequency broadening Δ _m (right column) as functions of time in units of ion gyroperiods for three distinct modes: fundamental ion Bernstein mode (IB1), its second harmonic (IB2), and oblique ion cyclotron mode (IC). The evolution profile for the energy partition and the frequency broadening is determined under the conditions of ion beta β _i=0.05 (dotted line), β _i=0.1 (dashed line), and β _i=0.2 (solid line in gray).

It is found that the peak time is shifted or delayed at higher values of beta. The peak time of the energy partition for the three modes (fundamental Bernstein mode, its second harmonic, and oblique ion cyclotron mode) is summarized in Table 2. The growth profile, however, differs between cases. Both the Bernstein modes exhibit rather monotonous growth of the energy partition toward the peak, while the oblique cyclotron mode exhibits the formation of a plateau followed by a sudden increase in the energy partition (e.g., beta 0.1 and 0.2).

Table 2 The peak time of the energy partition for the three modes

The evolution of the frequency broadening exhibits the same tendency: growth, peak, and decay. The right column in Figure 4 displays the time evolution of the frequency broadening for the three modes. The peak time of the frequency broadening for the three modes is summarized in Table 3. Again, the growth curve can be monotonous (Bernstein fundamental for beta 0.05 or 0.1) or can exhibit an abrupt increase toward the peak (cyclotron mode for beta 0.2).

Table 3 The peak time of the frequency broadening for the three modes

Transition time

The partition and sideband waves appear to peak almost simultaneously. This observation raises the question as to which effect comes first among the peak (or saturation) in the linear-mode growth and the enhancement of sideband waves. We test for the hypothesis that the linear mode energy reaches a peak first and subsequently the sideband wave growth saturates. The quantity Δ t=t _L−t _NL is used to verify this hypothesis. The symbol t _L denotes the transition time for linear waves, and is defined through

$$\begin{array}{@{}rcl@{}} \left. \frac{\mathrm{d} E_{m}}{\mathrm{d}t} \right|_{t=t_{\mathrm{L}}} &=& 0 \end{array} $$

((4))

$$\begin{array}{@{}rcl@{}} \left. \frac{\mathrm{d}^{2} E_{m}}{\mathrm{d}t^{2}} \right|_{t=t_{\mathrm{L}}} &<& 0. \end{array} $$

((5))

The symbol t _NL denotes the transition time for nonlinear (or sideband) waves, and is defined through

$$\begin{array}{@{}rcl@{}} \left. \frac{\mathrm{d} \Delta_{m}}{\mathrm{d}t} \right|_{t=t_{\mathrm{L}}} &=& 0 \end{array} $$

((6))

$$\begin{array}{@{}rcl@{}} \left. \frac{\mathrm{d}^{2} \Delta_{m}}{\mathrm{d}t^{2}} \right|_{t=t_{\mathrm{L}}} &<& 0. \end{array} $$

((7))

Here, t _L and t _NL are measured from the linear mode energy partition and the sideband waves, respectively. The evolution scenario is found to be dependent on the value of beta. The saturation or peak times are evaluated by averaging those of the three wave modes and are plotted as functions of beta in the top panel of Figure 5. Though only three data points are used in this panel, the delay of saturation time for the linear-mode fluctuations at higher values of beta can be clearly seen. A closer inspection yields the result that the effect of beta on the delay can well be fitted by a power law as $t_{\mathrm {L}} \propto \beta _{\mathrm {i}}^{\alpha }$ with the exponent α≃0.59.

Figure 5

Beta dependence and the time difference between the peak of partition and that of broadening. Beta dependence of the peak time (or saturation time) t _L of the energy partition for the linear mode fluctuations (top panel) and the time difference between the peak of partition and that of broadening Δ t (bottom panel). The peak times are averaged over the three wave modes shown in Figure 4 and are given in units of ion gyroperiods. Error bars include the variation among the three wave modes.

The transition time is also evaluated for the peak times in the broadening. As in the case for the saturation time, the broadening peak time is dependent on beta: the peak time is increasingly delayed with higher values of beta. The difference in the two times (one is the saturation time in linear-mode fluctuations and the other is the peak time in broadening) is of the order of 50 gyroperiods. The bottom panel in Figure 5 displays the time difference Δ t=t _L−t _NL (here the subscripts L and NL denote the linear-mode and the sideband or nonlinear waves, respectively). The saturation time data are again averaged over the three wave modes. At lower values of beta (β _i=0.05 and β _i=0.1), the broadening saturates first and then the linear-mode fluctuation saturates by about 50 gyroperiods. At a higher value of beta (β _i=0.2), the order is reversed and the broadening peak time arrives approximately 50 gyroperiods after the saturation of linear-mode fluctuations.

Discussion and conclusions

The dispersion relation plays a central role in wave evolution into turbulence. We conclude that waves evolve into a more turbulent state by exciting both normal-mode waves and sideband waves. Both fluctuation components serve as energy reservoirs. The existence of a peak in the energy partition is suggestive of a scenario in which a turbulent field cannot be constructed by normal modes only, as sideband waves arise naturally with normal-mode waves. It is also worth mentioning that beta influences or regulates the peak time or saturation time; the higher the value of beta is, the later the transition arrives. Beta dependence can be described in such a way that the transition is earlier and clearer when the value of beta is lower. In other words, cold plasmas show a clearer and sharper transition from the early stage to the late stage without the formation of a plateau (neither in the energy partition nor in the frequency broadening). Warm plasmas, in contrast, have delayed transition after the plateau formation.

In our simulation, the particle velocity distribution (of protons) is set to Maxwellian, and furthermore, the system is nearly homogeneous in that the mean magnetic field and plasma density can be regarded as nearly constant. The small-scale fluctuations on the kinetic scales (at wavenumbers higher than that of the ion inertial length) cannot be excited by linear instability processes - neither micro-instability caused by gradients in the velocity coordinate nor macro-instability caused by gradients in the spatial coordinate. To verify how exactly the fluctuation power is transferred to small-wavelength waves by wave-wave interactions (Bernstein and cyclotron modes), a higher-order statistical analysis is needed such as bispectral or bicoherence analysis for three-wave couplings. In the present study, a Maxwellian distribution is used for the initial velocity distribution of protons, and the evolution of wave spectra is associated with not instability but nonlinear wave coupling. This situation is different from that treated in the quasi-linear theory of wave-particle interactions, which is deeply related to the evolution of plasma wave instability.

It is true that ions are heated as waves evolve into turbulence, and in that sense, the value of beta is not strictly constant. In fact, our preceding work using the AIKEF code has already shown that the ion temperature is enhanced by a factor of approximately 1.7 preferentially in the direction perpendicular to the mean magnetic field (Verscharen et al. 2012). The possible effect of the temperature increase is in the determination of the exact scaling or relation between the value of beta and the saturation time. However, the qualitative nature of the relation still holds because all the input values of beta are merely shifted. Tracing the time evolution of ion velocity distribution is an important task that is not feasible with the current computation capacity. In our simulations, computational loads become increasingly demanding as we perform simulations at higher values of beta and for longer runs. To stabilize the simulation run, increasingly more particles must be put into the simulation box to compensate for the increasing mobility of particles (that can easily move from one mesh cell to another), which accordingly takes a huge amount of memory in a high-performance computer.

It is important to emphasize that particle distribution functions should be studied in detail to advance the knowledge derived from our analysis. A quasi-stationary state for wave-particle interactions is obtained in the quasi-linear theory by using the so-called plateau equation, addressing that the particle velocity distribution function exhibits pitch-angle plateau formation such that the energy is no longer transferred between electromagnetic fields and particles (e.g., ref. Marsch and Bourouaine (2011)). It might be possible that the plateau in the time history of the wave spectral energy is associated with the plateau formation in the distribution function.

How much does inhomogeneity affect our analysis? Qualitatively, there is always an effect of inhomogeneity in wave fields by considering the wave field as the source of inhomogeneity, which could introduce a certain amount of frequency deviation from the dispersion relation for the homogeneous medium. The dispersion relation in the weakly nonlinear treatment including the effect of finite amplitude A (e.g., magnetic field magnitude) is given by (cf. Equations 11.19 to 11.20 in the book of Treumann and Baumjohann (1997))

$$ \omega = \omega_{0} (k) + \left.\frac{\partial \omega}{\partial \left(A^{2}\right)} \right|_{0} A^{2} + {\mathcal O}\left(A^{4}\right) $$

((8))

We use the magnetic field magnitude for the amplitude A. Then, we estimate the frequency mismatch for the inhomogeneity (δ ω)_inh and compare it with the measured broadening, Δ _m ∼0.01 at most. The frequency deviation associated with the inhomogeneity is evaluated as

$$ (\delta \omega)_{\text{inh}} \sim \frac{\delta \omega}{\delta (A^{2})} A^{2}. $$

((9))

The field-amplitude variation is taken from the spatial distribution of the magnetic field magnitude, and it is of the order of 15%. The squared-amplitude variation is approximately δ(A ²)≃(1.15)²−(1.0)²≃0.33, while the squared amplitude itself is of the order of unity, A ²∼1. Hence, the presence of inhomogeneity amplifies the frequency modulation δ ω by a factor approximately 3.3. The modulation is estimated from the spatial scale convected by the propagating wave packets. We estimate the frequency modulation as δ ω∼v _gr δ k, where the typical group speed is estimated as v _gr∼0.01V _A and the wavenumber of inhomogeneity is estimated for the system size as $\frac {\delta k V_{\mathrm {A}}}{\Omega _{\mathrm {p}}} \sim 0.01$. The group speed is so low because the frequencies of Bernstein or cyclotron modes asymptote to harmonics or subharmonics of the cyclotron frequency. The frequency broadening due to the inhomogeneity is therefore evaluated as

$$ \frac{(\delta \omega)_{\text{inh}} }{ \Omega_{\mathrm{p}} } = \frac{v_{\text{gr}} }{ V_{\mathrm{A}} } \frac{\delta k V_{\mathrm{A}} }{\Omega_{\mathrm{p}} } \frac{A^{2} }{ \delta (A^{2}) } \sim 3.3 \times 10^{-4} $$

((10))

The measured frequency broadening in the simulation data is in the range

$$ 10^{-4} \le \Delta_{m} \le 10^{-2}. $$

((11))

The presence of inhomogeneity explains the lower limit of the measured broadening.

The detailed evolution profile is diverse. In the early evolution phase, the energy can be stored into the normal-mode waves until a peak is attained. This effect is suggestive of the existence of saturation in exciting the normal mode waves. We propose that dispersion relations should be considered as channels of the energy flux in the spectral domain. The late-stage spectra with enhanced broadening (e.g., bottom panels in Figure 2) are suggestive of a scenario in which fluctuations evolve first as a set of normal mode waves and then into sideband waves. However, this scenario should be regarded as the lowest-order picture in describing wave evolution into turbulence, since not only the partition but also the frequency broadening exhibits decay after saturation. If fluctuations were excited only as sideband waves after the peak, the magnitude of broadening should exhibit a monotonous increase, rather than decrease. A possible reason for the decay of broadening is the wave damping due to cyclotron or Landau resonance processes.

A recent study of the parametric decay instability in laboratory plasmas shows that the three wave coupling is well supported or mediated by dispersion relations, and furthermore, the existence of harmonic branches (such as in the Bernstein mode) even allows many open channels of the spectral energy transfer (Jenkins et al. 2013). In this picture, the dispersion relations of normal-mode waves serve as the stations of wave-wave coupling, and fluctuations are anticipated to live longer when excited as normal-mode waves. To construct a higher-order picture of wave evolution into turbulence, dispersion relations and wave damping need to be taken into account for wave-wave interactions and wave-particle interactions.

A potential application of our discoveries to spacecraft observations is the a remote sensing of wave excitation regions. If the energy partition and frequency broadening are evaluated from observational data, the time elapsed from pump wave excitation can be constrained. The transition time found in our numerically generated waves is of the order of 1000 ion gyroperiods. This time scale corresponds roughly to the convection spatial scale of 400,000 km or 63 R_e (in units of Earth radii) in interplanetary space (e.g., solar wind and foreshock) when using the typical values Ω _p∼1 rad/s for the gyrofrequency and V∼400 km/s for the flow speed. This spatial scale of approximately 63 R_e is by far larger than that of the foreshock or magnetosheath regions, but it may be of significance in the solar wind observation.