text:
# Contents

## 1 Introduction

Most of the available information about matter in our Universe, and in our Galaxy in particular, comes indirectly from the collection of the electromagnetic radiation (from meter waves to $`\gamma `$ rays) that was emitted or absorbed by this matter. A completely different information is provided by the cosmic ray nuclei, which constitute a genuine sample of galactic matter. Many different nuclei species are observed, in a wide range of energy, and with different origins. Some of them come unaltered from the sources (they are called primaries), others (secondaries) come from nuclear reactions between the primaries and the interstellar medium, or from the disintegration of unstable species. Moreover, the trajectories of these nuclei from creation to detection are rather erratic, due to the influence of the galactic magnetic field on all the charged particles, and it is generally not possible to follow the direction of an incoming nucleus back to the source. If we were able to understand clearly the processes by which all these nuclei are produced, accelerated and propagated in the Galaxy, the wealth of data available now or in the near future would yield most valuable information about the matter content and magneto-hydrodynamical properties of our Galaxy. In principle, it would even be possible to discover some evidence for new physics (e.g. supersymmetry) or new objects (e.g. primordial black holes or stars made of antimatter) as they can give rise to the emission of charged antinuclei and make an extra contribution to the observed cosmic ray fluxes. This review presents a summary of the work made in this direction by the authors from 1999 to 2002, in a lapth-isn-iap collaboration. As a first step, we tried to reach a quantitative understanding of the propagation of cosmic ray nuclei in the energy range 100 MeV/nuc-100 GeV/nuc. More precisely, we described propagation with a diffusion model, in which the free parameters are adjusted to account for the available data on cosmic rays. This provides the regions of the parameter space allowed by the data. As a second step, we took advantage of this model to investigate several points concerning astrophysics and astroparticle physics.

During this study, various aspects related to the “standard” or to more speculative processes were examined in detail. These may be summarized as follows:

* Standard cosmic rays

+ Diffusion parameters from secondary-to-primary ratio
    + The flux of standard secondary antiprotons and antideuterons
    + Spatial origin
    + Radioactive species and the local bubble
    + Evolution of composition with energy
* Exotic cosmic rays

+ Baryonic Dark Matter
    + Antimatter
    + Supersymmetric particles
    + Primordial Black Holes

The extraction of the diffusion parameters is the central goal since all conclusions follow from their values.

## 2 Propagation models

### 2.1 The context

The study of cosmic ray radiation raises a great number of questions. The most obvious ones are how the charged nuclei are accelerated and how do they propagate in chaotic magnetic fields. These questions are actually still actively studied and debated, but the nature of the sources and propagation media on which cosmic ray physicist focus has slightly changed in the last twenty years. The pioneering work of Parker , who first used the diffusion-convection equation, focused on Solar modulation. The same equation was used in the late seventies  to investigate the acceleration of Galactic Cosmic Rays (hereafter GCR, in opposition to Solar Cosmic Rays, SCR) in sources that were strongly suspected to be supernovæ (SN) remnants. Whereas SCR deal with the Sun magnetic field and plasma in the Solar cavity, GCR have to face the question of transport in the largely unknown galactic magnetic field. These studies reached an even larger scale, since the discovery of ultra high energy cosmic rays (UHECR) that are certainly extragalactic in origin. The exact nature of this radiation is not clearly established but if these are charged nuclei, the relevant medium in which these nuclei propagate is the extra-galactic magnetic field. From the UHECR source point of view, one has to imagine some powerful astrophysical sources where “standard” acceleration occurs, but an astroparticle solution (i.e. heavy meta-stable decaying new particles) could furnish the energetic particles as well. New particles could thus be discovered by the study of cosmic radiation, which is an interesting coming back to the very first concern of this field of research, when the muon and the positron were first observed in cosmic rays.

Notwithstanding the fact that UHECR is one of the driving subjects in the development of the so-called astroparticle field, the study of much less energetic particles, with energies ranging from GeV to PeV, may also put constraints on the existence of new particles. As discussed at length in the following, these low energy particles have a galactic origin, and by their study, many questions concerning the Galaxy may be addressed, such as LiBeB primordial abundances, dark matter content of our Galaxy, or the nature of the sources. It might be possible to reach a consistent picture of either conventional (SNs, Wolf-Rayet stars) or possible more exotic (micro-quasars , anti-globular clusters ) galactic sources. As a consequence, we underline that GCR nuclei in the GeV-PeV energy range provide a quite interesting laboratory for both astrophysicists and particle physicists. From now on, we will only consider cosmic rays in this energy range.

The great amount of data in various energy ranges has led to a quite good understanding of charged nuclei propagation (see e.g. ) along with the induced $`\gamma `$-ray production at low and high latitude in the galactic plane. As the gyration radius of charged particles in the galactic magnetic field is small, the propagation is intimately related to the detailed structure of this turbulent magnetic field. The latter is not observed directly and one would like to use the cosmic rays to infer its properties. It turns out to be a difficult task as, despite the success outlined above, several unknowns and inconsistencies remain at the quantitative level. Unlike the Solar case for which we have in situ observations of this turbulence – as early as in the mid-sixties (see e.g. ) – favoring a Kolmogorov spectrum, some recent MHD simulations  along with our secondary/primary studies point towards greater spectral index of turbulence. Even though a satisfying global picture emerges, some nuclei resist to a simple interpretation. If the first enduring problem – the depletion of the grammage distribution at low values for sub-Fe/Fe  – seems now to be solved thanks to new cross section measurements , another lasting discussion is related to acceleration and selection mechanisms that lead to the observed abundances. Is it a chemical selection, i.e. First Ionization Potential bias (see e.g. ), or a volatility bias related to grain destruction in SN explosions ? The accuracy of current data is not sufficient to conclude and the question remains in suspense, albeit the importance to elucidate the injection mechanism in acceleration models. It is at least known to a great certainty from radioactive primaries (57Co, 59Ni) synthesized in SN that $`10^5`$ years have past between synthesis and acceleration . The energy spectrum produced by the sources has also some indeterminacy. The acceleration models agree about a power-law dependence in rigidity $`^\alpha `$ (where the rigidity is defined as momentum per unit charge $`p/Z`$) but numerical estimates of the spectral index $`\alpha `$ can be rather different. Values $`\alpha 2.0`$ are preferred by acceleration theory (see  and in particular  for a short and readable introduction on diffusive shock acceleration), but if Kolmogorov spectrum for diffusion is retained, one is left with $`\alpha 2.5`$ from observed spectra that seems more problematic. Spectra may also differ from pure power laws. This depends also on the kind of sources involved (SN, explosion in wind bubbles, superbubbles) along with their plasma and magnetic states. There could be less explosive sources such a Wolf-Rayet stars which eject species through a powerful wind (this is the 22Ne abundance anomaly, see e.g. ).

This list is far from being exhaustive, even if we restrict ourselves to the case of $``$GeV/nuc charged nuclei. To conclude on what we call “standard” CR studies, we emphasize that the limiting factor is related to our poor knowledge of most production cross sections.

### 2.2 Physical motivations for diffusive propagation

The history of cosmic ray propagation theory traces back to the pioneering work of Fermi, which provides a statistical description of the way random magnetic irregularities scatter the particles and accelerate them. This leads to diffusion in energy. Earlier in the same decade, Chandrasekhar  had shown rigorously that diffusion could be equivalently described in terms of random walk, i.e. sequence of small erratic steps. By then, it was possible to imagine that the scattering mentioned above could lead to spatial diffusion. Actually, this hypothesis was first proposed on a phenomenological basis, as emphasized in Berezinskii et al. , in the context of Solar particles transport. It was then applied to cosmic ray transport in the Galaxy. The notion of diffusive motion was confirmed and refined by more fundamental approaches, based on relativistic Boltzmann equations . It now appears as a most valid description, that accounts for a wide range of observations (both in acceleration and propagation mechanisms). The reader is referred to  for a historical background about transport equations (see also  for a derivation from the kinetic equation level and  for the diffusion/convection equation).

The linearized kinetic theory approach provides grounds for a consistent derivation of the transport equation, which reads, neglecting spallations and energy losses for the sake of clarity

$`{\displaystyle \frac{f}{t}}\stackrel{}{}(K\stackrel{}{}f\stackrel{}{V_c}f){\displaystyle \frac{\stackrel{}{}.\stackrel{}{V_c}}{3}}{\displaystyle \frac{1}{p^2}}{\displaystyle \frac{}{p}}(p^3f)`$ (1)
$`={\displaystyle \frac{1}{p^2}}{\displaystyle \frac{}{p}}p^2K_{pp}{\displaystyle \frac{}{p}}f+{\displaystyle \frac{dQ}{dp}},`$
where $`ff(t,\stackrel{}{r},\stackrel{}{p})`$ is the phase space distribution. This equation contains most of the effects described below, like spatial diffusion, galactic wind, with the adiabatic energy loss associated, and diffusion in momentum space. In this equation, spatial diffusion has been assumed to be isotropic. Actually, the diffusion coefficient $`K`$ should be replaced by a tensor, with parallel and transverse components. As regards the first one, there is a strong consensus about a form

$$K_{}()=K_0\beta ^{2\kappa },$$
(2)
where $`\kappa `$ is the spectral index of the turbulence spectrum. The transverse component is still debated, and the two main propositions (given by quasi-linear or the Bohm conjecture) are probably wrong . All the results presented here are based on the usual assumption that diffusion is isotropic, with a diffusion coefficient

$$K()=K_0\beta ^\delta ,$$
(3)
where the normalization $`K_0`$ and the spectral index $`\delta `$ should ideally be related to the astrophysical properties of the interstellar medium. Unfortunately, our knowledge in this field is still demanding, and the value of the two parameters $`K_0`$ and $`\delta `$ can only be determined indirectly by the analysis of cosmic ray observations.

The different diffusion schemes mentioned above also lead to different forms for the energy dependence of the reacceleration term $`K_{pp}`$ . The reader is referred to , and we will not discuss this point further.

### 2.3 Overview of the effects affecting propagation

This section is devoted to a brief overview of the physical effects that play a role in the propagation of cosmic rays. The next section (see 2.4) will enter into more details and focus on the modelling of these effects.

#### 2.3.1 Geometry and content of the Galaxy

The propagation of cosmic rays ceases to be of diffusive nature beyond some surface where they can freely stream out of the diffusive volume. The density then drops to nearly zero, so that this surface may be considered as an absorbing boundary. The exact shape and dimensions of this boundary are not known, but direct observations of the radio halo of external galaxies suggest that it might radially follow the galactic disc, with a greater thickness. In this work, the diffusive halo will be modelled as a cylinder of radius $`R=20`$ kpc and half-height $`L`$ whose numerical value, still to be determined, is probably greater than a few kpc. The boundary thus imposes that the density satisfies $`N(r=R,z)=N(r,z=\pm L)=0`$. Embedded in this diffusive halo lies the galactic disc (hereafter “the disc”, in short) containing the stars and the gas . The gas is mostly made of hydrogen (90%), neutral and ionized, and helium (10%) (the heavier nuclei that may be present are of negligible importance for our concerns). The different components (stars and gas) have different half heights $`h_i`$, of the order of $`h100\text{pc}`$; they all satisfy $`hL`$, so that the disc will be considered as infinitely thin for all practical purposes. The density of interstellar matter is observed to be about $`n_{\mathrm{ISM}}1\text{part}\text{cm}^3`$ for all radii, so that we take $`n(r,z)=2h\delta (z)n_{\mathrm{ISM}}`$. Sources and interactions with matter are confined to the thin disc and diffusion which occurs throughout disc and halo with the same strength, is independent of space coordinates. The Solar System is located in the galactic disc ($`z=0`$) and at a Galactocentric distance $`R_{}=8\text{kpc}`$ . A schematic view of the galactic model is shown in Fig. 2.

#### 2.3.2 Spallations: the importance of cross sections

When a cosmic ray crosses the disc, it may interact with an interstellar hydrogen or helium nucleus and initiate a nuclear reaction (spallation). The importance of this effect is governed by the corresponding cross sections. Actually, it is important to know not only the reaction or total cross section, which determines the rate of destruction of a given CR species, but also the branching ratio into a particular channel, which gives the formation rate of new nuclei. The cross section for a given channel is often referred to as the spallation or fragmentation cross section. The determination of all these cross sections is a nuclear physics problem that shall be addressed in details in Sec. 2.4.5.

#### 2.3.3 Energy losses from interaction with the ISM

There are two types of energy losses which are relevant for nuclei: ionization losses in the ISM neutral matter ($`90`$% H and $`10`$% He), and Coulomb energy losses in a completely ionized plasma, dominated by scattering off the thermal electrons. The other effects like bremsstrahlung, synchrotron radiation and inverse Compton are negligible in the conditions considered here.

The Coulomb energy loss rate is given in

$$\left(\frac{dE}{dt}\right)_{\mathrm{Coul}}4\pi r_e^2cm_ec^2Z^2n_e\mathrm{ln}\mathrm{\Lambda }\frac{\beta ^2}{x_m^3+\beta ^3},$$

where

$$x_m(3\sqrt{\pi }/4)^{1/3}\sqrt{2kT_e/m_ec^2};$$

$$\mathrm{ln}\mathrm{\Lambda }\frac{1}{2}\mathrm{ln}(\frac{m_e^2c^4}{\pi r_e\mathrm{}^2c^2n_e}.\frac{M\gamma ^2\beta ^4}{M+2\gamma m_e}).$$

