# Dynamics of electron beams in the inhomogeneous solar corona plasma

## Abstract

Dynamics of an spatially limited electron beam in the inhomogeneous solar corona plasma is considered in the framework of weak turbulence theory when the temperature of the beam significantly exceeds that of surrounding plasma. The numerical solution of kinetic equations manifests that generally the beam accompanied by Langmuir waves propagates as a beam-plasma structure with a decreasing velocity. Unlike the uniform plasma case the structure propagates with the energy losses in the form of Langmuir waves. The results obtained are compared with the results of observations of type III bursts. It is shown that the deceleration of type III sources can be explained by the corona inhomogeneity. The frequency drift rates of the type III sources are found in a good agreement with the numerical results of beam dynamics.

## 1 Introduction

Type III bursts seem to be the most intensively studied type of solar sporadic radio emission. The theory of these bursts involves different aspects of solar physics, plasma physics, and nonlinear physics. Current understanding of the problem was summarized in the following reviews [Goldman (1983), Suzuki and Dulk (1985), Melrose (1990), Muschietti (1990)]. Thus, fast electrons propagating along open magnetic field lines become a source of Langmuir waves. In their turn, Langmuir waves are partly transformed into observable radio emission via nonlinear plasma processes. Propagation of electron beams accelerated during the solar flares is considered to be one of the key problems in the theory of type III solar radio bursts [Melrose (1990)]. The difficulty of the problem originates from the fact that electrons propagating in a collisionless plasma generate Langmuir waves, which influence on the electron propagation. For most type III events, the level of Langmuir waves generated by the beam is quite low [Grognard (1985)] and the beam dynamics is considered in the framework of weak turbulence theory [Vedenov, Velikhov and Sagdeev (1962), Drummond and Pines (1962)]. Even in the simplest, one-dimensional case the system of kinetic equations cannot be solved without simplifying assumptions [Ryutov and Sagdeev (1970)]. Therefore, there are many works treating this problem analytically [Ryutov and Sagdeev (1970), Zaitsev, Mityakov, and Rapoport (1972), Grognard (1975), Takakura and Shibahashi (1976), Mel’nik (1995)] as well as numerically [Magelssen and Smith (1977), Takakura and Shibahashi (1976), Grognard (1980), Takakura (1982), Grognard (1985), Mel’nik, Lapshin, and Kontar (1999)].

The presence of coronal plasma inhomogeneity is the important feature in the description of type III solar radio bursts. The electron beam propagating in the plasma with decreasing density (i.e. decreasing local plasma frequency) can emit at the plasma (fundamental) and double plasma (harmonic) frequencies. The relation between the drift rate measured at the local plasma frequency and the radial source velocity is given by

(1) |

where is the radial distance from the center of the Sun, is the local plasma frequency, is the plasma electron density. The presence of the local plasma frequency gradient leads to two physical effects in the kinetics of Langmuir waves:

i) The Langmuir wave propagating in the inhomogeneous plasma experiences a shift of wavenumber , due to the variation of the local refractive index.

ii) The time of beam-plasma interaction (quasilinear time) depends on the local density and therefore the resonance condition for the plasmons may itself change during the course of beam propagation.

The first effect is considered to have main impact on Langmuir wave kinetics and effect ii) can be neglected in comparison with effect i) [Coste et al. (1975)].

The influence of plasma inhomogeneity is usually considered to be a suppressing process for quasilinear relaxation [Ryutov (1969), Breizman and Ryutov (1969), Krasovskii (1978)]. It was shown that the inhomogeneity of plasma can significantly change the rate of quasilinear relaxation for low-density electron beams, when the length of quasilinear relaxation is comparable with the characteristic length of plasma inhomogeneity. Langmuir waves generated by the beam move out of the resonance velocities that decrease wave generation and prevent plateau formation at the electron distribution function.

