© ESO, 2012

1. Introduction

In the 1 GeV to 100 TeV energy range, cosmic ray (CR) nuclei that reach the Earth have a Galactic origin. They were accelerated in sources, the nature of which is still a subject of discussion and research, supernova-driven shock waves being a good candidate. The CR nuclei subsequently reach the Earth, after diffusing in the Galactic magnetic field. The exact locations and ages of the supernova explosions are not known, and properties of Galactic cosmic rays are often studied under the assumption that sources form a distribution that is continuous, in space as in time. In the standard model of CR propagation, sources are modelled as a jelly that extends inside the Galactic disk and steadily injects particles into the interstellar medium (ISM). This picture has proven to be quite successful so far. The CR spectra are indeed dominated by sources which are distant from Earth (for instance see Fig. 8 in Taillet et al. 2004) and the continuous hypothesis is expected to be valid. On the other hand, this approach may fail at the high-energy end of the CR spectra, which is mainly produced by young sources. The flux at high energy could well be dominated by a handful of individual sources, the positions and ages of which have a strong influence on the spectrum. This is of tremendous importance because the latest measurements, such as those of CREAM (Yoon et al. 2011) or PAMELA (Adriani et al. 2011), exhibit spectral features that are difficult to understand in a continuous and steady-state modelling of CR sources.

In this article, we study the relation between the mean-field approach and the so-called “myriad model” (Higdon & Lingenfelter 2003) according to which CR primary species, such as protons or helium nuclei, diffuse from point-like sources that are distributed randomly in space and time. In the myriad model, the flux of cosmic rays is a random variable. When averaged over many realisations of possible positions and ages of the sources, the myriad model gives the same results as the continuous model, as expected.

However, we expect fluctuations from the average flux, from one realization to the other. If these fluctuations are small, the myriad model predictions of the flux are robust and similar to that of the continuous model. In that case, the spectral features seen by CREAM and PAMELA should be explained within the continuous model. Invoking a distorted spectral injection spectrum at the sources (see for instance Yuan et al. 2011) or a non-standard behaviour of the Galactic diffusion coefficient with energy are two possible, although unlikely, solutions. On the other hand, if the statistical fluctuations of the myriad model are high, it could be easier, and to some extent more natural, to explain the CREAM and PAMELA anomalies in terms of the particular history and morphology of the population of sources within which we live. If so, there is little reason to trust the continuous model anymore, on which most of CR studies are nevertheless based. The “real” sources could generate CR fluxes at the Earth quite different from those predicted in the conventional approach. Unfortunately, we do not know where and when the “real” supernova explosions took place. We are reduced to rely on statistics to gauge how probable CR fluxes are that differ by much from the values predicted by the continuous model.

The amplitudes of the statistical fluctuations of the myriad model have recently been studied through the standard deviation of the random variable associated to the CR proton flux (see e.g. Blasi & Amato 2012). This quantity turns out te be infinite when computed without precautions. This unexpected result is, a priori, a severe blow for the mean-field approach on which most of the public codes of CR Galactic transport, such as galprop¹, dragon² or usine³, are based. These codes would have to be entirely modified. Abandoning the conventional CR model would also have drastic consequences on the problem of the astronomical dark matter (DM). The nature of this component, which contributes substantially to the mass budget of the universe, is still an unresolved issue. Weakly interacting massive particles (WIMP) have been invoked as a plausible solution. These putative DM species are expected to annihilate themselves inside the Galactic DM halo, producing in particular rare antimatter cosmic rays such as positrons and antiprotons. Distortions in the spectra of these particles would be an indirect probe of the presence of WIMPs inside the Milky Way. Theoretical predictions have been based so far on the continuous CR model. They need to be entirely revised should the myriad model replace the conventional scheme.

This article is devoted to the myriad model and to the problem arising from the infinite variance of the primary CR flux. We suggest two ways to solve this problem. First, we show that the statistical distribution of the flux is such that the confidence intervals are finite and well-defined, even though the variance is infinite. We study how the width of these intervals depends on the CR propagation model. Second, we notice that the divergence of the variance is due to the sources that are very close to the Earth and very young. In our Galaxy, theses objects are known actually, and the region in which we have catalogues of the sources should be excluded from the statistical ensemble of the myriad model, and should be treated as known sources. This paper explores the consequences of this approach.

2. Propagator for diffusive propagation from discrete sources

After their acceleration by supernova-driven shock waves, primary cosmic ray (CR) nuclei subsequently propagate through the Galactic magnetic fields and bounce off their irregularities – the Alfvén waves. The resulting particle transport is well described by space diffusion with a coefficient (1)which increases as a power law with the rigidity R ≡ p/Z of the particle. In addition, since these scattering centres move with a velocity V_a ~ 20 to 100 km s^-1, a second-order Fermi mechanism is responsible for some diffusive re-acceleration, which turns out to be mostly relevant at low energy, below a few GeV. Since we are interested here in the excess in CR protons and helium nuclei measured by CREAM, and more recently by PAMELA, we disregarded diffusive re-acceleration as well as energy losses, which are also negligible at high energy. In addition to space diffusion, Galactic convection wipes cosmic rays away from the Galactic disk with a velocity V_C ~ 5 to 15 km s^-1. The master equation for the space and energy number density ψ ≡ dn_p/dT_p of, say, CR protons with kinetic energy T_p takes into account space diffusion and Galactic convection (2)This equation applies to any nuclear species – protons and helium nuclei in particular – as long as the effective rate of production q_eff is properly accounted for. This general framework, summarized in Eq. (2), is generally implemented within an axisymmetric two-zone model that is extensively discussed in Maurin et al. (2001) or Donato et al. (2001), and whose salient features are briefly recalled here.