In these expressions $`r_e`$ and $`m_e`$ denote the classical radius and rest mass of the electron (Particle Data Group 111http://pdg.lbl.gov/), $`n_e0.033`$ cm-3 and $`T_e10^4`$ K denote the density and temperature of the interstellar electrons , $`Z`$ and $`M`$ are the charge and mass numbers of the incoming nucleus and $`\mathrm{ln}\mathrm{\Lambda }4050`$ is the Coulomb logarithm.

The relativistic expression giving the ionization losses is

$$\left(\frac{dE}{dt}\right)_{\mathrm{Ion}}(\beta \beta _0)\frac{2\pi r_e^2m_ec^3Z^2}{\beta }\underset{\mathrm{s}=\mathrm{H},\mathrm{He}}{}n_sB_s$$

where

$$B_s\mathrm{ln}\left(\frac{2m_ec^2\beta ^2\gamma ^2Q_{\mathrm{max}}}{I_s^2}\right)2\beta ^2;$$

$$Q_{\mathrm{max}}\frac{2m_ec^2\beta ^2\gamma ^2}{1+[2\gamma m_e/M]};$$

and $`\beta _0c0.01c`$ is the typical velocity of bound electrons in the hydrogen atom, $`I_s`$ is the geometrical mean of all the ionization and excitation potentials of the considered atom, ($`I_\mathrm{H}=19\text{eV}`$ and $`I_{\mathrm{He}}=44\text{eV}`$), $`Mm_e`$ is the incident nucleon mass, and $`n_s`$ is the density of the target atom in the ISM.

#### 2.3.4 Adiabatic losses from convective wind

Among the phenomena affecting the propagation of cosmic rays, the magnitude of those presented in this section and the following are more subject to debate because of the uncertainty associated to parameters $`V_c`$ and $`V_a`$.

It is very likely that the medium responsible for diffusion is moving away from the disc, with a velocity $`V_c`$. This is referred to as convective or galactic wind, in analogy with the Solar wind. One of the effects of this galactic wind is to dilute the energy of the particle located in the disc in a larger volume . This adiabatic expansion results in a third type of energy loss, depending on $`.V_c`$. Throughout our works, a very simple and tractable form for $`V_c`$ is adopted, following . It is considered to be perpendicular to the disc plane and to have a constant magnitude throughout the diffusive volume, so that $`dV_c/dz=0`$ except at $`z=0`$ where a discontinuity occurs, due to the opposite sign of the wind velocity above and below the galactic plane. In this case, the dilution effect on the energy is given by a term that can be expressed in the same form as ionization and coulomb losses

$$\left(\frac{dE}{dt}\right)_{\mathrm{Adiab}}=E_k\left(\frac{2m+E_k}{m+E_k}\right)\frac{V_c}{3h}.$$
(4)
$`E_k`$ stands for the total kinetic energy and it should not be confused with the kinetic energy per nucleus frequently used in this paper. We emphasize that this term corresponds to a process occurring only in the disc but not in the halo.

#### 2.3.5 Reacceleration

Along with a the spatial diffusion coefficient $`K`$, Eq. (1) also contains the momentum diffusion coefficient $`K_{pp}`$; both are related to the diffusive nature of the process. The latter coefficient $`K_{pp}`$ is related to the velocity of disturbances in the hydrodynamical plasma, called Alfvén velocity. From the quasi-linear theory, it is given by (see e.g. )

$$K_{pp}=\frac{h_{\mathrm{reac}}}{h}\times \frac{4}{3\delta (4\delta ^2)(4\delta )}V_{a}^{}{}_{}{}^{2}p^2/K(E).$$
(5)
In this expression, $`h_{\mathrm{reac}}`$ stands for the half-height of the cylinder in which reacceleration occurs. In our model, $`h_{\mathrm{reac}}h`$, but as $`K_{pp}`$ depends only on the combination $`V_a^2h_{\mathrm{reac}}/h`$, the same diffusion in momentum space is obtained for $`h_{\mathrm{reac}}h`$, provided that the true Alfvén velocity is given as a function of the parameter $`V_a`$ by $`V_a\times h/h_{\mathrm{reac}}`$.

#### 2.3.6 Full propagation equation

The transport equation (1) can be rewritten for each species $`j`$ using the cosmic ray differential density $`N^j(E)dn^j/dE`$. As the momentum distribution function is normalized to the total cosmic ray number density ($`n=4\pi 𝑑pp^2f`$), we have $`N^j(E)=(4\pi /\beta )p^2f^j`$ to finally obtain, assuming steady-state (see Sec. 2.4.6),

$`\stackrel{}{}\left[K\stackrel{}{}N^j(E)\stackrel{}{V_c}N^j(E)\right]\mathrm{\Gamma }^jN^j`$ (6)
$`{\displaystyle \frac{(\stackrel{}{}.\stackrel{}{V_c})}{3}}{\displaystyle \frac{}{E}}\left[{\displaystyle \frac{p^2}{E}}N^j(E)\right]=𝒬^j(E)+`$
$`{\displaystyle \frac{}{E}}\left[b_{\mathrm{tot}}(E)N^j(E)+\beta ^2K_{pp}{\displaystyle \frac{N^j(E)}{E}}\right];`$
where the following notation has been used for the total energy loss term $`b_{\mathrm{tot}}=b_{\mathrm{loss}}+b_{\mathrm{reac}}`$, with

$$b_{\mathrm{loss}}(E)=\left(\frac{dE}{dt}\right)_{\mathrm{Ion}}+\left(\frac{dE}{dt}\right)_{\mathrm{Coul}}+\left(\frac{dE}{dt}\right)_{\mathrm{Adiab}}$$

and the reacceleration drift term defined as

$$b_{\mathrm{reac}}(E)=\frac{(1+\beta ^2)}{E}K_{pp}.$$

We also use a compact notation to describe the most general form for a source term

$$𝒬^j(E)=q_0Q^j(E)+\underset{k}{\overset{m_k>m_j}{}}\mathrm{\Gamma }^{kj}N^k(0),$$
(7)
which includes primary sources – normalized abundance $`q_0`$, spectrum $`Q^j(E)`$ –, but also secondary sources, coming from spallations (see Sec. 2.4.6) or radioactive decay of a heavier species (see Sec. 2.4.8). The relative magnitude of all the effects affecting propagation can be estimated from the typical timescale associated with these effects, as displayed in Fig. 1.

Taking advantage of cylindrical symmetry and adding radioactive contributions localized in the disc and the halo, the previous equation may be rewritten as (making implicit the energy dependence)

$`0=\left(_{\mathrm{diff}}\mathrm{\Gamma }_{\mathrm{rad}}^j\right)N^j(r,z)+{\displaystyle \underset{k}{\overset{m_k>m_j}{}}}\mathrm{\Gamma }_{\mathrm{rad}}^{kj}(E)N^k`$ (8)
$`+2h\delta (z)\left(q_0^jQ(E)q(r)\mathrm{\Gamma }^j(E)N^j(r,0)\right)`$
$`+2h\delta (z){\displaystyle \underset{k}{\overset{m_k>m_j}{}}}\mathrm{\Gamma }^{kj}(E)N^k(r,0)`$
with

$$_{\mathrm{diff}}=V_c\frac{}{z}+K(E)\left(\frac{^2}{z^2}+\frac{1}{r}\frac{}{r}\left(r\frac{}{r}\right)\right).$$

One needs to solve a complete triangular-like set of coupled equations since a given nucleus can only be obtained by spallation of a heavier one and never by a lighter. Quantities in this equation are functions of spatial coordinates (not time, steady-state being assumed) and of kinetic energy per nucleon (energy for short) since this is the appropriate parameter to be used, as it is conserved in spallation reactions.

#### 2.3.7 Solar modulation

The cosmic rays that we detect on Earth had to penetrate the Solar cavity, a process by which they lose energy. This phenomenon is called Solar modulation (see  for a review), and may be pictured as follows. The Sun emits low energy particles in the form of a fully ionized plasma having $`v400`$ km s-1. This so-called Solar wind shields the Solar cavity from penetration of low energy GCR. It was first studied by Parker  who established the evolution of the flux in the Solar cavity. Analog to the propagation in the Galaxy, it is a diffusion equation in a quite different geometry (spherical). For practical purposes, this equation can be solved numerically, or one can use the force-field approximation.

##### Force-field

Perko  provided a useful and compact approximation to the full modulation equation. The final result for a nucleus with charge $`Z`$ and atomic number $`A`$ is a mere shift in the total energy

$$E^{\mathrm{TOA}}/A=E^{\mathrm{IS}}/A\left|Z\right|\varphi /A.$$
(9)
Here $`E^{\mathrm{TOA}}`$ and $`E^{\mathrm{IS}}`$ correspond to the top-of-atmosphere (modulated) and interstellar total energy, respectively. The Solar modulation parameter $`\varphi `$ has the dimension of a rigidity (or an electric potential), and its value varies according to the 11-years Solar cycle, being greater for a period of maximal Solar activity. Another equivalent quantity often used is $`\mathrm{\Phi }=\left|Z\right|\varphi /A\frac{1}{2}\varphi `$. Once the momenta at the Earth $`p^{\mathrm{TOA}}`$ and at the boundaries of the heliosphere $`p^{\mathrm{IS}}`$ are determined, the interstellar flux of the considered nucleus is related to the TOA flux according to the simple rule

$$\frac{\mathrm{\Phi }^{\mathrm{TOA}}\left(E^{\mathrm{TOA}}\right)}{\mathrm{\Phi }^{\mathrm{IS}}\left(E^{\mathrm{IS}}\right)}=\left\{\frac{p^{\mathrm{TOA}}}{p^{\mathrm{IS}}}\right\}^2.$$
(10)
The determination of the modulation parameter $`\mathrm{\Phi }`$ is totally phenomenological and suffers from some uncertainties. As explained in the following, we will deal with data taken around periods of minimal Solar activity, for which we fixed $`\mathrm{\Phi }=250\text{MV}`$ (or equivalently $`\varphi =500\text{MV}`$).

The effect of Solar modulation may be decoupled from the problem of interstellar propagation. Would a more careful treatment of Solar modulation be needed, e.g. , an interstellar flux can easily be obtained from force-field modulated fluxes by demodulation (the force-field approximation we used is reversible). This interstellar flux thus obtained could then be used as an input for any other preferred treatment of Solar modulation. Of course, Solar modulation induces some uncertainty, but this question is still debated, and a rigorous treatment of this effect is beyond the scope of this paper and moreover has not the same impact on antiprotons/protons and on other nuclei. For the first ones, the effect is greater but as one is mostly interested in setting upper limits and drawing exclusion plots on exotic sources, it is not worth using complex modulation schemes. The other nuclei are less sensitive to the exact modulation parameter and to the scheme used. Thus, our results should not be too sensitive to them.

### 2.4 Modelling the effects affecting propagation

We gave above the basic ingredients that enter in propagation models. We now come back to several aspects for a more detailed discussion, the purpose being twofold: first, some ingredients support enduring developments and for these, it is not always an easy task to choose among the variety of parameterizations and approaches available. Thus, taking a sort of step by step walk – leading sometimes to kind of zoological approach –, we will see how problems were tackled in the literature in the past, up to the point they reached today. Secondly, as modelling the galactic environment is a never ending story, because the description is richer and richer as one dives into fine structure subtleties, we find it advantageous to stress the level of detail we wish to include in our model. Finally, some choice will be made for the sake of tractability of the propagation problem.

#### 2.4.1 Matter content, ISM and LISM

The propagation of Galactic Cosmic Rays is influenced by the distribution of matter in the disc and the magnetic structure of the Galaxy. These two ingredients require a careful modelling. As regards the first point, the matter content is quite well known (see  for a review) from in situ observations of gas thanks to multi-wavelength surveys. As regards the second point, it is more difficult to determine because of our particular position in the galactic disc. Studies of other spirals are particularly useful in that respect. For instance, the presence of a diffusive halo or a galactic wind is more easily established in these other galaxies by radio observations.

##### Local interstellar medium (LISM)

As will be discussed in Sec. 6.3 some radioactive species of cosmic rays are sensitive to the very fine structure of the interstellar medium (ISM), in a region of about $``$ 100 pc around the Solar neighborhood, usually referred to as the local interstellar medium (LISM). A general review of our local environment can be found in . The LISM is defined as a region of radius $``$ 65-250 pc (the bubble) containing extremely hot ($`10^510^6`$ K) and low density ($`n`$ 0.005 cm-3) gas, surrounded by a dense neutral gas boundary (hydrogen wall. A finer description of this bubble tells that the Sun is located in a local fluff with $`N_{\mathrm{HI}}`$ 0.1 cm-3, $`T10^4`$ K and a typical extension of $`50`$ pc. To end with, it is of great importance for further modelling to realize that the local bubble is highly asymmetric , and that several cloudlets are present in the bubble (for a schematic representation, see Fig. 3). Various models have attempted to explain the formation of this local bubble , but this subject is far beyond our concern.

Several approximation levels may be considered to model this bubble. The first level is that for most species of cosmic rays it may be ignored, as the presence of the bubble has almost no effect on their propagation. A next level is necessary for some radioactive species which are very sensitive to the local environment. In the thin disc approximation, this bubble may then be represented by a lower density in some disc area surrounding the Sun. In our model, it is a circular hole of radius $`r_{\mathrm{hole}}`$ (to be determined) with a null density. It will be called “the hole” for short. A third level of approximation is the one proposed in , who use a three-layer model.

##### Interstellar medium (ISM)

Apart from the radioactive species, cosmic rays travel for hundreds of kpc before reaching us, so that they are not sensitive to the kind of structures mentioned above but only to the smoothed out properties of the ISM. From this coarse point of view, the ISM is a rather homogeneous mixture of neutral hydrogen, helium, molecular hydrogen, ionized gas and dust (see e.g. ) located in a narrow disc. The total density of matter varies with the distance $`r`$, but for the sake of analytical tractability, we consider a constant density of $`1\text{particle}\text{cm}^3`$. This is not a crucial assumption, as the cosmic rays detected on Earth have only probed a region of a few kpc of extension, so that this effective number can be considered as the mean density in this region.

To end with the ISM question, a last point deserves attention. Gas forms clouds, and stars form with gas, especially where dense clouds are. In , authors have studied the consequence of having most of these sources in dense clouds. Notwithstanding the question of a larger production of secondaries in these clouds, they could have different magnetic fields and then be more or less transparent to cosmic rays. As a matter of fact, we would be left with two more adjustable parameters (filling factor, transparency), and in a certain way, other usual propagation parameters fit to observations would have the meaning of averaged effective quantities, analogously to a uniform distribution of source and gas. Anyway, with our model, we do not intend precisely to provide better than effective parameters, useful and applicable for all species that seem to have the same propagation history. This notion of effectiveness and domain of validity of our model will often come back all the paper long.

#### 2.4.2 Diffusive halo, regular and erratic magnetic field

Many aspects of cosmic rays propagation call for the existence of an extended magnetic halo (see e.g. ). First, clear evidences were obtained from non–thermal radio emission in NGC 4631 . Two components of the galactic magnetic field coexist: a regular one (average value about a few $`\mu `$G, parallel to the galactic plane, responsible for confinement) and a stochastic one which is responsible for charged nuclei diffusion (as well as diffusive reacceleration), that has about the same strength. Tab. 1, taken from , summarizes several observational evidences.

In the last ten years, several studies tried to constrain the halo size using a combination of stable plus radioactive nuclei ratio measurements. They found $`L4`$ kpc , $`L4`$ kpc , $`L412`$ kpc  and $`L24`$ kpc . As all the authors used quite different diffusion parameters, with or without wind, reacceleration, and as they took quite different diffusion slopes $`\delta `$, the trend seems to favor small halos. However, the radial $`\gamma `$-ray distribution in our Galaxy is rather flat, which points towards a large halo (see e.g. ). To face the problem,  have recently proposed a more complex diffusion model, where some quantities (galactic wind) have a spatial dependence. We feel that such an approach becomes unavoidable with the present development of cosmic rays studies. However, we wish to mention two arguments against the claim of inconsistency advanced by these authors to motivate their model: first, we showed in  that considering B/C ratio in diffusion/convection/reacceleration models, all halo sizes are possible. Second, focusing on radioactive species, we showed  that if we agreed with usual results for a homogeneous LISM, the existence of a local underdensity leaves the choice of $`L`$ almost free (too large halo sizes, i.e. $`L>12`$ kpc are excluded). Thus, from our point of view, the discussion about the halo size is far from being closed.

Finally, it must be kept in mind that halos are dynamical objects. In particular, there seems to be a transition zone (thick disc) between the more quiet halo and the active gaseous disc. This is deduced from the clouds, hot ionized gas and stars kinematics. The halo could be substructured in chimneys, galactic fountains and galactic holes . The importance of this substructure for cosmic ray studies is an open question.

#### 2.4.3 Standard sources in the disc

One easily imagines that the location of sources plays an important role for the density of all species, primaries, secondaries and radioactive nuclei. Moreover, it is intimately related to the question of validity of the stationary hypothesis in cosmic ray propagation models. Nearby sources are for example invoked in  to explain the observed behavior at PeV energy, i.e. the knee. If one goes one step further into this line of thought, nothing would guarantee that even at low energy the steady state approximation is justified. It could be that propagation parameters provided by stationary models are just effective parameters, that would need to be reinterpreted in terms of time dependent propagation models. Such a view could seemingly solve some present problems of cosmic rays (Maurin et al., in preparation). Anyway, this paper focuses on stationary models.

##### Radial distribution of sources $`q(r)`$

Measurements of galactic $`\gamma `$ rays in the seventies have raised the question of the radial distribution of cosmic rays. This distribution is needed in order to evaluate the resulting $`\gamma `$ emissivity at different galactocentric locations. The first distribution used was that of Kodaira  following the radial distribution of supernovæ which is also close to that of pulsars. This is consistent with the present picture of cosmic rays where supernovæ provide the energetic budget and mechanism to accelerate nuclei. The description of the galactocentric distribution has been steadily improved  thanks to new observations of pulsars and supernovae. We take here the distribution of Case & Bhattacharya  which is an improvement of their earlier analysis  and happens to be closer to the distribution adopted by Strong & Moskalenko . All these distributions could be compared to flat distribution, i.e. $`q(r)=1`$, that is also widely used in one-dimensional propagation models . These three radial distributions are

* Model a: flat distribution

$$q(r)=1.$$
* Model b: Case & Bhattacharya

$$q(r)=\left(\frac{r}{8.5}\right)^{2.0}\mathrm{exp}\left(3.53\times \frac{(r8.5)}{8.5}\right).$$
* Model c: Strong & Moskalenko

$$q(r)=\left(\frac{r}{8.5}\right)^{0.5}\mathrm{exp}(1.\times \frac{(r8.5)}{8.5}).$$

They are displayed in Fig. 4 along with the metallicity gradient discussed below.

We emphasized in  that the difference between a flat and a more realistic $`q(r)`$ was a mere rescaling of the propagation parameters. However, this difference is not so important if one considers these parameters as effective and uses them in a consistent framework. To be more specific, one should make the same hypothesis on $`q(r)`$ when i) determining the propagation parameters from the observation of some cosmic species and ii) predicting the abundance of some species from these parameters. On the other hand, if one is interested in a physical interpretation of these parameters and wants to deduce some galactic property or to predict some quantity which does not depend only on the diffusion properties (such as proton-induced $`\gamma `$-ray production), one should be careful to consider the right radial distribution.

##### Metallicity gradient

It is also expected that the composition of SN should depend on its position in the Galaxy. Indeed, the inner Galaxy is richer in heavy stellar material than the outer Galaxy; this is the case for most spiral galaxies. Early studies showed a gradient for O/H from observations of ionized nebulae in galaxies like M33, M51, and M101, but later works observed this trend in our Galaxy for many other abundances (see  for a review and Sec. 2 of ). Several recent observations (see Tab. 1 of  for a compilation and  for latest results) lead to very similar conclusions for the metallicity gradient. It can roughly be parameterized, for all ion species X but He, as

$$\frac{d[\mathrm{X}/\mathrm{H}]}{dr}=0.05\text{ dex kpc}^1,$$
(11)
where $`[\mathrm{X}/\mathrm{H}]`$ is defined as $`log_{10}(\mathrm{X}/\mathrm{H})log_{10}(\mathrm{X}/\mathrm{H})_{}`$. Associated with the radial distribution, this gives an additional factor

$$q(r)10^{0.05(r8.)}\times q(r),$$
(12)
except for H and He.

##### Spiral structure

Because of the presence of spiral arms, there is also an angular dependence of the source distribution. As  showed, SNe II and Ib are likely to follow the gas distribution. As also confirmed in , SNe are seemingly correlated with the spiral arms position. The structure can be found in a recent meta-analysis of the Milky Way  where the structure is likely to resemble a 4 logarithmic spiral arms (see its Fig. 2). Notice that, despite their possible usefulness, such models already fall to describe the arms outside the range $`313`$ kpc. Anyway, it is shown in  that this structure actually has very little effect on the origin of sources that contribute to the flux detected at Earth. However, we are left once more with the question of the meaning of the propagation parameters we will derive.

Nevertheless, if a spiral structure is retained, the other ingredients should follow as well (magnetic field structure, gas distribution,…); such a modelling is beyond the scope of this study, though this is probably an important point. Joint progress in several research fields (diffusion mechanisms, better knowledge of galactic structure, etc) associated with numerical developments will certainly allow to cope with these questions in a couple of years. Anyway, we will not discuss further this point in this paper.

##### Source spectrum

The energy spectrum of the particles emitted by the sources is determined by the acceleration process at work. There is a strong belief, based on theoretical work, that this energy spectrum $`dQ/dp`$ is a power-law in rigidity . For a species $`j`$, the differential spectrum in energy $`dQ^j/dEQ^j()`$ is then given by

$$Q^j()=\frac{q_0^j}{\beta }\left(\frac{}{1\text{GV}}\right)^\alpha $$
(13)
in our Eq. (6) (see details in ), where the value of $`\alpha `$ is still debated (grossly from about 1.5 to 2.5 ). We will refer to this case as the pure power law spectrum. Among the various acceleration models, some of them are able to produce various spectra for different nuclei , but also a break in slopes around PeV energies (see references in ).

Because of collective effects, one can also consider a possible deviation from this power-law at low energy, and we choose

$$Q^j()=q_0^j\left(\frac{}{1\text{GV}}\right)^\alpha .$$
(14)
We refer to this as the modified spectrum.

#### 2.4.4 Exotic sources in the halo

Cosmic rays may also be created by the pair annihilation of exotic particles – such as the neutralino predicted by supersymmetry (SUSY) –, as discussed in Sec. 7.3, or by the evaporation of primordial black holes – hereafter PBH – as discussed in Sec. 7.4. Both populations are supposed either to fill (SUSY) or to follow (PBH) the dark matter profile, which can be assumed to be a spherical isothermal profile, i.e. for each species $`\rho _i(r,z)=\rho _i^{}\times f(r,z)`$ with

$$f(r,z)=\frac{R_c^2+R_{}^2}{R_c^2+r^2+z^2}.$$
(15)
For definiteness, the dark matter halo has been assumed to have a core radius $`R_c=3.5`$ kpc. The uncertainties on $`R_c`$ and the consequences of a possible flatness have been shown to be irrelevant at least in the PBH case .

The source density in the Solar neighborhood is about $`\rho _\chi ^{}=0.4`$ GeV cm-3. For PBH the normalization $`\rho _{\mathrm{PBH}}^{}`$ is very different (far less) and can be constrained through CR signatures (see Sec. 7.4). This is at variance with the SUSY case where the SUSY space parameter is constrained instead of source density (see Sec. 7.3).

##### Effective source term

Finally, the production rate, i.e. the effective source term is different in the two cases: it is proportional to $`\rho _{PBH}^{}\times f(r,z)`$ for cosmic rays evaporating from PBH, but to $`(\rho _\chi ^{}\times f(r,z))^2`$ for those coming from neutralino annihilations.

The resulting source profiles are shown in Fig. 5, compared to the typical radial extension of the Galaxy. One easily imagines that sources closer to the galactic center will contribute more to the detected flux, even more in the SUSY case than in the PBH case (see also discussion in Sec. 7.2).

#### 2.4.5 Catastrophic losses: destruction on interstellar matter

The overview in Sec. 2.3.2 emphasized the role of nuclear reactions and decay in the propagation of cosmic rays. The importance of these events is not the same for all species, and one should compare the typical propagation time $`20`$ Myr to the effective time of the process involved. The most important loss is from spallations, i.e. inelastic interaction of the nucleus with the matter content of the thin gaseous disc (as a crude estimation, H : He : CNO : Fe $`10:1:0.1:0.01`$, about the same proportions as for cosmic rays). This statement follows from the mere presence of secondary species: LiBeB and sub-Fe nuclei (Sc, Ti, V) are absent in sources but they are observed with the proportion B/C $``$ sub-Fe/Fe $`0.1`$ in cosmic ray radiations; they were formed by spallations. Of course, the heavier nuclei have larger cross sections and break more easily during their wandering through the Galaxy. This preferential destruction in flight is, as a matter of fact, responsible for part of the nuclear enrichment in heavies observed till PeV energies (, see also Sec. 6.4).

When a cosmic ray nucleus A of energy $`E_A`$ interacts with a nucleus H or He (for simplicity, we write the formulae for H only) at rest in the interstellar medium, a nuclear reactions may occur, producing some daughter species B, C, …of different energies $`E_B`$, $`E_C`$, …Actually, it turns out that the kinetic energy per nucleus $`E`$ is approximately conserved during the reactions, so that all the information we need is contained in the set of differential cross sections $`d\sigma _{A+HB+C+\mathrm{}}/dE`$. From these, the total inelastic cross section for the species $`A`$ may be obtained by summing over the final states $`B`$, $`C`$, …The production cross section for the species $`B`$ can also be found by summing over the initial states $`A`$.

On a purely theoretical point of view, the computation of the $`d\sigma _{A+HB+C+\mathrm{}}/dE`$ is a well-defined but hard to solve nuclear physics problem (see  for a review). A nucleus-nucleus collision is generally an out of equilibrium process that need to be dynamically treated. At a few hundreds of MeV/nuc, the system expands and the equation of state of nuclear matter enters in the spinodal region (region of great mechanical instability): small perturbations are amplified until they reach the size of the system. At GeV energies, one is left with several nuclei (multi-fragmentation); a statistical description can be used (the number of fragments does not depend on energy but only on the nucleus involved). More specific processes, such as single and double nucleon removal are indeed the dominant processes for many nuclei at these energies. Purely hadronic processes arise at a few hundreds of GeV; more and more resonances are excited until quark fragmentation becomes dominant. Among microscopical approaches , one could quote probabilist models (statistical multi-fragmentation), dynamical models (quantum molecular dynamics) and kinetic models (Boltzmann/Vlasov-Uehling-Uhlenbeck equations). However, none of these models enables a quick and precise evaluation of the desired cross section. At most, they can reproduce data with only a reasonable agreement. This approach is certainly promising as empirical formulae have now probably reached their highest point of refinement and development.

From the experimental point of view, the accurate determination of the cross sections $`d\sigma _{A+HB+C+\mathrm{}}/dE`$ is difficult, as it requires the identification of all the nuclei produced in the reaction. The total reaction cross section for a species $`A`$ is much easier to measure, as one has only to know that the reaction actually occurred. This is the reason why the question of total reaction cross section is usually separated from the problem of production cross section of secondaries. They will accordingly be considered separately in the next two parts.

##### Semi-empirical formulae

The first ingredient is the reaction (or inelastic) cross section $`\sigma _{\mathrm{inel}}(E)`$, giving the probability that a cosmic ray nucleus of energy $`E`$ undergoes a nuclear reaction with an interstellar nucleus. It should ideally be given by a fundamental approach, but as discussed above, this program has not been fulfilled yet. Simpler considerations lead to semi-empirical formulae which turn to be quite powerful in reproducing all existing data on nuclear reactions. In the simplest approach, the cross section is proportional to the geometrical area of the nucleus ($`\pi R_{eff}^2`$) that scales as $`A^{2/3}`$. In 1950, Bradt & Peters have proposed a first correction to take into account the overlap of the two nuclei wavefunctions,

$$\sigma _{\mathrm{inel}}=\pi r_0(A_{\mathrm{target}}^{1/3}+A_{\mathrm{proj}}^{1/3}b_0)^2,$$

where $`r_0`$ is the nucleon radius and $`b_0`$ the overlap parameter (or transparency).

In the eighties, Letaw and coworkers provided a refinement which is still widely used, in particular for $`p+A`$ reactions. At high energy, the cross section tends to the asymptotic energy independent form

$`\sigma _{\mathrm{inel}}^{\mathrm{HE}}\text{(mb)}=45A^{0.7}[1+0.016\mathrm{sin}(5.32.63\mathrm{ln}A)].`$
At lower energy, the behavior becomes energy dependent: one finds a slight hollow around 200 MeV/nuc, followed by a sharp enhancement and a sharp decrease till $`20`$MeV/nuc;

$`\sigma _{\mathrm{inel}}^{\mathrm{HE}}(E_k)=\sigma _{\mathrm{inel}}^{\mathrm{HE}}\text{(mb)}`$
$`\times [10.62\mathrm{exp}(E_k/200)\mathrm{sin}(10.9E_k^{0.28})].`$
In the above formula, $`E_k`$ denotes the kinetic energy per nucleon in MeV/nuc. More generally, there are some necessary corrections for light nuclei (due to their particular nuclear structure) and the reactions He+A are evaluated through an energy dependent rescaling factor $`[\sigma _{\mathrm{He}}/\sigma _\mathrm{p}](E_k)`$. Later on, in 1993, Sihver and coworkers, gave a formulation closer to the geometrical picture given above, to take also into account A+A reactions, thanks to new experimental data. Finally, in 1996, Wellish and Axen introduced refined functions in Sihver et al. formulæ in order to be able to describe the whole data set with a 2% accuracy from 6.8 MeV/nuc to 9 GeV/nuc (more than 1400 data points, but they only retain points whose accuracy is better than 4%).

##### Universal parameterization

As was shown in , the parameterization that yields the better fit to most measurements of nucleus-nucleus inelastic cross sections is currently the universal parameterization proposed by . It is given by

$`\sigma _{\mathrm{inel}}=\pi r_0^2\left(A_{\mathrm{proj}}^{1/3}+A_{\mathrm{target}}^{1/3}+\delta _E\right)^2`$ (16)
$`\times \left(1R_c{\displaystyle \frac{B}{E_{cm}}}\right).`$
Here, $`r_0=1\text{fm}`$ and $`E_{cm}`$ denotes the kinetic energy per nucleus in the rest mass frame. $`R_c`$ is a parameter introduced for light systems in order to keep the same formalism for all reactions. The term