If we assume that the electron beam generating type III bursts has low density then the beam can propagate in the state close to marginal stability giving rise to Langmuir waves due to plasma density fluctuations [Robinson (1992b), Robinson and Cairns (1993)]. On the contrary there are some evidences [Thejappa and MacDowall (1998)] that electron beams associated with the bursts are so dense that they can generate Langmuir waves in strong turbulence limit. Therefore, it is natural to take an intermediate beam density for consideration. Thus, the electron beam is dense enough to form a quasilinear steady state but the level of plasma turbulence is within the weak turbulence limit. For such electron beams the quasilinear length is much smaller than that of plasma inhomogeneity and the effects connected with the plasma inhomogeneity are usually not taken into account. However, at the starting frequencies of type III bursts MHz the drift rate is high (and consequently the plasma density gradient) and the exclusion of the plasma gradient cannot be justified. Therefore, to achieve the quantitative understanding in theory of type III bursts one has to include the influence of plasma density gradient.

For a theoretical descriptions of experiments in laboratory plasma [Ryutov (1969), Breizman and Ryutov (1969), Krasovskii (1978)] and some astrophysical applications [McClements (1989)] it is necessary to consider the stationary boundary-value problems involving the injection of a beam into a plasma half-space. However, for the interpretation and theoretical description of type III bursts the initial value problem seems to be more natural.

In this paper we investigate the dynamics of spatially limited fast electron beam in the inhomogeneous solar corona plasma with monotonically decreasing density. In our calculations we use quite realistic isothermal solar corona density model [Mann et al. (1999)]. The consideration of electron propagation is based on the system of quasilinear equation of weak turbulence theory where the influence of plasma density gradient is included. For the first time the dynamics of the electron cloud is investigated at large distances when the influence of plasma inhomogeneity is taken into account. We performed numerical computations of the beam dynamics for a wide range of initial beam parameters: initial velocity of the beam from to and initial beam density from cm to cm. These parameters allows us to obtain the source velocities in the wide range from to that covers the majority of observational data. The results obtained are compared with the results of observations.

## 2 Quasilinear equations for inhomogeneous plasma

Let us consider the propagation of the electron beam cloud when the energy density of excited Langmuir waves is much less than that of surrounding plasma

(2) |

where is the energy density of Langmuir waves, is the temperature of the surrounding plasma, is the wave number, and is the electron Debye length. Our analysis uses one-dimensional kinetic equations following Ryutov and Sagdeev (1970), Zaitsev, Mityakov, and Rapoport (1972), and Grognard (1985). Firstly, the one-dimensional character of beam propagation is supported by numerical solution of 3D equations Churaev and Agapov (1980). Secondly, in the application to III type solar bursts electrons propagate along the magnetic field and Suzuki and Dulk (1985) that ensures one-dimensional character of electron propagation.

In this case, as was shown in Vedenov, Velikhov and Sagdeev (1962); Drummond and Pines (1962) one can use equations of quasilinear theory

(3) |

(4) |

where is the electron distribution function, is the spectral energy density of Langmuir waves. plays the same role for waves as the electron distribution function does for particles. The system (3,4) describes the resonant interaction of electrons and Langmuir waves, i.e. an electron with the velocity can emit or absorb a Langmuir wave with the phase velocity . In the right-hand side of equations (3,4) we omitted spontaneous terms due to their smallness in comparison with induced effects Ryutov and Sagdeev (1970).

The evolution of the spectral function of the quasiparticles due to the spatial movement of Langmuir waves can be described with the aid of the Liouville equation Camac et al. (1962); Vedenov, Gordeev and Rudakov (1967)

(5) |

where the two last terms in right hand side of (5) describe propagation of Langmuir waves and the corresponding shift in wavenumber . Equation (5) describes the propagation of Langmuir waves in the approximation of geometrical optics or the WKB approximation Galeev and Rudakov (1963); Vedenov, Gordeev and Rudakov (1967). Clearly, this approximation imposes restrictions Ryutov (1969). The wavelength variation over the density fluctuations scale should be a small fraction of itself

(6) |

Using the fact that the frequency of a Langmuir wave does not change over its propagation that implies a wavenumber , which changes with position. Thus,

(7) |

or

(8) |

Using , , and eq. (8) we have from equation (6)

(9) |

where is the scale of the local inhomogeneity (1). If the condition (9) is fulfilled than propagation of Langmuir waves can be considered within the WKB approximation.

The clouds of fast electrons are formed in the spatially limited regions of solar corona where acceleration takes place. Therefore, spatially bounded beam is taken for consideration. The initial electron distribution function is

(10) |