The CR diffusive halo (DH) is pictured as a thick wheel that matches the circular structure of the Milky Way. The Galactic disk of stars and gas – where primary CR protons and helium nuclei are accelerated – lies in the middle. It extends radially 20 kpc from the centre and has a half-thickness h of 100 pc. Confinement layers, where cosmic rays are trapped by diffusion, lie above and beneath that disk. The interGalactic medium starts at the vertical boundaries z =  ±L as well as beyond a radius of r = R_gal =  20 kpc. Notice that the half-thickness L of the DH is not known and possible values range from 1 to 15 kpc. The diffusion coefficient D is taken to be the same everywhere while the convective velocity is exclusively vertical, with component (3)The Galactic wind – produced by the bulk of the disk stars like the Sun – drifts away from them along the vertical directions, hence the particular form assumed here for V_C. The effective production term q_eff takes into account the release of primary CR nuclei described by the positive term q_acc as well as a negative contribution that accounts for their interactions with the interstellar gas of the disk with rate q_col. In most propagation codes, the injection of cosmic rays is described by a smooth function of space and is constant in time. In this article, we treat supernova explosions as point-like events. The production rate of CR nuclei through acceleration is given by (4)where each source i belonging to the population P of supernovae contributes a factor q_i at position x_i and time t_i. The collisions of CR nuclei on the hydrogen and helium atoms of the interstellar medium (ISM) tend to deplete the high-energy regions of their spectra. For CR protons, which we will consider throughout this section as an illustration, they act as a negative source term for the energy bin located at T_p, with amplitude (5)The ISM is distributed within an infinitely thin disk, hence the presence of a δ(z_S) function in the previous expression. The collision rate may be expressed as (6)where the densities n_H and n_He have been averaged to 0.9 and 0.1 cm^-3, respectively. The total proton-proton cross section σ_pH was parameterized according to Nakamura et al. (2010) while σ_pHe is related to σ_pH by the Norbury & Townsend (2007) scaling factor 4^2.2/3. Similar scaling factors were used to derive the CR helium nuclei collision cross sections from the proton case.

2.1. Solution of the diffusion equation

Our description of the propagation of CR protons through the DH, and of any other CR nucleus for that matter, relies on the existence of the propagator G_p. This Green function translates the probability for a CR proton injected at position x_S ≡ (x_S,y_S,z_S) and time t_S to travel through the Galactic magnetic field until it reaches, at time t, an observer located at x ≡ (x,y,z). The CR proton density at the Earth may be expressed as the convolution over space and time of the Green function G_p with q_acc(7)where q_acc(x_S,t_S) is the CR proton injection rate at the source located at x_S and at time t_S. Notice that for a population P of point-like sources, the injection rate q_acc is given by Eq. (4) which, once inserted into the previous expression for ψ, leads to the proton flux at the Earth (8)The Green function G_p plays a crucial role in our analysis and is the ideal tool to study the effect of the discreteness of the sources. It is a solution of the CR transport equation (9)The resolution of this equation is presented in Appendix A.

2.2. Diffusion parameters

The computed flux depends on the diffusion parameters. The results will be presented for three benchmark sets of diffusion parameters, consistent with the energy dependance of the B/C ratio (Maurin et al. 2001). They are labelled “min”, “med” and “max”, according to the value of L. The values are indicated in Table 1.

3. Mean flux

The flux at solar position is obtained by summing the contributions of the point sources. The mean value can be computed once the distribution of sources is known. We first provide the result for a homogeneous distribution of sources, because an analytic expression can be found in this case. We then turn to a more realistic distribution.

3.1. Homogeneous distribution of sources

The average flux from one point source, drawn from a statistical ensemble with distances ranging from R_min to R_max and ages from t_min to t_max, is given by averaging G₁ ≡ G_p over distances and times. To obtain an order of magnitude estimate, we first consider that the sources are distributed evenly in space and time in the Galactic disk, which is assumed to be infinitely thin and to extend to R_max → ∞. We obtain Considering ν sources per unit time and unit area in the disk, the average flux from sources is given by By integrating over time (t_max → ∞) and using (Gradshteyn et al. 2007, 3.471.9) (10)where K_ν is the νth order modified Bessel function of the second kind, leads to where G_steady stands for the one source steady-state propagator, given by (11)Finally, using ∫x   K₀(x)   dx = −x   K₁(x), we have (12)for an infinite disk, R_min = 0 and R_max → ∞, so that, using xK₁(x) → 1 when x → 0, It can be shown that  ∑ _n1/(n + 1/2)² = π²/2, and finally for the infinite disk, which is what is also obtained by directly solving the diffusion equation in steady state. The mean value of the flux from randomly distributed point sources is equal to the steady-state flux obtained with a continuous source distribution. This is also true in more general cases such as thick disk, finite radius, with wind and spallation.