$$\delta _E=1.85S+\frac{0.16S}{E_{cm}^{1/3}}C_E+\frac{0.91(A_t2Z_t)Z_p}{A_tA_p}$$

describes several effects:

* $`S=A_p^{1/3}A_t^{1/3}/(A_p^{1/3}+A_p^{1/3})`$ is the asymmetrical mass term that is related to the overlapping volume of the colliding system;
* $`C_E=D\left[1\mathrm{exp}\left(E/T_1\right)\right]0.292\mathrm{exp}\left(E/792\right)\mathrm{cos}(0.229E^{0.453})`$ models the energy dependence at intermediate and large energies, mainly through transparency and Pauli blocking ($`E`$ is the colliding kinetic energy per nucleus in MeV/nuc). The number $`D`$ is related to the system density and as $`T_1`$, it is adjusted once for all for almost all reactions with a more careful treatment and specific adjustment for light systems ($`A4`$).
* Finally, the last term in $`\delta _E`$ takes into account isotopic dependence of cross sections.

The remaining terms in Eq. (16), respectively

$$B=\frac{1.44Z_pZ_t}{R}$$

and

$$R=r_p+r_t+\frac{1.2\left(A_p^{1/3}+A_t^{1/3}\right)}{E_{cm}^{1/3}},$$

correspond to the coulomb barrier (that depends on the energy) and the radius to evaluate the barrier height. In the latter term, $`r_i`$ is the equivalent radius of the hard sphere and is related to $`r_{rms,i}`$ thanks to $`r_i=1.29r_{rms,i}`$ (this last quantity is obtained directly through experiment ).

##### Experimental data and accuracy

As an illustration, Fig. 6 shows these inelastic cross sections for the reactions (p,He)+9Be, (p,He)+12C along with the experimental data points. We see that this parameterization reproduces well the non-trivial features of the cross section.

However, it must be stressed that no new measurement has been made since 1985 for the reaction of protons on light nuclei, and that many isotopes have even never been measured to date. Moreover, the helium induced reactions are even more difficult to probe experimentally. Some of the results concerning the abundance of some species of cosmic rays rely on the faith that the parameterization is valid for these nuclei. Considering that when data exist, the parameterization gives an accuracy of a few %, it is reasonable to consider that the above parameterization is $`25`$% accurate for proton induced reactions and $`1020`$% for He induced reactions.

#### 2.4.6 Secondary production from spallations on the interstellar matter

When a nucleus undergoes a spallation, it may generate a large variety of lighter nuclei, which are not all present in standard sources. These nuclei may typically be classified as 2H, 3He (p and He induced), LiBeB (CNO induced) and sub-Fe (Fe induced). They give important clues to understand the characteristics of galactic propagation. As another example of the importance to have a good knowledge of these production cross sections, one can also mention cosmogenic nuclei (due to the penetration of GeV CR protons in extraterrestrial bodies) which allow an estimation of the stability of cosmic ray fluxes in the past billion years . To study the nuclei population $`N_2`$ and to draw an inference for propagation models, we need to now, as accurately as possible, processes such as

$$N_1+(p,He)N_2.$$
(17)
In a more rigorous way, one should also take into account contributions from heavier species of the ISM (e.g. $`N_1`$+CNO). However, these cross sections are only poorly estimated, not to say not known; moreover, these have been correctly included for antiproton production but were shown to be negligible contributors.

##### Semi-empirical and empirical formulæ

In the literature, there exists basically two fast procedures to evaluate the spallation cross sections corresponding to reactions (17). They were provided by cosmic ray physicists because the energy range in which these reactions are needed is generally not of great interest for nuclear physicists. Silberberg & Tsao  approach has well-motivated theoretical grounds and takes advantage of some observed regularities in i) the mass difference between target and daughter nuclei ($`\mathrm{\Delta }AA_{\mathrm{target}}A_{\mathrm{produced}}`$); ii) the ratio between the number of neutrons and protons in the daughter nucleus.

In the nineties, many new data issued from various targets were obtained, and this semi-empirical parameterization gave only a $`\pm 35`$% accuracy for these new data. It led Webber and coworkers  to develop a new approach fully based upon experimental regularities. In this parameterization, having remarked that production of different isotopes has similar $`E`$ dependence, parameters were fit to one kinetic energy per nucleon, and then extended to other energies,

$`\sigma ^{p+N_i}(Z_i,A_i;Z_f,A_f;E)=\sigma _0(Z_f,Z_i)`$
$`\times f_1(Z_f,A_f,Z_i,A_i)f_2(E,Z_f,Z_i).`$
We do not enter into the details (the reader is referred to quoted papers), but only briefly comment the three terms that enters in this equation: the first one describes the dependence on the fragment charge (the production is exponentially suppressed far from the stability valley) and contains in particular two parameters ($`\sigma _{Z_f}`$ and $`\mathrm{\Delta }_{Z_f}`$) that are chosen to fit the data. Notice that these parameters, thanks to new measurements, were updated in a later version of their code; we implemented in our propagation code parameters as given in Tab. V of . The second term describes isotopic distribution of fragments for a given species. This is independent of energy as long as we are not too far from the stability valley. The third term is the energy dependence and is only related to charges involved. Added to these three terms, the particular single and double nucleon-removal reactions are also taken into account, and He-induced spallations are tabulated from  formulæ.

##### Some remarks about these parameterizations

Actually, for our purpose, the most suitable set of formulae is given by Webber and coworkers’ empirical formulae. It seems that the more general approach of Silberberg and coworkers is more useful to propagate $`Z>30`$ heavy nuclei (see also ) as the corresponding cross section have not been measured yet. It was found in  that for the existing data points, the dispersion obtained with Webber’s code is better than the dispersion using theirs. However, the two approaches have probably reached their latest development. As remarked recently in , any refinement of these formulae asks for complicate modifications of the present formulae. If one takes into account isospin effect, odd-even effect, neutron and proton stripping due to coulomb barrier, it would require an adjustment nucleus by nucleus. Webber’s code is only valid for $`Z>3`$: light systems require further modelling because they behave very differently from other nucleus (see e.g. ). To improve these codes, some authors have focused on peculiar single and double nucleon removal (coulomb fragmentation) that are dominant for most secondary productions . A further point of discussion is the reliability of the straight-ahead approximation, i.e.

$`{\displaystyle _0^{\mathrm{}}}\beta ^{}N^k(T^{})\sigma ^{kj}(T,T^{})𝑑T^{}=\beta N^k(T)\sigma ^{kj}(T).`$
In practice, this corresponds to the fact that each fragment has the same kinetic energy per nucleon than has the parent nucleus. If this approximation seems quite justified at high energy (above a few GeV/nuc), its relaxation may affect, for example, the B/C ratio around 1 GeV/nuc . In our propagation code, we use Webber et al.’s code that is not only very well suited, but also widely used by almost all cosmic ray physicists (for the nuclei we deal with).

##### Experimental data and accuracy

At the present time, it is estimated that roughly $`95`$% of the reactions involving cosmic ray nuclei have been measured at least at one energy. As an illustration, we give in Tab. 2 the cross sections for the reactions involving 12C. It illustrates the importance of single nucleon removal (12C $``$ 11C, 11B) as mentioned above. Double nucleon removal (12C $``$ 10B) is less marked but visible. Let us notice also a maximum around 6Li that corresponds to half of the initial nucleus. This table also emphasizes the multi-fragmentation character of spallations.

Tab. 3 illustrates another aspect of spallations in cosmic ray propagation: the output is from Webber’s spallation code and has been weighted by the relative abundances of sources to isolate the most important contributions. This table underlines the fact that, if most secondaries are created through spallation of the most abundant species, all the others come from small contributions of a great numbers of nuclei. Finally, the present data accuracy ranges from $`23`$% to $`20`$% depending on the produced fragment. The old activation measurements seem to be unusable, and despite numerous experiments since the beginning of the eighties, not all the nucleus have been at least measured twice. Recently, Flesh et al. (1999)  observed systematic deviations from Webber measurements. It is thus not an easy task to estimate the accuracy of the above parameterizations. For the case of B/C ratio that is particularly important in our studies, one can estimate the accuracy to lie in the range of about $`510`$% for p induced reactions and $`1020`$% for He induced reactions. It is worth noting that production cross sections will probably be the limiting factor of cosmic ray studies in the near future, when high statistic experiments (pamela, ams) will obtain new very accurate data on CR fluxes. This is already the case for antiproton predictions (see  and Sec. 6.1). For secondary nuclei such as B, sub-Fe, consequences of these errors are difficult to estimated, because of the complex reaction chain involved (see  for a discussion).

##### Ghost nuclei

There remains a last point which is only rarely addressed. As mentioned above, a complete grid of nuclear reactions is needed to implement the cross sections, including short-lived species which are not seen directly but that contribute to some other species through their decay products. As an example, let us consider the production of 9Be. One has to add all the contributions from all heavier nuclei, including the specific contribution from intermediate reactions such as

$$N^k+\text{(p,He)}$$

$$\{\begin{array}{ccc}{}_{}{}^{9}\text{Li}& \stackrel{\beta ^{}(t_{1/2}=178ms)}{}& {}_{}{}^{9}\text{Be}(r=0.49).\\ {}_{}{}^{11}\text{Li}& \stackrel{\beta ^{}+\mathrm{\hspace{0.25em}2}n(t_{1/2}=8.6ms)}{}& {}_{}{}^{9}\text{Be}(r=0.041).\end{array}$$

The nuclei 9Li and 11Li are called ghost nuclei: they contribute indirectly to the flux, but do not have to be propagated as their half-lives are negligible in the face of propagation time scale. This will be the case for all the nuclei having $`t_{1/2}20`$ Myr. It turns out that there is a gap in the lifetimes of nuclei between 10 kyr and a few Myr, so that all nuclei having $`t_{1/2}1`$ kyr could be considered as ghosts. Among them, those which have too short a lifetime are not seen in the experiments and are already accounted for in the measured cross sections. We emphasize that, would we be interested in supernovae explosions, the nuclei relevant would not be the same because of a much shorter typical evolution time.

In conclusion, even if it will be transparent in the rest of the paper – and also in almost all propagation papers –, one has to keep in mind that each time we evaluate production cross section, in the true computation intervenes the quantity ($`r(Xj)`$ is the branching ratio of $`X`$ into $`j`$)

$$\sigma _{ij}^{\mathrm{effective}}=\sigma _{ij}^{\mathrm{direct}}+\underset{\text{ghosts}X}{}\sigma _{iX}r(Xj).$$

Nuclear grids for all species (including very heavy nuclei) can be found in . Thanks to recent nuclear compilations , we reevaluated this data grid (and the associated ghost nuclei) for $`Z<30`$; we found many new ghost nuclei or modified branching ratio. However, these new nuclei are generally far from the stability valley, such that their production cross section is clearly very small. We checked that for most species, this completion have a negligible effect, so that the ’s table of ghost nuclei is sufficient and does not need to be further revised 222A complete table of ghost nuclei for $`Z<30`$ can be found in Maurin 2001 (PhD dissertation)..

##### Production of antiprotons by spallations

Among the nuclei that are created by spallations, the antiprotons require a specific treatment, as their creation relies on the conversion of collision energy into mass and not on the mere breaking of a heavier nucleus. In order to evaluate the contribution from the $`\mathrm{p}\mathrm{H}_{\mathrm{ISM}}`$ interactions – the dominant process at first approximation – the Tan & Ng parameterization based on experimental data is used . Collisions that involve heavier nuclei are more difficult to obtain. The contribution $`\mathrm{p}\mathrm{He}_{\mathrm{ISM}}`$ has been considered firstly by means of a simple geometrical approach . However, using a more sophisticated nuclear Monte Carlo treatment, it was noticed in  that this reaction does not only enhance the antiproton flux as a whole but also change its low energy spectrum, mostly for kinematical reasons. By the way, these authors also considered $`\mathrm{He}\mathrm{p}_{\mathrm{ISM}}`$, $`\mathrm{He}\mathrm{He}_{\mathrm{ISM}}`$ as well as CNO contributions, the latter being negligible (see their Fig. 7). Unfortunately, very few experimental data are available on antiproton production cross sections in nuclear collisions and a model-based evaluation – such as the dtunuc program  – seems to be unavoidable. This is what we will use in the following.

#### 2.4.7 Antideuteron production

When an antiproton and an antineutron are formed by a reaction, they may merge into an antideuteron nucleus. The calculation of the probability for the formation of an antideuteron proceeds in two steps. We first need to estimate the probability for the creation of an antiproton-antineutron pair. Then, there is some probability that those antinucleons merge together to yield an antinucleus of deuterium.

As regards the first step, the differential probability for the production of a single antiproton or antineutron (resp. $`𝒫_{\overline{p}}`$ or $`𝒫_{\overline{n}}`$) is known for each initial process – spallation, neutralino annihilation or black hole evaporation – and a first guess would be that the production of two antinucleons is proportional to their product. This so-called factorization hypothesis is fairly well established at high energies. For spallation reactions, however, the bulk of the antiproton production takes place for an energy $`\sqrt{s}10`$ GeV which turns out to be of the same order of magnitude as the antideuteron mass. Pure factorization should break in that case as a result of energy conservation and needs to be slightly adjusted. We have therefore assumed that the center of mass energy available for the production of the second antinucleon is reduced by twice the energy carried away by the first antinucleon. This yields the following factorization of the probability to form a $`\overline{p}`$-$`\overline{n}`$ pair (see  for a more detailed discussion)

$`𝒫_{\overline{p},\overline{n}}(\sqrt{s},\stackrel{}{k_{\overline{p}}},\stackrel{}{k_{\overline{n}}})=`$
$`{\displaystyle \frac{1}{2}}𝒫_{\overline{p}}(\sqrt{s},\stackrel{}{k_{\overline{p}}})𝒫_{\overline{n}}(\sqrt{s}2E_{\overline{p}},\stackrel{}{k_{\overline{n}}})`$
$`+(\stackrel{}{k_{\overline{p}}}\stackrel{}{k_{\overline{n}}}).`$
Once the antiproton and the antineutron are formed, there is a finite probability that they combine together to give an antideuteron. The coalescence function $`𝒞(\stackrel{}{k_{\overline{p}}}\stackrel{}{k_{\overline{n}}})`$ describes the probability for a $`\overline{p}\overline{n}`$ pair to yield by fusion an antideuteron, as a function of the difference of the initial momenta. An energy of $`3.7`$ GeV is required to form by spallation an antideuteron whereas the binding energy of the latter is $`B2.2`$ MeV. The coalescence function is therefore strongly peaked around $`\stackrel{}{k_{\overline{p}}}\stackrel{}{k_{\overline{n}}}=\stackrel{}{0}`$.

We considered that the antinucleons merge together if the momentum of the corresponding two-body reduced system is less than some critical value $`P_{\mathrm{coal}}`$. That coalescence momentum is the only free parameter of this factorization and coalescence scheme. As shown in , the resulting antideuteron production cross section in proton-proton collisions is well fitted by this simple one-parameter model. The result, for a typical value $`p_0=160`$ MeV , as well as the comparison with the antiproton yield is displayed in Fig. 7.

New calculations for such astrophysical applications of the antideuteron production are under progress in the diagrammatic approach to the coalescence model .

#### 2.4.8 $`\beta `$ and EC decay

A second catastrophic loss to add to inelastic reactions is the decay of unstable nuclei. One can distinguish two decay processes: $`\beta `$ decay which is spontaneous and electronic capture decay (hereafter denoted EC for short) which can happen only if an electron of the interstellar medium has been first attached in the K shell (so that the presence probability of the electron is finite in the region of the nucleus).

##### $`\beta `$ decay

The ubiquitous example of this case is 10Be$``$10B whose rest half-life is $`t_{1/2}=1.51`$ Myr. The ratio 10Be/9Be observed in cosmic rays is smaller than the ratio of two normal secondaries coming from the same primary progenitors, which is merely the ratio of the production cross sections. This has often been interpreted as an evidence for an extended diffusive halo where the cosmic rays would spend most of their time. This interpretation may be misleading and could lead to a wrong intuition . One of the main interest of these nuclei is that the short-lived unstable species propagate only very locally, and thus are very sensitive to the local environment. In particular, they may be used to give an independent evidence for the existence of a local underdense bubble, as shown in  (see Sec. 6.3 for more details).

We compile in Tab. 4 the nuclei that can be propagated as pure $`\beta `$ unstable. Those with half-time shorter than 14C are ignored (see Sec. 2.4.6).

In particular, we have checked that the EC mode can be neglected for 27Al and 36Cl but not for 54Mn and 56Ni (EC half-life determination for these two last species is particularly difficult, see resp.  and ).

##### Electronic capture

A second kind of unstable species is given by nuclei decaying under EC process. A nucleus such as 59Ni is formed during last stages of stellar nucleosynthesis; it decays as 59Ni$``$59Co in $`t_{1/2}=80`$ kyr once it attaches an electron in its K-shell. Studies of this species allows for example to show that there was a delay of at least $`10^5`$ yr between nucleosynthesis and acceleration. Attachment and the converse process, electron stripping, are crucial processes that determine the effective lifetime of the EC unstable nucleus. Notice that because of these processes, the energy dependence of lifetimes is more complex in the EC mode than in the $`\beta `$ mode.

There are two ways to attach an electron, radiative and non-radiative process: after attachment of a free electron, the momentum and energy balance are restored by the emission of a photon or the recoil of the nucleus respectively. Radiative capture is dominant around a few hundreds of MeV/nuc in hydrogen. Non radiative capture scales as $`Z_t^5`$ and dominates for the heavier species of the ISM. Only a few measurements are available for the latter process, and they suggest that this contribution is anyway less than a few %; it is neglected in this work. It must be said that these processes suffer from the same limitations that the production cross sections, i.e. they are not studied by atomic physicist, but by cosmic ray astrophysicists. The results widely used nowadays have been derived in the seventies and to our knowledge no updates have been made at GeV energies.

Fig. 8 displays the rate of attachment and stripping for several charges. Some general conclusions can be drawn from this single result: i) comparing $`\tau _{\mathrm{attach}}`$ to the general propagation time $`20`$ Myr, only low energy and large $`Z`$ nuclei are likely to attach an electron and subsequently decay; ii) comparing $`\tau _{\mathrm{attach}}`$ to the stripping time, once an electron is captured, decay is much more likely than stripping; iii) consequence of the previous point is that for $`Z<30`$ species, nuclei are always at most singly attached.

Actually, as attachment arises only when there are available electrons, one should take into account the local properties of the interstellar medium (LISM) in all studies involving EC-species (as for radioactive $`\beta `$-decay, for different reasons). These attachment and stripping cross sections are included in the propagation code used for the work presented here. We should mention that these effects may be neglected in most cases, e.g. for B/C studies.

#### 2.4.9 Convection and reacceleration

The ingredients discussed here were introduced with the first theoretical development of cosmic ray physics, but they were only taken into account quite recently in the propagation models.

##### Convective wind

It has been recognized for a long time that a thin disc configuration would be disrupted by cosmic ray pressure . It can be stabilized by the presence of a halo, but further considerations imply that this halo would not be static either. The stellar activity and the energetic phenomena associated to the late stage of stellar evolution may push the interstellar plasma and the magnetic field associated with it out of the galactic plane. The net result would be the presence of a convective current directed outwards from the galactic plane and called galactic wind, which adds a convective term to the diffusion equation. The properties of this wind may be studied by hydrodynamical  and magneto-hydrodynamical (MHD)  methods.

Consequences of a wind have been first investigated by Ipavitch . Since then, it has been observed in other galaxies through a flattening of the electron induced radio spectrum around 1 GHz in the halo  that could be neither predicted by a diffusion model, nor by a pure convection model. In our own Galaxy, this effect is more difficult to be clearly established , and its effects have been investigated in particular on nuclei in various models .

Tab. 5 compiles values obtained in several studies for very simple forms of galactic wind. As one can see, almost all the authors set only upper limits on this possible wind.

Actually, the exact form of galactic winds is not known. From a self-consistent analytical description including MHD calculations of the galactic wind flow, cosmic ray pressure and the thermal gas in a rotating Galaxy, a wind increasing linearly with $`z`$ up to $`z15`$ kpc, with a $`z=0`$ value of about 22.5 km s-1 was found in  (see also references therein). Following a completely different approach,  extract the velocity form able to reproduce data from a one-dimensional diffusion/convection model. They obtain a decrease from 35 km s-1 to 12 km s-1 for $`z`$ ranging from 40 pc to 1 kpc followed by an increase to 20 km s-1 at about 3 kpc. These two recent works raised a positive conclusion about presence of a wind. For reference, the values for our best B/C fits correspond to about 15 km s-1. The difficulty to compare constant wind values to other $`z`$-dependences is related to the fact that cosmic rays do not spend the same amount of time at all $`z`$, so that there cannot be a simple correspondence from one model to another. As a result, all the above-mentioned models are formally different, with, moreover, many different inputs (spectral index, diffusion slope) that complicate a possible comparison. Nevertheless, their values are roughly compatible, the grounds for a physical motivation for this wind being provided in . One can notice that other models (e.g. ) dislike convective wind in their fit, but we will come back to this point in Sec. 5.3. To conclude about galactic winds, a radial dependence is not excluded .

##### Reacceleration

The term acceleration is used for the process by which low energy nuclei gain a huge amount of energy in a short time, promoting them to the rank of Cosmic Ray Nuclei. Once the cosmic rays are accelerated, there are other processes that can lead to less sudden energy gains, referred to as reacceleration. There are two classes of such processes. For the first, called sporadic reacceleration , it is assumed that there are some well localized reaccelerating centers (e.g. supernovæ), on which the cosmic rays can scatter and gain a small amount of energy at each of these scatterings. For the second, referred to as continuous , reacceleration occurs during the wandering of the charged nuclei through the Galaxy, as an unavoidable consequence of spatial diffusion. The quasi-linear theory of Boltzmann equation, where all effects are averaged over statistical properties of the plasma (see e.g. ), indicates that the net effect is a diffusion in energy space . The typical times for diffusion in space and energy are related through the relation $`\tau _{\mathrm{spatial}}.\tau _{\mathrm{reac}}L^2/V_a^2`$ (see e.g. ), where $`V_a`$ is the Alfvén speed of the medium. Some basic features of these two modes of reacceleration can be found in Giler et al. (1989) . It is likely that both these modes are present (see as a good example, the treatment of ).

