Äîêóìåíò âçÿò èç êýøà ïîèñêîâîé ìàøèíû. Àäðåñ îðèãèíàëüíîãî äîêóìåíòà : http://star.arm.ac.uk/preprints/2010/571d.pdf Äàòà èçìåíåíèÿ: Wed Sep 26 13:34:55 2012 Äàòà èíäåêñèðîâàíèÿ: Tue Oct 2 03:06:02 2012 Êîäèðîâêà: ISO8859-5 Ïîèñêîâûå ñëîâà: supernova remnant

Fast Gyrosynchrotron Co des
Gregory D. Fleishman1,2 & Alexey A. Kuznetsov3,4 ABSTRACT Radiation produced by charged particles gyrating in a magnetic field is highly significant in the astrophysics context. Persistently increasing resolution of astrophysical observations calls for corresponding 3D modeling of the radiation. However, available exact equations are prohibitively slow to compute comprehensive table of high-resolution models required for many practical applications. To remedy this situation, we develop approximate gyrosynchrotron (GS) codes capable of fast calculating the GS emission (in non-quantum regime) from both isotropic and anisotropic electron distributions in non-relativistic, mildly relativistic, and ultrarelativistic energy domains applicable throughout a broad range of source parameters including dense or tenuous plasmas and weak or strong magnetic field. The computation time is reduced by several orders of magnitude compared with the exact GS algorithm. The new algorithm performance can gradually be adjusted to the user needs depending of what--precision or computation speed-- is to be optimized for the given modeling. The codes are made available for the users as a supplement to this paper. Subject headings: radiation mechanisms: non-thermal--methods: numerical-- Sun: radio radiation--Sun: flares--stars: flares--radio continuum: planetary systems

1.

Intro duction

Radiation produced by charged particles moving in magnetic fields (magnetobremsstrahlung) plays an exceptionally important role in astrophysics making dominant contribution to radio emission in most of the astrophysical ob jects, ma jor contribution to gamma-ray and
1 2 3 4

Center For Solar-Terrestrial Research, New Jersey Institute of Technology, Newark, NJ 07102 Ioffe Physico-Technical Institute, St. Petersburg 194021, Russia Armagh Observatory, Armagh BT61 9DG, Northern Ireland Institute of Solar-Terrestrial Physics, Irkutsk 664033, Russia


-2- X-ray emission (in compact ob jects, supernova remnants, and gamma ray burst sources), and important contribution to the IR/optical/UV (in active galactic nuclei or extragalactic jets), which implies necessity of its calculation in highly diverse conditions specific to the ob ject or phenomenon under study. A straightforward way of performing such computations is the use of exact formulae for the emissivities and absorption coefficients (Eidman 1958, 1959; Melrose 1968; Ramaty 1969) valid for arbitrary conditions (when the magnetic field changes in space only smoothly and no quantum effect is important). However, these exact formulae are cumbersome and extremely slow computationally, especially, when high harmonics of the gyrofrequency are involved (that is, when the emission frequency is much larger than the gyrofrequency). This is why a number of simplified approaches have been developed, whose applicability regions and/or accuracy are limited. For example, for ultrarelativistic emitting particles a so-called synchrotron approximation (Korchak & Terletsky 1952; Getmantsev 1952; Korchak 1957; Syrovatskii 1959; Razin 1960) is valid (see Ginzburg & Syrovatskii 1965, 1969, for a review), which is often sufficient for such astrophysical ob jects as supernova remnants or radio galaxies. For the gyrosynchrotron (GS) regime, when nonrelativistic and mildly relativistic particles are involved--a typical case, e.g., for solar and stellar flares--and the synchrotron approximation breaks down, a few more approximations have been suggested. The simplest one is an analytical Dulk & Marsh (1982) approximation derived for isotropic power-law energy electron distributions, for a limited range of harmonic numbers (20-100) and moderate viewing angles (30-80 ) relative to the magnetic field; no thermal plasma effect is included (Dulk 1985). In this range the accuracy of these analytical formulae is within dozens percent. Although this approximation can be used for rough estimates and parametric dependences, it is apparently insufficient for quantitative treatment or detailed modeling. A better (fast numerical) approximation was proposed by Petrosian (1981) to calculate GS radiation ignoring the plasma and then generalized by Klein (1987) to include the plasma effect. This (Petrosian-Klein, PK) approximation is valid for a broader range of the parameters and yields typically an accurate radiation intensity, although it is not so precise in polarization and it does not reproduce the GS harmonic structure at low frequencies, where it can often be important. But the most severe limitation of the PK approximation is that it is only valid for isotropic (and weakly anisotropic) angular distributions of the emitting electrons. However, there is currently ample evidence of anisotropic electron distributions in many cases, for example, in microwave bursts accompanying solar flares (Lee & Gary 2000; Melnikov et al. 2002; Fleishman et al. 2003; Altyntsev et al. 2008). The only way of GS calculation from the anisotropic distributions is currently the use of exact formulae (Fleishman


-3- & Melnikov 2003a), which fitting inversions (Altyntsev for the development of new, into account the anisotropy makes both the detailed GS modeling and, especially, forward et al. 2008; Fleishman et al. 2009) prohibitively slow. This calls computationally effective schemes of the GS calculations taking (Fleishman et al. 2009).

This paper describes new approximate GS codes we developed for both isotropic and anisotropic electron distributions. To derive these codes we critically re-evaluated (and improved where needed) the assumptions made to obtain the PK approximation. These new codes allow for a large flexibility in selection of the computational options; e.g., they can be optimized for accuracy or computation speed with the full range of the intermediate options between these two extremes. The code mode, which is precise for most practical applications, is two-three orders of magnitude faster than the exact code; the mode fully optimized for speed is even about one more order of magnitude faster. No competing approximate GS code is currently available for the anisotropic electron distributions; for the isotropic distributions our code provides better accuracy than the PK code with a comparable computation speed.

2. 2.1.

Metho dology

Exact GS equations

Exact equations for the GS emissivity and absorption coefficient for an eigen-mode in a magnetized plasma have the form (Eidman 1958, 1959; Melrose 1968; Ramaty 1969):
jf =

2 e2 n f 2 ç 2 c 1 + T


ç

s= -

T (cos - n Å) + L sin Js () + Js () n sin

2

1-Å

2

ç d3 p, (1a)

ç F (p) f (1 - n Å cos ) - 2 e2 ç 2 n (1 + T )


sfBe

= -

ç

s= -

T (cos - n Å) + L sin Js () + Js () n sin

2

1-Å

2

ç d3 p, (1b)

ç

1 F (p) n cos - Å F (p) sfBe + f (1 - n Å cos ) - p p Å