3.2. Realistic distribution of sources

We now consider the case of a general distribution of sources, with a radial distribution in the Galactic disk, as well as a distribution across the thickness of the disk. It is widely accepted that up to energies corresponding to the knee, cosmic ray accelerators are supernova remnants (SNR). Unfortunately, supernovae are pretty rare events and their spatial distribution is difficult to measure accurately. However, because pulsars are created in SNR and are easier to detect, it is a fair asumption that the cosmic ray sources follow the pulsar distribution. In this work we assume that the radial profile of sources follows the pulsar distribution given in Yusifov & Küçük (2004), where distances are expressed in kpc. The distribution along the z axis is given by where z₀ is set to the half-thickness of the Galactic disk, z₀ = h. In this section, z_s, r_s, θ_s and t_s refer to the position and age of sources. The solar system is located at r_⊙ = 8.5kpc, θ = 0 and z = 0. Moreover, we define ρ(r_s,θ_s) as the distance from a source to the solar system: We evaluate the scatter of the flux that arises when we consider point-like sources with positions and ages following a given probability distribution. The sources are supposed to be independent, so the mean value of the propagator coming from all the N sources is just given by (13)Assuming cylindrical symmetry and a uniform age distribution, up to a maximum age T_max, the normalization factor A is given by We find (the antisymmetric part of the propagator vanishes at z = 0) We separate the integration and define (14)The integration is easy using the exponential form of sinus and leads to (15)The integration over time (t_max → ∞) is made using Eq. (10): (16)Then we can write the mean value as (17)This result is integrated numerically over the two coordinates r_s and θ_s. The mean flux is then given by where g(E) is the source energy injection term. We checked that this calculation gives the same result as the steady state model. In the next section, we perfom the calculation of the variance associated to the flux.

4. Statistical analysis: variance of the flux

The flux obtained in the myriad model depends on the exact positions and ages of all sources. It is bound to be different from the mean values computed above, at some level. How different from the mean value are the typical fluxes in the myriad model likely to be? To address this question, the usual method is to compute the standard deviation associated to the flux, considered as a random variable.

4.1. Realistic distribution of sources

The variance of the propagator for one source is defined as The variance for N independent sources is given by We must compute the average value of , given by (18)At fixed m and n, we define (19)This is easily computed, and we can integrate the time distribution to obtain the variance from sources of all ages: (20)Using Eq. (10) and considering t_max → ∞ (21)we find (22)We are interested in the behaviour of this quantity at the lower bound in r. When , the property K₁(x) → 1/x as x → 0 can be used to show that the above integral diverges as ln(R_min) as the lower bound R_min of the spatial integral goes to zero. The spatial distribution probability is f_r ~ 1 when r → 0, so that (23)

4.2. Divergence of the variance: how strong is it?

The mean value, computed considering r_s from r_s = 0 to r_s = 20 kpc and t_s from t_s = 0 to t_s = ∞, is equal to the flux obtained in the steady-state model, as commonly obtained in propagation theories. But as we have just seen, the calculation related to the standard deviation of this quantity diverges when we allow the possibility of sources that are close (r_s → 0) and young (t_s → 0), as already noticed previously (Blasi & Amato 2012).

The standard deviation is commonly interpreted as the typical spread of the random values around the mean, and a high standard deviation could be interpreted as if the actual value of the flux had a disturbingly high probability to be very far from the mean value.

One could argue that the problem we considered is physically irrelevant, because we know for sure that there is no supernova remnant with zero age and null distance to the Earth. One can impose a lower cut-off in ages and distances, based on observations. However, even with reasonable values for the cut-off, the variance can be finite but large (see Sect. 4.3). One can also adjust the cut-off to eliminate the very rare events that make the standard deviation very high, without contributing significantly to the the mean value. This is the approach adopted by Blasi & Amato (2012). These authors used a cut-off given by . This condition gives a fair order of magnitude of the spread of the values around the mean. It is difficult though to interpret it in rigourous statistical terms. Indeed, the value of the variance depends quite strongly on the exact position of the cut-off, as featured in Fig. 1 where the evolution of the mean and the standard deviation of the total proton flux derived with the max model at 10 TeV is presented. Choosing t_min/2 or 2    ×    t_min rather than t_min, for instance, has a small effect on the mean but a drastic effect on the standard deviation. The value chosen by Blasi & Amato (2012) is indicated by an arrow in the figure. The mean obtained with this cut-off is about 10% lower than the true mean. A lower cut-off would give a more precise mean value, but a much larger variance.