The most direct observational evidence for reacceleration comes from cosmic ray species which are unstable via electronic capture. Their abundance relative to stable isotopes shows that they had some time to decay and therefore to attach an electron, which is more difficult at high energy. This means that these EC unstable species had a lower energy at the beginning of their propagation stage: the nuclei have gained about 100-200 MeV/nuc for a kinetic energy of a few hundreds of MeV/nuc . A less direct evidence is related to the question of isotropy of cosmic rays measurements, which favors small diffusion slope (e.g. Kolmogorov $`\delta =1/3`$) whereas the classical propagation models (no reacceleration) prefer larger values $`\delta 0.6`$ (see  and ). It turns out that reacceleration allows to redeem small $`\delta `$.

In our studies , we further investigated the allowed values of the reacceleration parameter, i.e. the Alfvén speed $`V_a`$. We tried to compare these values with those used by other authors . The difficulty is that the secondary to primary ratios are determined by an effective value:

$$V_a^{\mathrm{eff}}=V_a^{\mathrm{true}}\times \sqrt{\frac{h}{h_{\mathrm{reac}}}}\times \sqrt{\omega },$$
(18)
where $`h`$ and $`h_{\mathrm{reac}}`$ are respectively the height of the diffusive and reacceleration zone. The parameter $`\omega `$, which may depend on $`z`$, characterizes the level of turbulence and is often set to 1 . Our model, as others, uses

$`\omega (z)=\{\begin{array}{cc}1& \text{if }z<h_{\mathrm{reac}}\text{ ,}\hfill \\ 0& \text{otherwise ;}\hfill \end{array}`$ (21)
as a crude approximation of the more complex reality. The total reacceleration rate, at least in a first approximation, is given by the convolution of the time spent in the reacceleration zone with the corresponding true Alfvén speed in this zone. There are no direct observational clues about the size of the reacceleration zone, or about $`\omega (z)`$. As a result, even if the analysis of cosmic rays can give values of $`V_a^{\mathrm{eff}}`$, this does only give the value $`V_a^{\mathrm{true}}`$ up to an unknown factor $`h/h_{\mathrm{reac}}`$. The situation is even more complex if no assumption is made about $`\omega (z)`$. If $`\omega (z)`$ strongly depends on $`z`$ in a large reacceleration zone, then the expression above may still be used, with $`h_{\mathrm{reac}}`$ being an effective size depending on the vertical distribution of cosmic rays. However, this situation seems to be unfavored by MHD simulations .

