A&A 402, 971-983 (2003)
DOI: 10.1051/0004-6361:20030318

Spatial origin of Galactic cosmic rays in diffusion models

I. Standard sources in the Galactic disk[*]

R. Taillet1,2 - D. Maurin1,3

1 - Laboratoire de Physique Théorique LAPTH, 74941 Annecy-le-Vieux, France
2 - Université de Savoie, 73011 Chambéry, France
3 - Institut d'Astrophysique de Paris, 98bis Bld Arago, 75014 Paris, France

Received 6 December 2002 / Accepted 25 February 2003

The propagation of Galactic cosmic ray nuclei having energies between 100 MeV/nuc and several PeV/nuc is strongly believed to be of diffusive nature. The particles emitted by a source located in the disk do not pervade the whole Galaxy, but are rather confined to a smaller region whose spatial extension is related to the height of the diffusive halo, the Galactic wind and the spallation rate. Following the pioneering work of Jones (1978), this paper presents a general study on the spatial origin of cosmic rays, with a particular attention to the role of spallations and Galactic wind. This question is different, and to a certain extent disconnected, from that of the origin of cosmic rays. We find the regions of the disk from which a given fraction of cosmic rays detected in the solar neighborhood were emitted (f-surfaces). After a general study, we apply the results to a realistic source distribution, with the propagation parameters obtained in our previous systematic analysis of the observed secondary-to-primary ratios (Maurin et al. 2002a). The shape and size of these f-surfaces depend on the species as well as on the values of the propagation parameters. For some of the models preferred by our previous analysis (i.e. large diffusion slope ${\delta }$), these f-surfaces are small and in some extreme cases only a fraction of a percent of the whole Galactic sources actually contribute to the solar neighborhood cosmic ray flux. Moreover, a very small number of sources may be responsible for more than 15% of the flux detected in the solar neighborhood. This may point towards the necessity to go beyond the approximations of both homogeneity and stationarity. Finally, the observed primary composition is dominated by sources within a few kpc.

Key words: ISM: cosmic rays

1 Introduction

The propagation of charged cosmic ray nuclei, in the energy range going from a few 100 MeV/nuc and a few PeV/nuc, is strongly affected by the Galactic magnetic field. It is a diffusive process, so that the cosmic rays emitted by a single source spread out in time, pervade the whole Galaxy, and can escape the Galaxy when reaching its boundaries. Those coming from a source located far from the Sun have a larger probability of escaping than reaching the solar neighborhood. It is the opposite for nearby sources, so that the cosmic ray fluxes in the solar neighborhood are more sensitive to the properties of the local sources (as opposed to the remote sources). Other effects like spallations and Galactic wind further limit the distance cosmic rays travel before being detected. Some consequences of the Galactic wind were studied in Jones (1978) where convective escape was compared to escape through the top and bottom boundaries of the Galaxy.

The goal of this paper is to go one step beyond by providing a general study on the spatial origin of cosmic rays, i.e. to answer the question "from which region of the Galaxy were emitted the cosmic rays detected in the solar neighborhood?''. This question is different, and to a certain extent disconnected, from that of the origin of cosmic rays ("What are the astrophysical objects which are responsible for the acceleration of cosmic rays?'') which is still much debated. We believe that it is nevertheless an interesting question, for several reasons. First, we find that the answer may cast some doubt on the validity of the stationary model, upon which most studies on cosmic rays are based. Second, it gives some clues about the spatial range beyond which the cosmic ray studies are blind to the sources. Finally, this study may be of interest to optimize the propagation codes based on Monte-Carlo methods, by focusing the numerical effort on the sources that really contribute to the detected flux.

The reader who does not want to go through the pedagogical progression can go directly from the general presentation of the method in Sect. 2 to its application in realistic cases in Sect. 7. For the others, the effect of escape is studied in Sect. 3 and that of spallations and Galactic wind is studied in Sect. 5. Then, Sect. 6 studies the effect of a realistic source distribution. Finally, the fully realistic case is considered in Sect. 7. The results and the perspectives are discussed in the last section. For convenience, we will use the word f-surfaces to describe the surfaces in the thin disk within which the sources form the fraction f of cosmic rays detected at the observer location.

2 Description of the method

A stationary point source emits particles that diffuse in a given volume. At the boundaries of this volume, the particles are free to escape and the density drops to zero. After a sufficiently long time, the stationary regime is eventually reached and the density profile is established inside the diffusive volume. If several sources are present (or even a continuous distribution of sources), their contributions add linearly at each point.

The question we wish to answer is the following: a cosmic ray being detected at the position $\vec{r}_{\rm o}$ of an observer (in practice, this will be the position of the Sun, and we refer to this position as the solar neighborhood), what is the probability density

\begin{displaymath}\frac{{\rm d}{\cal P} \left\{{\rm emitted:} ~
\vec{r}_{\rm s...
...\rm s} \vert \vec{r}_{\rm o} \right\}}{{\rm d}\vec{r}_{\rm s}}
\end{displaymath} (1)

that this cosmic ray was emitted from a source located at the position $\vec{r}_{\rm s}$? Such a question falls among classical problems of statistics. A rigorous theoretical frame is provided by the Bayes approach that summarizes the proper use of conditional probabilities. A cruder but sufficient (and equivalent) treatment is given by the frequency interpretation. The probability written above is simply given by

\begin{displaymath}\frac{{\rm d}{\cal P} \left\{ \vec{r}_{\rm s} \vert \vec{r}_{...
...m s}}
{{\cal N} \left[ \rightarrow \vec{r}_{\rm o} \right] },
\end{displaymath} (2)

where ${\cal N} \left[ \rightarrow \vec{r}_{\rm o} \right]$ is the number of paths reaching $\vec{r}_{\rm o}$ and ${\rm d}{\cal N} \left[ \vec{r}_{\rm s}\rightarrow \vec{r}_{\rm o} \right]/{\rm d}\vec{r}_{\rm s} $ is the density of paths going from $\vec{r}_{\rm s}$ to $\vec{r}_{\rm o}$. We finally notice that the latter number determines the density of cosmic rays that reach the position $\vec{r}_{\rm o}$, when a source is placed at $\vec{r}_{\rm s}$. We can thus write

 \begin{displaymath}\frac{{\rm d}{\cal P}\left\{ \vec{r}_{\rm s} \vert \vec{r}_{\...
...d}\vec{r}_{\rm s}}
\equiv N_{\rm r_{\rm s}}(\vec{r}_{\rm o}),
\end{displaymath} (3)

where the density $N_{ r_{\rm s}}(\vec{r}_{\rm o})$ is the solution of the propagation equation for a point source located at $\vec{r}_{\rm s}$. The normalization factor in this relation is obtained by imposing that ${\rm d}{\cal P}/{\rm d}{\vec{r}_{\rm s}}$ actually is a probability density, i.e. is normalized to unity. We refer to the contours on which the probability density is constant as isodensity contours.

If the sources are distributed according to $w(\vec{r}_{\rm s})$, the probability that a cosmic ray detected at $\vec{r}_{\rm o}$ was emitted from a surface ${\cal S}$ is given by

 \begin{displaymath}{\cal P} \left\{ {\cal S} \vert \vec{r}_{\rm o} \right\}=
... N_{ r_{\rm s}}(\vec{r}_{\rm o})
{\rm d}\vec{r}_{\rm s}}\cdot
\end{displaymath} (4)

This probability contains all the physical information about the spatial origin of cosmic rays. We define the f-surfaces, inside which the sources contribute to the fraction f of the detected flux, by the relation ${\cal P} \left\{ {\cal S} \vert \vec{r}_{\rm o} \right\}=f$. Actually, even for a given value of f, there are many different surfaces, delimited by different closed contours, fulfilling this condition. We focus on the smallest of these surfaces, which is precisely delimited by an isodensity contour. We also use the term $r_{\rm lim}$-probability for the quantity ${\cal P} \left\{ r_{\rm s} < r_{\rm lim} \vert \vec{r}_{\rm o} \right\}$.

3 The escape through the diffusive volume boundaries

The region in which diffusion occurs is limited by surfaces (hereafter the boundaries) beyond which diffusion becomes inefficient at trapping the particles, so that they can freely escape at a velocity close to c. The density outside the diffusive volume is very small, and it is very reasonable to suppose that the boundaries are absorbers, i.e. they impose a null density (N=0).

It is well-known that the shape and location of the boundaries play a crucial role for diffusive propagation. This section shows that the cosmic rays emitted from standard sources in the disk are not sensitive to the radial extension of the Galaxy, but only to its top and bottom edge. To this aim, it is sufficient to concentrate on pure diffusion and to neglect spallations, the Galactic wind and reacceleration. Indeed this is a conservative case as these effects can only make the diffusion process even less sensitive to the presence of the boundaries (see below). Moreover, we consider the case of a homogeneous source distribution located in the disk $w(\vec{r}_{\rm s})\propto \delta(z)$, which also leads to a conservative result if compared to a realistic radial distribution of sources.

We first consider the pure diffusion equation with a Dirac source term

 \begin{displaymath}-K\triangle N(\vec{r})=\delta(\vec{\vec{r}-\vec{r}_{\rm s}}).
\end{displaymath} (5)

In unbounded space, the solution is given by $N_{r_{\rm s}}(r_{\rm o}) =1/4\pi K \vert\vert\vec{r_{\rm o}}-\vec{r_{\rm s}}\vert\vert$. The influence of the boundaries is estimated by solving this equation in three situations: first we consider only a side boundary, then only a top plus bottom boundary, and finally all the boundaries.

3.1 Boundaries influence

\par\includegraphics[width=4.5cm,clip]{3376f1.eps} \end{figure} Figure 1: Geometry of the diffusive volume.
Open with DEXTER

Our Galaxy can be represented as a cylindrical box with radial extension R and height L (see Appendix A for further details). The probability density ${\rm d}{\cal P}_{\rm cyl}(\vec{r}_{\rm s}\vert\vec{r}_{\rm o})/{\rm d}\vec{r}_{\rm s}$can be computed for arbitrary source and observer positions $\vec{r}_{\rm s}$ and $\vec{r}_{\rm o}$, using a Fourier-Bessel decomposition of the density. In our case, the observer is located near the Sun, at a Galactocentric distance $R_\odot \sim 8.5$ kpc. Unless the diffusive halo height is very large, the top and bottom boundaries located at $z=\pm L$ are nearer to us than the side boundary located at R=20 kpc. As a result, we expect the effect of the side boundary to be smaller. The first simplified situation we consider is that of an observer located at the center of the Galaxy. (In the case of an infinite disk, i.e. $R\rightarrow \infty $, this amounts to a mere redefinition of the origin of the disk).

With $\vec{r}_{\rm o} = \vec{0}$, the solution for a point source in this particular geometry is given in Appendix A. The probability density that a particle reaching the observer was emitted from a point located at a distance $r_{\rm s}$ from the center is thus given by (with $\rho_{\rm s}\equiv r_{\rm s}/R$)

$\displaystyle {\rm d}{\cal P}_{\rm cyl}(\vec{r}_{\rm s}\vert O)=\frac{{\rm d}^2...
...m_{i=1}^\infty \frac{\tanh(\zeta_iL/R)}
{\zeta_i^2 J_1(\zeta_i)}
\right\}^{-1},$     (6)

normalization being obtained by imposing $\int_{0}^R
{\rm d}{\cal P}(r_{\rm s}\vert O)= 1$. The $r_{\rm lim}$-probability is given by
$\displaystyle {\cal P}_{\rm cyl} (r_{\rm s} < r_{\rm lim}\vert O )$ = $\displaystyle \left\{
\sum_{i=1}^\infty \frac{J_1(\zeta_i r_{\rm lim}/R)}
...=1}^\infty \frac{\tanh(\zeta_iL/R)}
{\zeta_i^2 J_1(\zeta_i)}
\right\}^{-1}\cdot$ (7)

