# Double plasma resonance instability as a source of solar zebra emission

###### Key Words.:

Instabilities – Methods: numerical – Sun: radio radiation^{†}

^{†}offprints: ,

###### Abstract

Context:The double plasma resonance (DPR) instability plays a basic role in the generation of solar radio zebras. In the plasma, consisting of the loss-cone type distribution of hot electrons and much denser and colder background plasma, this instability generates the upper-hybrid waves, which are then transformed into the electromagnetic waves and observed as radio zebras.

Aims:In the present paper we numerically study the double plasma resonance instability from the point of view of the zebra interpretation.

Methods:We use a 3-dimensional electromagnetic particle-in-cell (3-D PIC) relativistic model. We use this model in two versions: a) a spatially extended ”multi-mode” model and b) a spatially limited ”specific-mode” model. While the multi-mode model is used for detailed computations and verifications of the results obtained by the ”specific-mode” model, the specific-mode model is used for computations in a broad range of model parameters, which considerably save computational time. For an analysis of the computational results, we developed software tools in Python.

Results:First using the multi-mode model, we study details of the double plasma resonance instability. We show how the distribution function of hot electrons changes during this instability. Then we show that there is a very good agreement between results obtained by the multi-mode and specific-mode models, which is caused by a dominance of the wave with the maximal growth rate. Therefore, for computations in a broad range of model parameters, we use the specific-mode model. We compute the maximal growth rates of the double plasma resonance instability with a dependence on the ratio between the upper-hybrid and electron-cyclotron frequency. We vary temperatures of both the hot and background plasma components and study their effects on the resulting growth rates. The results are compared with the analytical ones. We find a very good agreement between numerical and analytical growth rates. We also compute saturation energies of the upper-hybrid waves in a very broad range of parameters. We find that the saturation energies of the upper-hybrid waves show maxima and minima at almost the same values of as the growth rates, but with a higher contrast between them than the growth rate maxima and minima. The contrast between saturation energy maxima and minima increases when the temperature of hot electrons increases. Furthermore, we find that the saturation energy of the upper-hybrid waves is proportional to the density of hot electrons. The maximum saturated energy can be up to one percent of the kinetic energy of hot electrons. Finally we find that the saturation energy maxima in the interval of = 3-18 decrease according to the exponential function. All these findings can be used in the interpretation of solar radio zebras.

Conclusions:

## 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 = , where and , , and 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

(1) |

where is wave vector, , , and are the hot electron velocities perpendicular and parallel to the magnetic field; is the electron mass, is the relativistic Lorentz factor, is the gyro-harmonic number, and 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 in the form

(2) |

where 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 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 . 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 , 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 Section 2 we describe our numerical model and initial conditions for the studied double plasma resonance instability. In Section 3 we present the results for different model parameters. Discussion of these results and conclusions are in Section 4.

## 2 Model

We use a 3-D 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.^{1}^{1}1https://www.terrapub.co.jp/e-library/cspp/text/10.txt. 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 broad range of model parameters. The model size
in x-, y-, and z-directions is for the multi-mode model, and for the specific-mode model, respectively, where
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 , electron plasma frequency , initial magnetic field is in the direction, and electron cyclotron frequency varies from to approximating the gyro-harmonic numbers . For dependencies of the growth rate on temperatures, we use in the range of . Higher values of 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 c, corresponding to the temperature in the interval 5.4-14.8 MK, and hot electrons with the DGH distribution having the velocity 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 per cell and the ratio of background to hot electrons was taken as 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 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 , where is the initial electric energy, 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 in Figure 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.

Comp. no. | |||
---|---|---|---|

1M | 0.15 c | c | 6 |

2M | 0.2 c | c | 6 |

3M | 0.2 c | c | 3 |

4M | 0.3 c | c | 5 |

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 , the gyro-harmonic number , the maximal growth rate, and Model 2M are presented in Figure 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 depends on the plasma density of hot electrons.

### 3.1 DPR instability 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 for Models 2M and 4M (Table 1; see Figure 3). Comparing changes of the distribution in these models (left and right columns in Figure 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 . In Figure 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 (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 Figure 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 (Figure 4). As seen in this figure, the interval of k-vectors is relatively narrow and remains practically the same during the 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.

### 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 important for 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 Figure 5. While in the multi-mode model we use the model size , in the specific-mode model in this case we use the model size . This specific-mode model covers interval of above 5.23 (see Figure 4). This means that this specific-mode model covers the -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 Figure 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.

### 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.

Model no. | ||
---|---|---|

1S | 0.15 c | c |