To conclude with an order of magnitude, our best models indicate that $`\sqrt{h/h_{\mathrm{reac}}}`$ must be $`1/4`$ in order to give realistic values 333There is a misprint in  where 1/4 should be read instead of 4. for $`V_a`$ (with $`\omega 1`$.

### 2.5 Different approaches of propagation

The equations describing diffusion have been written down in Sec. 2.3.6, as Eq. (8). There are many approaches to solve the same diffusion equation, i.e. direct resolution (if possible), extraction of the Green function of the problem, or separation of astrophysical and nuclear aspects through the weighted slab technique. Most of the theoretical work on cosmic rays makes an extensive use of the weighted slab model and of the so-called Leaky Box model. We now present a synthetic overview of these models and give some clues about their successes and failures.

#### 2.5.1 Leaky Box

One of the simplest (though very simple) model is the so-called Leaky Box, in which the Galaxy is described as a finite propagation volume, delimited by a surface. Inside this volume, the densities of sources, interstellar matter and cosmic rays are homogeneous. Moreover, each nucleus has a probability per unit time $`1/\tau _{\mathrm{esc}}`$ to escape from the box. In the stationary regime, the densities are given by

$$\frac{N^j}{\tau _{\mathrm{esc}}}+nv\sigma ^jN^j=q^j+\underset{\mathrm{heavier}k}{}nv\sigma ^{kj}N^k.$$
(22)
This model has been successful to explain most of observed cosmic ray stable fluxes at different energies by a single function $`\tau _{\mathrm{esc}}(E)`$. This function can be either adjusted to the data, its physical interpretation being found afterwards, or extracted directly from more complete propagation equations . The second approach is a good way to understand why Leaky Boxes work so well, and is more easily seen in the framework of the weighted slab formulation. Thus we now turn to this point, introducing first the slab model.

#### 2.5.2 Slab model

Independently of the exact processes which are responsible for propagation, the first certainty we have about cosmic ray nuclei is that they have crossed some interstellar matter between their creation in the sources and their detection. This leads to spallation reactions that alter the initial (source) composition by destroying primary species and producing secondaries. The ratio of secondary to primary species fluxes gives information about the quantity of matter crossed. It is convenient to introduce the column density of matter crossed by a particle, also called the grammage, expressed in g cm-2. All the nuclei of a given species, with a given energy, do not have the same propagation history, in particular they have not crossed the same amount of matter, so that a distribution of grammages is associated with each species.

In a first step, though, one can assume that all nuclei of a given species have crossed the same grammage. This is called the slab model. Formally, the number $`\stackrel{~}{N}^j(x)`$ of nuclei $`j`$ that have crossed the grammage $`x`$ is related to the destruction rate (the inelastic cross section is denoted $`\sigma ^j`$) and the creation rate (i.e. the spallation rate of all heavier nuclei $`k`$ giving $`j`$, the cross section being denoted $`\sigma ^{kj}`$) by

$$\frac{d\stackrel{~}{N}^j(x)}{dx}+\frac{\sigma ^j}{\overline{m}}\stackrel{~}{N}^j(x)=\underset{\mathrm{heavier}k}{}\frac{\sigma ^{kj}}{\overline{m}}\stackrel{~}{N}^k(x),$$
(23)
with the initial condition $`\stackrel{~}{N}^j(x=0)=q^j`$. The resolution of this equation yields the secondary to primary ratios as functions of $`x`$. Comparison with observations give the value of $`x`$. This leads to a contradiction, as weighting correctly the production rate by all the parent nuclei, the observed LiBeB/CNO ratio $`0.28`$ gives $`x4.8`$ g cm-2, whereas the ratio sub-Fe/Fe$`1.5`$ gives $`x0.8`$ g cm-2 (see e.g. ).

#### 2.5.3 Weighted slab and other propagation models

The previous model is too simple and one should take into account that the nuclei of the same species and same energy may cross different amounts of matter. Introducing the probability $`G(x)`$ that a nucleus $`j`$ has crossed the grammage $`x`$, the density of a nucleus $`j`$ is given by

$$N^j=_0^{\mathrm{}}\stackrel{~}{N}^j(x)G(x)𝑑x.$$
(24)
The function $`G(x)`$ is called the Path Length Distribution (hereafter referred to as PLD). It is remarkable that the use of the same function $`G(x)`$ for all nuclei can solve the apparent contradiction between the grammage seen by CNO and Fe: the Fe nuclei have a larger destruction cross section, so that they are more easily destroyed. The associated function $`\stackrel{~}{N}^j(x)`$ decreases more rapidly and the sub-Fe/Fe ratio is more sensitive to the low $`x`$ part of the grammage distribution function $`G(x)`$.

##### Weighted slab as models or as a general technique

The probability distribution $`G(x)`$ can be determined by the choice of a propagation model. For example, the slab model is obtained by considering $`G(x)=\delta (xx_0)`$. Actually, the weighted slab approach has been widely used in literature, under two slightly different forms. In the weighted slab model, one tries empirically to modify the Path Length Distribution $`G(x)`$ to account for the data. The physical meaning of $`G(x)`$ may be explored in a second step. The weighted slab technique is more general, and the name refers to the possibility to introduce Eq. (24) containing the additional and meaningful variable $`x`$ (the grammage) in any propagation equation (such as the LB model, diffusion model, diffusion/convection model…). Example of direct extraction of PLD can be found in . Up to the end of this section, we will omit this subtility. The weighted slab approach allows to link Leaky Box models with more realistic diffusion models, explaining why these Leaky Boxes work so well.

##### Path Length Distribution of the Leaky Box model

As an alternative to the direct analytical resolution of Eq. (22), one can insert Eq. (24) in Eq. (22) which separates the nuclear part from the other effects. This gives

$`\{\begin{array}{cc}G^{\mathrm{LB}}(x)={\displaystyle \frac{1}{\lambda _{\mathrm{esc}}}}\mathrm{exp}\left({\displaystyle \frac{x}{\lambda _{\mathrm{esc}}}}\right);\hfill & \\ {\displaystyle \frac{d\stackrel{~}{N}^j(x)}{dx}}+{\displaystyle \frac{\sigma ^j}{\overline{m}}}\stackrel{~}{N}^j(x)={\displaystyle \underset{k>j}{}}{\displaystyle \frac{\sigma ^{kj}}{\overline{m}}}\stackrel{~}{N}^k(x).\hfill & \end{array}`$ (27)
In this expression, the PLD $`G(x)`$ is expressed as a function of $`\lambda _{\mathrm{esc}}(E)=\overline{m}nv\tau _{\mathrm{esc}}(E)`$. An interesting feature of this model is that the average grammage, given by $`x=xG(x)𝑑x`$ is exactly the escape length, i.e $`x=\lambda _{\mathrm{esc}}(E)`$. Actually, this quantity solely depends on energy if $`G^{\mathrm{LB}}(x)`$ is the same for all species. This PLD may be modified to give new weighted slab models. For example, some discrepancies have led authors in the past to introduce the double PLD  which is a mixture of two $`G^{\mathrm{LB}}(x)`$ with two different characteristic functions $`\lambda _{\mathrm{esc}}`$ and $`\lambda _{\mathrm{esc}}^{}`$. This can be interpreted as the presence of two different propagation zones (Simon et al., 1977, 1979 ), which supports the more realistic disc-halo model. Let us notice that the terms escape length and Path Length Distribution are sometimes mixed in the literature, which can generate confusion.

Finally, it is worth noting a specific point quite independent of the rest of the discussion. It has been noticed in  that a different but equivalent formulation could be obtained from the above set (27) of LB Weighted Slab representation, by the substitution

$$G^{\mathrm{LB}}(x)=(1/\lambda _{\mathrm{esc}})\mathrm{exp}(x/\lambda _{\mathrm{esc}})G(x)=1,$$

and the addition of a term $`N^j/\lambda _{\mathrm{esc}}`$ in the nuclear part of Eq. (27), i.e.

$$\frac{d\stackrel{~}{N}^j(x)}{dx}+\frac{N^j}{\lambda _{\mathrm{esc}}}+\frac{\sigma ^j}{\overline{m}}\stackrel{~}{N}^j(x)=\underset{k>j}{}\frac{\sigma ^{kj}}{\overline{m}}\stackrel{~}{N}^k(x).$$

The meaning of the new set of equations is clear: either one solves the stationary Eq. (22), or one solves the same but non-stationary equation integrating over all times, i.e. $`\stackrel{~}{N}^k(x)G(x)𝑑x\stackrel{~}{N}^k(x)𝑑x`$ (here, $`x`$ is a dummy variable, so that it can be called time). Such a procedure is very general. It is for example used in .

##### Weighted slab technique applied to diffusion models

The link between the diffusion equation and the escape time introduced in Leaky Box models was clarified by Jones in . He showed that diffusion models, energy losses included, can be reinterpreted in terms of a Leaky Box model (with some exceptions, such as $`e^{}`$ for which the synchrotron or inverse Compton losses are too important). This reinterpretation, called the leakage lifetime approximation can be understood by writing the diffusion equation, neglecting energy gains and losses and assuming steady state,

$$.(K^jN^j)+nv\sigma ^jN^j=q^j+\underset{k>j}{}nv\sigma ^{kj}N^k.$$
(28)
When applying the weighted slab technique (Eq. 24), the previous equation can be split into a purely diffusive and a purely spallative part

$`\{\begin{array}{cc}K\mathrm{}G(𝐫,x)=Q(𝐫,x);\hfill & \\ {\displaystyle \frac{d\stackrel{~}{N}^j(x)}{dx}}+{\displaystyle \frac{\sigma ^j}{\overline{m}}}\stackrel{~}{N}^j(x)={\displaystyle \underset{k>j}{}}{\displaystyle \frac{\sigma ^{kj}}{\overline{m}}}\stackrel{~}{N}^k(x).\hfill & \end{array}`$ (31)
The Path Length Distribution $`G(𝐫,x)`$ encodes all the propagation properties. This function has been extracted for various geometries (thin disc, spherical halo) and various forms of source spatial distribution $`Q(r)`$ . A general result is that as the operator $`\mathrm{}`$ is hermitic, the function $`G`$ can always be set under the form

$$G(𝐫,x)\underset{j=1}{\overset{\mathrm{}}{}}c_n(𝐫)\mathrm{exp}\left(\frac{x}{\lambda _n}\right).$$

Comparing to $`G^{\mathrm{LB}}(x)`$ derived above in Eq. (27), this shows that diffusion models can be equivalently expressed in terms of a sum of Leaky Boxes. Jones showed  that in most situations, only the first terms contribute in the sum above. This explains why the simplest Leaky Box models are so successful to describe diffusion, even in a complex geometry.

The specific case of primaries is noticeable. In this case, the nuclear part of Eq. (31) reduces to a very simple form to finally give

$`N^j(𝐫)={\displaystyle _0^{\mathrm{}}}G(𝐫,x)\mathrm{exp}(\sigma ^jx)𝑑x.`$ (32)
Mathematically, the solution is the Laplace transform of parameter $`\sigma ^j`$ of the function $`G(𝐫,x)`$. Several solutions corresponding to several diffusive geometries can be tabulated with this Laplace transform . This approach gives an alternative way to evaluate the average grammage $`x`$, i.e. the equivalent Leaky Box description, from $`N^j`$

$$x=\frac{_0^{\mathrm{}}xG(x)𝑑x}{_0^{\mathrm{}}G(x)𝑑x}=\left(\frac{d}{d\sigma ^j}\mathrm{ln}N^j\right)|_{\sigma ^j=0}.$$

To end with these relations between treatment of the same equations, the latter formulation is close to what one obtains using the random walk approach . This should be not too surprising since  showed a long time ago the equivalence of random walk description and diffusion processes.

#### 2.5.4 Limitation and usefulness of all these models

Actually, the fluxes derived with the weighted slab technique as depicted above differ slightly from the direct evaluation. Whereas Lezniak  tried to extend the simple scheme presented in Eq. (24) to include energy losses, Jones  demonstrated that the separation between the nuclear side and the propagation side is only approximately valid even in the high energy regime. The reason is that in the rewriting $`N^j(𝐫)=_0^{\mathrm{}}\stackrel{~}{N}^j(x)G(𝐫,x)𝑑x`$, the separation between $`\stackrel{~}{N}^j(x)`$ and $`G(𝐫,x)`$ can have only one parameter in common, otherwise, equations obtained from the initial propagation equation are generally inconsistent. Hence, one can think that energy appears only as a mute parameter, but this is not exactly the case: the astrophysical part, i.e. $`G(𝐫,x)`$ is the same for all species as long as rigidity is considered. Conversely, for the nuclear part, the equation obtained above is only valid for a fixed kinetic energy per nucleus.

However, redefining some parameters and under several conditions, the technique can be made exact , and this idea has been numerically and quantitatively validated in . Despite that, the leakage lifetime approximation is not valid in some cases a mentioned above. For some nuclei, a description in terms of Leaky Box may lead to wrong results. This is the case for radioactive species in a realistic Galaxy, as was shown by .

##### Concluding remarks

The Leaky Box models, due to their simplicity, are very well suited to the extraction of source abundances (elemental as well as isotopic). They can also be used to compute the secondary antiproton production, since the same processes as for secondary stable nuclei are at work. However, as emphasized in  (see also Sec. 2.4.6), they are not able to predict any primary contribution in the antiproton signal, since it requires the knowledge of the spatial distribution of primary progenitors. It is also well known that the Leaky Box parameters are just phenomenological with only a distant connection to physical quantities. However, even if apparently many existing effects cannot be correctly included, it is shown in  that the galactic wind can be accounted for in a Leaky Box description and emphasized that the phenomenological behavior of the escape length at low energy could be due to the presence of this galactic wind. This idea was investigated further in  with a generation of several equivalent phenomenological escape lengths from several possible physical configurations of a one-dimensional diffusion model. The relation between one-dimensional models and Leaky Box models is thus firmly established and very well understood. This relation also elucidates some of the physical contents of Leaky Box models.

Our model furnished the following step towards realistic description, because it is equivalent, up to several minor modifications, to the one-dimensional model of . Finally, in the Strong et al. model , all subtle effects can be studied and modelled, with the counterpart that the numerical approach makes the physical intuition of the results less straightforward and the computation slower. However, the physical input is almost the same as in our model, which explains why the results are very close.

Even if all these models are equivalent to describe the local observations of charged cosmic rays, they lead to very different conclusions and interpretations when the spatial variation of the cosmic ray density is considered. As an illustration of the poor current understanding of this global aspect, we mention the ever-lasting problem of the gamma ray excess about 1 GeV towards the galactic center  or the too flat radial $`\gamma `$-ray distribution observed in the disc .

### 2.6 Numerical implementation

From the above discussions, it should be clear that the evaluation of the fluxes involves the computation of a nuclear reactions grid and the resolution of differential equations describing diffusion.

##### Nuclear part

The flux of a given nucleus is the sum of the primary (source) contribution and all the secondary contributions, from the spallations of heavier nuclei. It may be convenient to write a matrix relation (easily diagonalized, see coefficients in ) between the fluxes $`\{N^j\}^{j=1\mathrm{}n}`$ and the sources $`\{q^j\}^{j=1\mathrm{}n}`$

$$[N]=[\alpha ][q]+[\beta ][N],$$

where the matrices $`[\alpha ]`$ and $`[\beta ]`$ contain the information about destructive spallations and secondary creation. This is particularly useful in the Leaky Box model, as this relation is straightforwardly inverted and can provide source abundances as well as associated errors without too much efforts . However, as soon as energy losses are considered, this becomes more complicated.

We chose a more direct method, i.e. the cascade method. The nuclei are classified according to their masses. The heaviest nucleus $`N_1`$ has no secondary contribution and is given directly by its source contribution $`q_1`$. Eq. (8) is solved for this first nucleus. The flux of the second heavier nucleus $`N_2`$ is a mere combination of the secondary contribution coming from $`N_1`$ and its own source contribution $`q_2`$. Then, Eq. (8) is solved for this second nucleus, using the result obtained for $`N_1`$. The flux of the third $`N_3`$ is a combination of its own source contribution $`q_3`$, plus spallative contributions from the two previous ones $`N_1`$ and $`N_2`$. This procedure is repeated for all the nuclei. These implicit contributions are displayed in Fig. 9: source abundances have been arbitrarily normalized to 100 (dashed right bars, Solar abundances  times first ionization potential, see e.g. ).

Left bars are for propagated fluxes at about 1 GeV/nuc in a simple Leaky Box model without energy losses. A first point is to notice that heavy nuclei are more subject to destruction during propagation than light nuclei. Secondary species also appear clearly on the plot. Coming back to the cascade procedure, three quantities are displayed: one has (left bars) the sole propagation of sources (dashed bars), i.e. primaries; faint shaded bars show the sole secondary contribution (direct contribution of heavier sources); the higher bars correspond to the summation of all intermediary steps. The effect of these higher order secondaries is particularly important for heavier species. An iterative method can also be used , implementing the previous step propagated abundances as sources for the next step, until convergence is reached.

##### Spatial part

The spatial part of the diffusion equation involves second order differential equations. In our semi-analytical approach, one deals with a second order equation in energy but no spatial derivative. The resolution is much simplified, compared for example to  where the full second order partial derivative transport equation is solved. As regards the energy part, the Runge-Kutta method is not suited and one must turn to more refined methods, such as the Crank-Nicholson scheme. The general method can be found in many numerical books . An alternative approach is random-walk Monte Carlo simulations. A diffusive process can be mimicked by a random walk, with step related to the diffusion coefficient . This approach is very time consuming but a few authors has used this technique since mid-seventies (e.g. Owens, 1976b ) to propagate self-consistently proton and electrons . Such an approach is useful when one wants to model very complex and inhomogeneous environments.

### 2.7 Quick survey of charged nuclei behaviors and their interest

We now have all the elements needed to compute the fluxes of all the cosmic ray nuclei. They all have different propagation histories, and they carry information of different nature (see  for a review). They can be classified as follows:

* Stable primaries and mixed species (i.e. secondary plus primary contributions): most elemental abundances point towards a “standard” origin (from a particle physicist point of view). Some nuclei present anomalous contributions that could be due to wind-induced enrichment from Wolf-Rayet stars surfaces. The correct abundance determination of nuclei belonging to this class also allows a better characterization of the acceleration processes.
* Pure stable secondary species: the three groups of secondary species ($`Z=12`$, $`Z=35`$ and $`Z=2123`$) give information on the propagation history for quite different charges. Their importance for the determination of the propagation parameters is extensively demonstrated in this review.
* $`\beta `$ unstable secondary species: though radioactive species were first used to show the existence of a large diffusive halo , we underlined in  the interest of 10Be, 26Al and 36Cl (54Mn is more difficult to tackle since it has a mixed decay mode, $`\beta `$ and EC) for the study the LISM, or as a possibility to derive the diffusion coefficient once LISM properties are fixed .
* Primary EC clocks: through Co/Ni, Co/Fe or the more precise 59Ni/60Ni ratios, these EC clocks indicate the time elapsed between synthesis and acceleration and are shown as circles in Fig. 10 (the latest ace results point towards $`t5\times 10^5`$ yr ).
* Secondary reacceleration indicators: other K-capture nuclei (open diamonds in Fig. 10) indicate that reacceleration during propagation could be as large as a few hundreds of MeV/nuc at low energy. The latest results from the ace experiment seem to confirm this reacceleration, but the propagation of EC unstable nuclei in a refined propagation model and a thorough analysis including peculiarities of the LISM has not yet been performed (Maurin et al., in preparation).

This classification is summarized in Fig. 10.

## 3 Solutions of the diffusion/convection equation for two-zone cylindrical models

### 3.1 General remarks

The density of cosmic rays of energy $`E`$ at the position $`(r,z)`$ is obtained by solving the energy-dependent diffusion equation (see Sec. 2.3.6). It turns out that it is possible to first focus on the spatial diffusion and, then to take the energy changes into account, as long as these energy changes take place in the disc. The density is then obtained by solving a Laplace equation in a cylindrical geometry. The standard method is to develop all the quantities over a suitable complete set of orthogonal functions involving the first Bessel function $`J_0`$. Introducing $`\rho =r/R`$, one can write

$$N(r,z)=\underset{i=1}{\overset{\mathrm{}}{}}N_i(z)J_0(\zeta _i\rho ),$$

with

$$N_i(z)=\frac{2}{J_1^2(\zeta _i)}_0^1\rho N(\rho R,z)J_0(\zeta _i\rho )𝑑\rho ;$$

the $`\zeta _i`$ are the successive zeros of the function $`J_0`$. This development is inserted into the diffusion equation, which can then be rewritten, using the properties of Bessel functions, as

$`N_i^{\prime \prime }(z)\left({\displaystyle \frac{2h\mathrm{\Gamma }_{\mathrm{inel}}}{K}}\delta (z)+{\displaystyle \frac{\mathrm{\Gamma }_{\mathrm{rad}}}{K}}+{\displaystyle \frac{\zeta _i^2}{R^2}}\right)N_i(z)`$ (33)
$`{\displaystyle \frac{V_c}{K}}N_i^{}(z)={\displaystyle \frac{q_i^{\mathrm{source}}(z)}{K}},`$
where $`\mathrm{\Gamma }_{\mathrm{inel}}`$ and $`\mathrm{\Gamma }_{\mathrm{rad}}`$ denote respectively the spallation rate over the interstellar medium and the radioactive decay rate. The Bessel expansion $`q_i^{\mathrm{source}}(z)`$ for a unspecified source term $`q(r,z)`$ has also been used.

### 3.2 Full solutions

We do not recall here the procedure used to solve Eq. (8) and we refer the suspicious reader to . Suffice to say that very general solutions can be written for each species. They consist in the sum of two terms, the first involving sources located in the disc (primary or secondary), and the second involving sources located in the whole halo (radioactive origin). For the sake of clarity, we present separately these two contributions.

#### 3.2.1 Contribution from progenitors in the disc

The galactic disc is a particular place for GCR, as it contains the standard primary sources $`q(r,z)=q^{\mathrm{disc}}(r)\delta (z)`$ (see Sec. 2.4.3), and it is the place where spallations occur. The contribution to the density is given by

$`N(r,z)=\mathrm{exp}\left({\displaystyle \frac{V_cz}{2K}}\right)`$ (34)
$`\times {\displaystyle \underset{i=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{𝒬_i}{A_i}}{\displaystyle \frac{\mathrm{sinh}\left[S_i(Lz)/2\right]}{\mathrm{sinh}\left[S_iL/2\right]}}J_0\left(\zeta _i{\displaystyle \frac{r}{R}}\right)`$
with

$`\{\begin{array}{cc}A_i=& 2h\mathrm{\Gamma }_{\mathrm{inel}}+V_c+KS_i\mathrm{coth}\left({\displaystyle \frac{S_iL}{2}}\right)\text{ ;}\hfill \\ S_i^2=& {\displaystyle \frac{4\zeta _i^2}{R^2}}+{\displaystyle \frac{V_c^2}{K^2}}+4{\displaystyle \frac{\mathrm{\Gamma }_{\mathrm{rad}}}{K}}\text{ .}\hfill \end{array}`$ (37)
For a pure primary $`𝒬_i=q_0Q(E)\times q_i^{\mathrm{disc}}`$ – where $`q_i^{\mathrm{disc}}`$ denotes the Bessel expansion of $`q(r)`$ – and for a pure secondary $`𝒬_i=_k^{m_k>m_j}\mathrm{\Gamma }^{kj}N_i^k(0)`$ – where the superscript $`j`$ denotes the nucleus evaluated, omitted in above expressions –.

When energy losses and diffusive reacceleration are taken into account, the differential equation is obtained, following the procedure described e.g. in , as

$`A_iN_i(0)=𝒬_i`$
$`2h{\displaystyle \frac{}{E}}\left\{b_{\mathrm{tot}}(E)N_i(0)\beta ^2K_{pp}{\displaystyle \frac{}{E}}N_i(0)\right\}`$
where the energy loss and reacceleration terms were discussed in Sec. 2.3.3 and 2.3.5. Subscript $`i`$ refers to Fourier-Bessel coefficients for the above equation expended on cylindrical Bessel basis $`J_0(\zeta _i\rho )`$.

#### 3.2.2 Contribution from the halo: radioactive $`𝜷`$ decay

All the nuclei treated here have at most one unstable progenitor. The contribution of these radioactive nuclei may be unimportant in some cases, but we should take it into account as it is the dominant process for some others. In the simple example of 10Be$`^{10}`$B, neglecting this channel would give an error of about 10% on the B flux, whereas considering that this term is only located in the disc would give an error of about 3% compared to the rigorous treatment given above. Due to lack of space, we do not reproduce the expression in this case, but the reader is referred to .

#### 3.2.3 Contribution from the halo: primary sources

We will consider the hypothesis that primary sources of cosmic rays may be present in the diffusive halo (see Sec. 2.4.4). The solution of the diffusion equation is somewhat different than in the standard case of sources located in the disc. Denoting $`q^{\mathrm{halo}}`$ the corresponding source term, the Bessel terms of the density are given by (see )

$`N_i(z)={\displaystyle \frac{y_i(z)}{KS_i}}+{\displaystyle \frac{y_i(L)e^{V_c(|z|L)/2K}}{A_i\mathrm{sinh}(S_iL/2)}}`$
$`\times \left[\mathrm{cosh}(S_iz/2)+{\displaystyle \frac{V_c+2h\mathrm{\Gamma }_{\mathrm{inel}}}{KS_i}}\mathrm{sinh}(S_iz/2)\right]`$
with

$`y_i(z)=2{\displaystyle _0^z}\mathrm{exp}\left({\displaystyle \frac{V_c(zz^{})}{2K}}\right)`$
$`\times \mathrm{sinh}\left({\displaystyle \frac{S_i(zz^{})}{2}}\right)q_i^{\mathrm{halo}}(z^{})dz^{}.`$
In particular, the density in the disc ($`z=0`$) is given by

$$N_i(0)=\mathrm{exp}\left(\frac{V_cL}{2K}\right)\frac{y_i(L)}{A_i\mathrm{sinh}(S_iL/2)}.$$

It can be checked that with a source term localized in the disc, i.e. $`q^{\mathrm{halo}}(r,z)=2h\delta (z)q^{\mathrm{disc}}(r)`$, the expression (34) is recovered.

#### 3.2.4 Modelling of the local bubble

As mentioned above, it is very likely that the Sun lies in an underdense bubble. The main effect on the cosmic rays is that spallations, and hence the secondary production and destruction, are locally strongly reduced. In the thin disc model, this can be taken into account assuming position dependent spallation terms. Actually, as will be discussed below (see Sec. 6.3), this effect is not so important for stable species, but it is crucial for the low energy radioactive species, which are only sensitive to the local propagation conditions. Once created at a given point, they diffuse in a small region surrounding this point before they decay. The density of their progenitors does not vary much in this region, and it can be considered as homogeneous. Moreover, the propagation of the radioactive is almost not altered by the boundaries of the diffusive volume, so that we can then put the center of the diffusive volume at the position of the Sun. The bubble can then be modelled as a hole in the disc. The rates $`\mathrm{\Gamma }_{\mathrm{spal}}(r)=n_{\mathrm{LISM}}(r).v.\sigma _{\mathrm{spal}}`$ and $`\mathrm{\Gamma }_{\mathrm{inel}}(r)=n_{\mathrm{LISM}}(r).v.\sigma _{\mathrm{inel}}`$ now depend explicitly on $`r`$ via the local interstellar density which reads

$$n_{\mathrm{LISM}}(r)=\mathrm{\Theta }(rr_{\mathrm{hole}})n_{\mathrm{ISM}}$$
(39)
where $`\mathrm{\Theta }`$ is the Heaviside distribution (see Fig. 11).

As for the no-hole model, a solution is found by expanding the density over the Bessel functions $`J_0(\zeta _ir/R)`$. The expression for the solution may be found in .

### 3.3 Physical interpretation of the formulæ

For the sake of clarity, we wish to further discuss two particular cases. To begin with, the case of primary stable species is of particular interest. Their density in the disc is given by

$$N(r,z=0)=\underset{i=1}{\overset{\mathrm{}}{}}\frac{q_i^{\mathrm{disc}}}{A_i}J_0(\zeta _i\rho )$$

which, when compared to the Bessel expansion of the source

$$q(r)=\underset{i=1}{\overset{\mathrm{}}{}}q_i^{\mathrm{disc}}J_0(\zeta _i\rho ),$$

contains basically all the physics of diffusion. Diffusion in Bessel space may be seen as a filtering process, with a filter $`1/A_i`$. The small scale features, which correspond to large values of the index $`i`$, are spread out and erased by diffusion as $`A_i2K\zeta _i\pi Ki`$ for large $`i`$. The large scale structure is determined by the small $`i`$ behavior. The effect of small $`L`$, spallations and galactic wind is to limit the decrease of $`A_i`$ for small values of $`i`$, i.e. they filter the large scale features out. In other words, they make diffusion a local process. More precisely for small $`L`$ and small $`i`$ we have $`A_i2h\mathrm{\Gamma }_{\mathrm{inel}}+V_c+2K/L`$. It appears that three characteristic lengths, although of different origin, play a very similar role in erasing the large scale features of the source distribution, namely $`L`$, $`r_{\mathrm{wind}}2K/V_c`$ and $`r_{\mathrm{spal}}2K/2h\mathrm{\Gamma }_{\mathrm{inel}}`$. Structures larger than these are destroyed by diffusion. This point is further discussed in Secs. 4.3 and 6.2.

## 4 Origin of cosmic rays

### 4.1 The question of the spatial origin

The cosmic rays emitted from distant sources are more likely to escape from the diffusive halo or to be destroyed by spallation before reaching Earth than those coming from nearby sources. Conversely, a CR eventually detected on Earth has a greater probability to come from a nearby source than a distant source. In this section, we investigate the question of the spatial origin in more details, and we find that under quite reasonable conditions, most of the cosmic rays detected in Earth were emitted by sources located in a radius less than $`3L`$.

### 4.2 Method

The question we wish to address is the following: a cosmic ray being detected at the position $`\stackrel{}{r}_o`$ of an observer (in practice, this will be the position of the Sun), what is the probability

$`d𝒫\{\mathrm{emitted}:\stackrel{}{r}_s,\stackrel{}{r}_s+d\stackrel{}{r}_s|\mathrm{observed}:\stackrel{}{r}_o\}`$
$`𝒫\left\{\stackrel{}{r}_s\right|\stackrel{}{r}_o\}d\stackrel{}{r}_s`$
that it was emitted from a source located at the position $`\stackrel{}{r}_s`$? Such a question falls among classical problems of statistics. A rigorous theoretical frame is provided by the Bayes approach that summarizes proper use of conditional probabilities. A cruder but sufficient (and equivalent) treatment is given through the frequency interpretation. The probability is simply given by the number of particular realizations, i.e. paths $`𝒩\left[r_sr_o\right]`$ (the cosmic rays travels from the source at position $`\stackrel{}{r}_s`$ to the observer at position $`\stackrel{}{r}_o`$) divided by the number of all realizations of $`𝒩\left[\stackrel{}{r}_o\right]`$ (the cosmic rays reach the observer position). Notice that the latter number is just a normalization so that the first quantity determines completely the problem. In the frequency interpretation, we thus have

$$𝒫\left\{\stackrel{}{r}_s\right|\stackrel{}{r}_o\}=\frac{𝒩\left[r_sr_o\right]}{𝒩\left[r_o\right]}.$$
(40)
We finally notice that the number of paths $`𝒩\left[r_sr_o\right]`$ determines the number of cosmic rays that reach position $`\stackrel{}{r}_o`$, when a source is placed at position $`\stackrel{}{r}_s`$. We can thus write

$$𝒫\left\{\stackrel{}{r}_s\right|\stackrel{}{r}_o\}N_{\mathrm{r}_\mathrm{s}}(\stackrel{}{r}_o),$$
(41)
where the density $`N_{\mathrm{r}_\mathrm{s}}(\stackrel{}{r}_o)`$ is the solution of the propagation equation for a point source located at $`\stackrel{}{r}_s`$. The normalization factor of this relation is obtained by imposing that $`𝒫`$ actually is a probability, i.e. is normalized to unity. If the sources are distributed according to $`q(\stackrel{}{r}_s)`$, the probability that a cosmic rays detected at $`\stackrel{}{r}_o`$ was emitted from a volume $`𝒱^{}`$ (or a disc surface) is given by

$$𝒫\left\{𝒱^{}\right|\stackrel{}{r}_o\}=\frac{_𝒱^{}q(\stackrel{}{r}_s)N_{\mathrm{r}_\mathrm{s}}(\stackrel{}{r}_o)d^3\stackrel{}{r}_s}{_{𝒱_{\mathrm{tot}}}q(\stackrel{}{r}_s)N_{\mathrm{r}_\mathrm{s}}(\stackrel{}{r}_o)d^3\stackrel{}{r}_s}.$$
(42)
This probability contains all the physical information about the spatial origin of cosmic rays and may be computed for different situations (see  for further details).

### 4.3 Pure diffusion case and key parameters for spallations and wind

The extent of the zone from which the cosmic rays detected on Earth actually come from is limited by all the phenomena that keep them from propagating over large distances: escape through the boundaries, spallations and galactic wind. It comes as no surprise, then, that we can define three typical lengths, beyond which the cosmic rays do not reach us

$$r_L=L;$$

$`r_{\mathrm{spal}}K/h\mathrm{\Gamma }_{\mathrm{inel}}`$
$`3.17\text{kpc}\times {\displaystyle \frac{K/\beta }{0.03\text{kpc}^2\text{Myr}^1}}{\displaystyle \frac{100\text{mb}}{\sigma _{inel}}};`$
$`r_{\mathrm{wind}}2K/V_c`$
$`5.87\text{kpc}\times {\displaystyle \frac{K}{0.03\text{kpc}^2\text{Myr}^1}}{\displaystyle \frac{10\text{km}\text{s}^1}{V_c}}.`$
The role of the length $`r_{\mathrm{wind}}`$ in the context of spatial origin was recognized in . For a given set of diffusion parameters, these lengths can be computed, as well as the map of the probability that a particle reaching Earth was emitted from each source located in the disc. The lengths $`r_{\mathrm{wind}}`$ and $`r_{\mathrm{spal}}`$ depend on $`K`$, and are larger at higher energy. This is of course not the case for $`r_L=L`$, so that at sufficiently high energy, $`r_{\mathrm{wind}}L`$ and $`r_{\mathrm{spal}}L`$, the spatial origin being solely dictated by the halo size $`L`$ (side boundary $`r=R`$ is unimportant). The purely diffusive case is illustrated in Fig. 12, which shows the regions from which a given fraction of the cosmic rays reaching the Earth were emitted from. We refer the reader to  for the details.

Conversely, it turns out that when realistic diffusion parameters are considered, these lengths may be rather small, in particular for heavy nuclei which are more sensitive to spallations, so that the observation of cosmic rays in the Solar neighborhood only gives information about the local conditions under which diffusion occurs. This will be discussed in more details in Sec. 6.2.

## 5 Experimental determination of the diffusion parameters

We wish to study the constraints on the diffusion parameters (the spectral index of sources $`\alpha `$, the normalization $`K_0`$ and the spectral index $`\delta `$ of the diffusion coefficient, the height of the diffusive halo $`L`$, the galactic convective wind speed $`V_c`$ and the Alfvénic speed $`V_a`$) that come from the measured fluxes of cosmic rays. For the aim of our analysis, we can consider different classes of flux ratios: primary-to-primary (e.g. C/O), secondary-to-primary (e.g. B/C or sub-Fe/Fe), secondary-to-secondary (e.g. Li/B or Be/B), ratios of either stable (e.g. 10B/11B) or unstable (e.g. 10Be/9Be) isotopes. Each of these may be an indicator of some dominant physical phenomenon and be particularly sensitive to the corresponding diffusion parameters (see also Sec. 2.7). The ratio of two primaries is practically insensitive to changes in all the parameters, since they have the same origin and undergo the same physical processes (but keep in mind that extreme nuclei such as p and Fe present drastically different destruction rates, see also Sec. 6.4). Similar conclusions, even if less strong, may be drawn for the ratio of two isotopes of the same species, such as 10B/11B. Indeed, at very low energy values this quantity is slightly affected by changes in the injection spectra, but the effect is too weak to constrain the diffusion parameters.

One of the most sensitive quantity is B/C, as B is purely secondary and its main progenitors C and O are primaries. The shape of this ratio is seriously modified by changes in the propagation coefficients. Moreover, it is also the quantity measured with the best accuracy, so that it is ideal to test models. Indeed, as a ratio of two nuclei with similar Z, it is less sensitive to systematic errors and to Solar modulation than single fluxes or other ratios of nuclei with more distant charges. For the same reasons, the sub-Fe/Fe may also be useful. Unfortunately, since existing data are still affected by sizeable experimental errors, we can only use them to cross-check the validity of B/C but not to further constrain the parameters under scrutiny.

### 5.1 The method

For a given set of parameters, the source abundances of all nuclei (i.e. primaries and mixed nuclei) are adjusted so that the propagated top of atmosphere fluxes agree with the data at 10.6 GeV/nuc (see ). The nuclear cascade is started from Sulfur, as we checked that the heavier nuclei do not contribute significantly to the B/C ratio.

The top of atmosphere fluxes are deduced from interstellar fluxes using the force field modulation scheme (see  and references therein). The resulting B/C spectrum is then compared to the data (see below) and the $`\chi ^2`$ is computed for the chosen set of parameters, in our case

$$\chi ^2=\underset{i}{}\frac{\left((\mathrm{B}/\mathrm{C})_{i,\mathrm{exp}}(\mathrm{B}/\mathrm{C})_{i,\mathrm{model}}\right)^2}{\sigma _{i,\mathrm{exp}}^2},$$
(45)
where the sum runs over 26 experimental values from heao-3 . In general, if the experimental set-up is such that the measured (experimental) values differ from the “real” values by a non biased quantity with a given probability distribution, then the $`\chi ^2`$ value gives a quantitative estimate of the probability that the model is appropriate to describe the data. However, this condition is probably not fulfilled for heao-3, as for some measured quantity, the quoted errors $`\sigma _{i,\mathrm{exp}}^2`$ are much smaller (e.g. oxygen fluxes) or much larger (e.g. sub-Fe/Fe ratio) than the dispersion of data themselves. For this reason, it is meaningless to associate a likelihood to given $`\chi ^2`$ values. As a consequence, in an ideal situation in which very good and consistent data on B/C and sub-Fe/Fe ratios were available, the best attitude would be to make a statistical analysis of the combined set of data. Unfortunately, this is not currently the case. One can follow two ways to extract information from the sub-Fe/Fe data. First, as a check, one can compare the sub-Fe/Fe ratio predicted by our model – using the parameters derived from our above B/C analysis – with data from the same experiment. Second, one can search directly the minimum $`\chi _{\mathrm{Fe}}^2`$ of the sub-Fe/Fe ratio, with no prior coming from B/C. This procedure is more hazardous since the statistical significance of the sub-Fe/Fe data is far from clear.

### 5.2 The available data sets

#### 5.2.1 B/C and sub-Fe/Fe

As emphasized above, the results presented here are mostly based on the data taken by heao-3 . They have been taken in 1979-80, for elements with charges from 4 to 28 and for energies ranging from 0.6 to 35 GeV/nuc around a minimal Solar activity. In the case of B/C, the quoted 1-$`\sigma `$ heao-3 relative errors are 2-3%. We also considered data from balloons  and from the isee-3 experiment , even if the relevant error bars are wider. The first one collected data in 1973-75 for energies spanning from around 1.7 to 7 GeV/nuc. The second experiment was operating during 1979-81 on board a spacecraft, in the energy range 100-200 MeV/nuc. In some of the figures presented below, we also plot – for purely illustrative goals – the data point from imp-8  and the voyager experiments  (we did not include ulysses data point , since it corresponds to a period of maximal Solar activity). Nevertheless, we have checked that the best $`\chi ^2`$ values were not significantly modified when these points were added.

As regards the sub-Fe(Sc+V+Ti)/Fe ratio, we used data from heao-3  and from balloons . In both cases the error bars, around 10%, are significantly larger than for B/C.

#### 5.2.2 Radioactive species

Several experiments in the last twenty-thirty years have measured radioactive isotopes in cosmic rays with increasing precision, at energies of a few hundreds of MeV/nuc. The first data – usually presented as the ratio of some radioactive isotope to its stable companion(s) – were affected by errors of around 25-30%. The latest published data have error bars reduced by a factor of two or three. In the following, we implicitly refer to three satellite experiments, namely voyager, ulysses and ace. Other experiments will sometimes be shown on figures but they will be purely illustrative since their accuracy is far smaller.

The best measured ratio is probably 10Be/9Be which corresponds to the lowest Z $`\beta `$-radioactive nucleus. Data from ulysses  and from ace  are consistent, the quoted error bars being smaller for ace. They are also consistent with the voyager data point  for which the quoted error is larger. We do not use the smili data, as the possibility that they are plagued by statistical fluctuations is not ruled out .

As regards the radioactive chlorine isotope 36Cl, results are usually provided as 36Cl to total Cl ratio. The only available data, to our knowledge, are those from ulysses  with a 1-$`\sigma `$ error of about 35%, and ace  whose errors (even taken at 3-$`\sigma `$) are completely included in the ulysses 1-$`\sigma `$ upper error band.

Finally, the measurement of the 26Al/27Al ratio is more problematic. Indeed, the data from ulysses  and ace  do not seem to be compatible (the ace central point is much lower than ulysses’ one). Even enlarging ace error bars (which are smaller than ulysses) to 3-$`\sigma `$ does not improve significantly the compatibility. On the other side, the ulysses data are fully compatible with 1-$`\sigma `$ voyager  ones, whose uncertainty is still much greater than for the other two experiments. The possible discrepancy between some of these data is addressed in .

#### 5.2.3 p and He fluxes

We do not study here the compatibility of p and He fluxes with the prediction from diffusion models. However, we need the observed values of these fluxes to compute the secondary antiproton fluxes. The contribution of heavier nuclei to the antiproton production is negligible. Until recently, the spectra of p and He were known with a modest accuracy and the data from different experiments were often incompatible at high energy. This induced an uncertainty of some tens of percents in the predicted antiproton spectrum. Recent measurements made by the balloon-borne spectrometer bess  and by the ams detector during the space shuttle flight  dramatically reduced the uncertainties both on proton and helium spectra. We fitted the high energy (T $`>`$ 20 GeV/nuc) part of these measured spectra with the power law:

$$\mathrm{\Phi }(\mathrm{T})=\mathrm{N}(\mathrm{T}/\text{GeV/nuc})^\gamma ,$$
(46)
where the kinetic energy per nucleon T is given in units of GeV/nuc and the normalization factor N in units of $`\text{m}^2\text{s}^1\text{sr}^1\text{(GeV/nuc)}^1`$. This provides a good description down to the threshold energy for the antiproton production.

We fitted the bess and ams data both separately and combined, obtaining very similar results. This is obvious since the data from the two experiments are now totally compatible. The best fit corresponds to $`N=13249\text{m}^2\text{s}^1\text{sr}^1\text{(GeV/nuc)}^1`$ and $`\gamma =2.72`$; spectra obtained from the best fits on the single bess and ams data completely overlap. We did the same for helium and the corresponding numbers are $`N=721\text{m}^2\text{s}^1\text{sr}^1\text{(GeV/nuc)}^1`$ and $`\gamma =2.74`$. The 1-$`\sigma `$ deviation from the best fit spectrum does not exceed 1% for both species. The uncertainty induced on the antiproton spectrum is smaller than the ones induced by the uncertainties on nuclear physics and diffusion parameters. The situation has significantly improved since 1998, when an error of $`\pm `$ 25% was quoted .

### 5.3 The diffusion parameters estimated from heao-3 B/C measurements

#### 5.3.1 Partial exploration of the parameter space

In a first step , we did not consider $`\alpha _j`$, the slope of the source spectrum of each species $`j`$ (see Sec. 2.4.3, p. 2.4.3), as a free parameter but instead we set it to the value obtained by subtracting $`\delta `$ from the spectra measured at high energy (see ). We display models that have $`\chi ^2<40`$.

##### Parameters are strongly degenerated

Many configurations for $`L`$, $`K_0`$, $`\delta `$, $`V_c`$ and $`V_a`$ lead to the same B/C ratio. This can be understood as meaning that they all lead to the same effective description in terms of the associated $`\lambda _{\mathrm{esc}}^{\mathrm{LB}}`$ and grammage (see discussion in Sec. 2.5.3). In Fig. 14, each value of $`\delta `$ gives a different contour plot in the $`K_0/LL`$ plane. It appears that they all can be superimposed to a single curve by a rescaling $`K_0/LK_0/L\times f(\delta )`$, where $`f`$ is a function of $`\delta `$ only. For the contours displayed in Fig. 14, it takes the values $`f(0.46)=0.51`$, $`f(0.5)=0.62`$, $`f(0.6)1`$, $`f(0.7)=1.54`$ and $`f(0.85)=2.78`$. In the $`V_a/\sqrt{K_0}V_c`$ plane (Fig. 15), the values of $`V_c`$ are shifted downward as $`\delta `$ is decreased but the allowed range of $`V_a/\sqrt{K_0}`$ does not significantly move. In particular, with the above criterion, we do not find any model having a good $`\chi ^2`$ without convection ($`V_c=0`$) or without reacceleration ($`V_a=0`$). We emphasize the fact that when $`\alpha `$ and $`\delta `$ are fixed, the four remaining parameters are strongly correlated so that, when one of them is determined, in principle the allowed ranges for the others will be narrower than what could naively appear from the figures.

##### Partial conclusion

Several comments about these results can be made. First, the diffusion slope is constrained and the particular value $`\delta =1/3`$ is excluded by our analysis (even if a less stringent criterion on $`\chi ^2`$ is adopted). We find that the power law index for the diffusion coefficient is restricted to the interval $`0.450.85`$, the best $`\chi ^2`$ being $`25.5`$ for $`\delta =0.70`$, leading to $`\alpha 2.1`$. For any $`\delta `$ in this interval, the good parameters in the $`K_0/LL`$ and $`V_a/\sqrt{K_0}V_c`$ planes can be straightforwardly deduced from the corresponding values for $`\delta =0.6`$ by a simple scaling law. We also exclude any model without a convective velocity or without reacceleration for any combination of the three other diffusion parameters. These conclusions could get more stringent by new measurements in the whole energy spectrum for all nuclei, as all our results were obtained using the best data, which are rather scarce and more than 20 year-old.

#### 5.3.2 Full exploration of the parameter space

This study was then refined  by adding four updates: i) a more conventional reacceleration term is used (i.e. similar to terms used by other authors); ii) effect of deviation from a pure power law of source spectra at low energy is inspected – see Sec. 2.4.3, p. 2.4.3 –; iii) different diffusion schemes (slab Alfvén wave turbulence, isotropic fast magnetosonic turbulence and a mixture of the two case) are tested and, iv) above all, the source slope $`\alpha _j\alpha `$ for all species is set as a free parameter. The net result, is that all the above trends are confirmed. As a bonus, one recovers that for small diffusion slope, convection is unfavored, as found in .