This probability is independent of the value of the diffusion coefficient K.

3.1.1 Side boundary

\par\includegraphics[width=3.6cm,clip]{3376f2.eps} \end{figure} Figure 2: Geometry of the diffusive volume in the limit $L\rightarrow \infty $.
Open with DEXTER

The escape from the side boundary (located at r=R) is disentangled from the escape from the $z=\pm L$ boundary by first considering the limit $L\rightarrow \infty $. For the sake of simplicity, we will, as above, only study the effect of this boundary on observations performed at the center of the Galaxy. In the limit $L\rightarrow \infty $, we have $\coth(\zeta_iL/R)\approx 1$ in expression (6). This gives for the $r_{\rm lim}$-probability,

\begin{displaymath}{\cal P}_{\rm R} (r_{\rm s} < r_{\rm lim}\vert O ) =
\sum_{i=1}^\infty \frac{1}
{\zeta_i^2 J_1(\zeta_i)}

3.1.2 Top and bottom boundaries

The influence of the $z=\pm L$ boundaries, in the case of an infinite disk ( $R\rightarrow \infty $) is now considered. In this limit, the sum over Bessel functions can be replaced by an integral and we obtain (see Appendix B.3)

 \begin{displaymath}{\rm d}{\cal P}_{\rm L}(\vec{r}_{\rm s}\vert O) \propto
~ \tanh \left( \frac{xL}{r_{\rm s}} \right) ~ {\rm d}x,
\end{displaymath} (8)

which allows to compute the $r_{\rm lim}$-probability ${\cal P}_{\rm L} (r_{\rm s} < r_{\rm lim}\vert O )$ as before, which is a function of $r_{\rm lim}/L$ only. These integrals are somewhat intricate to compute numerically, due to the very slow convergence. In this particular case, the accuracy of the numerical calculation can be checked for $r \gg L$, as a detailed study of the function (8) shows that in this limit

N^{\rm L}_{(r_{\rm s},0)}(O) \approx \frac{1}{4\pi K r_{\rm...
2\sqrt{\frac{r_{\rm s}}{L}}
~ {\rm e}^{-\pi r_{\rm s}/2L}.
\end{displaymath} (9)

It is also noticeable that the quantity

\begin{displaymath}f_{\rm esc}(r_{\rm s}) \equiv 1- \frac{N^L}{N^{L=\infty}}
= ...
...-\tanh \left( \frac{xL}{r_{\rm s}} \right)
\right\}~ {\rm d}x

gives the fraction of cosmic rays emitted from a distance $r_{\rm s}$ that has escaped the diffusive halo before reaching us.

\par\includegraphics[width=4.4cm,clip]{3376f3.eps} \end{figure} Figure 3: Geometry of the diffusive volume in the limit $R\rightarrow \infty $.
Open with DEXTER

3.2 Summary: The effect of boundaries on primary species

Figure 4 shows the probability density computed above as a function of $r_{\rm s}$ for unbounded space, for the cylindrical geometry with several halo sizes L, i.e. Eq. (7), and for the two limiting cases corresponding to $L\rightarrow \infty $ or $R\rightarrow \infty $.

\par\includegraphics[width=8.8cm,clip]{3376f4.eps} \end{figure} Figure 4: Cosmic ray probability density as a function of $r_{\rm s}$ (distance of the $\delta(\vec{r_{\rm s}})$ source in the disk), for several values of L and for a disk of radius R=20 kpc. Big stars are for unbounded model, dotted line is for a spherical boundary at radius R, small stars are for top and bottom boundaries, and solid lines are for cylindrical boundaries.
Open with DEXTER

We also show, in Table 1, the radii of the f-surfaces, in the two cases R= 20 kpc and $R\rightarrow \infty $. It can be noticed that even if the source distribution is infinite in extent, the finite size of the halo limits the quantity of cosmic rays that reach a given point. The mean distance from which the cosmic rays reach the center is given by $\langle r_{\rm s} \rangle = 1.4 ~ L$. This effect dominates over the leakage through the side boundary $r_{\rm s}=R$, and it will be even more negligible in realistic situations, as (i) the source density is small near the edge of the disk, (ii) when the spallations and Galactic wind are considered, most cosmic rays are destroyed or blown out of the disk before they have a chance to reach this side boundary.

Table 1: This table indicates the radius $r_{\rm lim}$ inside which a given fraction $f(r_{\rm lim})\equiv {\cal P}(r_{\rm s} < r_{\rm lim}\vert O )$ of cosmic rays reaching the center were emitted from, for several L and in the case of the infinite disk and R=20 kpc.

An important consequence is that as long as the observer and the sources are not too close to the side boundary, the density only depends on the relative distance to the source in the disk, so that it may be assumed, for numerical convenience, that the observer is either at the center of a finite disk, or in an infinite disk. In all the paper, i.e. for standard sources in the disk, we will consider the limit $R\rightarrow \infty $, i.e. we use the integral representation described in Appendix B.3.

4 Secondary and radioactive species

4.1 Progenitors of stable secondaries

As can be seen in Appendix A.3, the secondary distribution from point-like primary sources is related very simply to the primary distribution itself. One could find strange to speak about secondaries as we have not, for the moment, included spallations in the model. The right picture is the following: a primary emitted at $r_{\rm s}$ propagates and from time to time crosses the disk (mostly filled with hydrogen, density $n_{\rm ISM}$). During this crossing, there is a probability $n_{\rm ISM}.v.\sigma_{{\rm prim}\rightarrow{\rm sec}}$to create a secondary, that in turn propagates in the diffusive volume until it reaches (or not) the experimental setup. This will be taken into account properly in the next section. However, in order to have a compact expression, a crude estimation can be obtained by neglecting the influence of spallations on the primary and secondary component. This is obtained if one discards $\Gamma_{\rm inel}$ in the terms $A_i^{\rm prim}$ and $A_i^{\rm sec}$ of Eq. (A.6). The net result will be an overestimation of the distance the secondaries come from since their destruction is discarded two times; once under their primeval primary form and once in their secondary form.

We find, in the case $R\rightarrow \infty $ (see Appendix B.3), and for a homogeneous distribution of sources,

\begin{displaymath}\frac{{\rm d}{\cal P}_{\rm sec}(r_{\rm s}\vert R_\odot)}{{\rm...
...times \tanh^2 \left( \frac{x L}{r_{\rm s}} \right) ~ {\rm d}x.

The resulting integrated probabilities are shown in Table 2. The source of the primary that will give the secondaries observed at a given point is located farther away than the sources of the primary we detect (compare Tables 2 and 1). This may be of importance if for instance the source composition or the source intensity varies with position: in the ubiquitous secondary-to-primary ratio, the numerator is sensitive to sources located on a greater range than the denominator. Moreover, these secondaries set the size of an effective "local'' zone outside of which the particles reaching the solar neighborhood have never been. The local observations tell nothing about the propagation conditions outside of this zone. One could object that this conclusion is mainly based on the f-surfaces which refer to the sources contributing to observed CR, but that the cosmic rays reaching us from these sources actually sample (via random walk) a much larger volume. This is actually not the case, as a particle wandering too far has a very small probability to ever come back to us. This point can be made more quantitative, as a simple reasoning shows that the probability ${\cal P}[ACB]$ that a particle emitted in A and reaching B has passed through C is given by ${\cal P}[AC]{\cal P}[CB]$, which is closely related to $N_A(\vec{r}_C) \times N_B(\vec{r}_C)$. This later quantity is small as soon as C is too far from A or B.

Table 2: This table indicates the radius inside which a given fraction $f(r_{\rm lim})$ of secondary cosmic rays reaching the center were emitted from, for several L and in the case of a disk of radius R=20 kpc. The last line shows that for small L, the effect of the side boundary is completely negligible.

4.2 Radioactive secondaries

In the case of an unstable species with a lifetime $\tau$, formula (A.5) can be written as

 \begin{displaymath}{\rm d}{\cal P}_{\rm rad} \left\{ r_{\rm s} \vert O \right\} ...
\sqrt{R^2 \Gamma_{\rm rad}/K + \zeta_i^2} ~ J_1^2(\zeta_i)},
\end{displaymath} (10)

where $\Gamma_{\rm rad}= \tau^{-1} = \gamma^{-1}\tau_0^{-1} $. This expression can be transformed using the identity (Lebedev 1972)

\begin{displaymath}\frac{1}{\sqrt{\zeta_i^2 + \alpha^2}} =
\int_0^\infty \frac{...
... e}^{-\alpha \rho}}{\rho} \rho J_0(\zeta_i
\rho) {\rm d}\rho.

The approximation in the last step is valid if the exponential term decreases with $\rho$ fast enough (i.e. $\alpha$ is large so that the upper limit can be set to 1 in the integral). We then recognize in (10) the Fourier-Bessel transform of $\exp(-\alpha \rho)/\rho$, so that finally the normalized probability reads

 \begin{displaymath}{\rm d}{\cal P}_{\rm rad} \left\{ r_{\rm s} \vert O \right\} ...
...}{2\pi ~r_{\rm s} .~ l_{\rm rad}}
~ {\rm d}^2\vec{r}_{\rm s},
\end{displaymath} (11)