where is the characteristic size of the electron cloud, is the velocity of the electron cloud, and is the width of the beam in the velocity space. It is also implied that initially the spectral energy density of Langmuir waves

(11) |

is of the thermal level and uniformly distributed in space. The system of kinetic equations (3-5) is nonlinear with three characteristic time scales. The first is the quasilinear time that is determined by the interaction of particles and waves. The second is the characteristic time of Langmuir wave velocity shift due to plasma gradient. The third scale is the time length of the electron cloud . And we are interested in dynamics of the electron cloud at the time scale .

To solve kinetic equations of quasilinear theory we use finite difference method Smith (1985); Samarskii (1989); Thomas (1995). Each differential operator is replaced by finite-difference one. To achieve high accuracy of transport terms in space and velocity space we use monotonic scheme van Leer (1974, 1977a, 1977b) that is considered to be the best in such kind of problems Hawley, Smarr and Wilson (1984). To approximate the right hand terms we use implicit finite difference scheme Takakura (1982); Samarskii (1989). This scheme seems to have better results near plateau maximum velocity Kontar (2001).

## 3 Solar corona density model

In order to obtain the density dependency of the solar corona we follow the heliospheric density model Mann et al. (1999). On one hand this is quite a simple isothermal model. However, it gives results which agree with observations (see Mann et al. (1999) for details). These two features justify the choice of the model.

To derive the required heliospheric density model the magnetohydrodynamic equations

(12) |

(13) |

(14) |

are employed with the velocity of the flow, the mass density , the thermal pressure , the magnetic field , and the gravitational force , where - gravitational constant; - is the mass of the Sun, is the unit vector along the radial direction. The system of equation (12-14) is completed with the isothermal equation of state

(15) |

- electron number density, - Boltzmann’s constant, is the temperature. The mass density and the electron number density are related by ( is the proton mass). In the solar corona and the solar wind has a value of Priest (1982).

The stationary spherical symmetric solutions of the system (12-15) can be found. Assuming all the variables to be radially directed equation (12) can be integrated to

(16) |

with the constant . And equations (13,14) can be transformed to the second equation

(17) |

with the critical velocity and the critical radius , the distance where corona becomes supersonic.

The density distribution is found by numerical integration of the equations (16,17). The constant appearing in (16) is fixed by satellite measurements near the Earth orbit (at AU cm ). Thus following Mann et al. (1999), constant is determined to be s. For calculations we take the electron temperature to be K. To solve equations (16,17) numerically rugula falsi method is used Korn and Korn (1961). The density and corresponding local plasma frequency profiles are presented in fig. 1.

In the limit , (16) and (17) can be solved analytically in terms of a barometric height formula:

(18) |

with and cm.

Note that the main factor ruling the density distribution at the distances of our interest is the gravitational force and the density distribution must be close to the barometric formula. Therefore, the choice of the different model and improvements of equations (12-15) should not lead to significant change in density profile.

## 4 Results of numerical solution of kinetic equations

The initial location of the source is taken to be km Magelssen and Smith (1977). In the paper we refer to the attitude as the initial location of the electron cloud . Electron beam is assumed to be close to monoenergetic, and the spatial size of the beam is taken to be cm for all calculations presented in the paper.

### 4.1 Initial evolution of the electron beam

The value determining the time of beam - plasma interaction and consequently the character of beam propagation is the quasilinear time. In the case of inhomogeneous plasma quasilinear time depends on distance whereas in the uniform plasma it is a constant value. When the quasilinear time is smaller than the time of electron beam propagation we have quasilinear propagation and free propagation in the opposite. If the initial electron beam density so that then electrons initially move without interaction with plasma

(19) |

As a result of electron propagation in plasma with decreasing density quasilinear time

(20) |

becomes smaller. Thus, at the maximum of electron density equation (20) becomes

(21) |

where , , the plasma density and the local plasma frequency are defined by barometric formula (18). Thus, quasilinear relaxation appears for the first time when given by (21) becomes equal to - the time of electron propagation.

Numerical solution of kinetic equations (3-5) gives us that various electron beam densities lead to a different points where relaxation occurs for the first time. In fig. 2 the

