MONTE CARLO CALCULATIONS OF COSMIC RAY ELECTRON AND NUCLEI DIFFUSION IN THE GALAXY - A COMPARISON OF DATA AND PREDICTIONS

W.R. Webber1 and J.M. Rockstroh2

1. Astronomy Department, New Mexico State University, Las Cruces, NM.

2. 33 Fairway Drive, Hudson, N.H.

ABSTRACT

We have made detailed calculations of cosmic ray electron and nuclei spectra using a Monte Carlo program for diffusion and energy changes during propagation in the galaxy. A comparison of the results with data indicates that all of the secondary to primary ratios, 3He/4He, B/C and (Z=21-23)/Fe can be fit with the same set of diffusion parameters. These parameters also reproduce well the measured electron spectrum. Instead of specifying a path length in g/cm2 to fit the observation as required in leaky box propagation models, the Monte Carlo calculations require physical and observable parameters. For example the best fit parameters require the line integral of n, the total hydrogen column density perpendicular to the disk, to be þ15 X 1020 cm-2 which is close to the measured value.

INTRODUCTION

The currently most popular model to describe cosmic ray propagation in the galaxy is the so called leaky box model (LBM) (Englemann et al., 1990). This model assumes a uniform distribution of sources within a uniform volume. Diffusion and boundary effects are taken into account by means of simple loss terms. This model has a minimum of physical reality, but since it can be modeled accurately it provides a very useful set of predictions that are representative of a simple exponential distribution of path lengths.

Diffusion models are the next step in trying to describe this propagation and have over the years become ever more sophisticated. For example a recent diffusion model (Webber, Lee and Gupta, 1992) included a thin matter disk and a hypothetical source distribution along with convection effects to describe cosmic ray diffusion and escape perpendicular to the disk of the galaxy. But these analytical models are limited in the description of the many parameters and their spatial and energy dependence that are needed to fully describe cosmic ray propagation in the galaxy. We have now developed a Monte Carlo model (Webber, 1993a) which incorporates realistic parameters to describe this motion and have calculated the spectra of both electrons and primary and secondary nuclei along the disk and as a function of distance perpendicular to the disk. This paper gives a brief summary of this model and examples of results for both electrons and nuclei.

THE MODEL

If one considers a one-dimensional transport equation with an arbitrary diffusion coefficient K(Z,E) then þn/dt = S(Z) + þ/þz (K þN/þz) where n is the particle density, S is a source function taken here to be in the Z direction perpendicular to the galactic disk. This equation is identical to the equation obtained from a one dimensional random walk in position Z with a diffusive step size þ given by (þZ)d = (2K þt)« = þ. The distribution of sources is simulated by choosing starting positions Zs for each of the injected particles. The cosmic ray diffusion is described by a Monte Carlo process in which the particles position at time tn+1 is related to its position at tn by Zn+1 = Zn +aþ where a is a random sign. The particle is allowed to escape when it reaches a boundary, ZB = ñL.

The description of this motion through physical space and energy space is carried out as follows in our program. 104 particles are injected at each energy which is part of an energy grid of 50 energies covering 5 decades in energy. The diffusion mean free path or step size is taken in units of the boundary distance L. Typical values of L/þ may range from 103-104. The distance L and the step size þ are chosen to give a particular diffusion coefficient. For example, in our nominal case Ko = 2 X 1028 cm2 s-1 at þ1 GeV so that the step size þ = 2Ko/v = 1.5 X 1018 cm and if L is set equal to 6 X 1021 cm (2 Kpc) then there are 4 X 103 steps to the boundary. A particle is injected and allowed to random walk a certain number of steps (a big step) - typically 0.25 X the number of steps to the boundary - in this case 103 steps. At that point its location is evaluated and the probability of all of the possible things that could have happened to the particle is determined. Then the next particle is injected and both particles are followed to the next step where the process is repeated. This process is continued until 104 particles have been injected at each energy. The total time for this overall process is always taken to be > ç = L2/2K, the effective diffusion lifetime of the cosmic rays. In this way the particle distribution reaches an equilibrium situation at all energies. The diffusion coefficient (step size) is allowed to vary with both energy and distance from the plane.