-4- where f is the emitting frequency, fBe = eB /(2 me c) is the electron gyrofrequency, e and me are the electron charge and mass, B is the magnetic field, c is the speed of light, n , T , and L are the refraction index, and transverse and longitudinal (relative to the wave vector) components of the polarization vector respectively (see Appendix A), is the angle between the wave vector and the magnetic field vector, p and = v /c are the electron momentum and normalized (by c) velocity, Å = cos , is the electron pitch-angle (i.e., the angle between the electron momentum and the magnetic field vector), = (1 - 2 )-1/2 is the Lorentz-factor, Js () and Js () are the Bessel function and its derivative over the argument , f n sin 1 - Å2 . = (2) fBe The electron distribution function F (p) is normalized as follows (we assume that it is azimuthally symmetric, which results in 2 factor)
p
2

1

F (p) d3 p = 2
p
1

p 2 dp
-1

F (p, Å) dÅ = ne ,

(3)

where ne is the number density of electrons having the momentum between p1 and p2 . As has been said, equations (1) are computationally expensive, so only a very limited number of studies utilizing the exact GS formulae for anisotropic electron distributions (see, e.g., Fleishman & Melnikov 2003a,b; Fleishman 2006; Altyntsev et al. 2008; Tzatzakis et al. 2008; Reznikova et al. 2009) have been performed. Most of the GS studies, therefore, were performed for the isotropic electron distributions, although a number of 3D models considering either isotropic (Preka-Papadema & Alissandrakis 1992; Bastian et al. 1998; Fleishman et al. 2009) or anisotropic (Tzatzakis et al. 2008) electron distributions is yet very limited. The main reason for this lack of analysis is a long time required to compute the exact GS emissivity and absorption coefficient whose multiple calculations are needed to build a realistic 3D model. For example, in a recently developed "GS Simulator" (Fleishman et al. 2009), which is capable of modeling GS emission from a predefined magnetic structure populated with a thermal plasma and accelerated electrons, the user currently can select between fast PK code valid for the isotropic distributions only and slow exact code correct for arbitrary anisotropic distribution. However, while using the PK code yields the full model output over a few minutes (of a standard PC CPU), the exact code requires a day or so to do a comparable simulation, which apparently restricts practical applications of this tool. Thus, the fast GS codes applicable for a broad range of isotropic and anisotropic distributions and plasma parameter combinations are called for.


-5- 2.2. Petrosian-Klein approximation for the isotropic case

The most successful algorithm available now for fast and reasonably precise computation of the GS radiation was proposed by Petrosian (1981) and then developed by Klein (1987). Below we consider the approximations adopted in the PK algorithm and then generalize them as needed to accommodate the anisotropic case and improve the overall accuracy of the output. A key step of the PK algorithm is replacement of the summation over the discrete harmonic number s by integration over corresponding continuous parameter


s= -

U (s)
-

U (s) ds,

(4)

which is formally valid if (i) |s| 1 and (ii) contributions of adjacent terms, s and s + 1, are comparable, i.e., U (s) changes smoothly for large s. Within the PK approximation, this replacement is adopted for all the harmonic numbers; then, the integral over s can easily be taken analytically using the -function. Although replacement (4) is formally incorrect for low harmonic numbers and so the PK approximation cannot, in particular, reproduce a discrete (harmonic) structure of GS radiation at low frequencies, the approximate spectra do reproduce the mean level of the exact spectra remarkably well. The mentioned replacement implies that in place of integer-order Bessel functions, realorder Bessel functions enter the integrals, whose numerical computation can even be more demanding than that of the integer-order Bessel functions. This calls for use of reliable fast approximations for the real-order Bessel functions. Klein (1987) points out that the approximate expressions for the Bessel functions and their derivatives proposed by Wild & Hill (1971) are applicable; so the use of the Wild & Hill (1971) expressions represents the second important assumption within the PK algorithm. At this stage the GS emissivity and absorption coefficient represent double integrals over dp and dÅ, which is still a computationally expensive task. The third key simplification employed by Petrosian (1981) and Klein (1987) is the analytical Laplace estimation of the angular integral over dÅ. The Laplace integration method replaces the integrand by a best-fit gaussian function and then analytically takes the integral of this gaussian profile. To find this gaussian function requires to determine the gaussian peak Å0 and width Å0 . In a general case the peak position Å0 obeys an algebraic transcendental equation, although both Petrosian (1981) and Klein (1987) use a truncated simple analytical solution of this equation--the fourth approximation--approximately valid for the isotropic electron distributions, but insufficient to account for the pitch-angle anisotropy.


-6- And finally, because the standard PK algorithm is not precise at the low frequencies anyway, Klein (1987) adopts a fifth approximation--neglects the longitudinal component of the eigen-mode polarization vectors, which is indeed correct at high frequencies. Petrosian (1981) does not need this approximation since he considers the vacuum case, where no longitudinal wave component is present. Apparently, this approximation is minor and one can easily get rid of it by considering the full polarization vectors of the high-frequency eigen-modes in the magnetized plasma (see Appendix A).

2.3.

Strategy of Petrosian-Klein algorithm generalization

To obtain a better approximation within the Petrosian (1981) and Klein (1987) guidelines we have to reevaluate all five approximations described above in ? 2.2. In what follows we get rid of the fifth approximation and always use the full polarization vectors of the eigen-modes including the longitudinal component L , see Eq. (1). At the first stage of the generalization we apply the replacement (4) for all harmonics to obtain a more precise analogy of the PK expressions valid for both isotropic and anisotropic electron distributions. In doing so, like in Klein (1987), we employ both integration over ds with the -function and Wild & Hill (1971) approximate expressions for the Bessel functions. However, we more precisely evaluate the angular integral by more accurate fitting the corresponding integrand with a gaussian. The output of this step yields "continuous" (i.e., without the discrete low-frequency harmonic structure) spectra, which are reasonably precise at the high frequencies, while give correct mean level at the low frequencies even when the harmonic structure is prominent in the exact expressions. This output is often even more precise for the anisotropic than for isotropic distributions and computationally comparably fast compared with the original Klein's algorithm. We refer to these algorithm and code as "continuous" ones. At the second stage of the generalization we improve the algorithm towards reproducing the harmonic structure at low frequencies. We tried a number of approaches to this problem and found that the best results are obtained when we use the exact (i.e., composed of summation over integer harmonics, in contrast to "continuous" contribution discussed above) C expressions at some low frequencies, f fcr , and apply the continuous approximation at higher frequencies; we refer to the corresponding output as "hybrid" algorithm or "hybrid" code. At the third stage we apply a number of optimizations to the continuous and hybrid algorithms. For example, for the discrete contribution we use the exact Bessel functions at low frequencies (where the use of approximate expression was found to result in noticeable


-7- errors for some parameter combinations), typically at f /fBe < 12, while use the Wild & Hill (1971) approximations at higher frequencies. Also, we introduce correcting multiplicative factors to ensure perfect matches at all matching frequencies, where we switch from exact to approximate Bessel functions or from discrete to continuous contribution. As a result, we obtain a gradually tunable algorithm of GS calculations, which is overall fast and precise, and whose precision can further be increased at the expense of computation speed, and vice versa. The fastest (the least precise) version of the code is still precise enough for most practical applications (except those specifically interested in the low frequency harmonic structure of the GS emission), providing the spectrum accuracy typically better than 5% and better than 10% for all studied parameter combinations within the code applicability range, see ? 3.3. 3. 3.1. Continuous GS co de Analytical derivation