dependency on electron beam density of the distance where Langmuir waves start to grow for the first time and a plateau is formed in the electron distribution function is presented. Later, we see that the maximum of the Langmuir wave turbulence is located at this point. We see the strong tendency of the first relaxation point to decrease with density growth. Thus for beam densities cm relaxation and correspondingly the starting frequency of bursts takes place near the initial electron beam location . For beams of low density cm the distance between the initial location and relaxation point can be significant (about for cm) (fig.2).

### 4.2 The electron distribution function and the spectral energy density of Langmuir waves

The electron distribution function and the spectral energy density are presented in fig. 3,4. In every spatial point the electron distribution function is a plateau from almost zero velocity up to some maximum velocity. This maximum velocity is not a constant value as it takes place in uniform plasma but a decreasing function of distance. It is the direct result of plasma inhomogeneity.

We also see that the slow spatial movement of Langmuir waves leads to the shift of phase velocities toward smaller ones. In fig. 4 we see the result – there is no plasma waves with the phase velocity more than .

The spectral energy density consists of two parts. The first part is the Langmuir turbulence located near initial beam location (see fig. 3,4). This part of Langmuir waves drifts toward small phase velocities. Thus, the peak is placed close to at s, while at s the peak is close to (initially the peak appears at at ). In the homogeneous plasma the spectral energy density of Langmuir waves also consists of two parts but the waves located near injection point are time independent. Following Kontar, Lapshin and Melnik (1998) we have

(22) |

In the case of inhomogeneous plasma these waves (22) evolve in accordance with equation (5), where we can neglect the spatial movement of plasma waves due to the smallness of group velocity and the interaction with electrons due to the electron absence in this region

(23) |

where . Integrating (23) and using (22) as an initial condition we obtain the following solution

(24) |

where is the maximum phase velocity for these Langmuir waves. The equation (23) is mathematically correct and found in a good agreement with the numerical results. Indeed, since the group velocity of Langmuir waves is a small but a finite value , the spatial movement of waves is hardly observable. Thus, at s the peak is close to , and only after s the peak approaches to . However, the neglect of wave movement violates the limitations of geometrical optics (9), the velocity shift of Langmuir waves is only due to the spatial wave movement.

The second part of Langmuir turbulence accompanies the electron cloud. The electrons together with plasma waves present a beam-plasma structure Mel’nik (1995); Kontar, Lapshin and Mel’nik (1998). These plasma waves are also strongly influenced by the plasma density gradient. Thus electrons arriving at a given point excite Langmuir waves. The level of plasma waves is growing up to the moment when the maximum of electron density arrives to this point. After this moment electrons mainly with arrive to this point that leads to absorption of Langmuir waves (see Mel’nik (1995), Mel’nik, Lapshin and Kontar (1999), Kontar, Lapshin and Mel’nik (1998) for details). Negligible spatial movement of waves causes the plasma waves to move toward smaller phase velocities. As a result the part of Langmuir waves appeared to be out of resonance with particles. The back front has more low velocity plasmons that can be absorbed by the electrons, and the fastest electrons have no Langmuir waves with the corresponding phase velocity to absorb. Therefore, a part of beam-plasma structure energy is lost by the structure in the form of low velocity Langmuir waves (compare fig. 4). In the case of homogeneous plasma all the plasma waves generated at the front are absorbed at the back of the structure Mel’nik (1995); Kontar, Lapshin and Mel’nik (1998).

### 4.3 Energy density distribution in electron and plasma wave subsystems

The change of the spectrum of plasma waves influences the energy distribution in the beam-plasma system. In the homogeneous plasma beam-plasma structure propagates conserving its total energy but if the density gradient is different from zero the energy is pumped out of the structure. The balance between electron and Langmuir wave subsystems is not observed in the inhomogeneous plasma. The distribution of energy density is presented in fig. 5 at two different time moments.

The distribution of energy demonstrates that approximately a half of the initial beam energy is concentrated near the site of initial electron distribution as it takes place in the case of uniform plasma Mel’nik (1995); Kontar, Lapshin and Mel’nik (1998); Mel’nik, Lapshin, and Kontar (1999)

(25) |

where is the initial beam energy density. The high level of Langmuir waves concentrated near the initial beam location is due to initially nearly monoenergetic electron beam (10). If the initial distribution function (10) is proportional to Mel’nik, Lapshin, and Kontar (1999) then and the high level of Langmuir waves is not excited in the region of initial relaxation.

