© ESO 2018

1 Introduction

The loss-cone type of distribution of hot electrons superimposed on the distribution of much denser and colder background plasma is unstable due to the double plasma resonance instability (Zheleznyakov & Zlotnik 1975; Melrose & Dulk 1982; Zaitsev & Stepanov 1983; Winglee & Dulk 1986; Ledenev et al. 2001; Zlotnik 2013; Karlický & Yasnov 2015). This instability generates the upper-hybrid waves, which have their maxima close to the gyro-harmonic number s = ω_UH∕ω_ce, where ${\omega }_{\text{UH}}^2={\omega }_{\text{pe}}^2+{\omega }_{\text{ce}}^2$ and ω_UH, ω_pe, and ω_ce are the upper-hybrid, electron-plasma, and electron-cyclotron frequency, respectively. Owing to this property, the double plasma resonance (DPR) instability is used in models of solar radio zebras (Zheleznyakov & Zlotnik 1975; Winglee & Dulk 1986; Ledenev et al. 2001; Treumann et al. 2011; Zlotnik 2013; Karlický & Yasnov 2015; Yasnov & Karlický 2015; Benáček et al. 2017).

The resonance condition for the double plasma resonance instability in the relativistic case can be expressed as $${\omega }_{\text{UH}}^2={\omega }_{\text{pe}}^2+{\omega }_{\text{ce}}^2$$(1)

where k = (k_∥, k_⊥) is wave vector, u = (u_⊥, u_∥), u_⊥ = p_⊥∕m_e, and u_∥ = p_∥∕m_e are the hot electron velocities perpendicular and parallel to the magnetic field; m_e is the electron mass, $\gamma =\sqrt{1+u^2/c^2}$ is the relativistic Lorentz factor, s is the gyro-harmonic number, and c is the speed of light. In theoretical models studying the double plasma resonance instability, the distribution of hot electrons is usually taken as the Dory–Guest–Harris (DGH; Dory et al. 1965) distribution for j = 1 in the form $${\omega }_{\text{UH}}^2={\omega }_{\text{pe}}^2+{\omega }_{\text{ce}}^2$$(2) where v_t we call the thermal velocity of hot electrons, although the distribution function in the relation (2) is not Maxwellian. The distribution of the background plasma is taken as a Maxwellian one. The combination of both these distributions is considered to be the prototype distribution generating radio zebras (Winglee & Dulk 1986).

In our previous paper (Benáček et al. 2017) we studied the double plasma resonance instability analytically. We showed that the maxima of growth rates of the upper-hybrid waves are shifted to lower ratios of ω_UH∕ω_ce and the contrast between maxima and minima of the growth rate decreases as the temperature of hot electrons increases. On the other hand, when the temperature of the background plasma is increased, the contrast remains the same.

In studies of solar radio zebras (e.g. Ledenev et al. 2001), the frequencies of the zebra stripes and the contrast of these stripes to the background continuum are analyzed. The frequencies are used for the determination of the magnetic field in zebra radio sources. Therefore, it is important to know the relation between the zebra stripe frequencies and the gyro-harmonic number s ≈ ω_UH∕ω_ce. Furthermore, the contrast of the zebra stripes is believed to be connected with the temperature of hot electrons (Yasnov & Karlický 2004). Moreover, some zebras have many stripes and their frequencies correspond to high gyro-harmonic numbers (sometimes >20; Karlický & Yasnov 2015). All these facts require computations in a very broad range of the gyro-harmonic number s, which are difficult to make with the spatially extended multi-mode model. Fortunately, we found that the maximal growth rates and saturation energies computed in the multi-mode model agree very well with that computed in the spatially limited specific-mode model. This agreement is caused by the fact that the wave with the maximal growth rate becomes very soon dominant over other growing waves (due to its exponential growth). Using this specific-mode model, which considerably saves computational time, we compute the maximal growth rates and saturation energies of the DPR instability in a broad range of parameters. Because in the analytical analysis of growth rates of the DPR instability several assumptions were made, these numerical computations serve to verify the analytical results.