Now we turn to the actual implementation of our generalization strategy starting from more precise evaluation of the angular integrals. Making replacement (4) and performing integration over ds (using the -function) we get nf 2 4 2 e2 jf c fBe (1 + T 2 )
1 p
2

p 2 dp ç
p
1

ç
-1

T (cos - n Å) + L sin Js () + Js () n sin

2

1-Å

2

F (p, Å) dÅ, (5a)

4 2 e2 me c - nfBe (1 + T 2 )
1

p

2

2 dp ç
p
1

ç
-1

T (cos - n Å) + L sin Js () + Js () n sin çp

2

1-Å

2

ç

Fp (p, Å) FÅ (p, Å) F (p, Å) dÅ, (5b) + (n cos - Å) F (p, Å) F (p, Å)


-8- where we omitted the wave-mode indices for brevity, s is determined from the resonance condition and so has the form s= f (1 - n Å cos ). fBe (6)

For electromagnetic waves n < 1, thus, s > 0 and < s. Therefore, one can use the approximate expressions (Wild & Hill 1971; Klein 1987) Js (sx) 1 Z s (x) , 2 s a(s, x) (7a) (7b)

b(s, x) Z s (x) Js (sx) , 2 s x where sx = and x exp( 1 - x2 ) , Z (x) = 1 + 1 - x2 a(s, x) = (1 - x2 )3/2 + b(s, x) = (1 - x )
2 3/2

(8a) (8b) (8c)

A s

1/6

,

A = 0.503297, , B = 1.193000.

B + s

1/6

1-

1 5s2/3

Taking into account (6), parameter x receives the form x= n sin 1 - Å2 = . s 1 - n Å cos (9)

Substitution of the above expressions into formulae (5) for the emissivity and absorption coefficient yields f 2 e2 jf c n(1 + T 2 ) sin2 2 e2 me c - 3 n f (1 + T 2 ) sin2 where Q=
p
2

1

p 2 dp
p p
1

F (p, Å)Z 2s Q dÅ,

(10a)

-1 1

2

dp
p
1

F (p, Å)Z 2s QR dÅ,

(10b)

-1

[T (cos - n Å) + L sin + ab(1 - n Å cos )]2 , a2 (1 - n Å cos ) R=p
FÅ (p, Å) Fp (p, Å) + (n cos - Å) . F (p, Å) F (p, Å)

(11) (12)


-9- As has been explained, following Petrosian (1981) and Klein (1987) we estimate the integral over dÅ in (10) using the fastest descent (Laplace) method, which results in the following equations 2 e2 f jf c n(1 + T 2 ) sin2 2 e2 me c - 3 n f (1 + T 2 ) sin2
p
2

p2 F (p, Å0 )Z
p
1

2s 0

0

Q0

-

h

2 d p, (Å0 ) 2 d p, (Å0 )

(13a)

p

2

F (p, Å0 )Z
p
1

2s 0

0

Q0 R

0

-

h

(13b)

where index 0 means that the corresponding factor must be taken at (has yet to be calculated) value Å0 = Å0 (p) at which the function h(p, Å), h = ln[F (p, Å)Z 2s Q] = ln F (p, Å) + 2s ln Z + ln Q reaches its peak value as a function of Å for a given p value. Petrosian (1981) and Klein (1987) justified applicability of the Laplace method by the fact that, for isotropic or weakly anisotropic distributions, the variation of the angular integrands in (10) is dominated by the factor Z 2s whose dependence on Å is very sharp (and very similar to a gaussian, due to the nature of that factor, see Eq. 8a). Thus, one may use a gaussian fit instead the exact expression, and so expressions (13) are explicitly obtained by analytical integration of the corresponding gaussian fits. However, validity of the above approximation for the anisotropic distributions in a general case is not that obvious and requires further analysis. Figure 1 displays a representative set of angular integrands in (10a) for several pitch-angle distributions: (a) Isotropic distribution. (b) Gaussian distribution with the maximum at the pitch angle of 90 F (p, Å) exp - with Å = 0.1. (c) Beam-like distribution with the maximum at zero pitch angle F (p, Å) exp - with Å = 0.1. (Å - 1)2 Å 2 (16) Å2 Å
2

(14)

(15)


- 10 - The plots in Figure 1 correspond to the following parameters: X-mode, background plasma density n0 = 5ç109 cm-3 , magnetic field B = 270 G, emission frequency f = 10 GHz, emission direction = 75 , electron energy E = 1 MeV; the integrands are normalized to unity by their maximum values. Apparently, the gaussian fits to the integrands are comparably good for both isotropic and anisotropic cases, which makes the Laplace method of angular integral evaluation widely applicable. Note that the used value of Å = 0.1 describes the case of rather strong anisotropy, which has huge effect on the GS radiation. For each momentum p the peak position Å0 is determined as the root of equation h (Å) = 0, where h (Å) is the derivative of function h(p, Å) (14) over Å: h (Å) =
FÅ (p, Å) f n cos - Å 1 - x2 - n cos ln Z + (ln Q) , + 2 F (p, Å) fBe 1 - Å2

(17)

(ln Q) = 2

ab[(a /a + b /b)(1 - n Å cos ) - n Å cos ] - T n - T (cos - n Å) + L sin + ab(1 - n Å cos ) -2 a (A - ) = , a a6 (B - ) b2 - 1 b = , + 4n cos 6 b b1 b2 b1 = (1 - x2 )3/2 + f = , 6sfBe B s
1/6

n cos a + , (18a) a 1 - n Å cos (18b) (18c) b = b1 b2 , (18d) (18e)

,

b2 =

n cos , = s

1 - x2 = 3x (n cos - Å) . 1 - Å2
2

1-

1 5s2/3

,

Transcendental equation (17) has no closed solution for Å0 in a general case. For isotropic case the derivative of the distribution function vanishes, and then discarding terms ln Z and ln Q from (17) yields Å0,PK n cos , which is approximately valid for the isotropic case as is shown in Figure 1a. Petrosian (1981) and Klein (1987) applied this approximate root to derive their approximate GS codes. In the general case, however, when the pitch-angle anisotropy is not necessarily weak, this approximation breaks down. Inspection of Figure 1 suggests that the accuracy of the integral evaluation can be significantly improved if one uses the exact solution Å0 of transcendental equation (17) in place of the approximate one Å0,PK . Indeed, the solution of transcendental equations can easily be found numerically, e.g., with the bisection method (see e.g. Press et al. 1997). The vertical dotted lines in Figure 1 show the position of the approximate root Å0,PK used in the PK algorithm. This approximate root