##### Summary of the trends

In Fig. 16 we show the preferred values of the three diffusion parameters $`K_0`$, $`V_c`$ and $`V_a`$, for each best $`\chi ^2`$ in the $`\delta \gamma `$ plane, when $`L`$ has been fixed to 6 kpc (the behavior does not particularly depend on $`L`$).

The two upper panels show that the evolution of $`\alpha `$ does not affect $`K_0`$. On the other hand, we clearly see the (anti)correlation between the two parameters $`K_0`$ and $`\delta `$ entering the diffusion coefficient formula, because they need to give about the same normalization at high energy ($`K_0\times E_{\mathrm{thresh}}^\delta cte`$). Almost the same numbers are obtained for the pure power law and for the modified source spectra (these were defined in Sec. 2.4.3). $`K_0`$ spans the range 0.003 to 0.1 kpc2 Myr-1. The middle panels show the values for the convective velocity. Only very few configurations include $`V_c=0`$ km s-1, always when $`\delta =0.3`$, for both types of source spectra. Increasing $`\gamma `$ and $`\delta `$ at the same time makes $`V_c`$ change its trend. Finally, the Alfvén velocity $`V_a`$ doubles from $`\delta =1.0`$ to 0.3, whereas it is almost unchanged by a variation in the parameter $`\gamma `$ (or equivalently $`\alpha `$). B/C and sub-Fe ratios are displayed for several configurations in Fig.17.

The three parameters $`K_0`$, $`V_c`$ and $`V_a`$ behave very similarly with respect to a change in the source spectrum from “pure power law” to “modified”. It can be explained as the influence on the primary and secondary fluxes which can be factored out if energy changes are discarded (their effect is actually small on the derived parameters). The currently available data on B/C do not allow to discriminate clearly between these two shapes for the acceleration spectrum. This could be achieved by means of better data not only for B/C but also for primary nuclei (Donato et al., in preparation). We finally notice that pursuing the analysis with the sub-Fe/Fe ratio does not allow to go farther for reasons invoked above.

##### Conclusions

Thus, forgetting for a while some of our theoretical a priori about the diffusion power spectrum, a new picture of cosmic ray propagation seems to emerge, in which high values for the diffusion coefficient spectral index ($`\delta 0.60.7`$444It is noticeable that Strong et al. that used and preferred $`\delta =1/3`$ for the five past years took $`\delta =0.420.52`$ in their recent study (astroph/0210480). and source spectral indices $`\alpha 2.0`$ are favored. The conclusions of the full analysis can be summarized as follows: (i) the values $`\delta 0.70.9`$ and $`\alpha 2.0`$ are preferred; (ii) this preference holds whatever the specific form of the spectrum at low energy; the numerical values of the other parameters are also only slightly modified by this low energy dependence even though deviation from a power-law at low energy is preferred. The study of fluxes should give a more definite answer; (iii) $`K_0`$ scales logarithmically with $`\delta `$ and models with small halos tend to one-dimensional models with a simple relation between the surface mass density $`\mu `$, $`K_0`$, $`L`$ and $`V_c`$; (iv) several existing models are compared and the qualitative and quantitative differences between them are studied and partially explained (the reader is referred to  for a deeper discussion).

## 6 Diffusion parameters applied to “standard” CR physics

Armed with our set of parameters derived from the B/C ratio, it seems reasonable to suppose that other charged nuclei follow the same propagation history; this is already the case for sub-Fe/Fe nuclei as shown in Fig. 17. We therefore use the same diffusion model, with these sets of parameters, to study the propagation of all species. In this section we focus on the subjects pertaining to the astrophysics field. We start with secondary $`\overline{p}`$  (Sec. 6.1) that is a further cross-check of the validity of the propagation parameter derived. This nucleus is also particularly important to possibly constrain exotic primary sources in the diffusive halo  – this will be presented in Sec. 7 –. We devote Sec. 6.2 to the question of the spatial origin of all these stable nuclei . Then we focus on the study of radioactive species  in Sec. 6.3, where some clues about the local properties of the interstellar medium will be given. Finally, in Sec. 6.4 we study some consequences of our model on the composition at higher energies .

It should be reminded that a strong degeneracy between the propagation parameters was pointed out in Sec. 5.3. Consequently, much of the following work has to deal with inspecting carefully the consequences of this degeneracy, and if possible, to break it.

### 6.1 Secondary antiprotons and antideuterons

The secondary antiprotons and antideuterons are yielded by the spallation of cosmic ray nuclei over the interstellar medium, mostly protons and helium whose spectra $`\mathrm{\Phi }_\mathrm{p}`$ and $`\mathrm{\Phi }_{\mathrm{He}}`$ are by now well determined  as discussed in Sec. 5.2.3. The production cross sections have been discussed in Sec. 2.4.6.

#### 6.1.1 Secondary antiprotons

Unlike the secondary nuclei that are produced at fixed energy per nucleon, the source term for secondary antiprotons is obtained from the convolution of the energy spectra of the incident cosmic rays with the relevant differential production cross sections. In the case of impinging protons, this leads to

$`q_{\overline{p}}^{\mathrm{sec}}(r,E_{\overline{p}})=4\pi {\displaystyle _{E_{\mathrm{th}}}^+\mathrm{}}𝑑E_\mathrm{p}\mathrm{\Phi }_\mathrm{p}(r,E_\mathrm{p})`$
$`\times \left\{n_\mathrm{H}{\displaystyle \frac{d\sigma _{\mathrm{pH}\overline{\mathrm{p}}}}{dE_{\overline{p}}}}+n_{\mathrm{He}}{\displaystyle \frac{d\sigma _{\mathrm{pHe}\overline{\mathrm{p}}}}{dE_{\overline{p}}}}\right\}\left(E_\mathrm{p}E_{\overline{p}}\right).`$
A similar term arises from the cosmic ray helium radiation.

##### Tertiary production

Once they have been created, antiprotons may interact with the interstellar material in three different ways. First, they may undergo elastic scatterings on galactic hydrogen. The cross section for that reaction has been shown to peak in the forward direction  so that the corresponding antiproton energy loss is negligible. Antiprotons are not perturbed by these elastic scatterings as they survive them while their energy does not change. Then, they may also annihilate on interstellar protons. This process dominates at low energy, and its cross section $`\sigma _{\mathrm{ann}}^{\overline{p}p}`$ is given in . Last but not least, antiprotons may survive inelastic scatterings where the target proton is excited to a resonance. Antiprotons do not annihilate in that case but lose a significant amount of their kinetic energy. The cross section for these inelastic yet non-annihilating interactions is

$$\sigma _{\mathrm{non}\mathrm{ann}}^{\overline{p}p}=\sigma _{\mathrm{inel}}^{\overline{p}p}\sigma _{\mathrm{ann}}^{\overline{p}p}.$$
(47)
The energy distribution of antiprotons that have undergone such reactions has not been measured. It may be assumed to be similar to the proton energy distribution after p-p inelastic scattering. An impinging antiproton with kinetic energy $`T_{\overline{\mathrm{p}}}^{}`$ has then a differential probability

$$\frac{dN_{\overline{p}}}{dE_{\overline{p}}}=\frac{1}{T_{\overline{\mathrm{p}}}^{}}$$
(48)
to end up with the final energy $`E_{\overline{p}}`$, hence the differential cross section

$$\frac{d\sigma _{\overline{p}H\overline{p}X}}{dE_{\overline{p}}}=\frac{\sigma _{\mathrm{non}\mathrm{ann}}^{\overline{p}p}}{T_{\overline{p}}^{}}.$$
(49)
The corresponding source term for these so-called tertiary antiprotons may be expressed as

$`q_{\overline{p}}^{\mathrm{ter}}(r,E)=4\pi n_\mathrm{H}\sigma _{\mathrm{non}\mathrm{ann}}^{\overline{p}p}(E)\mathrm{\Phi }_{\overline{p}}(r,E)`$
$`+4\pi n_\mathrm{H}{\displaystyle _E^{\mathrm{}}}{\displaystyle \frac{\sigma _{\mathrm{non}\mathrm{ann}}^{\overline{p}p}(E^{})}{T^{}}}\mathrm{\Phi }_{\overline{p}}(r,E^{})𝑑E^{}.`$
The integral over the antiproton energy $`E`$ of $`q_{\overline{p}}^{\mathrm{ter}}(E)`$ vanishes. This mechanism does not actually create new antiprotons. It merely redistributes them towards lower energies and tends therefore to flatten their spectrum. Notice in that respect that the secondary antiproton spectrum that results from the interaction of cosmic ray protons impinging on interstellar helium is already fairly flat below a few GeV. Since it contributes a small fraction to the final result, the effect under scrutiny here may not be as large as previously thought . Helium is taken into account by replacing the hydrogen density in relation (6.1.1) by a geometrically-inspired factor of $`n_\mathrm{H}+4^{2/3}n_{\mathrm{He}}`$.

##### Conclusions

First, the values of all the inputs being either extracted from the analysis of nuclei (diffusion parameters $`\delta `$, $`L`$, $`K_0`$, $`V_c`$ and $`V_a`$) or measured (proton and helium fluxes), all the cosmic antiproton fluxes naturally coming out of the calculation are completely contained within the experimental error bars of bess data (see Fig. 18).

The major uncertainties come from nuclear physics and are already comparable to observational error bars. In Fig. 19, the two dotted lines feature the uncertainties related to nuclear physics. The upper curve is obtained with the set of maximal p-He, He-p, He-He cross sections while increasing the p-p cross section – as given by the Tan & Ng parameterization  – by a generous 10%. Similarly, the lower curve is obtained with the minimal values for those cross sections while decreasing the p-p cross section by 10%. Indeed, such a variation for p-p has been included for the sake of completeness even if it modifies the antiproton spectrum only by a few percents. The uncertainties which the upper and lower dotted curves exhibit are of the order of 22-25 % over the energy range 0.1-100 GeV.

Another source of uncertainty comes from the fact that the propagation parameters are degenerate, even if they are severely constrained by the analysis of B/C experimental results. This induces an error which is also displayed in Fig. 19 (solid lines).

As antiproton spectrum measurements should better in the near future, antiproton studies could be limited by nuclear indeterminacies. Further work and especially new measurements of antiproton production in the p-He channel would be of great interest.

#### 6.1.2 Secondary antideuterons

The production of secondary antideuterons also proceeds through the collision of cosmic ray high-energy particles – mostly protons and helium nuclei – on the atoms of the interstellar medium. The source term $`q_{\overline{D}}^{\mathrm{sec}}(r,E_{\overline{D}})`$ for secondary antideuterons can then be readily obtained, using the coalescence scheme described in Sec. 2.4.7.

The result is that the spectrum of secondary antideuterons sharply drops below a few GeV/nuc (see Fig. 31, p. 31), mostly because in the galactic frame, the production threshold is 17 $`m_\mathrm{p}`$. Energy losses tend to shift slightly the antideuteron spectrum towards lower energies with the effect of replenishing the low-energy tail with the more abundant species which, initially, had higher energies. That effect is nevertheless mild. In particular, there is no tertiary contribution $`q_{\overline{D}}^{\mathrm{ter}}(r,E_{\overline{D}})`$. An antideuteron that undergoes an inelastic scattering is likely to be destroyed unlike antiprotons for which inelastic yet non-annihilating interactions are possible (see Sec. 6.1.1).

As a consequence, the flux of secondary antideuterons is quite suppressed with respect to antiprotons. However, a preliminary study, which does not take advantage of the full set of propagation parameters allowed by B/C shows that a dozen events may be collected by the ams collaboration during its space station stage, in the energy range extending up to 100 GeV/nuc. These antideuterons correspond to IS energies in excess of $`3`$ GeV/nuc, a region free from the effects of Solar modulation. This result takes into account geomagnetic suppression as discussed in .

### 6.2 Spatial origin

The question of the spatial origin of cosmic rays in our stationary diffusion model was addressed in Sec. 4. It was shown that the origin is fully described by three parameters $`L`$, $`r_{\mathrm{wind}}`$ and $`r_{\mathrm{spal}}`$ (see also Sec. 4.3), which depend on the values of the propagation parameters. Having determined a definite range of realistic parameters in Sec. 5, we can go one step further in characterizing the origin of cosmic rays.

##### Realistic values for $`r_{\mathrm{spal}}`$ and $`r_{\mathrm{wind}}`$

In Fig. 20 are plotted $`\chi _\mathrm{w}(\delta ,L)2L/r_{\mathrm{wind}}`$ and $`\chi _{\mathrm{spal}}(\delta ,L)L/r_{\mathrm{spal}}`$ for several energies and species. Actually, $`\chi _\mathrm{w}(\delta ,L)`$ is a function of rigidity; it is more clever to use the latter choice instead of kinetic energy per nucleus, since $`\chi _\mathrm{w}(\delta ,L)`$ then does not depend on the species. The left panel of Fig. 20 displays $`\chi _\mathrm{w}(\delta ,L)`$ for three rigidities: 1 GV, 10 GV and 100 GV. The quantity $`\chi _\mathrm{w}(\delta ,L)`$ is an indicator of the competing role of convection and diffusion to keep the remote cosmic rays to reach Earth. The figure shows that above several tens of GV, diffusion has the main role, whereas convection may dominate below this value, at least for large values of $`\delta `$. Different values of the halo size $`L`$ yield the same conclusions. The right panel of Fig. 20 displays $`\chi _{\mathrm{spal}}(\delta ,L)`$ which is related to the spallation efficiency. It appears that at low energy, heavy nuclei such as Fe are preferentially destructed rather than swept out by the convective wind or by the escape through the boundaries.

##### Origin of primaries and secondaries

With the knowledge of $`r_{\mathrm{wind}}`$, $`r_{\mathrm{spal}}`$ and $`L`$ at 1 GeV/nuc for all species one can now answer the question of the origin of CR: results for p, CNO and Fe nuclei are shown in Fig. 21. The peculiar behavior of secondaries is fully explained in .

To summarize the situation, most cosmic rays detected on Earth were emitted from a limited zone of the disc. This is even more true for heavy species such as Fe, which are very sensitive to spallations and thus are unlikely to travel long distances. This implies that the information on the diffusive processes inferred from the study of the ratio sub-Fe/Fe are only valid on a very local region, all the more that the diffusion slope $`\delta `$ is large. Even if the propagation conditions were very different outside of this region, the observations made in the Solar neighborhood would almost not be affected, pointing preferentially once more towards large $`\delta `$ values. Otherwise, as several anomalies in some CR radiations indicate, CR fluxes could be variable with location, making all stationary models meaningless.

### 6.3 The radioactive species

The radioactive cosmic ray nuclei deserve a specific treatment. From their creation, they diffuse on a typical distance $`l_{\mathrm{rad}}\sqrt{K\gamma \tau _0}`$ before decaying. In this expression, not only the diffusion coefficient $`K`$, but also the lifetime $`\gamma \tau _0`$, depend on energy, due to the relativistic time stretch. The following table gives some values of this distance for three species (typical values $`K_0=0.033\text{kpc}^2\text{Myr}^1`$ and $`\delta =0.6`$ have been assumed in Tab. 6).

These nuclei are therefore very sensitive to the presence of the local bubble which has an extent of about $`50200`$ pc (see Sec. 2.4.1). As regards the diffusion process itself, as described by the coefficient $`K(E)`$, radio and $`\gamma `$-ray observations, which can test in situ the spectrum and density of cosmic rays , indicate that it is not affected by the presence of the bubble, i.e. diffusion is homogeneous.

The bubble has nevertheless an effect on the propagation. First it leads to a decrease in the spallation source term of the radioactive species. Second, it also leads to a local decrease of destructive spallations. Third, as there is less interstellar matter to interact with, the energy losses are also lowered. Because the typical propagation time of radioactive nuclei is short, the energy redistributions are negligible, as well as destruction, and the first effect is dominant. As mentioned above, the bubble is modelled as a hole, in the thin disc approximation (see Fig. 11, p. 11). The radius of this hole is considered as an unknown parameter in the analysis.

Fig. 22 shows, for $`r_{\mathrm{hole}}=200`$ pc, how the fluxes both for stable and radioactive species are affected in the neighborhood of the hole. The stables remain grossly unaffected by its presence, whereas radioactive are strongly suppressed for the reason just mentioned.

##### Effect and size of the hole

The major result we found  is that at the center of the bubble, the radioactive fluxes are decreased by a factor which can be approximated as

$$\frac{N^{r_{\mathrm{hole}}}}{N^{r_{\mathrm{hole}}=0}}\mathrm{exp}(r_{\mathrm{hole}}/l_{\mathrm{rad}}).$$

When all the sets of diffusion parameters allowed by the B/C data are used to compute the 10Be/9Be, 36Cl/Cl ratios, we find that each of the radioactive nuclei independently points towards a bubble of radius $`100\text{pc}`$, in relatively good agreement with direct observations. If these nuclei are considered simultaneously, only models with a bubble radius $`r_{\mathrm{hole}}60100\text{pc}`$ are consistent with the data. In particular, the standard case $`r_{\mathrm{hole}}=0\text{pc}`$ is disfavored. This is shown in Fig. 23, which is a projection of the parameter subspace allowed by B/C, 10Be/9Be and 36Cl/Cl on the $`L\delta `$ plane (left panel, no hole) or $`r_{\mathrm{hole}}\delta `$ plane (right panel, hole $`r_{\mathrm{hole}}`$). When the ratio 27Al/27Al is considered along with the two others mentioned above, the results become less clear , and it is suspected that the data (nuclear or astrophysical) on which they rely should not be trusted.

##### Whole energy spectrum

As 10Be/9Be is probably the best measured radioactive ratio (see Sec. 5.2.2), it is interesting to see how ratios depend on energy, given the constraint that parameters must fit both B/C and 10Be/9Be data. This is shown in Fig. 24 for the homogeneous model and a hole $`r_{\mathrm{hole}}=80`$ pc.

These curves show that it would be very important to have accurate measurements of the energy spectra of radioactive species such as 10Be, 36Cl and 27Al.

##### Conclusions

Most studies use the radioactive nuclei to constrain the halo size $`L`$ (see discussion in Sec. 2.4.2). The meaning of the results are very much dependent on the the treatment of the LISM. In particular, the presence of the local bubble leads to an exponential attenuation of the radioactive fluxes. This is of special importance for short-lived nuclei such as 14C, which is attenuated by a factor $`1`$, unless a local source is present. This may explain the fact that ace did not detect any such nucleus. Finally, if one takes $`r_{\mathrm{hole}}`$ as an input fixed by observations, the radioactive progenitor/parent may give further constraints on the diffusion parameters.

### 6.4 PeV fluxes

Recent measurements of the cosmic ray average logarithmic mass and all-particle spectrum around $`10^{15}`$ eV  gave new clues to understand the origin of the cosmic rays and in particular the puzzle of the knee in the energy spectrum. The highest energy particles are almost certainly extragalactic. A similar origin is not excluded near the knee, but it is difficult to account for the observed softening of the spectrum in this region. As a consequence, the intermediate region between $`10^{15}`$ and $`10^{19}`$ eV should be analyzed in terms of the same physical mechanisms than lower energy particles.

At these energies, fluxes are too low to be directly measured ($``$ m-2 sr-1 yr-1). With present techniques, only two quantities can be extracted in large ground array detectors (e.g. ); namely all-particle spectrum and $`\mathrm{ln}A`$. They are given by linear combination of the individual fluxes with different weights;

$$\mathrm{\Phi }^{\mathrm{all}}\underset{j}{}\mathrm{\Phi }_j\text{and}\mathrm{ln}A\frac{_j\mathrm{ln}A_j\mathrm{\Phi }_j}{_j\mathrm{\Phi }_j}.$$
(51)
As the experimental data are given in total energy rather than in kinetic energy per nucleus, we adopt the same presentation for the results above 100 GeV.

#### 6.4.1 Separation of key ingredients

We first estimate the evolution of $`\mathrm{ln}A`$ in the Leaky Box frame. This allows to apprehend more easily some basic ingredients that drive the mass evolution, but not all of them. In particular, the geometrical effects cannot be considered. Trends are also more easily understood with only one light and one heavy nucleus (e.g. p and Fe). If energy gains and losses are discarded, the Leaky Box equation for primary species, such as p and Fe, reads

$$\frac{N^j}{\tau _{\mathrm{esc}}}+\overline{q}_0^jQ^j(E)\mathrm{\Gamma }^jN^j=0.$$
(52)
The proportion of light and heavy species evolves with energy because of two effects. The source spectra $`\overline{q}_0^jQ^j(E)`$ may differ, and the spallative destruction rate $`\mathrm{\Gamma }^j`$ are definitely not the same. The escape time $`\tau _{\mathrm{esc}}=\tau _0^\delta `$ is supposed to be independent of species.

##### Source spectrum effect

The first effect is studied by neglecting the spallation term in Eq. (52). We suppose that the two species p and Fe are injected with different slopes related by $`\alpha _{\mathrm{Fe}}\alpha _\mathrm{p}=0.1`$ . If we start from the experimental value $`\left(\mathrm{Fe}/\mathrm{p}\right)_{100\mathrm{GeV}}=1/20`$ (e.g. ), Eq. (51) gives the following evolution, independently of $`\delta `$,

$`\mathrm{ln}A(E=100\mathrm{GeV}/\mathrm{\hspace{0.25em}10}\mathrm{TeV}/\mathrm{\hspace{0.25em}1}\mathrm{PeV}/\mathrm{})`$
$`0.19/\mathrm{\hspace{0.25em}0.30}/\mathrm{\hspace{0.25em}0.50}/\mathrm{\hspace{0.25em}4.03}.`$
##### Escape plus selective destruction

If the spallations are taken into account, Eq. (51) can be written as

$`\mathrm{ln}A(E)=\mathrm{ln}(A_{\mathrm{Fe}})\times [1+(q_\mathrm{p}^0/q_{\mathrm{Fe}}^0)\times \left({\displaystyle \frac{Z_\mathrm{p}}{Z_{\mathrm{Fe}}}}\right)^\alpha `$ (53)
$`\times {\displaystyle \frac{(\sigma _0(E/Z_{\mathrm{Fe}})^\delta +\sigma _{\mathrm{Fe}}^{\mathrm{inel}})}{(\sigma _0(E/Z_\mathrm{p})^\delta +\sigma _\mathrm{p}^{\mathrm{inel}})}}]^1.`$
where the cross sections are $`\sigma _\mathrm{p}^{\mathrm{inel}}30`$ mb and $`\sigma _{\mathrm{Fe}}^{\mathrm{inel}}710`$ mb). The quantity $`\tau _{\mathrm{esc}}`$ has been converted in mb, with a typical value $`\sigma _0=47`$ mb and $`\delta =0.61`$ . We are not interested here in the effect of injection spectra and the slopes have now been assumed to be the same for all species ($`\alpha _\mathrm{p}=\alpha _{\mathrm{Fe}}\alpha `$=2.2). With $`\mathrm{ln}A(100\mathrm{GeV})=0.19`$, this gives