This confusing situation, in which some rare events have a very small contribution to the mean, but give rise to a very high standard deviation, is not uncommon in physics (Levy flights, Cauchy distribution). The total flux φ is the sum of a myriad of N individual contributions ϕ_i from single sources. Assuming that these contributions are not correlated with each other allows one to relate simply the means and variances of the total and single source fluxes through (24)We show in the Appendix that the probability p(ϕ) of measuring a single source flux ϕ is given by p(ϕ) ∝ ϕ^−8/3 for objects located in a thick disk and p(ϕ) ∝ ϕ^−7/3 for sources located in a thin disk, in the limit where ϕ → ∞. The standard deviation σ_φ of the total flux is related to the integral , which diverges as ϕ → ∞. For this distribution with an infinite second moment, the central limit theorem does not hold, at least in its usual form. In this case, the standard deviation is a poor estimator of the typical spread of the values around the mean (in particular the probability distribution function p(φ) is no longer Gaussian), and there is no use trying to regularize this quantity by applying cut-offs. Notice that the average value  ⟨ φ ⟩  is still well-defined, because  ⟨ ϕ ⟩  = ∫ϕ   p(ϕ)   dφ is convergent, and more to the point, the confidence intervals of the global probability distribution p(φ) are well-defined too, since the integral of p(ϕ) is well-behaved. The spread of the random values of the total flux around its mean  ⟨ φ ⟩  can then be studied by computing the quantiles associated to the probability distribution p(φ), rather than using the standard deviation, which has no clear statistical relevance in this case.

All this discussion arises from the fact that there is a non-vanishing probability in the myriad model to find sources that are arbitrarily close and young. Although the variance σ_φ cannot be defined in this case, we can still infer the global distribution function p(φ) and delineate intervals inside which the total fux φ is mostly expected. Another possibility is to take advantage of the solid information that has been collected on these young and nearby sources, whose statistical properties are disconcerting. We can remove a space-time region around the Earth containing the very close and very young sources of our Galaxy, and replace it with a catalogue, built from observations.

Below, we investigate these points in more details. First we study how the standard deviation σ_φ can be regularised by applying some cut-offs related to astronomical information about the sources. We then focus on the confidence intervals of p(φ) and show that they are well-defined, with no need of a particular cut-off. Finally, we take into account the astronomical knowledge about the local sources of cosmic rays and explore in greater detail how introducing such a catalogue improves our statistical analysis.

4.3. Regularisation by cutting off the ages

From observations of the solar neighbourhood, we have a good idea of the distribution of sources that are young and close. Given the age t_min of the youngest local supernova remnant, one can compute the mean value and standard deviation of the total flux φ, by applying that cut-off t_min to the age distribution.

For the sake of illustration, Fig. 2 features the mean value  ⟨ φ ⟩  and the standard deviation σ_φ of the total proton flux using a lower cut-off of t_min = 100yr in the source distribution, for the min set of propagation parameters. In the same figure we have also plotted the data points from the CREAM and PAMELA experiments. With the chosen value for the cut-off, the standard deviation at high energies remains of the same order of magnitude as the flux. The relative value of the standard deviation, σ_φ/ ⟨ φ ⟩ , is fairly independent of the energy. This trend can be shown analytically to hold for sources located in a thin disk. In the framework of a purely diffusive model (no wind/spallation), the relative dispersion can be approximated by (25)This ratio does not depend on the diffusion coefficient D, hence it does not depend on energy. It is of the order of unity for t_min ~ 100yr.

For sources distributed in a disk with a finite thickness h, one obtains the same result as before as long as . In the opposite limit, one finds (26)In both cases, the relative standard deviation does not vary much with energy.

The value of σ_φ/ ⟨ φ ⟩  depends sensitively on the chosen cut-off t_min. For a very low value of t_min, the standard deviation is very large. Conversely, for higher values of t_min, the standard deviation decreases, and the average value  ⟨ φ ⟩  starts to be affected by the cut-off, as shown in Fig. 1.

4.4. The variance is infinite but the confidence levels are finite

The variance of the flux being infinite does not necessarily imply that the random values are typically very different from the mean. To illustrate this affirmation, consider the case of a unique point-like steady-state source located in the Galactic disk, with cosmic ray diffusion taking place in a boundless space. The solution of the diffusion equation is given by ϕ = a/r where a is a constant. Assuming that this source is uniformly distributed inside a disk of radius R leads to the probability distribution function (27)We can readily infer the mean flux (28)and the average value of the flux squared (29)where we have introduced a cut-off value ϵ at the lower end of the radial distribution to exhibit the divergence of  ⟨ ϕ² ⟩ . The variance of ϕ goes to infinity as ϵ → 0.

However, the distribution of ϕ (which is what we are really interested in) is well-behaved. From the relation between r and ϕ, we can write dr = a   dϕ/ϕ² so that (30)The probability that the flux is lower than a given value Φ may be expressed as (31)provided that Φ > Φ₀ ≡ a/R. Introducing the Heavyside distribution Θ leads to (32)The probability that Φ > 10    ⟨ ϕ ⟩  is only 1/400, even though the variance is infinite. Indeed, the flux is more likely to be lower than the mean value, whereas one might have guessed the opposite, considering the divergence of the variance.

When N sources are considered, the mean flux value and the variance are both just multiplied by N. The probability distribution p_N(φ) for the flux can be obtained by recurrence from (33)These are displayed in Fig. 4. The variance still diverges. In the high-φ region, the flux is dominated by the contribution of a single source and the probability distribution is given by (34)For N = 100 sources, the probability that φ > 2    ⟨ φ ⟩  is 2.5    ×    10^-3 and P(φ > 10    ⟨ φ ⟩ ) is vanishingly small.