- 11 - noticeably deviates from the true solution even for the isotropic electron distribution; thus, in addition to accounting of the anisotropic distributions, using the exact numeric solution for Å0 also improves the accuracy of GS computation from the isotropic distributions. The second derivative of h(p, Å) over Å has the form
FÅ (p, Å) FÅÅ (p, Å) - h (Å) = F (p, Å) F (p, Å)

1 - x2 (n sin )2 (n cos - Å)2 f 1+ - - 2 fBe 1 - Å2 (1 - n Å cos )3 (1 - x2 ) 2Å(n cos - Å) n cos (n cos - Å) - + (ln Q) . (19) + 1 - Å2 1 - n Å cos
2

The analytical expression for the term (ln Q) is very cumbersome. As a rule, the contribution of this term is minor. However, we found that for some combinations of parameters (including, e.g., isotropic distribution), the account of (ln Q) noticeably increases the computation accuracy. In such cases, it is more convenient to calculate (ln Q) by explicit numerical differentiation of (ln Q) , Eq. (18a). Electron distributions G(E , Å) over energy E (in place of momentum p) are often used in GS calculations. One can easily adjust equations (10-13) to this form of the distribution function by making simple variable transformation to integration over energy E . The emissivity and absorption coefficient then take the form f 2 e2 jf c n(1 + T 2 ) sin2
E
2

1

dE
E
1

Z 2s QG(E , Å) dÅ
-1 E
2

f 2 e2 c n(1 + T 2 ) sin2

2 Z0 s0 Q0 G(E , Å0 ) - E
1

2 dE , (20a) h (Å0 )

2 e2 c - 3 n f (1 + T 2 ) sin2

E

2

1

dE
E
1

Z 2s QR
-1

(E )

dÅ
E
2

2 e2 c - 3 n f (1 + T 2 ) sin2 where R
(E )

Z
E
1

2s 0

0

Q0 R

(E ) 0

-

h

2 dE , (20b) (Å0 )

=

G(E , Å) 1 + 2 n cos - Å G(E , Å) - G(E , Å) + E cp cp Å

(21)


- 12 - and the distribution function G(E , Å) is normalized as follows
E
2

1

2
E
1

dE
-1

G(E , Å) dÅ = ne .

(22)

Then, the logarithmic derivatives of the function F (p, Å) in (17) and (19) must be replaced by the corresponding derivatives of the distribution function G(E , Å)
GÅ (E , Å) FÅ (p, Å) = , F (p, Å) G(E , Å) FÅÅ (p, Å) GÅÅ (E , Å) = F (p, Å) G(E , Å)

(23)

at a given energy E .

3.2.

Numeric Implementation

The described algorithm was implemented in two (basically identical) numeric source codes--in FORTRAN and C++. We adopted a factorized form of the distribution function G(E , Å) = u(E )g (Å), (24)

where the type of the factors can be selected from a predefined list of options, similar to those used typically in X-ray modeling and data analysis. For example, the distribution over energy u(E ) can be single or double power-law over energy, momentum, or Lorentz-factor, include thermal component matched to a nonthermal tail (so-called thermal-nonthermal distribution, TNT, see Holman & Benka (1992); Benka & Holman (1994)), or have a form of kappa distribution. The angular distribution can be a gaussian or modified (with fourth-order term in the exponent) gaussian function, or a loss-cone distribution; isotropic distribution (g (Å)=const) is also a possibility. For the full list of the built-in energy and pitch-angle distributions, see Appendix B. For a factorized form, Eq. (24), the logarithmic derivatives of G(E , Å) are equal to corresponding logarithmic derivatives of g (Å), so we made the replacements GÅ (E , Å) g (Å) = , G(E , Å) g (Å)
GÅÅ (E , Å) g (Å) = . G(E , Å) g (Å)

(25)

As has been explained, finding precise solution for Å0 is very important to ensure the highest accuracy of the Laplace estimate of the angular integral, which is in turn a key of providing the overall accuracy of the algorithm. This root is a solution of transcendental equation h (Å) = 0. We found that this equation has always a unique solution if we discard the term (ln Q) from Eq. (17), so a simple bisection method of the root finding


- 13 - can efficiently be used. This approximation yields typically a very good approximation to the exact Å0 value, so the term (ln Q) is a small correction; however, its inclusion in the equation improves the accuracy of the algorithm for some parameter combinations, especially for O-mode radiation. A disadvantage of having this term in the equation is that the function h (Å) can become discontinuous (for some parameter combinations) and so the simple bisection method can converge to the discontinuity point instead of the true root. We found two equally good approaches to correctly account for the (ln Q) contribution. In the first of them, the bisection method is first used to find Å00 when the term (ln Q) is discarded and a correction to it, Å is determined analytically by the perturbation theory: Å = -(ln Q(Å00 )) /h (Å00 ). Then, this correction is used to set up the restricted interval of Å around Å00 in which the true solution of the full equation h = 0 (i.e., with (ln Q) included) is being specified by the bisection method. In the second one, the approximate root Å00 is also determined first and then it is improved with the secant method (Press et al. 1997). Both methods proved to give the same results with comparable computational expenses (the second one is marginally longer). The last step of the numeric implementation is integration over energy/momentum with the adopted distribution function. Since we use a number of different distribution functions, we also use a few different methods of the numeric integration. For the thermal distribution and the thermal-like (low-energy) parts of TNT and kappa distributions, we use gaussian integration or trapezoidal method with equally spaced nodes. For the power-law distributions and the power-law-like (high-energy) tails of TNT and kappa distributions, the trapezoidal or extended midpoint methods with logarithmic node placement (which gives an exact solution for true power-law integrand) are applied (Bastian 2006). We also used the trapezoidal integration method with step-by-step accuracy improvement (Romberg integration, Press et al. 1997). For all the methods, we tried to find an optimal balance between the accuracy and the computation speed, which was achieved by appropriate choice of either the number of nodes in the trapezoidal or extended midpoint methods or the desired accuracy in the Romberg method. In all cases the codes include also the free-free contribution to the radiation, which is extremely fast to compute. Figures 2-4 demonstrate a few representative examples of the GS spectra computed with (1) exact equations, (2) discrete approximation, i.e., exact equations in which the approximate expressions of the Bessel functions (Wild & Hill 1971) are used in place of exact ones, and (3) our new continuous approximation. The spectra were obtained under the assumption of a homogeneous source located at the Sun, i.e. the emission intensity (observed at the Earth) of the magnetoionic mode was calculated using the formula
S jf I = 2 (1 - e- R f L

),

(26)


- 14 - where S and L are the visible source area and its depth along the line-of-sight, respectively, X O and R 1.49ç1013 cm is the astronomical unit. Total emission intensity equals If = If +If , and the degree of polarization is defined as =
X O If - If . O X If + I f

(27)