where the following typical length has been introduced

 \begin{displaymath}l_{\rm rad}= \sqrt{\frac{K}{\Gamma_{\rm rad}}} =
0.17 \; \mb...
...\mbox{Myr}^{-1}} }
\sqrt{\frac{\tau}{ 1 \; \mbox{Myr}} }\cdot
\end{displaymath} (12)

Indeed, this result can be derived much more straightforwardly starting from the stationary equation $ -K\Delta_{\rm r} N(r) +\Gamma_{\rm rad}N(r)=0$ (with a source at the origin) in unbounded space. This is also in full agreement with the expression given in Appendix B (see also Sect. 4.1) of Donato et al. (2002), where we found the same expression starting from the propagator of the non-stationary diffusion equation in unbounded space.

To sum up, Eq. (11) is valid as long as $l_{\rm rad} \ll R$ and $l_{\rm rad} \ll L$: the propagation of the unstable species can be then considered as local, with a typical scale $l_{\rm rad}$. This is no longer the case if the lifetime $\tau = \gamma \tau_0$ is large, which is the case at high energy because of the relativistic factor $\gamma$, even if the proper lifetime $\tau_0$ is short. The $r_{\rm lim}$-probability is straightforwardly derived. As on these typical scales, the source distribution can safely taken to be constant, the distance $r_{\rm lim}$ is expressed as

 \begin{displaymath}r_{\rm lim}=-l_{\rm rad}\times \ln (1-f).
\end{displaymath} (13)

It means that the sources that contribute to the fraction f=(50-90-99)% of the radioactive species measured flux are located inside the disk of radius $r_{\rm lim}=(0.7{-}2.3{-}4.6)\times l_{\rm rad}$. The effect of a local underdensity around the Sun is discussed later.

4.3 Electrons and positrons

Cosmic ray sources also emit electrons and positrons. In contrast with the nuclei, these particles are light, so that they are subject to much stronger energy losses, due to synchrotron radiation and inverse Compton. This results in an effective lifetime given by (e.g. Aharonian et al. 1995) $\tau_{\rm loss} \sim 300 \; \mbox{Myr} \times (1 \; \mbox{GeV}/E)$. The results given in the previous section on radioactive species can be applied to this case, with a scale length

\begin{displaymath}r_{\rm loss} \sim 1 \; \mbox{kpc} \times \sqrt{\frac{1 \; \mb...
\sqrt{\frac{K}{0.03 \; \mbox{kpc}^2\; \mbox{Myr}^{-1}}}\cdot

Formulae (11) and (13) can be used with $l_{\rm
rad}\leftrightarrow r_{\rm loss}$. This effect is discussed by Aharonian et al. (1995) to show that a nearby source may be necessary to explain the high energy electron flux observed on the solar neighborhood.

4.4 Summary: Pure diffusive regime, an upper limit

The important conclusions at this point are that i) most of the stable primary cosmic rays that reach the solar neighborhood were emitted from disk sources located within a distance of the order of L, such that the $R=20 \; \mbox{kpc}$ boundary can reasonably be discarded ii) the secondary species composition is determined by sources located farther away than those determining the primary composition; iii) radioactive species may come from very close if their lifetime is so short that $\sqrt{K \gamma \tau_0} < L$, high energy electrons and positrons definitely do.

These conclusions are expected to be stronger when spallations, Galactic wind and a realistic source distribution are taken into account. All these effects will limit even more the range that the particles can travel before reaching the solar neighborhood.

5 The effects of spallation and convection

5.1 Pure convection

The diffusion of cosmic rays may be disturbed by the presence of a convective wind of magnitude $V_{\rm c}$, directed outwards from the disk. For numerical convenience, a constant wind has been considered, although other possibilities (especially a linear dependence) are probably more justified on theoretical grounds (see discussion in Maurin et al. 2002a). The effect is to blow the particles away from the disk, so that those detected in the solar neighborhood come from closer sources (compared to the no-wind case). With an infinite halo, the probability density in the disk is given by

                 $\displaystyle \frac{{\rm d}{\cal P}^{\rm wind}_{L \rightarrow
\infty} \left\{r_{\rm s} \vert 0 \right\}}{{\rm d}^2\vec{r}_{\rm s}}$ $\textstyle \propto$ $\displaystyle \int_0^\infty \frac{k J_0(kr_{\rm s})~ {\rm d}k}{V_{\rm c} + K
V_{\rm c}^2/K^2 + 4 k^2}}$  
  $\textstyle \propto$ $\displaystyle \frac{1}{r_{\rm s}} \int_0^\infty \frac{x J_0(x)~
{\rm d}x}{r_{\rm s}/r_{\rm w} +
(r_{\rm s}/r_{\rm w})^{2} + x^2}},$ (14)

where the characteristic radius $r_{\rm w}\equiv 2K/V_{\rm c}$has been defined. The expression in Eq. (14) is a function of $r_{\rm s}/r_{\rm w}$ only. The deviation from a pure $1/r_{\rm s}$ law, as well as deviations due to escape, radioactive decay and spallation (see next section), is shown in Fig. 5.
\par\includegraphics[width=8.8cm,clip]{3376f5.eps} \end{figure} Figure 5: Deviation from the pure $1/r_{\rm s}$ density profile, $N(r_{\rm s})/(1/4\pi K r_{\rm s})$, due to the various effects studied here: escape from the $z=\pm L$ boundaries, spallations and Galactic wind. In this latter case, the choice $r_{\rm scale}=2r_{\rm w}$ has been made to show the similar behavior at large $r_{\rm s}$. The case of a radioactive species has also been shown. It should be noticed, however, that in most interesting cases, the scale length $l_{\rm rad}$ is much smaller than the others, so that in this case the propagation is dominated by radioactive decay and spallations and Galactic wind can be safely discarded.
Open with DEXTER

The $r_{\rm lim}$-probability is given by
$\displaystyle {\cal P}^{\rm wind}_{L \rightarrow\infty}(<r_{\rm lim})
= \int_0^...
...\rm w} +
\sqrt{r_{\rm lim}^2/r^2_{\rm w} +
x^2}}{x}\right] \right\} ~ {\rm d}x.$     (15)

Some values are indicated in Table 4 and plotted in Fig. 6.

It is interesting to note that the effect of $2r_{\rm w}$ is similar (though not rigorously identical) to the effect of L (see Fig. 5). As a matter of fact, this was noticed by Jones (1978) who studied the propagation properties in a dynamical halo and provided a very simple picture (along with a rigorous derivation) of the effect of the wind. Consider a particle initially located at a distance z from the disk. It takes a time $t_{\rm diff}\approx z^2/K$ to diffuse back in the disk. In the meantime, convection sweeps the particle in a distance $z_{\rm w}\equiv
V_{\rm c} t_{\rm diff}\approx V_{\rm c} z^2/K$. Both processes are in competition and the particle will not reach the disk if $z_{\rm w}>z$. This define an effective halo size $L^*\approx K/V_{\rm c}$. This is our parameter $r_{\rm w}$ up to a factor 2.