The paper is structured as follows. In Sect. 2, we describe our numerical model and initial conditions for the studied double plasma resonance instability. In Sect. 3, we present the results for different model parameters. Discussion of these results and conclusions are in Sect. 4.

2 Model

We use a 3D particle-in-cell (PIC) relativistic model (Buneman & Storey 1985; Matsumoto & Omura 1993; Karlický & Bárta 2008) with multi-core message passing interface (MPI) parallelization. Further details can be found in Matsumoto & Omura (1993, p. 67–84) and on the link below¹ . In the present article we use this model in two versions: (a) the spatially extended model with many wave modes (multi-mode model), and (b) spatially limited model with specific wave mode (specific-mode model), which is used for relatively fast computing in the broadrange of model parameters. The model size in x-, y-, and z-directions is 128 × 60Δ × 128Δ for the multi-mode model, and λΔ × λΔ × 32Δ for the specific-mode model, respectively, where Δ = 1 is the grid size and λ is the wavelength of the specific upper-hybrid wave. For chosen gyro-harmonic numbers and plasma temperatures, we fit the size of the specific-mode model to find the wave mode with the maximal growth rate. In all models we use periodic boundary conditions in all three directions.

The model time step is d t = 1, electron plasma frequency ω_pe = 0.05, initial magnetic field is in the z-direction, and electron cyclotron frequency ω_ce varies from 0.38 ω_pe to 0.056 ω_pe approximating the gyro-harmonic numbers s = 3–18. For dependencies of the growth rate on temperatures, we use s in the range of s = 3–7. Higher values of s up to 18 are used for calculating the saturation energies of the generated upper-hybrid waves, which is needed for zebras with many zebra stripes. In the model we use two groups of electrons: cold background Maxwellian electrons with the thermal velocity v_tb = 0.03–0.05 c, corresponding to the temperature in the interval 5.4–14.8 MK, and hot electrons with the DGH distribution having the velocity v_t = 0.15–0.3 c. Higher values of the background plasma temperatures are given by the requirements of the PIC model. However, as known from the previous analytical study (Benáček et al. 2017) and shown in the following, variations of the background plasma temperature have little effect.

The electron density of background electrons was taken as n_e = 960 per cell and the ratio of background to hot electrons was taken as n_e ∕n_h = 8 with some exceptions mentioned in the following. The number of protons is the same as the number of electrons and their temperature is always the same as that of background electrons. The proton-to-electron mass ratio m_p ∕m_e is 1836.

To find the maximal growth rate for specific physical parameters, we changed the model size λΔ in the perpendicular direction to the magnetic field. Namely, we fit the wavelength of the upper-hybrid wave with the maximal growth rate to a distance of grids used in the PIC model. Thus, we get the maximal simulated growth rate. We found that the optimal model size does not change when computing with different ratios of the electron plasma and electron cyclotron frequencies, but is changed by changing the temperature.

At the beginning of DPR instability, the electric wave energy grows exponentially. We fit this part of the electric energy evolution by the function W_UH = A₀ exp(2Γt), where A₀ is the initial electric energy, t is time, and Γ is the growth rate in model units. Finally we expressed the growth rate of the upper-hybrid waves in the ratio to the upper-hybrid frequency.

3 Results

Firstly we compute an evolution of the DPR instability and generation of the upper-hybrid waves in the multi-mode models with the parameters summarized in Table 1. We show an example of the time evolution of the electric energy density in the parallel and perpendicular directions to the magnetic field lines in Model 2M normalized to the kinetic energy density of hot electrons for the maximal growth rate and the gyro-harmonic number s = 6 in Fig. 1. As can be seen here, the electric energy density in the perpendicular direction dominates over the parallel energy density by at least two orders of magnitude, which indicates the generation of the upper-hybrid waves.