The plots in Figures 2-4 correspond to the following parameters: background plasma density n0 = 2 ç 109 cm-3 and temperature T0 = 4 MK, magnetic field B = 370 G, accelerated electrons have a power-law energy spectrum F E - with = 3.5 between Emin = 0.01 MeV and Emax = 10 MeV, the number density of fast electrons nb = 9.3 ç 107 cm-3 , visible source area S = 1.8 ç 1018 cm2 , and source depth L = 6 ç 108 cm. In Figure 2, the electron pitchangle distribution is isotropic and the emission direction = 30 . In Figure 3, the electrons have loss-cone distribution (15) with Å = 0.2, and = 60 . In Figure 4, the electrons have beam-like distribution (16) with Å = 0.2, and = 60 . The integrals and series in the exact and approximate discrete expressions were calculated with relative accuracy 10-5 . Very high overall accuracy of radiation intensity, degree of polarization, and spectral index of the new approximation is evident, especially--at high frequencies, where the slope of the true spectrum is especially well reproduced in all considered cases. In some cases, however, a small systematic offset of the continuous approximation relative to the exact one is noticeable, which is partly related to the use of approximate Bessel function expressions (note the corresponding offsets between the exact and approximate discrete curves). Thus, using better approximation for the Bessel functions can further improve the accuracy of our continuous approximation. At low frequencies the continuous approximation is intrinsically unable to reproduce the harmonic structure of the GS emission, although the mean level of the spectra is reproduced remarkably well. Due to the same reasons, the continuous approximation cannot account for the cyclotron maser instability (negative absorption coefficient which results in coherent wave amplification). This instability can occur if the emission frequency is close to low harmonics of the cyclotron frequency, and typically requires an anisotropic electron distribution with R > 0, where R is specified by Eq. (12). One can see an example of the maser instability at Fig. 3c (loss-cone electron distribution, O-mode): the graphs for the exact equations and the discrete approximation have a very sharp peak at f 2fBe (2 GHz) while the continuous approximation misses this feature. We emphasize that the continuous approximation works equally well for both the isotropic distribution and for anisotropic distributions having strong effect on GS emission, which we refer to as "strong anisotropy" in this context. To show this explicitly, we added the degree of polarization and spectral index curves corresponding to the case of isotropic electron distribution (all other simulation parameters being the same as for the anisotropic distributions)


- 15 - in Figures 3g-h and 4g-h. One can see that the anisotropy affects both the spectral index and the degree of polarization quite significantly; in particular, for the beam-like distribution, the polarization even changes its sign compared to the isotropic case. This explicitly confirms quantitative applicability of our continuous codes to strongly anisotropic electron distributions. We performed all possible measures to optimize our codes, and the programs written in FORTRAN and C++ are checked to be comparably fast for comparable numeric settings. The continuous codes are much faster than the exact algorithm implementations. Figure 5 displays how the computation speed changes with the emission frequency (that, in turn, determines the maximal harmonic number used in the computations). For both exact and approximate discrete implementations the computation time increases nearly exponentially with the harmonic number, while for the continuous approximation this time does not depend on the harmonic number at all. Thus, the achieved improvement in the computation speed is remarkably great.

3.3.

Applicability region of the continuous co de

Perhaps no approximation can entirely substitute its exact prototype for all possible parameter combinations, thus, we specifically looked for those cases when the accuracy of the continuous codes decreases noticeably. As has been shown (see Figure 1), the basic reason of the overall success of the approximation is the possibility to precisely fit the angular integrand by a gaussian, whose parameters are determined based on analytical properties of the integrand. Therefore, we can expect, that the code performance will drop when the integrand cannot be well fitted by a gaussian profile, for example, when the integrand is essentially asymmetric. We found that for all analytical functions (continuous functions with all continuous derivatives) the approximation works well regardless the function sharpness. However, if the angular distribution function or its first derivative have a discontinuity at a certain c , then the angular integrand cannot be precisely fitted by a gaussian for some range of the emission angle around the angle c and the accuracy of the GS spectrum computation drops noticeably here: the approximate intensity can differ from the exact one by up to factor of 2. Such behavior of the approximate algorithm is in fact expected for the discontinuous angular electron distributions. We point out that even for an angular distribution function with a discontinuity of the second derivative at some c the accuracy of the continuous code


- 16 - can decrease around c . As an example, we consider a gaussian loss-cone distribution g (Å) = const for Å cos c , and g (Å) exp[-(Å - cos c )2 /Å2 ] for Å > cos c , whose second derivative has a discontinuity at c , while the first derivative and the function itself are both continuous. Figures 6-7 display an example when the angular integrand becomes noticeably asymmetric around the loss-cone boundary c and the gaussian function fails to provide a good fit to the integrand. This happens, however, only when the distribution function decreases very fast at Å > cos c ; in the given example we use c = 80 and Å = 0.05 (viewing angle = 80 , other simulation parameters are the same as at Fig. 3, and the integrands at Fig. 6 are plotted for the electron energy E = 1 MeV). One can see that the angular integrand deviates from a gaussian both for X- and O-modes, however, the integrand asymmetry is more significant for the O-mode. As a result, the continuous code fails to precisely reproduce the emission intensity, and the error is larger for the O-mode. We found that for smoother angular gradients, Å 0.1, which are still sharp enough for most practical applications, the integrand asymmetry remains small and the algorithm works well. We conclude that although the algorithm can fail for some models of the angular distribution, it has to work well for any natural distribution function since it (being a solution of a transport equation) is supposed to be an analytical function continuous with all the derivatives.

4.

Hybrid GS co des

Although the continuous approximation is sufficient for many practical applications, the low-frequency harmonic structure can also often be of interest, especially for the anisotropic case (Fleishman & Melnikov 2003b). Therefore, in this section we describe a number of improvements of the code towards recovery of the harmonic structure and overall higher accuracy of the codes.

4.1.

Harmonic structure recovery

We tried a few approaches to recover the GS harmonic structure at low frequencies and found that the simplest and the most straightforward way ensures the best results. In fact, C to recover the harmonic structure at frequencies f < fcr it is sufficient to use exact equations C at f < fcr , while the continuous approximation at higher frequencies. Since at low harmonics the computation speeds of both exact and approximate discrete codes are almost identical


- 17 - (Figure 5), while the error related to the use of approximate Bessel function expression is the biggest at these frequencies, we use the exact code (with exact Bessel functions) at the low frequencies. Our numeric experiments show that for full recovery of noticeable low-frequency harC C monics it is typically sufficient to use the exact code at f < fcr with fcr 12fBe and the continuous codes above. From Figures 2-4 we know that at high frequencies the continuous contribution displays the same slope as the exact one, although the approximate spectrum can systematically be shifted by a small value relative to the exact spectrum. Thus, simple transition from the exact contribution to the continuous one results in a matching jump or matching residual. To remove this jump we renormalize the continuous contribution as follows. First, we find the ratio of the values calculated by the exact and approximate forC C mulae at f = fcr , and then we apply this factor to all values calculated for f > fcr with the approximate code, i.e. j (f ) = j
f C f C jf E (fcr ) (f ) C C , jf (fcr )