If we now consider time-dependent sources spread homogeneously inside an infinite DH with pure diffusion, the variance is given by the integral (35)which diverges with the lower cut-off in ages as . For a 3D homogeneous distribution of steady-state sources, σ_ϕ does not diverge (see Table 2). For the sake of illustration, Fig. 5 presents histograms for the flux obtained with the propagator discussed in Sect. 2, at several energies.

4.5. Quantiles of the total flux

From now on, we show the quantiles associated to the fluxes. These are defined as regions with a given probability to find the flux. In Fig. 6, we present “deciles” (10% quantiles), as well as 68% confidence intervals (CI), which are more familiar, because they correspond to 1-σ intervals for Gaussian distributions.

To compute these figures, we ran a Monte Carlo simulation over  ~3000 populations. We also adjusted the injection rate so that the mean flux follows the PAMELA data below 200 GeV: (36)in cm^-2   s^-1   sr^-1   GeV^-1.

The spread of the fluxes around the mean is not sensitive to this choice. In spite of its infinite variance, the flux has a well-defined probability distribution as featured in Fig. 6. The possible flux lies on a band whose thickness depends sensitively on the propagation model. In particular, the width increases when the thickness L of the diffusive halo decreases.

At high energy, the CREAM and PAMELA data cannot be explained in the med and max cases. For the min case, taking ν = 1century^-1 would enlarge the uncertainty band, so that the explanation in terms of sources becomes more probable.

5. Known local sources – the catalogue

5.1. Regularisation using a catalogue

The situation described above occurs as long as we are ignorant of the positions and ages of the CR sources. The young and nearby objects are responsible for the divergence of the flux variance and potentially lead to the problems encountered above in the statistical analysis. However, we do have data concerning the distribution of nearby sources, for which catalogues are available. A natural way then to regularise the variance is to separate the sources into two lots. The first set contains the young and local sources, which can be extracted from the catalogues. The second group, ofwe have which little information, comprises the old or distant sources and will be treated according to the statistical analysis of Sect. 4. This procedure allows us to regularise the variance of the flux in the most natural way while reducing its uncertainties. Following Delahaye et al. (2010), we have used two catalogues.

• (i)

  The Green survey (Green 2009) compiles various informations on supernova remnants, but fails to systematically provide their ages or the precision with which their distances from the Earth have been determined. A quite thorough bibliographic work has been summarised in the Appendix of Delahaye et al. (2010), which we have borrowed as a complement to the Green catalogue. In total, we have collected 27 local SNR with their ages, distances from the Sun and, when possible, the corresponding observational uncertainties.

• (ii)

  Pulsars are not expected to be sources of primary CR nuclei. As residues of supernova explosions, they are nevertheless a good tracer of old SNR that are too old to be directly detected in radio waves. Moreover, because they are point-like objects, their distances from the Sun is much easier to measure. Their ages can also be estimated precisely through spin-down. After removing millisecond pulsars from the ATNF catalogue (Manchester et al. 2005) (by selecting Ṗ > 5 × 10^-18) and objects associated with known SNR, we are left with 157 objects with their ages and distances.

  Table 3

  Characteristics of nearby SNRs, from Delahaye et al. (2010).

It is fair to ask whether these 27+157 sources are representative of the local environment. As featured in Fig. 7, the number of objects found in the catalogues agrees well with what can be inferred from various Galactic distributions found in the literature, provided that the supernova explosion rate is, on average, approximately equal to 3 per century. At least, this is true within the 2 kpc nearby the Sun, for sources younger than 30 000 years. At later ages, SNR become too dim to be detected and the Green survey cannot be trusted anymore. However, the catalogues agree fairly well with theoretical expectations and do not suffer from major biases, at least not more than the theoretical models. In this work, we define the “local” region as the domain extending 2 kpc around the Sun, with sources younger than 50 000 years. A rate of one supernova explosion per century in our Galaxy is also quoted in the literature (Delahaye et al. 2010). If this is the case, we would be in a locally high-density region of sources. We have plotted the flux considering this low rate, but will discuss the validity of this assumption in a forthcoming Letter. The CR proton flux produced by the SNR and pulsars of our catalogue is presented in Fig. 8 for the “min” CR propagation benchmark model. The sources mainly contributing to the flux at high energy are given in Table 3, which was taken from Delahaye et al. (2010). In the PAMELA energy region, the flux is dominated by pulsars, whereas SNR come into play at the energies of CREAM.

5.2. Is the catalogue plausible?

To check the plausibility of our catalogue, we compare in Fig. 9 the flux φ_cat yielded by the sources that the catalogue contains (blue curve) to the flux φ_loc from a set of populations drawn randomly inside the local region (black curve). As explained above, we did not compute the variance of these random populations, but derived the confidence intervals through Monte-Carlo simulation. Assuming an explosion rate of three events per century in the Galaxy, we find that the flux φ_cat yielded by the objects of our catalogue lies within the 1-σ confidence interval surrounding the mean value φ_loc. This is true for the min, med, and max benchmark sets of CR propagation parameters. We conclude that the theoretical source distribution and explosion rate, which we have chosen here, agree fairly well with local realistic sources. They are not necessarily representative of the entire Galaxy though. How potential differences would impact the relative importance of φ_cat with respect to the total flux will be detailed in a forthcoming Letter.