2S | 0.2 c | c |

3S | 0.3 c | c |

4S | 0.3 c | c |

Growth rates have the maxima that are shifted to frequencies lower than those given by the gyro-harmonic number (Figure 6, upper part and Figure 7) and this shift increases with increasing temperature. The ratio between maximal and minimal growth rates for each temperature is 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 = 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 -vector wave, in numerical computations -vectors in some interval are always in operation. The effect of the temperature of background plasma electrons is small (Figure 6, bottom part), which agrees with the analytical theory (Benáček et al., 2017). We also determined the shifts of maxima of growth rate (Figure 7). They increase with increasing and also depends on temperature of hot electrons.

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

As shown in Figure 8, the frequency of the growth rate maximum also depends on the ratio of the background and hot electron densities . This dependence is exponential as shown by its fit and converges to for high values of the ratio. Namely, low values of the 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 ().

In most of our simulations the density ratio 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).

### 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 on Figure 9. It was made as follows. We chose the upper-hybrid frequency and hot electron density as s ( GHz) and cm. We changed the growth rate plot expressed in dependence on (Benáček et al., 2017) to that in dependence on . Then, we transformed the simulated growth rates, computed for , to those with . The transformation was made in two steps. We shifted the simulated growth rate maxima toward lower values of according to the relation expressed in Figure 6. Finally, the growth rates were multiplied by the factor using the linear relation between the growth rate and density of hot electrons (Yasnov & Karlický, 2004). As seen in Figure 9, there is good agreement between the numerically and analytically computed growth rates.

### 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 . 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 .

We computed one set of models changing the gyro-harmonic number from up to and keeping the same temperatures as in Model 2S. As shown in Figure 10, the saturation energy level in local maxima decreases with increasing and converges to the value (see the exponential fit as presented in Figure 10). This exponential function can be used for an extrapolation of the saturation energies for even higher parameter , which is useful for zebras with many stripes.

Maximal values of the saturated energy levels of the upper-hybrid waves for low values of 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 (Figure 11) have the maxima and minima shifted to lower ratios of than those given by the integer values of the gyro-harmonic number . 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 (Figure 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 (Figure 11, bottom part).

## 4 Discussion and conclusions

Using the 3-D 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 . This frequency shift increases with an increase of the temperature of 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 . We used this dependence to extrapolate our numerical results, made mostly for the ratio = 8, to those with the ratio = 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 = 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 -vector space, not only with one -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 ().

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 . 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 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 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

- Benáček et al. (2017) Benáček, J., Karlický, M., & Yasnov, L. 2017, A&A, 555, A1
- Buneman & Storey (1985) Buneman, O. & Storey, L. R. O. 1985, Simulations of fusion plasmas by A 3-D, E-M particle code, Tech. rep.
- Dory et al. (1965) Dory, R. A., Guest, G. E., & Harris, E. G. 1965, Physical Review Letters, 14, 131
- Karlický & Bárta (2008) Karlický, M. & Bárta, M. 2008, Sol. Phys., 247, 335
- Karlický & Yasnov (2015) Karlický, M. & Yasnov, L. V. 2015, A&A, 581, A115
- Ledenev et al. (2001) Ledenev, V. G., Karlický, M., Yan, Y., & Fu, Q. 2001, Sol. Phys., 202, 71
- Matsumoto & Omura (1993) Matsumoto, H. & Omura, Y. 1993, Computer space plasma physics: simulation techniques and software, Terra Scientific Pub. Co, p.305
- Melrose & Dulk (1982) Melrose, D. B. & Dulk, G. A. 1982, ApJ, 259, 844
- Treumann et al. (2011) Treumann, R. A., Nakamura, R., & Baumjohann, W. 2011, Annales Geophysicae, 29, 1673
- Winglee & Dulk (1986) Winglee, R. M. & Dulk, G. A. 1986, ApJ, 307, 808
- Yasnov & Karlický (2004) Yasnov, L. V. & Karlický, M. 2004, Sol. Phys., 219, 289
- Yasnov & Karlický (2015) Yasnov, L. V. & Karlický, M. 2015, Sol. Phys., 290, 2001
- Zaitsev & Stepanov (1983) Zaitsev, V. V. & Stepanov, A. V. 1983, Sol. Phys., 88, 297
- Zheleznyakov & Zlotnik (1975) Zheleznyakov, V. V. & Zlotnik, E. Y. 1975, Sol. Phys., 44, 461
- Zlotnik (2013) Zlotnik, E. Y. 2013, Sol. Phys., 284, 579