Examples of the time evolution of the energy density of the upper-hybrid waves normalized to the kinetic energy density of hot electrons for three values of n_h ∕n_e, the gyro-harmonic number s = 6, the maximal growth rate, and Model 2M are presented in Fig. 2. The growth rates, estimated in the early stages of the evolution, are written in this figure. As seen here, these growth rates are proportional to the density of hot electrons, in agreement with the paper by Yasnov & Karlický (2004). Furthermore, the normalized energy density of the upper-hybrid waves in all three models converges to the same saturation energy, which indicates that the saturation energy of the upper-hybrid waves is also proportional to the density of hot electrons. The kinetic energy density of hot electrons E_k,hot depends on the plasma density of hot electrons.

Fig. 1 Time evolution of the electric energy density in the parallel and perpendicular directions to the magnetic field lines in Model 2M (Table 1), normalized to the kinetic energy density of hot electrons for the maximal growth rate and the gyro-harmonic number s = 6.

Fig. 2 Time evolution of the energy density of the upper-hybrid waves normalized to the kinetic energy density of hot electrons for three values of nh ∕ne and for the maximal growth rate and s = 6. Other parameters are the same as in Model 2M.

3.1 DPRinstability in detail

To see the processes during the DPR instability in detail, we make a comparison of the distributions in the initial state and at ω_pe t = 1000 for Models 2M and 4M (Table 1; see Fig. 3). Comparing changes of the distribution in these models (left and right columns in Fig. 3), we can see that the instability with different model parameters causes different changes of the distribution function. The distribution function does not change only in one point of the phase space, but in the subspace defined by resonance ellipses corresponding to the range of k_∥ . In Fig. 3, in the third row, we show these resonance ellipses in the phase space, where the changes are dominant. We also show how the resonance ellipses shift across the distribution function in dependance on k_∥ (see the arrows). In both models, the loss-cones of the distributions are step by step fulfilled by electrons (see the red regions in the third row of Fig. 3) and thus the distributions become closer to the Maxwellian distribution. These changes are due to an increasing energy level of the upper-hybrid waves.

We analyzed time evolution of energies in waves with the k-wave vectors perpendicular to the magnetic field in Model 2M (Fig. 4). As seen in this figure, the interval of k-vectors is relatively narrow and remains practically the same duringthe instability evolution. It shows that the energy of the upper-hybrid waves is concentrated in the relatively narrow interval of k-vectors. As will be shown in the following, it enables us to use the specific-mode models.

Fig. 3 Changes of the electron distribution functions of hot electrons during the DPR instability. Left column: Model 2M with s = 6. Right column: Model 4M with resonance s = 5. First row: the distribution at the initial state. Second row: the distribution at ωpe t = 1000. Distributions are normalized to their maximal value. Third row: the difference of the distributions between the initial state and at ωpe t = 1000. Red regions are with enhanced densities and blue ones are with reduced densities. Elliptical lines show the resonance ellipses for a given k-wave vector along magnetic field. Black arrows show that the ellipses shift with increasing value of k∥ .

Fig. 4 Time evolution of energies in waves with the k-wave vectors in perpendicular direction to the magnetic field lines in Model 2M.

3.2 Comparison of specific-mode and multi-mode model saturation

Because our main objective is to determine the growth rates and saturation energies in a broad range of parameters (which is importantfor the interpretation of observed zebras), we search for ways to accelerate computations. We decided to use the specific-mode models that are much faster than the multi-mode ones. To justify their use, we compare the growth rate and saturation energies for all multi-mode models shown in Table 1 with the specific-mode models with the same physical parameters. An example of this comparison for Model 2M is shown in Fig. 5. While in the multi-mode model we use the model size 128Δ× 60Δ × 128Δ, in the specific-mode model in this case we use the model size 12Δ × 12Δ × 32Δ. This specific-mode model covers interval of k_⊥ c∕ω_pe above 5.23 (see Fig. 4). This means that this specific-mode model covers the k-vector waves that are important for the DPR instability. This explains why the results from the multi-mode and specific-mode models are very similar (see the following).

As seen in Fig. 5, the growth rate and saturation energy in both the models agree very well. The same result is found for all other models according to Table 1. This agreement is caused by a dominance of the wave with the maximal growth rate (exponential increase) during an evolution of the double plasma resonance instability. We utilize this agreement in the following computations of the growth rates and saturation energies in a broad range of the model parameters. The agreement between numerically and analytically computed growth rates, shown in the following, also justifies the use of the specific-mode models in our case.