\par\includegraphics[width=8.8cm,clip]{f6.eps} \end{figure} Figure 6: Integrated probability that a particle detected at the origin was emitted inside the ring of radius $r_{\rm s}$, in the three situations considered. The solid dark line is obtained when only the leakage through the $z=\pm L$ boundaries is considered, in which case the radii scale as $r/r_{\rm scale}$ with $r_{\rm scale} = L$. The dotted, respectively dashed, line is obtained when only the spallations, respectively only the convective wind, are considered. The solid grey line indicates the probability that the primary progenitor of a secondary detected in the solar neighborhood was emitted from within a given distance.
Open with DEXTER

5.2 Pure spallation

The Galactic disk contains interstellar gas mostly made of hydrogen. When cosmic rays cross the disk, they can interact with this gas. This interaction may result in a nuclear reaction (spallation), leading to the destruction of the incoming particle and to the creation of a different outgoing particle (secondary). We present two approaches to the problem of diffusion in presence of a spallative disk. When the halo is infinite in extent, the solution may be obtained by using the interpretation of diffusion in terms of random walks. This will be treated in Appendix C. In the general case, the Bessel developments can be used as before. Starting from Eq. (A.5), the expression for the probability density is readily obtained. The limit $L\rightarrow \infty $ is noteworthy, as the resulting expression isolates the influence of spallations:

$\displaystyle \frac{{\rm d}{\cal P}^{\rm spal}_{L \rightarrow\infty} \left\{ r_{\rm s} \vert
0 \right\}}{{\rm d}^2\vec{r}_{\rm s}}$ $\textstyle \propto$ $\displaystyle \int_0^\infty \frac{k J_0(kr_{\rm s})}{2h \Gamma_{\rm inel}
+ 2kK} ~ {\rm d}k$  
  = $\displaystyle \frac{1}{4\pi K r_{\rm s}}
\int_0^\infty \frac{x J_0(x)
{\rm d}x}{r_{\rm s}/r_{\rm sp} + x}$  

where the quantity $r_{\rm sp}\equiv K/(h \Gamma_{\rm inel})$has been defined. Would there be no spallation, the $1/r_{\rm s}$ behavior would be recovered. The term $2h \Gamma_{\rm inel}$ has the effect to kill the contributions of $k \la k_{\rm sp}$ in the integral, with $k_{\rm sp} \equiv h \Gamma_{\rm inel}/K$. It leads to a decrease of the integral on scales $r>r_{\rm sp} =
1/k_{\rm sp}$. Some typical values, for $K=\beta K_0 {\cal R}^\delta$(see Sect. 5.5) with $K_0 = 0.03 \; \mbox{kpc}^2 \; \mbox{Myr}^{-1}$and $\delta =0.6$ are given below at 1 GeV/nuc and 100 GeV/nuc. The heavy species are more sensitive to spallations, so that they come from a shorter distance. This could in principle affect the mean atomic weight of cosmic rays if the composition of the sources is not homogeneous (see e.g. Maurin et al. 2003a). See Sect. 7.1) for the results with realistic propagation parameters.

For small values of $r_{\rm s}/r_{\rm sp}$, the convergence of the previous integral is slow, and other forms obtained by integration by parts, as developed in the Appendix B.3, might be preferred. However, in this particular case, the identity

\begin{displaymath}\int_0^\infty x {\rm d}x J_0(x)/(x+\alpha)
= \int_0^\infty y {\rm d}y {\rm e}^{-\alpha y}/(1+y^2)^{3/2}

yields the more useful form

 \begin{displaymath}\frac{{\rm d}{\cal P}^{\rm spal}_{L \rightarrow\infty} \left\...
...\rm e}^{- y r_{\rm s}/r_{\rm sp}}}{(1+y^2)^{3/2}} \; {\rm d}y.
\end{displaymath} (16)

This expression is in full agreement with Eq. (C.3) obtained with the random walk approach (see Appendix C). For large values of $r_{\rm s}/r_{\rm sp}$, the convergence can be checked by comparing the results to the asymptotic development

\begin{displaymath}\int_0^\infty \frac{x J_0(x)}{x+\alpha} \; {\rm d}x
\approx ...
...\alpha^2} - \frac{9}{\alpha^4}+
\frac{225}{\alpha^6} + \ldots

Finally, the $r_{\rm lim}$- probability can be computed as before

\begin{displaymath}{\cal P}^{\rm spal}_{L \rightarrow \infty}(<r_{\rm lim})
= \...
...\left\{-x\frac{r_{\rm lim}}{r_{\rm sp}} \right\}

Some values are indicated in Table 4.

Table 3: Some values of the inelastic cross section and the associated spallation scale length.
  p O Fe
$\sigma~(\mbox{mb})$ 44 309 760
$r_{\rm sp}~(\mbox{kpc})$, 1 GeV/nuc 10.2 1.45 0.59
$r_{\rm sp}~(\mbox{kpc})$, 100 GeV/nuc 115 16.4 6.7

5.3 Comparison and combination of the different effects

To summarize, the effect of spallation and Galactic wind depends on the two parameters:

r_{\rm w} \equiv \f...
\frac{ 100 \; \mbox{mb}}{\sigma}\cdot
\end{array} \right.
\end{displaymath} (17)

The $r_{\rm lim}$-probability is displayed in Fig. 6 as a function of $r_{\rm lim}$. The effect of the Galactic wind is very similar to that of the top and bottom boundaries, whereas the effect of spallations is quite different. In the latter case, the cutoff in the density is a power law in $r_{\rm s}$ and decreases much more slowly than the exponential cutoff due to the wind or to escape. As a result, the 99%-surfaces are much larger than the 90%-surfaces. This can also be seen in the first three lines of Table 4.

When all the effects above are considered, Eq. (A.5) gives

\begin{displaymath}\frac{{\rm d}{\cal P} \left\{ r_{\rm s} \vert 0 \right\}}{{\r...
... \sqrt{\rho_{\rm w}^2 +
x^2}}{\displaystyle \rho_L}\right\} }

where $\rho_{\rm sp} = r/r_{\rm sp}$, $\rho_{\rm w} = r/r_{\rm w}$ et $\rho_L = r/L$. The smallest of these three numbers indicates the dominant effect. Various $r_{\rm lim}$-probabilities are shown in Table 4.

Table 4: This table indicates the radius of several f-surfaces, for several values of L, $V_{\rm c}$ and $\sigma _{\rm inel}$. We have introduced $L_{10}\equiv L/10$ kpc, $V_{10} \equiv V_{\rm c}/10$ km s-1 and $\sigma _{100}\equiv \sigma /100$ mb.

For a radioactive species, the spallations and the Galactic wind have a negligible effect on propagation as long as $l_{\rm rad}$(see Sect. 4.2) is smaller than L, $r_{\rm sp}$ and $r_{\rm w}$.

5.4 The number of disk-crossings in the general case

Several properties (energy losses, amount of reacceleration, secondary-to-primary ratio) of the cosmic ray flux detected in the solar neighborhood are determined by the number of times a given cosmic ray has crossed the disk since it was created. The distribution of disk-crossings is computed in Appendix C in the case of an infinite diffusive volume and in the absence of Galactic wind. In the most general situation, the mean number of crossings (though not the entire distribution of crossing numbers) can be computed as follows. Each time a particle crosses the disk, it has a probability $p=2h
\sigma_{\rm inel} n_{\rm ISM}$ of being destroyed by a spallation. The number N(r) of surviving particles can thus be obtained from N0(r), the number of particles diffusing without spallations, as

\begin{displaymath}N(r) = N_0(r) \times \left( 1- p\right)^{n_{\rm cross}},

so that the number of crossing is readily obtained from the densities with and without spallations as

 \begin{displaymath}n_{\rm cross}(r) = \frac{\ln (N(r)/N_0(r))}
{\ln( 1 - p)}\cdot
\end{displaymath} (18)

Notice that this expression applied to Eq. (16) leads to Eq. (C.2) when $L\rightarrow \infty $, $V_{\rm c} = 0$, and when p is small. As the surface density of the disk is $2hn_{\rm ISN} \sim 10^{-3} \; \mbox{g}
\; \mbox{cm}^{-2}$, the mean column density crossed by the particle (called grammage) is given by

\begin{displaymath}\Sigma (r_{\rm s}) = n_{\rm cross} (r_{\rm s})\times 2hn_{\rm...
\times \frac{n_{\rm cross} (r_{\rm s})}{10^4}\cdot

The evolution of the grammage with the distance of the source is displayed in Fig. 7. The effect of escape, spallations and Galactic wind is shown.
\par\includegraphics[width=8.8cm,clip]{f7.eps} \end{figure} Figure 7: Grammage crossed as a function of the origin, for some of the models discussed in the text and for a typical value of K=0.03 kpc2 Myr-1.
Open with DEXTER

As a cross-check, it can be noticed that in this approach, the mean grammage

\begin{displaymath}\langle \Sigma \rangle_{\rm spatial}=\int \Sigma (r_{\rm s})
...0 \right\}}{{\rm d}^2\vec{r}_{\rm s}} {\rm d}^2\vec{r}_{\rm s}

yields the right order of magnitude for the usual grammage derived from leaky box analysis ($\sim $9 g cm-2). Moreover, the knowledge of $n_{\rm cross} (r_{\rm s})$ allows to estimate the magnitude of energy losses and reacceleration rates.

5.5 The energy dependence

The diffusion coefficient actually depends on energy. A commonly used form (see Maurin et al. 2002b for a discussion) is

\begin{displaymath}K = K_0 \beta \left( \frac{\cal R}{\mbox{1 GV}} \right)^\delta