The second part of the wave energy is the result of beam-plasma structure propagation in the inhomogeneous plasma. The decrease in the number of waves resonant with the beam causes the energy losses. In fig. 5 we see nonzero level of Langmuir waves in the region where the structure has passed and the slow spatial movement of plasma waves can also be seen. Despite of the considerable energy losses the beam remains the source of a high turbulence level (fig. 5).

The dissipation of energy by the beam-plasma structure changes the spatial size of the beam-plasma plasma structure. The fast electrons, which appeared to be out of resonance with waves can overlap the structure electrons and pass some distance before they relax toward a steady state (i.e. a plateau in the electron distribution function). This effect increases the spatial size of the structure.

### 4.4 Velocity of a beam-plasma structure

Recently it was shown that the beam-plasma structure propagates in a uniform plasma with a constant velocity and the upper plateau border (maximum plateau velocity) remains equal to the initial value Mel’nik (1995); Kontar, Lapshin and Mel’nik (1998). The maximum velocity of the plateau in the inhomogeneous plasma is not a constant. As we see in the previous subsection the drift of the Langmuir waves leads to the decrease of the maximum wave phase velocity. The electrons with the velocity larger than the maximum wave velocity have no resonance waves with the same phase velocity and therefore propagate freely overtaking the interacting electrons. Thus the fastest electrons leave behind the electrons of the plateau and decrease the maximum velocity of the plateau i.e. decreasing the average velocity of the electrons at a given point. The dependency of the maximum plateau velocity measured at the electron density maximum (this velocity corresponds to the velocity of a beam-plasma structure) on distance is a decreasing function (fig. 6). The speed of the structure can be approximated as , where is the maximum velocity of the plateau near the maximum of electron density. Thus we obtain that the source of Langmuir waves propagates with a decreasing velocity.

## 5 Discussion and main results

From the physical point of view it is interesting to mention that the main physical process leading to the change of the electron distribution function and the spectrum of Langmuir waves is the shift of the wavenumber (or phase velocity) due to Langmuir wave spatial movement with group velocity. Despite the fact that the increment of beam-plasma instability depends on distance it can be neglected in comparison with the first effect. Indeed, as soon as we omit the terms responsible for the phase velocity shift in kinetic equations the numerical solution becomes very close to the homogeneous plasma case Kontar, Lapshin and Mel’nik (1998); Mel’nik, Lapshin, and Kontar (1999). This numerical result is consistent with the qualitative conclusions in Coste et al. (1975).

The results obtained in the previous section also bring us new understanding for the theory of the radio bursts. Unfortunately we do not directly observe the electron beams or Langmuir waves (the exception is the satellite observations, which are far from the initial beam acceleration site). However, following the standard model Ginzburg and Zheleznyakov (1958) that Langmuir waves are transformed into observable radio emission via nonlinear plasma processes one infers that the maximum of the spectral energy density of Langmuir waves corresponds to the maximum of radio emission.

The numerical results show that electron beams with different electron beam densities have initial relaxation point at various distances. Electron beams with low density initially propagate freely. However, the quasilinear time fastly decreases with distance due to the decreasing plasma density. Therefore, at some distnace, which is a function of beam density, quasilinear relaxation appears. The point where relaxation occurs for the first time give us the starting frequency of the type III bursts. Normal type III bursts usually start somewhere between 250 and 500 MHz Suzuki and Dulk (1985) and we can see from fig. 2 that specifically this range of frequencies seems to be the most probable for typical beam densities.

We also see that the electron beam propagates as a beam-plasma structure, which is the source of Langmuir waves and consequently radio emission. Observations give us that the sources of type III bursts propagate with almost a constant velocity, typically in the range Suzuki and Dulk (1985). Direct measurements of solar electrons by space probes imply that the velocity of the sources are slowing down Fainberg and Stone (1970). Radial variations of the exciter velocity are also required to satisfy consistency with other observational parameters Robinson (1992a). The numerical results demonstrate that the velocity of beam-plasma structure becomes less as the structure propagates into a plasma whereas in the uniform plasma the speed of the structure is a constant value Kontar, Lapshin and Mel’nik (1998); Mel’nik, Lapshin, and Kontar (1999).