$`\mathrm{ln}A(E=100\mathrm{GeV}/\mathrm{\hspace{0.25em}10}\mathrm{TeV}/\mathrm{\hspace{0.25em}1}\mathrm{PeV}/\mathrm{})`$
$`0.19/\mathrm{\hspace{0.25em}0.83}/\mathrm{\hspace{0.25em}1.05}/\mathrm{\hspace{0.25em}1.08}.`$
##### Geometrical source effect

As emphasized above, the study of the influence of geometry is not straightforward in Leaky Box models. From the previous discussion in the framework of diffusive propagation (see Sec. 4), we know that cosmic rays detected on Earth come from an average distance $`r(E)`$ which is larger at higher energy. As a metallicity gradient is present, at the level of $`0.05`$ dex kpc-1, i.e. $`[\mathrm{Fe}/\mathrm{H}]=0.05(rR_{})\mathrm{log}(q_{\mathrm{Fe}}/q_\mathrm{p})^{\mathrm{eff}}\mathrm{log}(\mathrm{Fe}/\mathrm{H})_{}`$ (see Sec. 2.4.3), the nuclei with different energies were emitted from sources having different properties. This leads to (the two previous effects are ignored)

$$\mathrm{ln}A(E)=\frac{\mathrm{ln}(A_{\mathrm{Fe}})}{1+\left(q_\mathrm{p}/q_{\mathrm{Fe}}\right)^{\mathrm{eff}}};$$

where

$$\left(\frac{q_\mathrm{p}}{q_{\mathrm{Fe}}}\right)^{\mathrm{eff}}=\left(\frac{q_\mathrm{p}}{q_{\mathrm{Fe}}}\right)_{}\times 10^{0.05(rR_{})}.$$

An upper bound is obtained for this effect by assuming that $`r(100\mathrm{GeV})=R_{}`$ and $`r(\mathrm{})=0`$ (galactic center), which yields

$$\mathrm{ln}A(E=100\mathrm{GeV}/\mathrm{})0.19/\mathrm{\hspace{0.25em}0.54}.$$
(54)
Even in this optimal case, this effect is seen to be negligible compared to the two others (see below).

#### 6.4.2 Results from our propagation model

The quantity $`\mathrm{ln}A`$ can also be computed with our propagation model, considering the different sets of propagation parameters derived from the B/C analysis (see Sec. 5.3).

##### Geometrical effects are subdominant

Fig. 25 shows spallations effects and geometrical effects. The first strong conclusion is that spallation effects affect dramatically the composition of cosmic rays.

Furthermore, at sufficiently high energy, as expected with the Leaky Box toy model, the asymptotic regime is reached and propagation ceases to affect $`\mathrm{ln}A`$. The geometrical effects (source distribution plus metallicity gradient) are less important than the others. They induce a change of at most 5% in the results. This would be even more true for lower $`L`$. For $`L=3`$ kpc, these geometrical effects are completely negligible.

##### Evolution above PeV energies

The evolution of $`\mathrm{ln}A`$ obtained with different sets of parameters compatible with B/C is presented in Fig. 26.

The right panel focuses on the source effect and shows that $`\alpha _\mathrm{p}>\alpha `$ is preferred. As suggested in our toy model, only source effects enable evolution of $`\mathrm{ln}A`$ above 1 PeV. It can also be seen that large diffusion slopes $`\delta `$ are preferred, as was already hinted in Sec. 5.3.

##### Inclusion of the knee

An important experimental fact, referred to as the knee, has not be taken into account in the analysis presented above. It is a drop of the slope of $`\mathrm{\Delta }\gamma 0.4`$ in the observed spectrum of cosmic rays at an energy of a few PeV. This may be due to a change in the source spectra or in the diffusion coefficient. We now want to study the effect of the knee on the evolution of $`\mathrm{ln}A`$. To this aim, we propose to model the knee in a very naive way, by a break in the source spectral indices at an energy $`E_{\mathrm{knee}}`$ that depends on the species (notice that a change in diffusive regime, e.g. as recently re-inspected in , is also possible). We consider two cases, a break at a given rigidity $`E_{\mathrm{knee}}=Z\times 4`$ PeV and a break at given energy per nucleon $`E_{\mathrm{knee}}=A\times 4`$ PeV. The resulting average mass $`\mathrm{ln}A`$ is shown in Fig. 27 (the all particle flux $`\mathrm{\Phi }^{\mathrm{all}}`$ can be found in ). Several features are visible. They are first generated by protons changing their spectral index, followed by others: the helium component, the CNO group, and Fe. At high energy, all species have the same post-knee spectral index. The transition is even more gradual when secondary species (upper curve) are taken into account (they also change the normalization). The choice of the break (in rigidity or in energy per nucleon) affects mostly the all-particle spectrum. The conclusion at this stage is that the scatter in the data is too large to provide a clear test of our theoretical predictions (see Fig. 27).

## 7 Propagation models applied to astroparticle physics

Connections between high-energy cosmic ray physics and cosmology are explored in these three sub-sections. We start by a general introduction and illustrate the above-mentioned ties by a discussion (i) on the limits on the cold invisible gas , (ii) on the baryon asymmetry in the universe (see e.g. ) and finally (iii) on the presence of putative neutralinos that could make up the astronomical dark matter  or primordial black holes that could have been formed in the early Universe. Then, in Sec. 7.2 we explore the question of the spatial origin of these exotic nuclei. Finally, Sec. 7.3 and Sec. 7.4 respectively deal with SUSY  and PBH  signatures in antiprotons and antideuterons.

### 7.1 Introduction

Ties are strong between high-energy cosmic ray physics and cosmology. In this section, we explore a few examples that illustrate these connections and show that a deep understanding on how cosmic rays propagate turns out to be crucial.

A long standing problem in astronomy lies in the existence of large amounts of unseen material whose gravitational effects are nevertheless large. To commence, the Milky Way is surrounded by an extended halo that induces a flat rotation curve in its plane . This trend has also been observed in many spiral systems. On larger scales, the presence of dark matter inside galactic clusters has been noticed since many decades . Finally, the observations of the Cosmic Microwave Background (CMB) anisotropies , combined either with the determination of the relation between the distance of luminosity and the redshift of type Ia supernovae , or with the large scale structure (LSS) information from galaxy and cluster surveys , give independent evidence for a cosmological dark matter density in the range $`\mathrm{\Omega }_{\mathrm{cdm}}h^2=0.13\pm 0.05`$ , to be compared to a baryon density of $`\mathrm{\Omega }_\mathrm{b}h^2=0.019\pm 0.002`$ as indicated by nucleosynthesis  and the relative heights of the first acoustic peaks in the CMB data, $`\mathrm{\Omega }_\mathrm{b}h^2=0.022\pm 0.003`$ . The nature of the astronomical dark matter is still unresolved and most of it could be made of non-baryonic species or primordial black holes. But observations also point towards the presence of dark baryons insofar as the amount of luminous material is significantly smaller than $`\mathrm{\Omega }_\mathrm{b}`$.

##### Limits on the Dark Gas inside the Milky Way

As a dark matter of fact, De Paolis et al.  have outlined a scenario in which dark clusters of compact objects pervade the halo of the Milky Way together with clouds of molecular hydrogen, at distances larger than 10 to 20 kpc. Pfenniger et al.  have also suggested that a flat rotation curve beyond the Solar circle could be explained by a thin disc of cold molecular hydrogen, widening at large distances from the center. If unseen gas was present in the galactic halo or in the ridge of the Milky Way, it would be impacted by cosmic rays originating from the disc. This would lead to a strong $`\gamma `$-ray signal showing up as a new component in the galactic diffuse radiation. By using the cosB results at high galactic latitude, Gilmore  has already inferred a conservative limit $``$ 15 % on the amount of hidden gas in the halo. That analysis can be significantly refined because a key ingredient – the propagation of cosmic rays – is now better understood. A reliable calculation of their density as a function of galactocentric radius $`r`$ and height $`z`$ above the plane is definitely possible, as discussed in the previous sections. It allows for the determination of the gamma ray hydrogen emissivity $`I_H`$ per hydrogen atom as a function of location

$$I_H(E_\gamma ,r,z)=\left\{\frac{N(r,z)}{N(r_{},0)}\right\}I_H(E_\gamma ,),$$

where the spectral index of cosmic ray protons – the dominant contribution – is found to remain roughly constant throughout the Galaxy. In the local range, the emissivity is

$`I_H(E_\gamma >100\text{MeV},)=`$
$`(1.84\pm 0.10)\times 10^{26}\mathrm{photons}\mathrm{H}^1\mathrm{s}^1\mathrm{sr}^1`$
The spallation of cosmic rays with the gas potentially concealed in the Milky Way produces an extra $`\gamma `$-ray diffuse emission whose flux obtained from the convolution along the line of sight of the density $`n_\mathrm{H}^{\mathrm{dark}}`$ of dark gas with the local $`\gamma `$-ray emissivity $`I_H`$ is

$$\mathrm{\Phi }_\gamma ^{\mathrm{extra}}(E_\gamma )=I_H(E_\gamma ,s)n_\mathrm{H}^{\mathrm{dark}}(s)𝑑s.$$
(55)
Observations collected by the egret instrument on board cgro have been compared to the $`\gamma `$-ray diffuse radiation that may be modeled with both HI and H2 contributions. An isotropic component is found in the residuals of the signal. Its integral flux between 500 MeV and 1 GeV is found to be $`2.5\times 10^6`$ photons cm-2 s-1 sr-1 in the direction of the Perseus arm (b = 0 , l = 115). By requiring that $`\mathrm{\Phi }_\gamma ^{\mathrm{extra}}`$ does not exceed that value, a constraint may be derived on the amount of cold molecular clouds – so far untraced by CO – which the Milky Way ridge could conceal. A Mestel’s disc with local surface mass density $`\mathrm{\Sigma }()`$ 200 $`\mathrm{M}_{}`$ pc-2 is assumed. The spatial density of that distribution is given by

$$\rho (r)=\left(\frac{\mathrm{\Sigma }()}{a_{\mathrm{disc}}}\right)\left(\frac{r_{}}{r}\right),$$
(56)
where $`a_{\mathrm{disc}}`$ denotes the disc thickness. Fig. 28 features the upper limit $`\eta _{\mathrm{max}}`$ on the fraction of cold gas which the outer parts of the disc ($`r>r_{\mathrm{min}}`$) may contain. In particular, for $`a_{\mathrm{disc}}=400`$ pc, the limit relaxes to 1 only if $`r_{\mathrm{min}}`$ exceeds 19.5 kpc. We therefore conclude that molecular hydrogen, sufficiently cold to have escaped detection and so far untraced in CO, should actually be sparse up to $``$ 20 kpc from the galactic center and that it cannot account for the rotation curve of the Milky Way.

##### Neutralino as a Dark Matter Candidate, PBH as a cosmological remnant

One of the favorite candidates to the astronomical missing mass is a neutral weakly interacting particle. Such a species is predicted in particular by supersymmetry, a theory that is actively tested at accelerators. There are indirect clues for supersymmetry, noticeably the existence of a single high-energy unification scale for supersymmetric grand-unified models alone. It is conceivable therefore that most of the dark matter in the halo of the Milky Way is made of such neutral particles. The relic abundance of these so-called neutralinos is relevant to cosmology. If present in the galactic halo, they should still annihilate mutually to yield, among a few other indirect signatures, a flux of antiprotons whose spectrum is discussed in Sec. 7.3.

The same line of argument applies if black holes make up a fraction of the galactic dark matter. Primordial black holes could have formed in the early universe and could have concealed baryons that are not accounted for by primordial nucleosynthesis since they would not have participated in the nuclear fusion of helium. The actual value for $`\mathrm{\Omega }_\mathrm{b}`$ could therefore be much larger than $``$ 5% and even reach $`\mathrm{\Omega }_{\mathrm{cdm}}`$, sparing us the need for non-baryonic particles. If the black holes are light, they could evaporate today and yield also antiprotons as is discussed in Sec. 7.4.

The primary antiprotons produced by neutralino annihilations or black hole evaporations will get mixed with a conventional population of secondary antiprotons that are produced by the spallation of cosmic ray nuclei on the interstellar gas of the Milky Way ridge. It is crucial to investigate in some detail this later mechanism in order to ascertain which of the primary or secondary species is dominant. This has been done in Sec. 6.1.

##### Antinuclei in the Galaxy

Another crucial problem of cosmology is the existence of a possible asymmetry between matter and antimatter. The current wisdom is that we live in a matter dominated universe and that antimatter islands seem to be excluded. The ongoing annihilations that would take place at the frontiers would create an intense gamma-ray emission that is not seen. The mechanism which would be required to account for a separation of matter and antimatter domains is not known. However, one should not subscribe to any dogmatic point of view. The amount of antimatter in cosmic rays is about to be measured with unequalled accuracy by the space station borne spectrometer of the ams collaboration . One of the most exciting goals of the experiment is the possible detection of antinuclei in the cosmic radiation. It is generally believed that the observation of a single anti-helium or anti-carbon would undoubtedly signal the presence of antistars because the conventional production of antinuclei through spallation is negligible. This point has been examined in Sec. 6.1 where we mostly concentrated on secondary antideuterons. Suffice it to say that the creation of any additional antinucleon during a single spallation event is suppressed by a factor of $`10^510^4`$ . As featured by Fig. 29, the $`\overline{D}/\overline{p}`$ ratio is $`3\times 10^5`$ while $`\overline{\mathrm{He}}/\overline{D}7\times 10^5`$ for a momentum per nucleon in excess of 10 GeV. Secondary antinuclei with atomic number $`A3`$ are so much suppressed that they are out of reach of the near future instruments. Another possibility would be that antiglobular clusters lies in the halo of the Galaxy and could make a substantial contribution to the anti-helium flux . There is more hope for antideuterons. Their primary production will be discussed in Sec. 7.3 for supersymmetric neutralinos and briefly mentioned in Sec. 7.4 in the case of evaporating black holes.

### 7.2 Spatial origin of the PBH or SUSY exotic primaries

As for species created in the disc (see Sec. 6.2), one can infer the closed surfaces from which most of the cosmic rays detected on Earth  were emitted from. We are now confronted with a volume distribution of sources (see Sec. 2.4.4). As for standard sources in the disc, geometrical effects are very important (see Sec. 4.3). Apart from these geometrical restrictions, spallations and convective wind greatly affect the shape of these isodensity surfaces and, e.g., location of the particularly interesting contour defined by $`𝒫(𝒱(x)|O)=99`$% (contour from where 99% of exotic cosmic ray originate). As a consequence, not only the characteristic size $`L`$ is important to estimate the $`x=99`$% surface, but also spallations and convection typical extension, i.e. parameters $`r_{\mathrm{spal}}`$ and $`r_{\mathrm{wind}}`$ such as defined in Eqs. (4.3) and (4.3). These parameters strongly depend on $`\delta `$ (the diffusion coefficient slope) when realistic configurations are retained.