After all 104 particles have been injected at each energy the number of particles in each energy bin is evaluated and weighted by a chosen input spectrum. The total number of particles in various Z bins and the time distribution of particles (in effect the path length distribution) are also evaluated.

For each big step all the various energy loss processes that could have happened to the particle are evaluated. For electrons these include synchrotron loss, inverse Compton loss, bremstrahlung loss and ionization loss as well as possible re-acceleration gain, all with their proper energy and Z dependence. For nuclei ionization loss and re-acceleration gain are included. These effects may move the particle to a new energy interval for the next step. For both electrons and nuclei a convection velocity V(Z) is also an option. And for nuclei both nuclear interactions and radioactive decay probabilities are also included. A primary nucleus may produce one secondary nucleus in the interaction (eg. CB) and this secondary nucleus is also followed throughout its lifetime including interactions and possible radioactive decay.

THE CALCULATIONS

Initially we take a source distribution located at Z = 0 and a spectral index S = -2.3 for both nuclei and electrons. The various parameters are specified at Z = 0 and their Z dependence is specified. For example for the B field, Bo = 6uG and has a Z dependance þexp -Z/0.5ZB. For the matter density no = 1.2 atoms/cm3, with a Z dependance þexp -Z/0.2ZB.

The first series of calculations involved keeping these basic parameters fixed and varying the diffusion coefficient K, and boundary distance L so as to get diffusion lifetimes ranging from þ12 M years to 150 M years (all at the nominal energy þ1 GeV/nuc). This corresponds to effective path lengths ranging from þ6g/cm2 to 30g/cm2 at this energy although we should note that in terms of the model the representative quantity is the line integral no x þn(z)dz, of hydrogen - a measurable quantity. The diffusion coefficient is taken to be þP0.6 above P = 3 GV and constant below that rigidity.

In Figures 1, 2 and 3 we show respectively the 3He/4He ratio, the (Z = 21-23)/Fe ratio and the electron spectrum for the four lifetimes (path lengths) considered. For the 3He/4He ratio, lifetime (2) corresponding to þ25 M years (þ10-12g/cm2) at þ1 GeV/nuc gives a good fit. This lifetime also gives a good fit to the B/C ratio (not shown) and the (Z = 21-23)/Fe ratio and also to the measured radioactive isotope 10Be surviving fraction of 0.20ñ0.05 (not shown). This is similar to taking a path length dependence þ = 31.6þR-0.60 > 3GV in the leaky box model which fits the B/C ratio but not the (Z = 21-23)/Fe ratio. The need to take different path lengths to fit different secondary to primary ratios (according to the interaction M.F.P. of the primaries) in the leaky box model seems to be obviated in the M.C. picture. Instead one sees that the production of Fe sec depends only weakly on the total matter traversed whereas the production of 3He (and also B), with a much longer M.F.P., depends strongly on the total matter traversed. These nuclei thus probe the extent of the Z distribution of matter whereas the Fe sec probe the more local distribution of matter.

The four different calculations imply quite different line integrals of matter in the Z direction because of the different values of L chosen. Some of the parameters appropriate to the four calculations are shown in Table 1. The model that fits the nuclei data best implies a matter line integral = 1.4 x 1021 cm2 which compares with values in the range of 1.0-1.5 x 1021 cm2 for the measured values of the sum of HI, H2 and HH summarized by Strong and Youssefi, 1995.

In order to compare electrons, since the absolute intensities can be arbitrarily normalized by the source function, we have normalized the different spectra at 1 GeV in Figure 3. At the high energy end the effect of the longer lifetimes is clearly seen in the steepening of the spectrum. Here the choice of the B field and the resulting synchrotron loss plays an important role - but the same set of diffusion parameters that fit the nuclei data also gives the best fit to the high energy electron spectrum. At the lower energies around þ10 MeV the different lifetimes make a factor >3 difference in the electron intensity. This should be observable directly from the intensity of bremstrahlung gamma rays at this energy as measured by COMPTEL.

TABLE 1

Diffusion parameters used in M.C. Calculations

K(cm2 s-1) L(Kpc) ç(106 years) þn-dþ (X 1020/cm2)