Fig. 5 Time evolution of the upper-hybrid wave energies for the multi-mode model 2M and specific-mode model with the same model parameters.

3.3 Effects of temperatures of hot and background plasma electrons on the growth rate

In the following computations with the specific-mode models, we use the model parameters as shown in Table 2. Figure 6 presents the effects of the temperature of hot and background electrons on the growth rate. Every point in this figure is computed for at least three model sizes of the specific-mode models. Error bars are estimated from their fit and the probability that the same initial parameters give the same result. Models with higher temperatures have lower errors because generated waves are more dominant over the background noise.

Growth rates have the maxima that are shifted to frequencies lower than those given by the gyro-harmonic number s (Fig. 6, upper part and Fig. 7) and this shift increases with increasing temperature. The ratio between maximal and minimal growth rates for each temperatureis up to one order and increases with the decreasing temperature of hot electrons. These effects are in agreement with the analytical results (Benáček et al. 2017). However, contrary to the analytical results, the growth rate for v_t = 0.3 is higher than that for lower temperatures. We think that this difference is caused by differences in the analytical and numerical approach. While in the analytical approach we calculate the growth rate only for one specific k-vector wave, in numerical computations k-vectors in some interval are always in operation. The effect of the temperature of background plasma electrons is small (Fig. 6, bottom part), which agrees with the analytical theory (Benáček et al. 2017). We also determined the shifts of maxima of growth rate (Fig. 7). They increase with increasing s and also depends on temperature of hot electrons.

Fig. 6 Growth rates in dependence on ωUH ∕ωce for Models 1S–4S. Top: plots for three different temperatures of hot electrons vt . Bottom: plots for two different temperatures of background plasma electrons vtb .

Fig. 7 Shifts of frequency of growth rate maxima from frequency given by gyro-harmonic number s.

3.4 Effect of the hot-cold plasma density ratio on the frequency ofthe growth rate maximum

As shown in Fig. 8, the frequency of the growth rate maximum also depends on the ratio of the background and hot electron densities n_e ∕n_h . This dependence is exponential as shown by its fit and converges to ω_UH ∕ω_ce = 5.1 for high values of the n_e∕n_h ratio. Namely, low values of the n_e∕n_h ratio shift the resonance of the DPR instability and thus the frequency of the growth rate maximum. In the analytical analysis it is supposed that the density of cold electrons is much greater than that of hot electrons (n_e ≫ n_h).

In most of our simulations the density ratio n_h ∕n_e is 1:8 in order to keep a low numerical noise. However, such dependencies enable us to extrapolate our results to much lower density ratios, for example to 1:100 as usually considered in the zebra interpretation (see the following).

Fig. 8 Frequency of growth rate maximum ${({\omega }_{\text{UH}}/{\omega }_{\text{ce}})}_{\text{Max}}$ as function of ratio of background and hot electron densities ne ∕nh for parameters in Model 2S and s = 6.

3.5 Comparison of the numerical and analytical growth rates

We compared the simulated growth rates with the analytical ones presented in the paper by Benáček et al. (2017). The comparison is shown in Fig. 9. It was made as follows. We chose the upper-hybrid frequency and hot electron density as ω_UH = 7 × 10⁹ s⁻¹ (f_UH= 1.11 GHz) and n_h =10⁸ cm⁻³. We changed the growth rate plot expressed in dependence on ω_ce (Benáček et al. 2017) to that in dependence on ω_UH . Then, we transformed the simulated growth rates, computed for n_e ∕n_h = 8, to those with n_e∕n_h = 100. The transformation was made in two steps. We shifted the simulated growth rate maxima toward lower values of ω_UH ∕ω_ce accordingtothe relation expressed in Fig. 6. Finally, the growth rates were multiplied by the factor 8∕100 using the linear relation between the growth rate and density of hot electrons (Yasnov & Karlický 2004). As shown in Fig. 9, there is good agreement between the numerically and analytically computed growth rates.

Fig. 9 Comparison of the simulated and analytically derived growth rates for vt = 0.2 c and vtb = 0.03 c.