where ${\cal R}$ stands for the rigidity, $K_0 \sim 0.01 {-} 0.1
\; \mbox{kpc}^2 \; \mbox{Myr}^{-1}$ and $\delta \sim 0.3{-}1$. The previous results were given for $K = 0.03 \; \mbox{kpc}^2 \; \mbox{Myr}^{-1}$, typical for a proton with an energy of 1 GeV. This implies that the parameters $r_{\rm w}$, $r_{\rm sp}$ are larger at higher energy. They eventually become larger than L, so that at high energy escape dominates. At low energy, the relative importance of spallation and convection can be evaluated by comparing $r_{\rm w}$ and  $r_{\rm sp}$. However, it must be noticed that even when $r_{\rm w}$ is greater than $r_{\rm sp}$, the Galactic wind may have a non negligible effect on the cosmic ray spatial origin because the cutoff due to $r_{\rm w}$ is much sharper (see Fig. 6). Moreover, the influence of the Galactic wind on the spectra is important because of the induced energy changes (adiabatic losses).

6 Realistic source distribution

For the sake of definiteness, we will consider from now on that the cosmic ray sources for stable primaries are located in the disk and that their radial distribution $w(r_{\rm s})$ follows that of the pulsars and supernovae remnants, given by

 \begin{displaymath}w_{\rm SN}(r_{\rm s})=\left(\frac{r_{\rm s}}{R_\odot}\right)^...
...\left(-\beta \times\frac{(r_{\rm s}-R_\odot)}{R_\odot}\right),
\end{displaymath} (19)

with $R_\odot = 8.5$ kpc, $\alpha = 2$, $\beta = 3.53$ for Case & Bhattacharya (1998). This distribution is now closer to the distribution adopted by Strong & Moskalenko (1998) ( $\alpha = 0.5$ and $\beta = 1$), a flatter distribution designed to reproduced radial $\gamma$-ray observations (see Fig. 16). This distribution can be inserted in Eq. (4), which is then used to compute the f-surfaces. These surfaces are displayed in Fig. 8 for three cases (L=2 kpc, L=5 kpc and L=10 kpc). For large halos, the source distribution acts as a cutoff and greatly limits the contributions from peripheric Galactic sources.

The results are not much affected by taking an angular dependence into account. Considering for example the spiral arms modelling of Vallée (2002), Fig. 8 shows that the extension of the f-surfaces is almost not affected by these small scale structures. In the rest of this paper, the purely radial distribution (19) is assumed.

7 Application to the propagation parameters deduced from the observed B/C ratio

 The previous sections present a complete description of the origin of cosmic rays in a stationary diffusion model (energy losses and gains are discarded). To each process by which a cosmic ray may disappear before it reaches the solar neighborhood is associated a parameter: L (escape through the top and bottom boundaries), $r_{\rm w}$ (convection), $r_{\rm sp}$ (destructive spallation). The relative importance of these parameters may be measured by the two quantities $\chi_{\rm w} \equiv L/r_{\rm w}$and $\chi_{\rm sp} \equiv L/r_{\rm sp}$. One can distinguish three regimes which determine the diffusion properties of the system: i) the escape through the boundaries dominates for $\chi_{\rm w}\ll1$ and $\chi_{\rm sp}\ll1$; ii) convection dominates for $\chi_{\rm w}\gg1$ and $\chi_{\rm w}\ga \chi_{\rm sp}$; iii) spallations dominate for $\chi_{\rm sp}\gg1$ and $\chi_{\rm w}\ll\chi_{\rm sp}$.

We now use the sets of diffusion parameters consistent with the B/C data given in Maurin et al. (2002a) (hereafter MTD02) to evaluate realistic values for these quantities.

7.1 Evolution of ${\chi _{w}}$ and ${\chi _{sp}}$ with ${\delta }$

In MTD02, we provide for each configuration $\alpha$(source spectral index), ${\delta }$ (diffusion spectral index) and L(diffusive halo size) the corresponding K0, $V_{\rm c}$ and $V_{\rm a}$(Alfvénic wind responsible for reacceleration) that fit best the ratio B/C. In this study, $V_{\rm a}$ is not very important since it only changes the energy of the particles: a cosmic ray emitted at 1 GeV/nuc and gaining a few hundreds of MeV/nuc during propagation will be detected at a slightly larger energy, for which the results given here will not be very different. This becomes even more true beyond a few GeV. Reacceleration will be ignored throughout this study, as well as energy losses, for the same reason. Moreover, the values of K0, $V_{\rm c}$ and $V_{\rm a}$ do not depend much on $\alpha$ (see Fig. 9 of MTD02), so that  $\chi _{\rm w}$and  $\chi _{\rm sp}$ depend mainly on ${\delta }$ and L. They depend on rigidity, through K(E), as can be seen in Fig. 9 where  $\chi _{\rm w}$ and  $\chi _{\rm sp}$are displayed as a function of ${\delta }$ for several species, several values of L and several rigidities.
\par\includegraphics[width=8.8cm,clip]{f8.eps}\end{figure} Figure 8: 99%-surfaces for R =20 kpc and three cases, L =2 kpc, L  = 5 kpc and L=10 kpc.
Open with DEXTER

The left panel displays $\chi_{\rm w}(\delta,L)$ for three rigidities: 1 GV, 10 GV and 100 GV. Up to several tens of GV, convection is in competition with escape; afterwards escape dominates. The noticeable fact is that models corresponding to $\delta\la 0.45$are escape-dominated, whereas convection dominates only for large ${\delta }$ at low energy. It appears that all other parameters being constant, $\chi_{\rm
wind}$ is fairly independent of L(indicating a similar relative importance of convection and escape for the models reproducing the B/C ratio, see MTD02). However, the spatial origin does depend on L and $r_{\rm w}$ and not only on their ratio.

The right panel of Fig. 9 plots $\chi_{\rm sp}(\delta,L)$for ${\cal R}=100$ GV and 1 GV for various nuclei. Protons are the most abundant species in cosmic rays. Boron and CNO family are important because they allow to constrain the value of the propagation parameters, e.g. through the B/C ratio. Last, the Fe group provides another test of the secondary production via the sub-Fe/Fe ratio. The evolution of $\chi _{\rm sp}$ for these species is conform to what is very well known from earlier leaky box inspired studies: for heavier nuclei, spallation dominates over escape and for this reason, the induced secondary production is particularly sensitive to the low end of the grammage distribution.

To summarize, the left panel shows the evolution from convection-domination to escape-domination as a function of ${\cal R}$ and ${\delta }$, the effect of the wind being negligible above $\sim $100 GeV whatever ${\delta }$ ( $\chi_{\rm w}\ga 10$). The right panel gives the evolution from spallation-domination to escape-domination as a function of ${\cal R}$${\delta }$ and the species under consideration. The effect of spallation is more important for heavy than for light nuclei, but this difference is too small to produce an evolution of the average logarithmic mass for high energy ($\sim $TeV) cosmic rays (Maurin et al. 2003a).

7.2 Spatial origin in realistic diffusion models at 1 GeV/nuc

From the previous discussion, it appears that spallations and Galactic wind play a role at low energy. The results will be shown for the particular value 1 GeV/nuc which is interesting for various astrophysical problems. First, once modulated, it corresponds to about the very lowest energy at which experimental set-ups have measured Galactic cosmic rays. Second, the low energy domain is the most favorable window to observe $\bar{p}$(resp. $\bar{d}$) from exotic sources (see companion paper Maurin & Taillet 2003), as the background corresponding to secondaries $\bar{p}$ (resp. $\bar{d}$) is reduced. Last, these energies correspond to that of the enduring problem of the diffuse GeV $\gamma$-ray radial distribution. This was first quoted by Stecker & Jones (1977) and further investigated by Jones (1979) taking into account the effect of a Galactic wind.

\par\includegraphics[width=8.8cm,clip]{f9.eps} \end{figure} Figure 9: Left panel: $\chi _{\rm w}(\delta ,{\cal R})$ as a function of the diffusion spectral index ${\delta }$ for different rigidities ${\cal R}$; from top to bottom, ${\cal R}=100$ GV, ${\cal R}=10$ GV and ${\cal R}=1$ GV. The parameter $\chi _{\rm w}$, as well as $\chi _{\rm sp}$, is not very sensitive to the halo size L. Right panel: $\chi _{\rm sp}(\delta ,{\cal R})$ as a function of ${\delta }$ for ${\cal R}=100$ GV (upper curves) and ${\cal R}=1$ GV (lower curves) for four species: p ( $\sigma \sim 40$ mb), $\bar{d}$ ( $\sigma \sim 100$ mb), B-CNO ( $\sigma \sim 250$ mb) and Fe ( $\sigma \sim 700$ mb). For the latter species we plotted the same three L values as in left panel. The behavior for other species is similar so that we only plotted the case L=6 kpc.
Open with DEXTER

From the sets of diffusion parameters that fit the B/C ratio, the values of the parameters $r_{\rm w}$ and $r_{\rm sp}$are computed (see Table 5) for the four nuclei shown in Fig. 9 and for three values of $\delta=[0.35,0.6,0.85]$.

Table 5: Values of $r_{\rm w}$ and $r_{\rm sp}$ for the sets of parameters that, for a given ${\delta }$, give the best fit to the observed B/C ratio. The mean square value  $\sqrt{\langle r^2 \rangle}$ of the distance to the sources is also shown.