(f ) =

C

(f )

C E (fcr ) , C C (fcr )

C for f fcr ,

(28)

where the upper indices E and C denote the values obtained using the exact and continuous codes, respectively. This renormalization provides both smooth transition from the exact to continuous contribution and also improves the accuracy of the spectra calculated at the high frequencies. Since this code involves both exact and continuous contributions, we call it the "hybrid" code. Having the matching frequency to be an adjustable parameter of the code allows the C C full range of options from purely continuous code (if fcr = 0) to the exact code (if fcr = or very large).

4.2.

Optimization of the hybrid co de

C For vast ma jority of applications the hybrid code with fcr 15fBe or even the purely continuous code seems to be sufficient. Nevertheless, some problems or settings requiring even more precise computations can exist in some cases (for example, various testing probC lems). For such cases we further optimize the code as follows. If fcr 10fBe is needed, the use of the exact Bessel functions becomes more and more computationally demanding. We found that we can significantly increase the computation speed if we use Wild & Hill (1971) expression for the Bessel functions, and so the discrete approximation in place of the exact W W W equations, at f > fcr H with fcr H 10fBe ; we use fcr H = 12fBe as a default value.

Apparently, the transition from exact algorithm to the discrete approximation results


- 18 -
W in a matching residual at f = fcr H , which is removed by the same renormalization as one C described in ? 4.1. As the frequency increases and reaches the next matching frequency fcr (if C W indeed fcr > fcr H ), the discrete approximation gives a way to the continuous approximation. This additional matching residual is removed by the same renormalization as before. ApC W parently, if fcr fcr H then the discrete approximation is not applied, and so no additional renormalization is needed.

4.3.

Optimization of the continuous co de

To round up the list of options and optimizations we discuss a number of different versions of the continuous code. First of all, if we neglect (ln Q) and (ln Q) terms in Eqns. (17) and (19) respectively, we obtain the code fully optimized for the computation speed. This code version produces the full GS spectrum at 100 frequencies and computes the degree of polarization in 15 - 30 ms depending on input distribution functions and the parameter combination. Typically, this version of the code produces spectra within the error better than 10%. The effect of neglecting (ln Q) and (ln Q) is illustrated in Figure 8, where "optimized" and "fast" refer to the codes with and without account of the mentioned terms, respectively (the electrons are assumed to have loss-cone distribution with Å = 0.3 and other parameters are the same as at Figure 3). Complementary, the account of the terms (ln Q) and (ln Q) reduces the error by a factor of two to the values better than 5% at the expense of increasing the computation time by around 70%. In many cases this code yields highly precise spectra (the error smaller than 2 - 3%) at those frequencies where no harmonic structure is pronounced. Nevertheless, we implemented one more optimization, which further improves the accuracy, when the code with (ln Q) and (ln Q) is not perfectly precise (the error is about 3 - 5%). The idea behind this improvement is very simple: the entire spectrum is still calculated with the approximate continuous algorithm, while the exact calculations are performed at a single frequency (we found that the case with f = 12fBe works well). Then, the continuous spectrum is renormalized as a whole to match the exact spectrum at this single frequency. This renormalization can reduce the errors of the continuous approximation by a factor of two or even better.


- 19 - 4.4. Nomenclature of the co des

Since we have developed and described many versions of the GS code, it is reasonable to summarize and compare them all, which is done in Table 1. The codes are ordered according to increase of the computation time. The fastest code is the continuous code (C) which typically provides appropriate accuracy but does not reproduce the low-frequency harmonic structure. The optimized continuous codes (those with either account of (ln Q), Q-optimized, or/and renormalization, R-optimized), which are basically the same continuous code but with better normalization, are only slightly slower and overall more precise than the continuous code (in particular, the polarization accuracy is improved greatly). No harmonic structure is reproduced either. Hybrid codes (H) composed of exact and/or approximate discrete contributions at low frequencies and of continuous contribution at high frequencies recover the harmonic structure. The first of them, the hybrid code, uses one matching frequency for transition from the exact to continuous code. This version of the code is very practical as it is both fast (only about five times slower than the continuous code) and highly precise if one selects C C fcr 12fBe . However, this code is much slower for fcr 10fBe , which makes it impractical for such settings. Optimized hybrid code, which uses the exact algorithm at f < 12fBe , while W C W the discrete approximation between fcr H = 12fBe and fcr > fcr H allows one to increase the C frequency fcr with less computational expenses still providing very high accuracy. However, C an increase of the computation time with fcr is noticeable in this code as well (see Figure 5). The approximate discrete code (D) uses the exact GS equations with the only one approximation--the approximate Bessel function expressions are used. This approximation basically yields very accurate results (with some exceptions at low frequencies and sometimes--a small systematic off-set of the optically thin part of the O-mode spectrum). However, since the summation over the harmonics is performed for arbitrarily high frequency in this version of the GS code, the computation time can be orders of magnitude longer than for the hybrid and continuous codes. Finally, the exact code (E) uses no approximation and so capable of producing results with arbitrarily high precision at the expense of the computation time, which can become much longer than for the approximate discrete code because the time required to compute the exact Bessel function Js increases with its order s. Overall, this set of the code implementations offers a complete list of options between fast and accurate continuous code through slower and more precise hybrid codes to the discrete and exact codes with approximate or exact Bessel functions respectively, so each potential user of the codes can make his own selection of the code option depending of the problem studied.


- 20 - 5. Co de Accessibility and Application

To ensure the widest practical application of these new codes including the perspective users who work more with IDL than with the computing languages like FORTRAN and C++, the above described codes were implemented as the Windows dynamic link libraries (DLLs) and Linux shared ob jects (SOs) callable from IDL. A brief description of these libraries is given in Appendix B. The libraries themselves with the complete description (including the calling sequence, input and output parameters, keys to switch between the different code versions, and the built-in distribution functions), examples of the calling IDL programs, and sample outputs are given in the online supplements to the paper. Although our GS codes can be applied for various astrophysical ob jects (e.g., active stars, planetary magnetospheres, cosmic jets, etc.) and even for laboratory settings, we have been developing them having in mind the solar flares as a primary target. In particular, this is why the calculated emission intensity is normalized to the distance from the Sun to the Earth. The online supplement contains a few sample spectra employing parameters typical for solar flares. Here we illustrate how the code application can be used to address key physics of solar flares. We note that using oversimplified approximate GS expressions has often led to a conclusion that the numbers of accelerated electrons deduced from hard X-ray (HXR) observations were insufficient to produce the observed microwave emission with the difference up to two-three orders of magnitude. Although the simulations using exact GS expressions (e.g., Altyntsev et al. 2008) do not confirm that conclusion, the use of the exact codes is highly time-consuming and so impractical. Let us consider specifically an occulted solar flare of 31 December 2007 recently reported by Krucker et al. (2010). This flare is particularly interesting because Krucker et al. (2010) argued that the observed coronal HXR source represented the very acceleration site of the flare; probing the acceleration regions is of tremendous importance for the modern astrophysics. Because of strong occultation, only the very top of the flaring loop was observed in the radio range, which implies that the observed source is reasonably uniform, so we can apply the uniform source model without explicit considering the radiation transfer. We took the source parameters derived from HXR by Krucker et al. (2010) (and selected a few of underconstrained parameters to have some typical values) and sent them to our GS code calculating the radio emission. As is apparent from Fig. 9, the model radio emission is in excellent agreement with the observed one. This means that the same electron population as needed to produce the HXR with the same spectral index extended to radio-producing energies (a few MeV) is fully consistent with the source radio properties.