5.3. Systematic errors from the catalogue

We have also plotted the uncertainty on φ_cat that arises from the ages and positions of the SNR of our catalogue. These have been varied within the ranges allowed by observations. We did not consider the uncertainties on the ages and positions of pulsars. The former are well determined. The rotation periods P and drifts in time Ṗ of pulsars are measured with good accuracy. The uncertainty on φ_cat corresponds to the grey bands of Figs. 9 and 11. The shape of these curves is explained by the fact that when L is lower (“min” model), the effective range of diffusion in the plane of the disk is reduced and the contribution of remote sources is lower. The relative contribution of the catalogue is then higher for lower values of L.

6. A mixed approach

The cosmic ray flux from the Galactic sources and the associated uncertainties can be computed as where φ_cat is the flux coming from the sources that belong to the catalogue (closer than 2 kpc and younger than 3    ×    10⁵ yr) and φ_ext(E) is the flux coming from sources that are more distant than 2 kpc or older than 3    ×    10⁵ yr. The uncertainties on φ_cat(E) are obtained by considering the observational errors on the source parameters (see Fig. 11). The uncertainty on φ_ext(E) can be evaluated by computing the confidence intervals, as above. In some conditions, this is equivalent to computing the usual standard deviation of the flux, because the variance of the flux located outside of the region covered by the catalogue is finite. This is shown in Fig. 10.

7. Conclusions

In the conventional model of Galactic CR propagation, the sources of primary nuclei are described as a jelly that spans the disk of the Milky Way and continuously injects particles into the ISM. The actual distribution is lacunary and consists of point-like objects that release cosmic rays in a very short time. Because measurements have become very accurate, taking into account the discreteness of the sources in our description of CR propagation has become timely. In particular, the PAMELA (Adriani et al. 2011) and CREAM (Yoon et al. 2011) observations point towards an excess in the CR proton and helium fluxes with respect to the pure power-laws predicted with the continuous model. Very few analyses have been devoted to the myriad model so far. A challenging problem lies in the divergence of the variance of the CR primary flux. As shown recently by Blasi & Amato (2012), the second moment of the flux probability distribution function (PDF) is infinite. Although the variance can always be regularised in some way or another, this intriguing result threatens the conventional model of CR propagation, so far extremely successful. That is why we have thoroughly reinvestigated here how the discreteness of the sources may modify our vision of Galactic CR propagation.

To commence, the CR flux of primary cosmic rays obtained in the myriad model has a well-defined mean, identical to what is derived in the conventional approach by averaging the source distribution over space and time. This good point gives credit to the continuous model. We have then concentrated on the problem of the flux variance, which diverges in the myriad model. Although the central limit theorem cannot be used in that case, at least in its ordinary form, the PDF of the flux is well defined everywhere inside the Galactic magnetic halo. We have run Monte Carlo simulations with several thousands of different populations of point-like sources. For the first time, we have derived the quantiles of the CR proton flux as a function of proton energy. The quantiles turn out to be well-defined and can be used to compute 68% confidence intervals. In spite of an infinite variance, meaningful error bars for the fluxes of primary CR nuclei can be defined. These depend on the propagation parameters and tend to decrease with the number of sources implied in the signal. That is why the uncertainty bands of Fig. 6 widen as the vertical extension L of the CR diffusive halo shrinks.

We have so far focused on the flux at the Earth, but the same procedure can be applied throughout the Galaxy. It is extremely time consuming, however, and we defer the analytical derivation of the flux PDF to a subsequent publication. Such an investigation is crucial insofar as secondary species are produced by the interactions of primary CR nuclei on the ISM and the actual distribution of the latter matters. It would be interesting to examine the effect of the myriad hypothesis on the fluxes of antiprotons or positrons at the Earth as well as on the Galactic gamma ray diffuse emission. This emission has so far been calculated in the framework of the continuous model. An essential prediction is a gamma ray power law spectrum that traces the CR proton and helium energy distributions. According to the continuous hypothesis, we expect a gamma ray spectral index of  ~−2.8 ± 0.1 everywhere. Measurements by the H.E.S.S. collaboration (Aharonian et al. 2006) of the diffuse emission from the Galactic centre indicate a photon index of 2.29    ±    0.07    ±    0.02, significantly below the prediction of the continuous model. Fluctuations of the CR proton and helium spectra around their mean values have thus already been observed.