In observations the frequency drift rate is normally measured as a function of frequency. Alvarez and Haddock (1973) give

(26) |

(where in Mhz, and in MHz/s) the least-square fit to the drift rates reported by various authors for the frequency range from MHz to kHz. Indeed, assuming that the velocity of the source is a constant and using (1) we cannot explain the observational drift rate (26). To obtain the observational properties of drift rate one should take into account the deaceleration of electron beams. The plasma inhomogeneity seems to be the physical process explaining the drift rate observed in the range from 412 MHz to 70 MHz (fig. 7). From fig. 7 it could be concluded that the ”typical” type III burst is generated by the electron beam with the initial velocity cm s (approximately ) and density cm.

Our numerical results demonstrate that the main factor defining the frequency drift rate is the initial electron density of the beam (fig. 7). The initial velocity of the beam does play a role but the impact on frequency dependency seems to be weaker (see fig. 7).

Note, that as the velocity of the structure decreases with structure propagation the decrease rate also declines. The highest velocity decrease rate is near the starting frequencies of type III bursts, where the plasma gradient has the biggest value (27). Thus, the gradient of local plasma frequency

(27) |

declines with the distance and respectively the influence of plasma inhomogeneity becomes smaller. Therefore, one can expect that the influence of plasma inhomogeneity on beam dynamics become smaller further away from the initial beam location.

The second part of Langmuir turbulence generated by the electron beam is located at the region of the initial electron beam relaxation. Beam-plasma structure propagating in the inhomogeneous plasma also leave behind quite high level of Langmuir turbulence. These Langmuir waves propagate slowly in a plasma with the group velocity cm s. These waves can easy transform into observable radio emission and might be responsible for the type V bursts. Indeed, if we note that Langmuir waves drift in velocity space due to inhomogeneity toward smaller phase velocities than we might expect the decrease of the transformation rate of Langmuir waves into observable radio emission and consequently the long duration of the emission. Moreover, the Langmuir waves drift very slowly and can be observed only in the range from the starting frequencies of type III bursts to the distance of a few solar radii, where plasma density gradient becomes much smaller. Finally we should note that these plasma waves are determined by the initial distribution function of the beam electrons. The distribution, which is proportional to the velocity Kontar, Lapshin and Mel’nik (1998) gives no waves in the initial beam location whereas monoenergetic beam leaves a half of its kinetic energy in this site Mel’nik (1995).

## 6 Summary

Dynamics of electron beam in the decreasing plasma density demonstrate rather strong influence of plasma inhomogeneity. The strongest influence is observed near starting frequencies of type III solar radio bursts (MHz). The beam of fast electrons propagates in the corona plasma as a beam-plasma structure Mel’nik (1995); Kontar, Lapshin and Mel’nik (1998) that explains on one hand the ability of beams to propagate over large distances being a source of plasma waves and on the other hand the decrease of velocity of the burst source. The plasma with declining density plays a role of dissipative media, where beam-plasma structure in the course of its propagation leaves a part of its energy in the form of Langmuir waves. The comparison of the obtained drift rate with observational data manifests a good agreement. The results obtained allows us to quantitatively connect the some observable properties of bursts (drift rate and starting frequency) with the basic plasma an beam parameters (plasma density, density gradient, beam density, and beam velocity).

However, more theoretical and numerical work is needed to enable detailed comparison of the observational properties and theoretical results. The brightness temperature, the size of the source, radiation flux can not be determined without involving the nonlinear plasma processes responsible for radio emission. Moreover, solar corona plasma is not a monotonic function of distance but rather inhomogeneous with different length scales. All these and probably other processes should be included to achieve quantitative understanding of electron beam dynamics.

### References