- 21 - However, the radio data analysis is not only a cross-check for the HXR diagnostics: in addition to it the radio diagnostics put further constraints to such source parameters as the magnetic field, the viewing angle, the highest energy of accelerated electrons, and on the angular distribution of the fast electrons (e.g., Fleishman et al. 2009). Note that to obtain a reasonable "by eye" fit of the data by the model spectrum presented in this Figure, we selected the mentioned parameters (see the full parameter list on-line in the IDL program Flare071231a.pro) from reasonable ranges and adopted a moderate loss-cone anisotropy similar to that observed for a number of radio loop-top sources (Melnikov et al. 2002) of extended flaring loops; an efficient account of the anisotropy, as has already been said, would not have been possible without the codes developed. Small excess of the observed flux compared with the model spectrum at 1 GHz (if not related to an observation error or a radio model uncertainty) could easily be provided by a minor source nonuniformity which is known to broaden the GS spectrum.

6.

Discussion and Conclusions

In this paper we present a set of GS codes, which allows a smooth transition from exact (but often computationally demanding) GS code to approximate (but much faster and still highly precise) codes, where the computation time can be substantially reduced at the expense of very modest (if any) reduction of the computation accuracy applicable to both isotropic and anisotropic angular distributions of fast electrons. Specifically, the computation time is reduced by orders of magnitude, while the computation error remains within 1 - 10% depending on the fast code option and parameter combination. The importance of this development is difficult to overestimate: the reduction of computation time makes it widely applicable to create and analyze sophisticated models of GS radiation produced in realistic 3D configurations. It was impossible until now because that long time required to compute a single GS spectrum from an anisotropic electron distribution made the solution of the whole problem prohibitively slow. As a result, very few realistic 3D GS models have yet been developed. Having the fast GS codes readily available enables analyzing many more interesting geometries and parameter combinations and so better understanding how the GS spectra and images relate to the underlying source parameter and their evolution. The developed software implementation of the fast GS codes is performed in the way complementary to the standard HXR tools. In particular, the same list of widely used model distribution functions is built in our codes, which allows a user to choose the same electron distribution function for both HXR and radio modeling from the set of predefined model


- 22 - functions and then vary the distribution parameters of interest. A combined analysis of the hard X-ray and radio observations (as has been shown in the example considered in the previous section) is a straightforward way of constraining the source parameters and testing the models of solar flares and other phenomena involving energetic particles. Then, the fast codes are suitable for inclusion into the forward fitting inversion codes as the source functions (Fleishman et al. 2009). This may have a science transforming effect, because the possibility to use an arbitrary pitch-angle anisotropy in the trial fitting function paves the way to diagnose both energy and angular distributions of the fast electrons accelerated in solar flares along with the magnetic fields (Fleishman et al. 2009) and so develop the remote sensing method for solar flares and other phenomena revealing themselves via GS emission. This work was supported in part by NSF grants ATM-0707319, AST-0908344, and AGS0961867 and NASA grant NNX10AF27G to New Jersey Institute of Technology, and by the Russian Foundation for Basic Research, grants Nos. 08-02-92204, 08-02-92228, 09-02-00226, 09-02-00624. We thank Dr. Tim Bastian, NRAO, for stimulating discussions on the sub ject and also for valuable numeric input. We are highly grateful to the Nobeyama Observatory, Japan, for invitation and support of our short-term visits to the Observatory, where key steps of the reported work were performed. We thank Professor Kiyoto Shibasaki for his persistent encouragement of this work and Dr. Stephen White for helpful discussion of Nobeyama radio observations. We have made use of NASA's Astrophysics Data System Abstract Service.

REFERENCES Altyntsev, A. T., Fleishman, G. D., Huang, G.-L., & Melnikov, V. F. 2008, ApJ, 677, 1367 Bastian, T. S. 2006, in Astronomical Society of the Pacific Conference Series, Vol. 358, Astronomical Society of the Pacific Conference Series, ed. R. Casini & B. W. Lites, 173-176 Bastian, T. S., Benz, A. O., & Gary, D. E. 1998, ARA&A, 36, 131 Benka, S. G. & Holman, G. D. 1994, ApJ, 435, 469 Dulk, G. A. 1985, ARA&A, 23, 169 Dulk, G. A. & Marsh, K. A. 1982, ApJ, 259, 350 Eidman, V. Y. 1958, Soviet Phys. JETP, 7, 91


- 23 - --. 1959, Soviet Phys. JETP, 9, 947 Fleishman, G. D. 2006, in Solar Physics with the Nobeyama Radioheliograph, 51-62 Fleishman, G. D., Gary, D. E., & Nita, G. M. 2003, ApJ, 593, 571 Fleishman, G. D. & Melnikov, V. F. 2003a, ApJ, 587, 823 --. 2003b, ApJ, 584, 1071 Fleishman, G. D., Nita, G. M., & Gary, D. E. 2009, ApJ, 698, L183 Getmantsev, G. G. 1952, Dokl. Akad. Nauk. SSSR, 83, 557 Ginzburg, V. L. & Syrovatskii, S. I. 1965, ARA&A, 3, 297 --. 1969, ARA&A, 7, 375 Holman, G. D. & Benka, S. G. 1992, ApJ, 400, L79 Klein, K.-L. 1987, A&A, 183, 341 Korchak, A. A. 1957, Soviet Astronomy, 1, 360 Korchak, A. A. & Terletsky, Y. P. 1952, Zh. Eks. Teor. Fiz., 22, 507 Krucker, S., Hudson, H. S., Glesener, L., White, S. M., Masuda, S., Wuelser, J.-P., & Lin, R. P. 2010, ApJ, 714, 1108 Lee, J. & Gary, D. E. 2000, ApJ, 543, 457 Melnikov, V. F., Shibasaki, K., & Reznikova, V. E. 2002, ApJ, 580, L185 Melrose, D. B. 1968, Ap&SS, 2, 171 Petrosian, V. 1981, ApJ, 251, 727 Preka-Papadema, P. & Alissandrakis, C. E. 1992, A&A, 257, 307 Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1997, Numerical Recipes in C. The Art of Scientific Computing. Second Edition. (Cambridge University Press) Ramaty, R. 1969, ApJ, 158, 753 Razin, V. A. 1960, Izv. VUZov Radiofizika, 3, 584