As regards the primary CR fluxes at the Earth, another important result is the use of a catalogue for the nearby and young sources that spans the so-called local region. These objects have been extracted from the Green and ATNF surveys and are based on the latest astronomical observations. The catalogue yields a contribution φ_cat to the total flux φ, whereas the sources extending beyond the local region contribute the complement φ_ext. This procedure provides a natural regularisation of the flux variance since it is based on observations. We have found that below a few hundreds of GeV, φ_cat is always negligible with respect to φ_ext, whatever the CR propagation parameters. Since the flux is dominated by distant or old sources, the continuous hypothesis applies and using the conventional CR propagation model is fully justified. We have succeeded in reconciling the myriad model with the continuous approach. Both should give identical results fo the low-energy CR fluxes at the Earth. At the TeV scale and beyond, the situation becomes more complicated. If the total number of objects that produce the flux is large with respect to the catalogue, the contribution φ_ext dominates and the continuous hypothesis is still valid at high energies. On the other hand, if the supernova explosion rate ν becomes lower than the canonical value of three events per century, or if the half-thickness L of the DH decreases, the catalogue becomes relatively important and may even dominate the flux above a few TeV. This is indeed the case for the black curve of Fig. 11 which corresponds to the min model and ν = 1century^-1. The CREAM data points lie inside the catalogue’s 68% error band. Notice how well the CR proton excess is naturally explained by local and young sources that have actually been detected. There is no need for a break in the injection spectrum, nor for a peculiar behaviour of the diffusion coefficient D with energy. It should even be possible to find a particuliar set of CR propagation parameters that would make the catalogue a natural explanation of both PAMELA and CREAM data. Since the local region alone is implied, these parameters ought not to be the same as for the bulk of the Galactic magnetic halo.

---

Acknowledgments

We thank Pasquale Blasi for very useful discussions. This work was supported by the Spanish MICINN’s Consolider-Ingenio 2010 Programme under grant CPAN CSD2007-00042. We also acknowledge the support of the MICINN under grant FPA2009-08958, the Community of Madrid under grant HEPHACOS S2009/ESP-1473, and the European Union under the Marie Curie-ITN program PITN-GA-2009-237920

References

Appendix A: Computation of the Green function

To solve Eq. (9) analytically, we need to simplify our description of the Milky Way DH and replace it by an infinite slab of half-thickness L with a gaseous disk in the middle at z = 0. Radial boundary conditions at r = R_gal are no longer implemented in the propagator. This simplification of the setup for CR propagation could be a problem if we were interested in the CR densities ψ close to the radial boundaries, at a distance of 20 kpc from the Galactic center. But our aim is to calculate these densities at the Earth, at a galactocentric distance of r_⊙ = 8.5 kpc, i.e. far from the radial boundaries. Furthermore, even though the propagator G_p is derived within the framework of an infinite diffusive slab, integrals on the sources of cosmic rays, such as relation (7), are still performed up to the radial boundaries at r = R_gal. We have checked that this procedure does not introduce any significant error on the CR fluxes at the Earth. For of a source term q_acc that is continuous in space and time, this approach yields results very close to those obtained with a method based on radial Bessel functions. With this simplified setup, the propagation of CR species becomes invariant under a translation along the horizontal axes x and y. The master Eq. (9) still needs to be solved along the vertical direction z, with the condition that G_p vanishes at the boundaries z =  ± L. The construction of the Green function for CR nuclei is inspired by the solution to the heat diffusion problem and has been given by many authors. The horizontal and vertical dependencies in G_p can be factored out by setting (A.1)where τ = t − t_S and ρ² = (x − x_S)²    +    (y − y_S)². The vertical function V_p is given by the expansion (A.2)where θ(τ) is the Heaviside function. The sum runs over even and odd eigenfunctions. The former can be expressed as (A.3)The corresponding even eigenwavevectors k_n are solutions of the equation (A.4)where k_w ≡ V_C/2D and k_s ≡ hΓ_p/D are typical wavevectors that account for the effects of Galactic convection and proton collisions on the ISM. Notice that at high energy, the dimensionless parameters k_wL and k_sL become very small as D increases. We then expect convection and spallations to have little effect with respect to space diffusion at high energy, as illustrated in Fig. A.1. The odd eigenfunctions are given by (A.5)where the odd eigenwavevectors are such that . This condition ensures that vanishes at z = 0 and is an odd function of the height z. The decay rates α_n and are related to the eigenwavevectors k_n and through (A.6)The eigenfunctions E_n(z) and may be understood as orthogonal vectors of a basis over which any vertical function that vanishes at z =  ± L can be expanded. They are normalized in such a way that (A.7)and (A.8)where the normalization factors C_n and are (A.9)The dimensionless parameter p is defined by 1/p = (k_w + k_s)L. It becomes infinite in the absence of Galactic convection and collisions on the ISM. In that limit, the eigenwavevectors k_n are equal to (n − 1/2)   π/L.