##### Realistic values of $`r_{\mathrm{spal}}`$ and $`r_{\mathrm{wind}}`$ for antideuterons

Antideuteron signal seems to be the most promising species to consider. We choose the interstellar energy 1 GeV/nuc; nuclei that reach the detector are Solar modulated so that the final energy corresponds roughly to $`200800`$ MeV/nuc depending of the modulation parameter. This is the window where the exotic signal becomes interesting. Tab. 7 summarizes the values of $`r_{\mathrm{wind}}`$ and $`r_{\mathrm{spal}}`$ at this energy for antideuterons for three halo sizes and three values of $`\delta `$.

We notice that the situation is very different for small or large $`\delta `$. For Kolmogorov power spectrum, only spallations act on the propagation and only weakly for this light nucleus: for large $`\delta `$ – the value $`\delta =0.85`$ is the one preferred in our B/C analysis (see Sec. 5.3) – models are convection/spallation dominated with spallation and convection acting at about the same footing.

##### Results

The production of particles by PBH is related to the dark matter profile whose weight is $`w_{\mathrm{PBH}}(r,z)=f(r,z)`$, whereas for supersymmetric annihilating particles, the production is related to the weight $`w_{\mathrm{SUSY}}(r,z)=f(r,z)^2`$ (see Sec. 2.4.4). Weighting the elementary probability by this production rate leads to the contours displayed in Fig. 30 for $`𝒫_{\mathrm{cyl}}(𝒱(x)|O)=99`$% (parameters are taken from Tab. 7). For $`\delta =0.35`$ (external contours), we recover basically contours that one would obtain in the high energy limit, i.e. contours driven by geometrical effects, pure diffusion (see also ). However, for larger $`\delta `$ (internal contours), these contours shrink naturally and all surfaces are distorted towards the galactic center, where the maximum of production occurs.

One has also to keep in mind that whatever these surfaces, most cosmic rays emitted from inside these volumes actually do not reach us. Escape preferentially occurs and it is obvious that the closer the source, a smaller part of the emitted flux escapes. The fraction that reaches us can be estimated to be about 0.01%-0.1% of the total emission in volume defined by $`x=99`$%, whatever $`\delta `$ (it slightly depends on the halo size $`L`$ for PBHs as well as for SUSY particles). Another meaningful number is the fraction of particles created in the volume $`x=99`$% with respect to the total number of particles created in the entire dark halo. This quantity is strongly dependent on the diffusive halo size and, for $`L=10`$ kpc compared to $`L=2`$ kpc, numbers are 2.% vs 0.2% (PBH) and 80% vs 20% (SUSY). This is a mere consequence of the different production terms.

### 7.3 Primary antiprotons and antideuterons from supersymmetric sources

The neutralinos that could be concealed in the halo of the Milky Way – and be responsible for the flatness of its rotation curve – should be steadily annihilating and produce antiprotons together with antideuterons. A key difference with respect to the production mechanisms that have already been discussed lies in the fact that this new yet putative source of antinuclei is not confined to the galactic disc. It spreads all over the halo far above and beneath the ridge. Its contribution to the local flux may nevertheless be easily derived with the formalism presented in Sec. 3.2.3. The propagation of primary species from the remote regions of the Milky Way neighborhood to the Earth has actually been treated in .

##### SUSY model

The neutralino naturally appears in the framework of the Minimal Supersymmetric extension of the Standard Model (MSSM)  as the lowest-mass linear superposition of the photino $`\stackrel{~}{\gamma }`$, the zino $`\stackrel{~}{Z}`$ and the two higgsino $`\stackrel{~}{H}_1^{}`$ and $`\stackrel{~}{H}_2^{}`$ states

$$\chi a_1\stackrel{~}{\gamma }+a_2\stackrel{~}{Z}+a_3\stackrel{~}{H}_1^{}+a_4\stackrel{~}{H}_2^{}.$$
(57)
This particle is neutral and interacts weakly. It generally turns out to be the lightest supersymmetric state and is therefore stable if R-parity is conserved.

##### Antiproton multiplicity

The differential multiplicity for antiproton production in a neutralino pair annihilation may be expressed as

$$\frac{dN_{\overline{p}}}{dE_{\overline{p}}}=\underset{\mathrm{F},\mathrm{h}}{}B_{\chi \mathrm{h}}^{(\mathrm{F})}\frac{dN_{\overline{p}}^\mathrm{h}}{dE_{\overline{p}}}.$$

The annihilation proceeds – through the various final states F – towards the quark or the gluon h with the branching ratio $`B_{\chi \mathrm{h}}^{(\mathrm{F})}`$. Quarks or gluons may be directly produced. They may alternatively result from the intermediate production of a Higgs or a gauge boson as well as of a top quark. Each quark or gluon h generates in turn a jet whose subsequent fragmentation and hadronization yields the antiproton energy spectrum $`dN_{\overline{p}}^\mathrm{h}/dE_{\overline{p}}`$. The source term for supersymmetric antiprotons

$$q_{\overline{p}}^{\mathrm{susy}}(r,z,E_{\overline{p}})=\sigma _{\mathrm{ann}}v\frac{dN_{\overline{p}}}{dE_{\overline{p}}}\left\{\frac{\rho _\chi (r,z)}{m_\chi }\right\}^2$$

supplements the spallation contributions $`q_{\overline{p}}^{\mathrm{sec}}`$ and $`q_{\overline{p}}^{\mathrm{ter}}`$.

Because the annihilation of a neutralino pair $`\chi \chi `$ occurs at rest with respect to the Galaxy, most of the antiprotons – and antideuterons for that matter – are produced at low energies. The resulting spectrum is fairly flat below a few GeV. This has important observational consequences insofar as the secondary antiproton spectrum is already fairly flat as shown in Fig. 19, p. 19. It may therefore prove difficult to disentangle a primary component from the secondary antiproton radiation. The primary flux nevertheless modifies the magnitude – if not the spectrum – of the cosmic ray antiproton radiation and could contribute significantly as shown in Tab. 8. This effect can be used to constrain the supersymmetric parameter space.

##### Supersymmetric $`\overline{D}`$ signal

As regards the antideuteron production, the factorization-coalescence scheme discussed above leads to the antideuteron differential multiplicity

$`{\displaystyle \frac{dN_{\overline{D}}}{dE_{\overline{D}}}}={\displaystyle \frac{4P_{\mathrm{coal}}^3}{3k_{\overline{D}}}}{\displaystyle \frac{m_{\overline{D}}}{m_{\overline{p}}m_{\overline{n}}}}`$
$`\times {\displaystyle \underset{\mathrm{F},\mathrm{h}}{}}B_{\chi \mathrm{h}}^{(\mathrm{F})}\left\{{\displaystyle \frac{dN_{\overline{p}}^\mathrm{h}}{dE_{\overline{p}}}}(E_{\overline{p}}={\displaystyle \frac{E_{\overline{D}}}{2}})\right\}^2.`$
It may be expressed as a sum – extending over the various quarks and gluons h as well as over the different annihilation channels F – of the square of the antiproton differential multiplicity. That sum is weighted by the relevant branching ratios. The antineutron and antiproton differential distributions have been assumed to be identical. This readily leads to the source term for supersymmetric antideuterons

$$q_{\overline{D}}^{\mathrm{susy}}(r,z,E_{\overline{D}})=\sigma _{\mathrm{ann}}v\frac{dN_{\overline{D}}}{dE_{\overline{D}}}\left\{\frac{\rho _\chi (r,z)}{m_\chi }\right\}^2.$$

The four supersymmetric configurations of Tab. 8 are presented in Fig. 31 together with the calculated secondary spectrum (see Sect.6.1). The antideuteron spectra are fairly flat at low energy. They yield a few events for the ams observatory on board ISS below an IS energy of 3 GeV/nuc – a region from which secondary antideuterons are absent –.

##### Consequence for SUSY parameter space

For each configuration of the whole supersymmetric parameter space, the $`\overline{D}`$ flux has been integrated over that low-energy range. The resulting yield $`N_{\overline{D}}`$ which ams may collect is presented as a function of the neutralino mass $`m_\chi `$ in the scatter plot of Fig. 32. During the ams shuttle mission, the Solar cycle was close to maximum.

Most of the configurations are gaugino like (crosses) or mixed combinations of gaugino and higgsino states (dots). A significant portion of the parameter space is associated to a signal exceeding one antideuteron – horizontal dashed line –. In a few cases, ams may even collect more than a dozen of low-energy $`\overline{D}`$ nuclei. Notice finally that numerical simulations have shown that neutralinos – and more generally cold dark matter – should cluster in very dense and numerous clumps . Because neutralino annihilations proceed through a two-body reaction, their antiproton and antideuteron signatures would be enhanced by a large factor that could even reach up a few hundred in the case of density profiles à la Moore. If so, the entire supersymmetric constellation in Fig. 32 would be shifted upwards by at least two orders of magnitude and detection would become crystal clear for the most optimistic configurations. Another important question is related to the clumpiness of the dark halo. If the dark matter is concentrated into clumps, then the cosmic ray signals on Earth may be sensitive to the spatial distribution of the clumps. This aspect can be studied with the formalism used to determine the spatial origin of cosmic rays, as described above (Bottino et al., in preparation).

### 7.4 Primary antiprotons and antideuterons from Primordial Black Holes

Very small black holes should have formed in the early Universe from initial density inhomogeneities (Hawking ). They should now evaporate intensely through the Hawking mechanism (Hawking ) if their initial masses were around $`M_{}5\times 10^{14}`$ g. Detecting such objects is a great challenge of modern physics and cosmology as it would both allow to give experimental grounds to the Hawking radiation which is one of the only tentative achievement of semi-classical quantum gravity, and to probe the very small scales of the early Universe that remains totally inaccessible to other observations. Although the standard cosmological model of structure formation, assuming a pure scale-invariant Harisson-Zeldovitch power spectrum normalized to CMB amplitudes, would lead to a very small amount of PBH dark-matter in the present Universe, several realistic inflationary scenarios and phase-transition phenomena can produce a significant amount of PBHs.

#### 7.4.1 Hawking evaporation and antiproton source term

The Hawking black hole evaporation process can be intuitively understood as a quantum creation of particles from the vacuum by an external field. The basic characteristics can be easily seen through a simplified model, and the interested reader is referred to  for more details.

##### Elementary energy spectrum

The accurate emission process was derived by Hawking, using the usual quantum mechanical wave equation for a collapsing object with a post-collapse classical curved metric. He found that the emission spectrum for particles of energy $`Q`$ per unit of time $`t`$ actually mimics a Planck law and is given by, for each degree of freedom:

$$\frac{\mathrm{d}^2N}{\mathrm{d}Q\mathrm{d}t}=\frac{\mathrm{\Gamma }_s}{h\left(\mathrm{exp}\left(\frac{Q}{h\kappa /4\pi ^2c}\right)(1)^{2s}\right)},$$
(58)
where the effect of angular velocity and electric potential have been neglected since the black hole discharges and finishes its rotation much faster than it evaporates . In this expression, $`\kappa `$ denotes the surface gravity, $`s`$ is the spin of the emitted species and $`\mathrm{\Gamma }_s`$ is the absorption probability, given by

$$\mathrm{\Gamma }_s=\frac{4\pi \sigma _s(Q,M,\mu )}{h^2c^2}(Q^2\mu ^2),$$
(59)
where $`\sigma _s`$ is the absorption cross section computed numerically  and $`\mu `$ the rest mass of the emitted particle. It is also convenient to introduce the Hawking temperature, defined by

$$T=\frac{hc^3}{16\pi kGM}\frac{10^{13}\mathrm{g}}{M}\mathrm{GeV}.$$
(60)
The expression (58) may then be written simply as a function of $`Q/kT`$.

##### Antiproton multiplicity

As it was shown in , when the black hole temperature is greater than the quantum chromodynamics confinement scale $`\mathrm{\Lambda }_{\mathrm{QCD}}`$, quarks and gluons jets are emitted instead of composite hadrons. To evaluate the number of emitted antiprotons $`\overline{p}`$ , one therefore needs to perform the following convolution:

$`{\displaystyle \frac{\mathrm{d}^2N_{\overline{p}}}{\mathrm{d}E\mathrm{d}t}}={\displaystyle \underset{j}{}}{\displaystyle _{Q=E}^{\mathrm{}}}{\displaystyle \frac{\alpha _j}{e^{Q/kT}(1)^{2s_j}}}{\displaystyle \frac{\mathrm{\Gamma }_{s_j}(Q,T)}{h}}`$
$`\times {\displaystyle \frac{\mathrm{d}g_{j\overline{p}}(Q,E)}{\mathrm{d}E}}\mathrm{d}Q,`$
where $`\alpha _j`$ is the number of degrees of freedom, $`E`$ is the antiproton energy and $`\mathrm{d}g_{j\overline{p}}(Q,E)/\mathrm{d}E`$ is the normalized differential fragmentation function, i.e. the number of antiprotons created with an energy between $`E`$ and $`E+dE`$ by a parton jet of type $`j`$ and energy $`Q`$. The fragmentation functions have been evaluated with the high-energy physics frequently used event generator pythia/jetset .

##### Whole spatial and spectral source distribution

Once the spectrum of emitted antiprotons is known for a single PBH of given mass, the source term used for propagation is given by

$$\frac{d^3N_{\overline{p}}(E)}{dEdtdV}=_0^{\mathrm{}}\frac{d^2N_{\overline{p}}}{dEdt}(M,t_0)\frac{d^2n}{dMdV}𝑑M,$$
(61)
where $`d^2n/dMdV`$ is the PBH mass spectrum today. It can be deduced from the initial mass spectrum through the mass loss rate which reads as $`dM/dt=\alpha (M)/M^2`$ (by a simple integration of the Hawking spectrum multiplied by the energy of the emitted quantum) where $`\alpha (M)`$ accounts for the available degrees of freedom at a given mass. In the assumption $`\alpha (M)const`$ it leads to:

$$\frac{d^2n(M)}{dMdV_i}=\frac{M^2}{(3\alpha t+M^3)^{2/3}}\frac{d^2n_i}{dM_idV_i}|_{M_i=M_i^0}$$
(62)
where the derivative is evaluated at $`M_i^0=(3\alpha t+M^3)^{1/3}`$. The initial values $`d^2n_i/dM_idV_i`$ were shown  to scale as $`M_i^{5/2}`$ for PBHs formed, as expected, in a radiation dominated Universe from a scale invariant power spectrum. The nowadays spectrum is therefore mostly identical to the initial one above $`M_{}3\alpha t_05\times 10^{14}`$g and proportional to $`M^2`$ below.

The last problem is that the spatial distribution of these PBHs is basically unknown. However, these objects should have formed in the very early stages of the history of the Universe: their initial mass is very close to the horizon mass at the formation epoch ($`M_HM_{\mathrm{Pl}}\times t_{\mathrm{form}}/t_{\mathrm{Pl}}`$ which means that the formation time $`t_{\mathrm{form}}`$ is of the order of $`10^{23}`$ s for $`M=M_{}`$). Once formed, they only interact through gravitation, so that they behave as cold dark matter, and their spatial distribution should be similar. As a consequence, we use the same profile for the PBHs distribution as for the neutralinos, as given in Eq. (15).

#### 7.4.2 Resulting upper limit and cosmological consequences

As shown in  the main uncertainties on the estimation of the antiproton flux from PBHs are associated with astrophysical parameters. Several possible effects due to the quite large size of the horizon mass at the end of the inflation period (mostly imposed by gravitinos and moduli fields constraints) and to the possible photosphere that could form around hot PBHs were studied in  and show that the present estimates could be substantially revised in the future. The emitted antiproton spectrum is mostly resulting from PBHs with masses between $`10^{12}`$ g and $`10^{14}`$ g: the lightest ones are not numerous enough because of the $`M^2`$ shape of the nowadays mass spectrum and the heaviest ones do not emit much because of their low temperature.

##### Density exclusion criterion

Superimposed with bess, caprice and ams data, the full antiproton flux, including the secondary component and the primary component, is shown in Fig. 33 for 20 values of $`\rho _{}^{\mathrm{PBH}}`$ logarithmically spaced between $`5\times 10^{35}`$ and $`10^{32}\mathrm{g}\mathrm{cm}^3`$ with fixed astrophysical parameters. The lowest curves are clearly in agreement with data whereas the upper ones

contradict experimental results. To derive a reliable upper limit, and to account for asymmetric error bars in data, we define a generalized $`\chi ^2`$ as

$`\chi ^2={\displaystyle \underset{i}{}}{\displaystyle \frac{(\mathrm{\Phi }^{\mathrm{th}}(Q_i)\mathrm{\Phi }_i^\mathrm{e})^2}{(\sigma _i^{\mathrm{e}+}+\sigma ^{\mathrm{th}+}(Q_i))^2}}\mathrm{\Theta }(\mathrm{\Phi }^{\mathrm{th}}(Q_i)\mathrm{\Phi }_i^\mathrm{e})`$
$`+{\displaystyle \underset{i}{}}{\displaystyle \frac{(\mathrm{\Phi }^{\mathrm{th}}(Q_i)\mathrm{\Phi }_i^\mathrm{e})^2}{(\sigma _i^\mathrm{e}+\sigma ^{\mathrm{th}}(Q_i))^2}}\mathrm{\Theta }(\mathrm{\Phi }_i^\mathrm{e}\mathrm{\Phi }^{\mathrm{th}}(Q_i)).`$
where $`\sigma ^{\mathrm{th}+}`$ and $`\sigma ^{\mathrm{e}+}`$ ($`\sigma ^{\mathrm{th}}`$ and $`\sigma ^\mathrm{e}`$) are the theoretical and experimental positive (negative) uncertainties. An upper limit on the primary flux for each value of the magnetic halo thickness is computed. The theoretical errors included in the $`\chi ^2`$ function come from nuclear physics ($`p+He\overline{p}+X`$ and $`He+He\overline{p}+X`$ ), and from the astrophysical parameters analysis which were added linearly, in order to remain conservative. The resulting $`\chi ^2`$ leads to very safe results as it assumes that limits on the parameters correspond to 1 sigma.

##### Resulting constraints

Fig. 34 gives the upper limits on the local density of PBHs as a function of $`L`$. It is a decreasing function of the halo size because a larger diffusion region means a higher number of PBHs inside the magnetic zone for a given local density. Between $`L=1`$ and $`L=15`$ kpc (extreme astrophysical values), the 99% confidence level upper limit goes from $`5.310^{33}`$ g cm-3 to $`5.110^{34}`$ g cm-3.

Such limits are of great cosmological relevance. They are a unique probe of the early Universe on very small scales that could be equaled only, in a far future, by the detection of primordial gravitational waves. Based on complementary, but nearly equivalent, limits coming from the 100 MeV gamma-ray background, numerous constraints were derived on the spectral index of scalar perturbations (see, e.g. ): $`n<1.27`$. A too “blue” spectrum would lead to an overproduction of PBHs in conflict with observations.

This upper limit on the local PBH density obtained with antiprotons can also be used to severely constrain models with Broken Scale Invariance (BSI). A jump in the first derivative of the inflaton potential  should lead to a huge increase of the PBHs formation probability when the corresponding scale re-enters the horizon . The antiproton limit can, therefore, be directly translated into an exclusion area in the parameter space of BSI inflationary models . Finally, it allows interesting prospects for PBH dark matter  or Planck relics investigations .

#### 7.4.3 A new window for detection: antideuterons

To go beyond an upper limit and try to detect PBHs it seems very interesting to look for antideuterons. Below a few GeV, there is nearly no background for kinematical reasons  and the possible signal due to PBHs evaporation could be easy to detect. The emission scheme is nearly the same but the probability that an antiproton and an antineutron merge into an antideuteron is taken into account. The possible detection range for the ams experiment  can be evaluated. It is shown on Fig. 35 as a function of the three unknown parameters: $`L`$, the height of the magnetic halo, $`p_0`$, the coalescence momentum to form an antideuteron, and $`\rho _{}^{\mathrm{PBH}}`$, the local density of primordial black holes. The sensitivity of the experiment should allow, for averaged parameters, an improvement in the current best upper limit by a factor of six, if not a positive detection. The situation is very different than for antiproton: the main limitation is due to the experimental sensitivity and not to the unavoidable physical background. Great improvements can, therefore, be expected in the future and this investigation means seems very promising.

## 8 Summary, conclusions and perspectives

A consistent framework to understand the propagation of CR nuclei in the energy range 100 MeV–100 GeV was presented in this paper. The observed fluxes of most species can be explained by assuming that once emitted from some sources located in the galactic disc, these nuclei undergo a diffusive propagation altered by escape through the boundaries, spallations, reacceleration, energy losses and galactic wind. The magnitude of these effects has been constrained using the B/C data, and the consistency of the model has been tested against the observed antiproton flux and by the study of radioactive species.

This well-tested model has then been used to study the propagation of cosmic rays of a more hypothetical origin, such as light antinuclei produced by SUSY galactic Dark Matter or Primordial Black Holes. In particular, limits on the abundance of primordial black holes in the galactic halo could be set.

There are several directions in which this work may be extended. First, the constraints on the propagation parameters could be refined by considering other species, stable or secondary. However, this approach is currently limited by the accuracy of the available data on cosmic ray fluxes and on the nuclear cross sections. Second, a specific study of the EC unstable species could provide valuable information about the processes responsible for the acceleration of cosmic rays. Third, the existing limits on the SUSY induced antiproton signal can be bettered by using the constraints on the propagation parameters in a fully consistent way. Finally, the propagation code we use should ultimately be able to yield the flux of all cosmic ray species (including gamma rays, electrons and positrons) at every position in the Galaxy.

## Acknowledgments

This work has benefited from the support of the Programme National de Cosmologie and PICS 1076, CNRS. The authors are grateful to E. Vangioni-Flam & M. Cassé for enduring support and fruitful collaboration on several aspects of CR sources and propagation. We acknowledge G. Pelletier for clarifying discussion on diffusion coefficients. One of the authors, D.M., is deeply indebted to A. Soutoul, E. Vangioni-Flam & M. Cassé for illuminating discussions and for their constant support.