From these values, the 50-90-99%-surfaces are derived and displayed in Fig. 10, for protons and Fe nuclei. The effect of ${\delta }$ (Fig. 12), of L(Fig. 13) and of the species (Fig. 11), are considered separately. As regards the first two effects, the f-surfaces are smaller: (i) for greater values of ${\delta }$, mainly because the effect of the wind is then greater, and (ii) for small values of L, as in this case escape is more important.

\par\includegraphics[width=4.3cm,clip]{f10a.eps}\includegraphics[width=4.3cm,clip]{f10b.eps} \end{figure} Figure 10: (50-90-99)%-surfaces (protons and Fe nuclei are considered), in a typical diffusion model with L=6 kpc and $\delta =0.6$.
Open with DEXTER

As regards the last effect, it can first be seen from Fig. 13 that the heavier species come from a shorter distance (because the spallations are more important). The secondary species can be treated simply by using a source function obtained by multiplying the primary density by the gas density. It would be straightforward to apply the previous techniques to a realistic gas distribution (taking into account, in addition to the fairly flat HI distribution, that of molecular H2 and ionized HII which are more strongly peaked in the inner parts, see e.g. Strong & Moskalenko 1998 for a summary and references) and to infer the contours inside which the secondaries are created. The corresponding f-surfaces are not shown here, as they would be quite similar to those of the primaries (see left panel). What we do display in the right panel are the f-surfaces of the primaries that lead to given secondaries, as these progenitors determine the secondary-to-primary ratios (see Sect. 4.1).

7.3 Effective number of sources

From the previous results, it appears that only a fraction of the sources present in the disk actually contribute to the flux in the solar neighborhood. In this paragraph we present the fraction $f_{\rm s}$ of the sources which are located inside given f-surfaces.

\par\includegraphics[width=4.3cm,clip]{f11a.eps}\includegraphics[width=4.3cm,clip]{f11b.eps}\end{figure} Figure 11: 99%-surfaces for several species. The left panel corresponds to primary species (protons, CNO and Fe) while the right panel corresponds to the progenitors of secondary species (B and sub-Fe), for L=6 kpc and $\delta =0.6$.
Open with DEXTER

\par\includegraphics[width=4.3cm,clip]{f12a.eps}\includegraphics[width=4.3cm,clip]{f12b.eps}\end{figure} Figure 12: 99%-surfaces for several ${\delta }$, in the case L=6 kpc. The left panel corresponds to protons while the right panel corresponds to Fe nuclei.
Open with DEXTER

\par\includegraphics[width=4.3cm,clip]{f13a.eps}\includegraphics[width=4.3cm,clip]{f13b.eps}\end{figure} Figure 13: 99%-surfaces for several L, in the case $\delta =0.6$. The left panel corresponds to protons while the right panel corresponds Fe nuclei.
Open with DEXTER

\par\includegraphics[width=8.8cm,clip]{f14.eps} \end{figure} Figure 14: Fraction $f_{\rm s}$ (in %) of the Galactic sources contributing to the fraction f of the cosmic ray flux at the solar position, for protons and Fe nuclei, for the particular diffusion model L=6 kpc, $\delta =0.6$.
Open with DEXTER

This fraction is presented in Fig. 14 for the particular model L=6 kpc, $\delta =0.6$, for protons and Fe nuclei. It can be read that for example, it takes 7.6% (resp. 1.5%) of the Galactic sources to make 90% of the protons (resp. Fe nuclei) reaching the solar neighborhood.

Figure 14 also shows that a very small fraction of sources may contribute to a non negligible fraction of the fluxes. For example, the sphere of radius $r\sim 100$ pc centered on the solar neighborhood contains only $2.5 \times 10^{-5}$ of the sources but for L=6 kpc and $\delta =0.6$, it is responsible for about 5% of the proton flux and 18% of the Fe flux. The mean age of the cosmic rays is given by $\langle t \rangle \sim \langle
r^2\rangle/2K \sim$ 7-400 Myr (see Table 7). For a supernova rate of three per century, the total number of sources responsible for the flux is ${\sim} 2 \times 10^5 {-} 10^7$. This tells us that in models with the largest ${\delta }$, 18% of the Fe flux can be due to only 5 sources. The approximation of stationarity and of continuous source distribution is likely to break down with such a small number of sources. Conversely, for small values of ${\delta }$(preferred by many authors), this approximation is probably better.

Table 6: Fraction $f_{\rm s}$ (in %) of the Galactic sources contributing to a given fraction (90% and 99%) of the cosmic ray flux at the solar position, for protons and Fe nuclei, for the diffusion models studied before.

Table 7: Mean cosmic ray age $\langle t \rangle \sim \langle r^2
\rangle/2K$ in Myr for the models studied in this paper.

7.4 Radioactive species, e+ and e- and the local bubble

Donato et al. (2002) emphasized that the existence of a local underdensity ( $n\la 0.005$ cm-3) around the solar neighborhood greatly affects the interpretation of the flux of radioactive species at low energy (we refer the interested reader to this paper for a deeper discussion and references on the local interstellar medium). The most important effect of this hole is that it exponentially decreases the flux by a factor $\exp(-a/l_{\rm rad})$ ( $a\la
65{-}250$ pc is the radius of the local underdense bubble and $l_{\rm rad}$ is given by Eq. (12)). This can be easily understood as there is almost no gas in this region, hence no spallations, leading to no secondary production. The local bubble is obviously not spherical, but this approximation is sufficient at this level. This attenuation factor is straightforwardly recovered starting from the probability density as given in Sect. 4.2, if correctly normalized to unity. To this end, the sources (here spallation of primaries on the interstellar medium) are considered to be uniformly distributed in the disk, except in the empty region r<a. The probability density is zero in the hole whereas outside, it is given by

\begin{displaymath}{\rm d}{\cal P}^{\rm hole}_{\rm rad} =
...}{2\pi ~r_{\rm s} .~ l_{\rm rad}}
~ {\rm d}^2\vec{r}_{\rm s}.
\end{displaymath} (20)

The quantity ${\cal P}^{\rm hole}_{\rm rad} (r_{\rm s} < r_{\rm lim}\vert O )$is obtained directly from the no hole case (see Sect. 4.2) by replacing $r_{\rm lim}$ by $r^a_{\rm lim}=r_{\rm lim}+a$. It means that the sources that contribute to the fraction f=(50-90-99)% of the flux of the radioactive species are located between a and $r^a_{\rm lim}=(0.7{-}2.3{-}4.6)\times l_{\rm rad}+a$. Hence, the hole only plays a marginal role for the origin of a radioactive species (unless $l_{\rm rad}\la a$), whereas the result for the flux is dramatically different.

We saw in a previous section that the high energy e+ and e-behave like unstable species. Their typical length $r_{\rm loss}$ can be compared to $l_{\rm rad}$

\frac{l_{\rm rad} \...
... \; \mbox{kpc}^2\; \mbox{Myr}^{-1}}}\cdot
\end{array} \right.
\end{displaymath} (21)

The dependence in the propagation model is similar for both expressions and is contained in the last term. There is a big difference, though, as the typical distances travelled by radioactive nuclei scale as $\sqrt{\cal R}$, whereas they scale as $1/\sqrt{E}$ for electrons and positrons.

Realistic values for $l_{\rm rad}$ and $r_{\rm loss}$ are presented in Fig. 15. At high energy, the Lorentz factor enhances the lifetime of radioactive nuclei, making their origin less local, whereas the energy losses are increased for electrons and positrons, making their origin more local (99-90-50% of 100 GeV $\rm e^+e^-$come from sources located in a thin disk with radius $r^{\rm e^+e^-}_{\rm lim}
\approx 1.1{-}0.55{-}0.38$ kpc).

\par\includegraphics[width=8.8cm,clip]{f15.eps} \end{figure} Figure 15: Realistic values of $l_{\rm rad}/\sqrt{\tau_0/{\rm
1~Myr}}$ and $r_{\rm loss}$ for two extreme halo sizes L and diffusion slope ${\delta }$. As all results in this section, propagation parameters fit B/C and are taken from MTD02.
Open with DEXTER

For 100 GeV e+ and e-, all models fitting B/C have the same value for $K_0{ \cal R}^\delta$, because at 100 GeV/nuc, spallations and convection are negligible (Maurin et al. 2002a). As a consequence, for GeV energies, $r_{\rm loss}$ increases more rapidly for small ${\delta }$ than it does for larger ${\delta }$. To study the effect of local contributions to the spectra of e+ and e-, Aharonian et al. (1995) used the value $\delta =0.6$ and compared to other works with $\delta=0$. As these authors noticed, the modelling in the whole energy spectrum is ${\delta }$ dependent, but our study gives the range compatible with B/C studies.