We still need to address a technical problem related to the behaviour of the vertical propagator (A.2) when t is close to t_S. We find that the variance of the CR proton flux is dominated by local and recent supernova explosions. In that case, the period of time τ that separates the injection of protons from their detection at the Earth becomes short. The expansion (A.2) needs to be pushed a long way up in terms of the number of eigenfunctions involved. When τ is short, the exponentially decreasing functions exp(−α_nτ) and are still close to unity unless n becomes very large. The numerical convergence of V_p requires then to sum expression (A.2) over an exceedingly large number of terms, hence a potential problem of CPU time. Because in that regime most of the expansion (A.2) is provided by high-order terms, Galactic convection and CR spallations become negligible with respect to diffusion. For these terms, the even eigenwavevectors k_n are very close to their pure diffusion values of (n − 1/2)   π/L, even if p is not very large. In this regime of basically pure diffusion, another solution to CR propagation is provided by the method of electrical images. To commence, if the half-thickness L of the slab is made infinite, we recover pure diffusion in infinite 3D space. In that case, the propagator is well-known, with a vertical contribution expressed as (A.10)As already discussed by Baltz & Edsjö (1999), the vertical boundaries of the DH can now be implemented by associating to each point-like source lying at height z_S the infinite series of its multiple images through the planes z =  ± L. The boundaries act as mirrors and give an infinite series of images of the source. The nth image is located at (A.11)and has a positive or negative contribution depending on whether n is an even or odd number. When the diffusion time τ = t − t_S is short, the 1D solution (A.10) is quite a good approximation. The relevant parameter is actually (A.12)and, in the regime where it is much larger than 1, i.e., for low values of τ, the propagation is insensitive to the vertical boundaries. When the diffusion parameter η decreases and the diffusion length becomes comparable to the half-thickness L of the DH, images of the point-like source need to be taken into account in the sum (A.13)When η is much smaller than 1, many terms must be taken into account in the previous sum. However, this regime corresponds to high values of the diffusion time τ for which expansion (A.2) converges very rapidly.

We have devised two complementary methods to calculate the propagator G_p. Depending on the value of the diffusion parameter η, we can use the relation (A.13) of electrical images (η ≥ 1) or the expansion (A.2) (η ≤ 1). Since the CR sources and the Earth are located inside the Galactic disk, we have set z_S and z equal to 0 in Fig. A.1. The vertical part of the Green function depends only on the time delay τ = t − t_S with V_p(τ) ≡ V_p(0,t    ←    0,t_S) = V_p(0,τ    ←    0,0). At low energy, Galactic convection and CR collisions on the ISM become important with respect to space diffusion. The method of electrical images is reliable only for low values of the diffusion time τ. Expansion (A.13) can still be used to calculate V_p in the regime where η becomes large. How large depends on the relative strength of the various CR propagation mechanisms. To obtain an impression of the range of validity over which electrical images can be used even if k_s and k_w are larger than 1/L, we have borrowed the MED set of CR propagation parameters from Donato et al. (2004). This benchmark configuration provides the best fit to the B/C measurements (see Table 1). The black short dashed line of Fig. A.1 corresponds to pure diffusion in infinite 3D space. The method of electrical images can only be used for pure CR diffusion and yields the red solid line. Expansion (A.13) has been calculated with n_sum = 100 images and is valid down to the low value of η ~ 0.016. Expansion (A.2) has also been pushed up to the 100-th term and converges even for η as large as 10³. It leads to the blue curves, which each correspond to a different proton kinetic energy T_p. Because diffusion takes over the other processes at very high-energy, the blue short dashed curve is completely superimposed on the red solid line. In that regime, expansions (A.2) and (A.13) yield the same result. As T_p decreases, convection and spallations come into play and the blue curves depart from their high-energy limit. Notice that the 100 GeV configuration is still fairly close to the high energy case. At 10 GeV, the blue dotted-long dashed curve differs noticeably from the red solid line. The pure diffusive regime is nevertheless obtained for η ≥ 100. Finally, the 1 GeV blue dotted-short dashed curve is significantly shifted towards higher values of η with respect to the pure diffusive case. To reliably calculate the proton propagator below  ~1 GeV, many terms need to be taken into account in expansion (A.2). This sum should be used up to high values of η before the electrical images provide an accurate result.

Appendix B: Probability distribution for the flux

In this section, we compute the high-φ behaviour of the probability density P(φ) for the flux generated by a point source drawn from a uniform spatial and temporal distribution. This high-φ tail of the distribution is the part that leads to the divergence of the variance. It is caused by the sources that are very close and very young, for which the effects of Galactic wind, escape and reacceleration are very small. We can neglect these effects, as long as we are only interested in the asymptotic behaviour of P(φ) at high φ.

The flux at distance r and at time t from a point source that instantly emits all its particles at r = 0 and t = 0 is given by

2D distribution of sources (thin disk)

Let us first consider sources with a given age t. The probability density that a source lies at a distance r is given by where R stands for the radius of the region containing the sources. We have so that For an age t, fluxes are in the interval which can be written as where W stands for the window function, which is equal to 1 in the interval and 0 outside. It is easily checked that The probability distribution for φ is obtained by where the upper bound is the lower bound t_min(φ) is a solution of If the age distribution p(t) is uniform between 0 and T, we have p(t) = 1/T so that This can be written as At high φ, we have

3D distribution of sources (thick disk)

If the sources are distributed in a volume instead of a surface, we still have but As we have As before, we can check that the total probability is 1. Let us compute We set x = φt^3/2/a, The integral is given by so that Finally, as φ_maxt^3/2 = a and φ_mint^3/2 = ae ^− R2/4Kt, The probability distribution is obtained by As before, p(t) = 1/T, where the bounds t_min and t_max are the same as before. We define y = tφ^2/3/a^2/3, with y_max = 1 and y_min(φ) solution of y = exp(− R²aφ^2/3/6Ky). For high values of φ, the lower bound vanishes and the integral does not depend on φ, which yields the final result

All Tables

All Figures