3.6 Saturation energy levels of the upper-hybrid waves in a broad range of gyro-harmonic numbers

In all models, we determined the saturation energy of the upper-hybrid waves. The saturation energy is reached when the growth rate is balanced by nonlinear effects. The time when the saturation begins depends on the model parameters. In our cases, the saturation time was in the time interval ω_pe t_satur = 300–1200. Mostly, the maximal growth rate leads to maximal saturation energy, but this is not a rule. In some cases, the saturated energy is higher in the model with lower growth rate. In other cases, the saturation energies have recognizable maxima, while the growth rate is nearly constant in the broad range of ratios of ω_UH ∕ω_ce.

We computed one set of models changing the gyro-harmonic number s from s =3 up to s = 18 and keepingthe same temperatures as in Model 2S. As shown in Fig. 10, the saturation energy level in local maxima decreases with increasing s and converges to the value W_UH ≈ 1.6 × 10⁻³E_h,kin (see the exponential fit as presented in Fig. 10). This exponential function can be used for an extrapolation of the saturation energies for even higher parameter s, which is useful for zebras with many stripes.

Maximal values of the saturated energy levels of theupper-hybrid waves for low values of s are about one percent of the kinetic energy of hot electrons. The saturated energy levels increase with an increase in the kinetic energy (temperature) of hot electrons.

Similarly to the growth rates, the saturation energies (Fig. 11) have the maxima and minima shifted to lower ratios of ω_UH∕ω_ce than those given by the integer values of the gyro-harmonic number s. Figure 12 shows these shifts. The frequency of the growth rate and saturation energy maximum is not generally the same. The saturation energy maxima are usually less shifted to lower frequencies than the growth rate maxima. However, contrary to the growth rates, the contrast between the maxima and minima increases with the increasing temperature of hot electrons (Fig. 11, upper part). Similarly, in the models that vary the temperature of background plasma electrons, the contrast between the saturation energy maxima and minima is higher than that for the growth rates (Fig. 11, bottom part).

Fig. 10 Maximal saturation energies of the upper-hybrid waves in dependence on the gyro-harmonic number s for Model 2S with the exponential fit.

Fig. 11 Saturation energies as a function of ωUH ∕ωce . Top: plots for three different temperatures of hot electrons vt . Saturation maxima are shifted to lower ratios ωUH ∕ωce . Bottom: plots for two different temperatures of background plasma electrons vtb .

Fig. 12 Shifts of the frequency of the saturated energy maxima from the frequency given by the gyro-harmonic number s.

4 Discussion and conclusions

Using the 3D PIC model in two versions (multi-mode and specific-mode models) we computed not only the growth rates of the double plasma resonance instability, but also saturation energies of the generated upper-hybrid waves. We described details of the DPR instability by showing how the distribution function of hot electrons changes during the DPR instability. We found that the growth rate as well as the saturation energy are proportional to the density of hot electrons.

Owing to many assumptions made in our previous analytical study, we checked the analytical results using the present numerical models. We found a very good agreement between the numerical and analytical results. This justifies a use of the specific-mode models in this case. We confirmed that the growth rate maxima are shifted to lower frequencies in comparison with those given by the gyro-harmonic number s. This frequency shift increases with an increase of the temperatureof hot electrons, in agreement with the analytical results. We confirmed that the contrast between maxima and minima of the growth rate decreases with the increasing of the hot electron temperature. On the other hand, the temperature of the background plasma has only a small effect on the growth rates.

We found that the frequency of the growth rate maximum depends on the ratio of the background and hot electron densities n_e ∕n_h. We used this dependence to extrapolate our numerical results, made mostly for the ratio n_e ∕n_h = 8, to those with the ratio n_e∕n_h = 100 considered usually in the analytical studies of zebras. Using this dependence in detailed comparison of the analytical and numerical growth rates, made in the interval of ω_UH∕ω_ce = 3–7 and the parameters considered in zebras, we found very good agreement.