Finally, radioactive nuclei are a very important tool for cosmic ray physics. They come from a few hundreds of parsec, and their fluxes are very sensitive to the presence of a local underdense bubble, through the attenuation factor $\kappa \equiv \exp(-a/l_{\rm rad})$. For example, for a typical bubble of size a=100 pc and an energy 800 MeV/nuc (interstellar energy), $\kappa_{0.35}\approx\exp(-0.33 \sqrt{1 \; \mbox{Myr}/\tau_0})$ if $\delta =0.35$, whereas $\kappa_{0.85}\approx\exp(-\sqrt{1 \; \mbox{Myr}/\tau_0})$. With $\tau_0^{\rm 10Be}=2.17$ Myr, $\tau_0^{\rm 26Al}=1.31$ Myr and $\tau_0^{\rm 36Cl}=0.443$ Myr, it leads to $\kappa_{0.35}^{\rm Be,~Al,~Cl}\approx0.80{-}0.75{-}0.61$ and $\kappa_{0.85}^{\rm Be,~Al,~Cl}\approx0.51{-}0.42{-}0.22$. For 14C, the attenuation is $\kappa^{14\rm C}\ll 1$ around 1 GeV/nuc, so that this species is heavily suppressed. However, it should be present around 10-100 GeV/nuc (as $\kappa^{14\rm C}\sim 1$ at these energies), with the advantage that solar modulation is less important at these energies.

The flux of radioactive species directly characterizes the local diffusion coefficient K0 if the local environment is specified. This would in turn allow to break the degeneracy in propagation parameters that one can not avoid at present. Last, even though the surviving fraction of a radioactive does depend on the halo size L, we emphasize that it is a very indirect way to derive the propagation parameters. In the forthcoming years, new measurements of radioactive species that do not depend on L (e.g. by PAMELA and AMS) should provide a promising path to update our vision of cosmic ray propagation.

8 Summary, conclusions and perspectives

The question of the source distribution is very present in cosmic ray physics. With the occurrence of the old problem of short pathlengths distribution in leaky box models (see for example Webber et al. 1998), Lezniak & Webber (1979) studied a diffusion model with no-near source in the solar neighborhood. Later, Webber (1993a,b) propagated ${\delta }$-like sources with diffusion generated by a Monte Carlo random walk for the same purpose. Brunetti & Codino (2000) follow this line but they introduce random walks in a more realistic environment, i.e. circular, elliptical and spiral magnetic field configurations. In a more formal context, Lee (1979) used a statistical treatment of means and fluctuations (see references therein) to characterize amongst others the possibility that nearby recent sources may dominate the flux of primaries. Finally, it is known that the present cosmic ray models are not able to reproduce accurately for example proton-induced $\gamma$-rays measurements. To illustrate this point, we plot in Fig. 16 the radial distribution of protons obtained with the same diffusion parameters as used above. None of the models shown match the data.

\par\includegraphics[width=8.8cm,clip]{f16.eps} \end{figure} Figure 16: Radial distribution of the proton flux for the models discussed in this study, compared to the source radial distribution of Case & Bhattacharya (1998) given Eq. (19). For each of the values L=2 kpc and L=10 kpc, the three values $\delta =0.35$, 0.6 and 0.85 are presented, the flatter distribution corresponding to the lower ${\delta }$. Also shown is the gamma-ray emissivity per gas atom ( COS-B Bloemen 1989), which is proportional to the proton flux, as given by COS-B (open circles, Bloemen 1989) and EGRET (triangles, Strong & Mattox 1996), along with the proton flux obtained with the Strong & Moskalenko (1998) distribution (see Sect. 6).
Open with DEXTER

One is left with two alternatives: either modifying the source distribution (for example, the distribution of Strong & Moskalenko 1998 yields a better agreement), or giving up the assumption that the diffusion parameters apply to the Galaxy as a whole (Breitschwerdt et al. 2002). It is thus of importance to understand to what extent the cosmic rays detected on Earth are representative of the distribution of the sources in the whole Galaxy.

We provide an answer to this question under the two important hypotheses that the source distribution is continuous and that we have reached a stationary regime: most of the cosmic rays that reach the solar neighborhood were emitted from sources located in a rather small region of the Galactic disk, centered on our position. The quantitative meaning of "rather small'' depends on the species as well as on the values of the diffusion parameters. For the generic values $\delta =0.6$ and L=6 kpc chosen among the preferred values fitting B/C (see Sect. 7), half of the protons come from sources nearer than 2 kpc, while half of the Fe nuclei come from sources nearer than 500 pc. Another way to present this result is to say that the fraction of the whole Galactic source distribution that actually contributes to the solar neighborhood cosmic ray flux can be rather small. For the generic model just considered, 8% (resp. 1.5%) of the sources are required to account for 90% of the proton (resp. Fe) flux. These fractions are smaller for higher ${\delta }$ and smaller L. To summarize, the observed cosmic ray primary composition may be dominated by sources within a few kpc, so that a particular care should be taken to model these source, spatially as well as temporally (Maurin et al. 2003b).

Independently of all the results, this study could be used as a check for more sophisticated Monte Carlo simulations that will certainly be developed in the future to explore inhomogeneous situations. Several other consequences deserve attention. First, the results may point towards the necessity to go beyondthe approximations of both continuity and stationarity. In particular, it could be that only a dynamical model, with an accurate spatio-temporal description of the nearby sources, provides a correct framework to understand the propagation of Galactic cosmic rays. The contribution from nearby sources would be very different in the low energy ($\sim $GeV/nuc) or in the high energy regime ($\sim $PeV) compared to the stationary background. Second, as discussed in Sect. 4.1, the diffusion parameters derived from the observed B/C ratio have only a local validity, and one should be careful before applying them to the whole Galaxy, since the cosmic rays are blind to most of it.

This work has benefited from the support of PICS 1076, CNRS and of the PNC (Programme National de Cosmologie). We warmly thank Eric Pilon for his expertise on asymptotic developments. We also thank the anonymous referee for his pertinent suggestions.


Online Material

Appendix A: General solutions of the diffusion equation in cylindrical geometry

A.1 General solution for our diffusion-convection model

For a primary species, the differential density (in energy) N(r,z) is a solution of the equation (see for example Maurin et al. 2002b and references therein)

$\displaystyle {\cal L}_{\rm diff}N(r,z) +\Gamma_{\rm rad} N(r,z)
+2h\delta(z)\; n_{\rm ISM}.\sigma_{\rm inel}.v \;N(r,z)$     (A.1)
        $\displaystyle = 2h\delta(z)w(r) ,$    (A.2)


\begin{displaymath}{\cal L}_{\rm diff} =V_{c} \frac{\partial}{\partial z}
...ial r}
\left(r\frac{\partial}{\partial r}\right)\right) \cdot
\end{displaymath} (A.3)

The various terms in Eq. (A.1) correspond respectively to (i) a differential operator ${\cal L}_{\rm diff}$ describing convection $V_{\rm c}$ out from the Galactic plane and isotropic diffusion K throughout the confinement volume; (ii) radioactive decay of the unstable nucleus in the whole Galaxy; (iii) destruction in flight $n_{\rm ISM}.\sigma_{\rm inel}.v$ when crossing the gaseous thin disk of constant height 2h and constant density $n_{\rm ISM}$; (iv) a source term in the thin disk. The solution is determined by the boundary conditions of the problem where cosmic rays freely escape, i.e. we demand N(r=R,z)=N(r,|z|=L)=0(L is the half-height of the diffusive halo, R the radial extension of the Galaxy). The solution in the disk (z=0) is given by

 \begin{displaymath}N(r,0)= \sum_{i=0}^{\infty}~ {\frac{{\cal Q}_i}{A_i}} ~
~ J_0 \left(\zeta_i\frac{r}{R} \right)
\end{displaymath} (A.4)


\begin{displaymath}A_i= 2h \Gamma_{\rm inel} + V_{\rm c} +
KS_i ~ \coth \left(\frac{S_iL}{2} \right)


\begin{displaymath}S_i^2=\frac{4 \zeta_i^2}{R^2} + \frac{V_{\rm c}^2}{K^2}
+ 4\frac{\Gamma_{\rm rad}}{K}

and ${\cal Q}_i$ is the Bessel transform of the source distribution (which may depend on the density of another species, in particular for secondary species).

For a primary point source, $w(r) = \delta(r)/(2\pi r)$and we find in the disk (z=0)

 \begin{displaymath}N_\delta^{\rm cyl}(r,0)=\sum_{i=1}^{\infty}
\frac{J_0\left(\zeta_i\frac{r}{R}\right)}{\pi J_1^2(\zeta_i)R^2A_i} \cdot
\end{displaymath} (A.5)

The generic solution for secondaries can be straightforwardly derived from that of primaries (e.g. Maurin et al. 2001),

N_\delta^{\rm cyl} (r,0)=2h\Gamma_{{\rm prim}\rightarrow{\r...
...)}{\pi J_1^2(\zeta_i)R^2
A^{\rm prim}_i A^{\rm sec}_i}
\end{displaymath} (A.6)

We use $\Gamma_{{\rm prim}\rightarrow{\rm sec}}=
n_{\rm ISM}.v.\sigma_{{\rm prim}\rightarrow{\rm sec}}$. The distinction between  $A^{\rm prim}_i$ and $A^{\rm sec}_i$ is necessary since both species have different destruction rates and rigidities.

Appendix B: Numerical evaluation of the point source solution in Bessel basis

In practice, the infinite sums above are truncated to some order $n_{\rm tronc}$, chosen as a compromise between accuracy (good convergence of the series) and computer time. In the case of a point source $\delta(\vec{r})$, the profiles are singular near the source and the convergence of the series appears to be very slow. A few methods are presented to speed up this convergence.

B.1 Softening of the source term

First, the source term may be spread out on a radius a, by replacing the ${\delta }$-function by $w(r) = \theta(a-r)/(\pi a^2)$for which an extra $2 J_1(\zeta_i a/R)/(\zeta_i a/R)$ term appears in the Bessel transform. With a judicious choice of the parameter a, the solution is very close to the original for $r\gg a$, but convergence is much faster due to the extra $1/\zeta_i$ factor.