- Alvarez, H., and Haddock, F.T.: 1973, Solar Phys. 29, 197.
- Breizman, B.N. and Ryutov, D.D.: 1969, JETP 57, 1401.
- Camac, M., Kantrovitz, A.R., Litvak, M.M., Patrick, R.M., and Petschek, H.E.: 1962, Nucl. Fusion, Suppl. 2 423
- Churaev, R.S. and Agapov, A.V.: 1969, Sov. J. Plasma Physics 6, 232.
- Coste, J.,Reinisich, G., Silevitch, M.B., and Montes, C.: 1975, Phys. Fluids 18, 679.
- Drummond, W.E., and Pines, D.: 1962, Nucl. Fusion, Suppl. 3 1049
- Fainberg, J., and Stone, R.G.: 1970, Solar Phys. 15, 433.
- Galeev, A.A., and Rudakov, L.I.: 1963, JETP 45, 647.
- Ginzburg, V.L., and Zheleznyakov, V.V.: 1958, Sov. Astron. J. 2, 653.
- Goldman, M.V.: 1983, Solar Phys. 89, 403.
- Grognard, R.J.-M.: 1975, Australian J. Phys. 28, 731.
- Grognard, R.J.-M.: 1980, In. Radio Physics of the Sun ed. M.R. Kundu, T.E. Gergely, D. Reidel Publishing Company, 303.
- Grognard, R.J.-M.: 1985, Solar Radiophysics ed. McLean, N.J.,Labrum, N.R., Cambridge Univ. Press, 289.
- Hawley, J.F., Smarr, L.L., and Wilson, J.R.: 1984, Ap. J. Suppl. 55, 211.
- Kontar, E.P., Lapshin, V.I., Mel’nik, V.N.: 1998, Plasma Phys. Rep. 24, 772.
- Kontar, E.P.: 2001, Computer Physics Communications, in press.
- Korn, G.A., and Korn, T.M.:1961, Mathematical Handbook for scientists and engineers, New York:cMGRAW-HILL BOOK COMPANY.
- Krasovskii, V.L.: 1978, Plasma Phys. Rep. 4, 1267.
- Magelssen, G.R., Smith, D.F.: 1977, Solar Phys. 55, 211.
- Mann, G., Jansen, F., MacDowall, R.J., Kaiser, M.L., and Stone, R.G.: 1999 Astron. Astrophys. 348, 614.
- McClements, K.G.: 1989 Astron. Astrophys. 208, 279.
- Mel’nik, V.N.: 1995, Plasma Phys. Rep. 21, 94.
- Mel’nik, V.N. Lapshin, V.I., Kontar, E.P.: 1999, Solar Phys. 184, 353.
- Melrose, D.B.: 1990, Solar Phys. 130, 3.
- Muschietti, S.: 1990, Solar Phys. 130, 201.
- Priest, E.R.: 1982, Solar Magnetohydrodynamics Reidel: Dordrecht
- Robinson, P.A.: 1992a, Solar Phys. 137, 307.
- Robinson, P.A.: 1992b, Solar Phys. 139, 147.
- Robinson, P.A. and Cairns, I.H.: 1993, Astrophys. J. 418, 506.
- Ryutov, D.D.: 1969, JETP 57, 232.
- Ryutov, D.D. and Sagdeev, R.Z.: 1970, JETP 58, 739.
- Samarski, A.A.: 1989, Theory of difference schemes, Moscow: Science, 616.
- Smith, G.D.: 1985, Numerical Solution of Partial Differential Equations: Finite Difference Methods, Oxford: Clarendon Press.
- Suzuki, S., and Dulk, G.A.: 1985, Bursts of Type III and Type V, In Solar Radiophysics, ed. McLean N.J.,Labrum N.R., Cambridge University Press, 289.
- Takakura, T. and Shibahashi, H.: 1976, Solar Phys. 46, 323.
- Takakura, T.: 1982, Solar Phys. 78, 141.
- Thejappa, G.T. and MacDowall, R.J.: 1998, Astrophys. J. 498, 465.
- Thomas, J.W.: 1995, Numerical Partial Differential Equations: Finite Difference Methods, New York: Springer.
- van Leer, B.: 1974, J. Comput. Phys. 14, 361
- van Leer, B.: 1977a J. Comput. Phys. 23, 263
- van Leer, B.: 1977b J. Comput. Phys. 23, 276
- Vedenov,A.A., Velikhov, E.P., and Sagdeev, R.Z.: 1962, Nucl. Fusion, Suppl. 2, 465
- Vedenov, A.A., Gordeev, A.V., and Rudakov, L.I.: 1967, Plasma Phys. 9, 719
- Zaitsev, V.V., Mityakov, N.A. and Rapoport, V.O.: 1972, Solar Phys. 24, 444.