We think that some small differences between the numerical and analytical results are caused by differences in the two methods. In numerical simulations, as in reality, the DPR instability works in some regions of the k-vector space, not only with one k-vector as assumed in the analytical theory. Furthermore, in some of the present simulations we found deviations from the assumption made in the analytical approach (k_⊥≫ k_∥).

We computed the saturation energies of the generated upper-hybrid waves, which is beyond the possibilities of the analytical theory. We compared the results of the multi-mode and specific-mode models considering the same physical parameters. We found very good agreement between the results from both types of model. This agreement is caused by a dominance of the wave with the maximal growth rate during an evolution of the double plasma resonance instability. Therefore we used the specific-mode models, which considerably save computational time, for the computation of the saturation energies in a broad range of the parameter s. A use of the specific-mode models is also justified by a very good agreement between the growth rates computed numerically and analytically. We found that the saturation energy decreases with increasing s and this decrease is exponential. This exponential dependance of the saturation energy can be used for an extrapolation of the saturation energies for the parameter s even greater than 18, which is useful for the interpretation of zebras with many zebra stripes (Karlický & Yasnov 2015; Yasnov & Karlický 2015).

Acknowledgements

We acknowledge support from Grants 16-13277S and 17-16447S of the Grant Agency of the Czech Republic. Computational resources were provided by the CESNET LM2015042 and the CERIT Scientific Cloud LM2015085, provided under the programme “Projects of Large Research, Development, and Innovations Infrastructures”.

References

All Tables

All Figures

	Fig. 1 Time evolution of the electric energy density in the parallel and perpendicular directions to the magnetic field lines in Model 2M (Table 1), normalized to the kinetic energy density of hot electrons for the maximal growth rate and the gyro-harmonic number s = 6.
In the text

	Fig. 2 Time evolution of the energy density of the upper-hybrid waves normalized to the kinetic energy density of hot electrons for three values of nh ∕ne and for the maximal growth rate and s = 6. Other parameters are the same as in Model 2M.
In the text

Fig. 3 Changes of the electron distribution functions of hot electrons during the DPR instability. Left column: Model 2M with s = 6. Right column: Model 4M with resonance s = 5. First row: the distribution at the initial state. Second row: the distribution at ωpe t = 1000. Distributions are normalized to their maximal value. Third row: the difference of the distributions between the initial state and at ωpe t = 1000. Red regions are with enhanced densities and blue ones are with reduced densities. Elliptical lines show the resonance ellipses for a given k-wave vector along magnetic field. Black arrows show that the ellipses shift with increasing value of k∥ .

In the text

	Fig. 4 Time evolution of energies in waves with the k-wave vectors in perpendicular direction to the magnetic field lines in Model 2M.
In the text

	Fig. 5 Time evolution of the upper-hybrid wave energies for the multi-mode model 2M and specific-mode model with the same model parameters.
In the text

	Fig. 6 Growth rates in dependence on ωUH ∕ωce for Models 1S–4S. Top: plots for three different temperatures of hot electrons vt . Bottom: plots for two different temperatures of background plasma electrons vtb .
In the text

	Fig. 7 Shifts of frequency of growth rate maxima from frequency given by gyro-harmonic number s.
In the text

	Fig. 8 Frequency of growth rate maximum ${({\omega }_{\text{UH}}/{\omega }_{\text{ce}})}_{\text{Max}}$ as function of ratio of background and hot electron densities ne ∕nh for parameters in Model 2S and s = 6.
In the text

	Fig. 9 Comparison of the simulated and analytically derived growth rates for vt = 0.2 c and vtb = 0.03 c.
In the text

	Fig. 10 Maximal saturation energies of the upper-hybrid waves in dependence on the gyro-harmonic number s for Model 2S with the exponential fit.
In the text

	Fig. 11 Saturation energies as a function of ωUH ∕ωce . Top: plots for three different temperatures of hot electrons vt . Saturation maxima are shifted to lower ratios ωUH ∕ωce . Bottom: plots for two different temperatures of background plasma electrons vtb .
In the text

	Fig. 12 Shifts of the frequency of the saturated energy maxima from the frequency given by the gyro-harmonic number s.
In the text