- 24 - Reznikova, V. E., Melnikov, V. F., Shibasaki, K., Gorbikov, S. P., Pyatakov, N. P., Myagkova, I. N., & Ji, H. 2009, ApJ, 697, 735 Syrovatskii, S. I. 1959, AZh, 36, 17 Tzatzakis, V., Nindos, A., & Alissandrakis, C. E. 2008, Sol. Phys., 253, 79 Wild, J. P. & Hill, E. R. 1971, Australian Journal of Physics, 24, 43

A.

Disp ersion of the magnetoionic mo des

The refraction index of the electromagnetic waves in a plasma satisfies the dispersion equation 2v (1 - v ) , (A1) n2 = 1 - 2(1 - v ) - u sin2 + D where D = u2 sin4 + 4u(1 - v )2 cos2 , (A2) u= fBe f
2

,

v=

fp f

e

2

,

(A3)

fpe = e N/( me ) is the electron plasma frequency, and N is the total number density of the plasma electrons. For X-mode, = -1; for O-mode, = +1. The polarization vectors of the waves in the reference frame with z -axes along the magnetic field and the wave-vector k in the (xz )-plane, has the form e = (T cos + L sin , i, -T sin + L cos ) 2 1 + T + L2 (A4)

with the parameters T and L defined as 2 u(1 - v ) cos , T = u sin2 - D v u sin + T uv sin cos L = . 1 - u - v + uv cos2

(A5) (A6)

Electromagnetic waves can propagate in plasma if their frequency exceeds the cutoff frequency, f > fc , where fcO = fpe , fcX = fBe + 2
2 fpe + 2 fBe . 4

(A7)

A This preprint was prepared with the AAS L TEX macros v5.0.


- 25 -

Fig. 1.-- The angular integrand (normalized) in expression (10a) for different electron distributions: a) isotropic; b) transverse gaussian; c) beam-like gaussian. The integrands themselves are shown by thick dashed lines, while solid lines are the corresponding gaussian fits.


- 26 -

Fig. 2.-- Calculated parameters of the gyrosynchrotron emission for the isotropic electron distribution; adopted parameters are given in the text.


- 27 -

Fig. 3.-- Calculated parameters of the gyrosynchrotron emission for the transverse gaussian electron angular distribution. In two lower panels, the dash-triple-dotted line shows the corresponding parameters for the isotropic distribution calculated using continuous code.


- 28 -

Fig. 4.-- Calculated parameters of the gyrosynchrotron emission for the beam-like electron distribution. In two lower panels, the dash-triple-dotted line shows the corresponding parameters for the isotropic distribution calculated using continuous code.


- 29 -

Fig. 5.-- Time required to calculate the intensity and polarization of the gyrosynchrotron X O emission (that is, the parameters jf , jf , X , O ) at a given frequency. The simulation parameters are the same as at Figure 3. 2 GHz Intel Pentium processor was used for the calculations.

Fig. 6.-- The angular integrand (normalized) in expression (10a) for the loss-cone distribution with c = 80 and Å = 0.05 (for the X- and O-modes).


- 30 -

Fig. 7.-- Calculated spectrum for the least favorable parameter combination of gyrosynchrotron emission (X- and O-modes) for the loss-cone distribution with discontinuous second angular derivative of the distribution function.

Fig. 8.-- Gyrosynchrotron emission for the transverse gaussian electron distribution calculated using different variants of the continuous code.


- 31 -

Table 1. Summary of the GS codes Code title Continuous Q-optimized continuousb R-optimized continuousb Optimized hybridc Hybridc Approximate discrete Exact
a

Approx. useda ds, WH, LI, ln Q ds, WH, LI ds, WH, LI, ln Q No / WH / ds, WH, LId No / ds, WH, LId WH No

Rel. time 1 1.6-1.8 2.0-2.5 10e 10e long longest

Accuracy < 5 - 10%, no harmonics < 3 - 5%, no harmonics < 3 - 5%, no harmonics exact / as exact / < 2 - 4% exact / < 2 - 4% as exact at f 10fBe exact

The following approximations can be used: "ds"--integration over harmonics instead of summation; "WH"--approximated Bessel functions; "LI"--Laplace integration; "ln Q"--the terms (ln Q) and (ln Q) are neglected. R-optimization and Q-optimization can be applied simultaneously which increases both the accuracy and the computation time.
c b

In these codes, different approximations are applied at different emission frequencies.

In the continuous (high-frequency) part of spectrum, one can choose to neglect the derivatives of ln Q as well. However, this makes only a minor increase of the computation speed while decreasing the accuracy.
C C This speed estimation was made for fcr = 12fBe . With an increase of fcr , the computation time for the optimized hybrid and hybrid codes increases and approaches the computation time for the approximate discrete and exact codes, respectively. e

d


- 32 - B. Co de implementation

Fast GS codes were implemented as the Windows dynamic link libraries and Linux shared ob jects callable from IDL. Basically, we developed two variants of the libraries with the following file names: Ç libGS Std HomSrc C Ç libGS Std HomSrc CEH The file extension depends on the target operation system. In the above names, Std refers to "standard" or pre-defined distribution functions (we are working on developing similar codes for arbitrary distributions given e.g. by arrays of values), and HomSrc refers to "homogeneous source" (we are working on developing the codes with numerical integration of the radiation transfer equation for inhomogeneous sources). The last letters in the file names refer to the supported code modes: continuous (C), exact + approximate discrete (E), and hybrid (H). Currently, the set of pre-defined distributions includes nine types of distributions over energy and five types of distributions over pitch-angle; any combination of the energy and angular distribution is possible. The following energy distributions are implemented: Ç thermal (THM); Ç single power-law over kinetic energy (PLW); Ç double power-law over kinetic energy (DPL); Ç thermal/nonthermal over kinetic energy (TNT); Ç kappa (KAP); Ç power-law over the absolute value of momentum (PLP); Ç power-law over Lorentz factor (PLG); Ç thermal/nonthermal over the absolute value of momentum (TNP); Ç thermal/nonthermal over Lorentz factor (TNG). The following pitch-angle distributions are implemented: Ç isotropic (ISO);


- 33 - Ç exponential loss-cone (ELC); Ç gaussian loss-cone (GLC); Ç gaussian (GAU); Ç supergaussian (SGA). The detailed description The calling sequence for keys to switch between t above libraries are given of the pre-defined distributions is given in the online supplement I. our GS codes (for IDL users), input and output parameters, and he different code versions as well as minor differences between the in the online supplement II.


- 34 -

Fig. 9.-- Model (solid) and observed (crosses) radio spectra for the solar flare of 31 December 2007. The spectra are plotted by IDL program Flare071231.pro available as on-line supplement; the list of simulation parameters is given in this program. The data point at 34 GHz is kindly provided by Dr. S. M. White.