B.2 Sum representation: Comparison to a known function

Part of the difficulty to evaluate numerically the Bessel expansions comes from the fact that the resulting functions are singular at the source position. If we know a reference function $f^{\rm ref}(r)$ which exhibits the same singularity and for which the Bessel coefficients $\mu_i^{\rm ref}$ are known, it is then judicious to write the density ( $\rho \equiv
r/R$) as

 \begin{displaymath}N^{\rm cyl}_{\delta}(r,0)=
\left[ \sum_{i=1}^{n_{\rm tronc}}...
...{\rm ref}\right)
J_0(\zeta_i \rho) \right] +f^{\rm ref}(\rho)
\end{displaymath} (B.1)

where the singularity is entirely contained in the f term, so that the Bessel expansion has been regularized. The choice of f may be guided by the behaviors observed in Sect. 3. For $L\sim R$, the solutions should be quite similar to the solution in spherical geometry, given by $f^{\rm
ref}(r)=(1-\rho)/(4\pi \rho K R)$. For very small halo size (L< 1 kpc), the effects of the top and bottom boundaries are dominant and a modified function $f^{\rm ref}(r)=(1-\rho)/(4\rho \pi KR) \times \exp \left( -\rho R/L \right)$is more adapted (see Eq. (9)). This method yields a very good and rapid convergence as long as sources are located in the thin disk z=0.

B.3 Integral representation for infinite radius disk

When the disk has an infinite radius, the Bessel sum can be replaced by an integral, and the end result is obtained from the Bessel sum by the substitution $\zeta_i/R \rightarrow k$ and $1/J_1^2(\zeta_i) \approx
\pi \zeta_i/2 \rightarrow k \pi R/2$, so that in the general case - see Eqs. (A.5) and Eq. (A.4) -,

\begin{displaymath}N^L(r, z) ={\rm e}^{-V_{\rm c}z/2K}
\int_0^\infty \frac{k J_...
...z)/2 \right\}
}{\sinh\left\{ S(k)L/2 \right\} }
~ {\rm d}k ,


\begin{displaymath}A(k) = 2h \Gamma_{\rm inel} + V_{\rm c} + KS(k) \coth \left(
\frac{S(k)L}{2} \right)


\begin{displaymath}S^2(k) = \frac{V_{\rm c}^2}{K^2} + \frac{4\Gamma_{\rm rad}}{K} +4 k^2.

The integrals of the form

\begin{displaymath}I[f] \equiv \int_0^\infty J_0(x) ~ f(x) {\rm d}x

where f(x) is a function such that $f(\infty) = 1$ are quite tricky to compute numerically, due to the very slow decrease of the oscillations in the integrand. Two remarks are of great help. First, as $\int_0^\infty J_0(x) {\rm d}x
=1$, we have

\begin{displaymath}I[f] = 1 - \int_0^\infty J_0(x) ~ (1-f(x)) {\rm d}x.

The convergence is faster, as $1-f \rightarrow 0$ when $x
\rightarrow \infty$. Second, using the identity (xJ1)' = xJ0 and integrating by parts, one has

= \int_0^\infty J_1(x) ~
\left( \frac{f(x)}{x} - f'(x)\right) {\rm d}x.

This expression is meaningful only if x f(x) is a bounded function near the origin x=0. Using the identity J0' = -J1 and integrating by parts again,
$\displaystyle I[f] = \left[J_0(x)\left(\frac{f(x)}{x}-f'(x)\right)
...ty J_0(x) ~
\left(f''(x) - \frac{f'(x)}{x} + \frac{f(x)}{x^2} \right)
{\rm d}x.$     (B.2)

The latter expression is meaningful only if x2 f(x) is a bounded function near the origin x=0. These expressions provide several efficient alternatives to evaluate I[f].

Appendix C: Alternative description of spallations: Random walk approach

A cosmic ray crossing the Galactic disk has a probability p to disappear in a nuclear reaction with interstellar matter. This probability is related to the reaction cross section $\sigma$ by

\begin{displaymath}p= \kappa_1 \sigma n_{\rm ISM} h = 6 \times 10^{-5} ~ \kappa_1 ~
\; \mbox{mb}},

where $\kappa_1 \sim 1$ contains the dependence on the incidence angle of the particle with the Galactic plane. The propagation in the z axis is a one-dimensional random walk, so that for a cosmic ray emitted in the disk and reaching again the disk after a number $t^\star$ of elementary random steps, the probability distribution of disk-crossing numbers n is given by (Papoulis 2002)

\begin{displaymath}{\rm d}{\cal P}_{\rm d} (n \vert t^\star) \equiv
{\rm d}{\ca...
...exp \left(
- \frac{n^2}{\kappa_2 t^\star} \right) ~ {\rm d}n.

In this expression, $t^\star$ is the number of steps of the walk $z=\sum_{i=1}^{t^\star} z_i$and $\kappa_2 \sim 1$ depends on its statistical properties (for instance, $\kappa_2 = 2$ for elementary steps $z_i = \pm \lambda$ and $\kappa_2 \approx 1.43$ for zi uniformly distributed in the interval $[-\lambda,\lambda]$). The diffusion coefficient is defined as

\begin{displaymath}K \equiv \frac{\langle z^2 \rangle}{2t}
= \frac{\kappa_3 \la...
...\kappa_3 \lambda^2}{2\tau} = \frac{1}{2} ~ \kappa_3 \tau

where $\kappa_3\equiv \langle z_i^2 \rangle/\lambda^2$ is the variance of the elementary random step (in units of $\lambda$), so that the physical time is related to $t^\star$ by $t = t^\star \times \tau = t^\star \times 2K/\kappa_3 v^2$, where $\lambda$ is the mean free path and v the velocity. We thus finally have,

\begin{displaymath}{\rm d}{\cal P}_{\rm d} (n \vert t) = \frac{4Kn}{\kappa_2 \ka...
- \frac{2Kn^2}{\kappa_2 \kappa_3 v^2 t}\right) ~ {\rm d}n.

We are now able to compute the probability distribution of disk crossings for cosmic rays emitted from a distance r in the disk as

 \begin{displaymath}\frac{{\rm d}{\cal P}_{\rm d} (n \vert r)}{{\rm d}n} = \int_0...
... \vert t)}{{\rm d}n} ~ {\cal P}_{\rm d} (t \vert r) ~{\rm d}t,
\end{displaymath} (C.1)

where the probability that a CR reaching distance r in the disk was emitted at time t is

\begin{displaymath}{\cal P}_{\rm d} (t \vert r) \propto \frac{1}{(Kt)^{3/2}}
\exp \left( - \frac{r^2}{4K t}\right)\cdot

The above integral C.1 can be performed, yielding the final result

 \begin{displaymath}\frac{{\rm d}{\cal P}_{\rm d} (n \vert r)}{{\rm d}n} = \frac{...
...\left( \displaystyle 1 + \frac{r_0^2 n^2}{r^2} \right)^{-3/2},
\end{displaymath} (C.2)

with $r_0^2 \equiv 8 K^2/\kappa_2 \kappa_3 v^2$. The average number of disk crossings is readily obtained:

\begin{displaymath}\langle n \rangle \equiv \int_0^\infty n \times \frac{{\rm d}...
...P}_{\rm d}
(n \vert r)}{{\rm d}n} ~ {\rm d}n = \frac{r}{r_0},

and the associated variance tends to infinity. We can also compute the integrated probability, that more that n0 crossings have occurred, as

\begin{displaymath}{\cal P}_{\rm d} (n>n_0 \vert r) =
\left( \displaystyle 1 + \frac{r_0^2 n_0^2}{r^2}
\right)^{-3/2} \cdot

A particle having crossed n times the disk has the probability $p_n = (1-p)^n \sim \exp(-np)$ of surviving, so that the survival probability at distance r is given by
                        $\displaystyle {\cal P}_{\rm surv} (r)$ = $\displaystyle \int_0^\infty \frac{{\rm d}{\cal P} (n \vert r,
z(t) = 0 )}{{\rm d}n} ~ {\rm e}^{-np} ~ {\rm d}n$  
  = $\displaystyle \int_0^\infty \frac{x ~ {\rm d}x}{\left(1+x^2\right) ^{3/2}} ~
{\rm e}^{-xrp/r_0}.$  

This can be written as

\begin{displaymath}{\cal P}_{\rm surv} (r) =
\int_0^\infty \frac{x ~ {\rm d}x}{\left(1+x^2\right) ^{3/2}} ~
{\rm e}^{-\alpha x r/r_{\rm sp}},

with $r_{\rm sp}$ defined in Eq. (17) and $\alpha \equiv
\kappa_1 \sqrt{\kappa_2 \kappa_3/2}$. The density of cosmic rays in the disk is then given by

 \begin{displaymath}N(r) = \frac{{\cal P}_{\rm surv} (r)}{4\pi K r}
= \frac{1}{4...
...left(1+x^2\right) ^{3/2}} ~
{\rm e}^{-\alpha x r/r_{\rm sp}}.
\end{displaymath} (C.3)

The quantity $\alpha$ seems to be related to the detailed statistical properties of the random walk under consideration, through $\kappa_1$, $\kappa_2$ and $ \kappa_3$. However, direct comparison with the alternative expression (16) obtained above indicates that these expressions are indeed equivalent, with $\alpha = 1$.

Copyright ESO 2003