CONSTRAINTS ON COSMIC-RAY PROPAGATION MODELS FROM A GLOBAL BAYESIAN ANALYSIS

The goal of this paper is to determine constraints on the propagation model parameters (introduced below) from observed CR spectra and we adopt a Bayesian approach to parameter inference (see, e.g., Trotta 2008 for further details). Bayesian inference is based on the posterior probability distribution function (pdf) for the parameters, which updates our state of knowledge from the prior by taking into account the information contained in the likelihood. A recent application to CR propagation models is given in Maurin et al. (2010) and Putze et al. (2010). Denoting by Θ the vector of parameters one is interested in constraining and by D the available observations, Bayes Theorem reads

$$\begin{equation} P(\Theta |{\bf D}) = \frac{P({\bf D}| \Theta)P(\Theta)}{P({\bf D})}, \end{equation} \tag{ 3 }$$

where P(Θ|D) is the posterior distribution on the parameters (after the observations have been taken into account), $P({\bf D}| \Theta) = {\mathcal L}(\Theta)$ is the likelihood function (when considered as a function of Θ for the observed data D), and P(Θ) is the prior distribution, which encompasses our state of knowledge about the value of the parameters before we have seen the data. Finally, the quantity in the denominator of Equation (3) is the Bayesian evidence (or model likelihood), a normalizing constant that does not depend on Θ and can be neglected when interested in parameter inference. The evidence is obtained by computing the average of the likelihood under the prior (so that the right-hand side of Equation (3) is properly normalized)

$$\begin{equation} P({\bf D}) = \int P({\bf D}| \Theta) P(\Theta) d\Theta. \end{equation} \tag{ 4 }$$

The evidence is the prime quantity for Bayesian model comparison, which aims at establishing which of the available models is the "best" one, i.e., the one that fits the data best while being the most economical in terms of parameters, thus giving a quantitative implementation of Occam's razor (see, e.g., Trotta 2007). The evidence can also be used to assess the constraining power of the data (Trotta et al. 2008) and to carry out consistency checks between observables (Feroz et al. 2009).

Together with the model, the priors for the parameters which enter Bayes' theorem, Equation (3), must be specified. Priors should summarize our state of knowledge and/or our theoretical prejudice about the parameters before we consider the new data, and for the parameter inference step the prior for a new observation might be taken to be the posterior from a previous measurement (for model comparison issues the prior is better understood in a different way; see Trotta 2008).

The problem is then fully specified once we give the likelihood function for the observations (see Section 4.4 below). The posterior distribution P(Θ|D) is determined numerically by drawing samples from it and Markov Chain Monte Carlo (MCMC) techniques can be used for this purpose. In this paper, we use both Metropolis–Hastings MCMC and the MultiNest algorithm, which implements nested sampling and provides a higher efficiency, guarantees a better exploration of degeneracies and multimodal posteriors, and computes the Bayesian evidence as well (which is difficult to extract from MCMC methods).

As a test case for this study, we choose the diffusion– reacceleration model, which has been used in a number of studies utilizing the GALPROP code (e.g., Moskalenko et al. 2002; Strong et al. 2004; Ptuskin et al. 2006; Abdo et al. 2009a, and references therein). The source distribution is specified in Section 3.

In this model, the spatial diffusion coefficient is given by

$$\begin{equation} D_{xx} = \beta D_{0} \left(\frac{\rho }{\rho _0}\right)^\delta, \end{equation} \tag{ 5 }$$

where D₀ is a free normalization at the fixed rigidity ρ₀ = 4 × 10³ MV. The power-law index is δ = 1/3 for Kolmogorov diffusion (see Section 3), but we take it as a free parameter for the purposes of this study. Fitting the B/C ratio below 1 GeV in reacceleration models is known to require large values of v_Alf. In these models, a break in the injection spectra is required to compensate for the large bump in the propagated spectra at low energies/nucleon. Therefore, the CR injection spectrum is modeled as a broken power law, with index below (−ν₁) and above (−ν₂) the break as free parameters, but with the location of the break fixed at a rigidity 10⁴ MV. The other free model parameters are v_Alf, the halo size z_h, and the normalization of the propagated CR proton spectrum at 100 GeV N_p, for a total of seven free model parameters, as summarized in Table 1. Other models discussed in the literature may be able to reproduce the B/C ratio without a break in the injection spectra, but the present paper is mainly intended as a presentation of the method and we defer a comprehensive study of different possibilities to a forthcoming paper.

Table 1. Summary of Input Parameters and Prior Ranges

Quantity	Symbol	Prior Range	Prior Type
Diffusion model parameters Θ
Diffusion coefficienta (1028 cm2 s−1)	D0	[1, 12]	Uniform
Rigidity power-law index	δ	[0.1, 1.0]	Uniform
Alfvén speed (km s−1)	vAlf	[0, 50]	Uniform
Diffusion zone height (kpc)	zh	[1.0, 10.0]	Uniform
Nucleus injection index below 104 MV	ν1	[1.50, 2.20]	Uniform
Nucleus injection index above 104 MV	ν2	[2.05, 2.50]	Uniform
Proton normalization (10−9 cm2 sr−1 s−1 MeV−1)	Np	[2, 8]	Uniform
Experimental nuisance parameters
Modulation parameters ϕ (MV)			Gaussian priorb
HEAO-3	$m_{\rm {\rm \it HEAO{\scriptsize\hbox{-}}3}}$	[420, 780]	${\mathcal N}(600, 60)$
ACE	mACE	[175, 475]	${\mathcal N}(325, 50)$
CREAM	mCREAM	[420, 780]	${\mathcal N}(600, 50)$
ISOMAX	$m_{\rm {\rm \it ISOMAX}}$	[370, 490]	${\mathcal N}(430, 20)$
ATIC-2	$m_{\rm {\rm \it ATIC{\scriptsize\hbox{-}}2}}$	0	Fixed (no modulation)
Variance rescaling parameters (j = 1, ..., 5)	log τj	[ − 1.5, 0.0]	Uniform on log τj

Notes. ^aAt ρ = 4 × 10³ MV. ^bWe use the notation $\mathcal {N}(\mu, \sigma)$ to indicate a Gaussian distribution of mean μ and standard deviation σ.

Download table as:  ASCIITypeset image

The nuclear chain used starts at ²⁸Si and proceeds down to protons. The source abundances of nuclei 6 ⩾ Z ⩾ 14 have an important influence on the B/C and ¹⁰Be/⁹Be ratios used in this study. In our analysis, they are fixed at values determined for ACE data at a few 100 MeV nucleon⁻¹ (Moskalenko et al. 2008), but the values are assumed to hold also at the GALPROP normalization energy of 100 GeV nucleon⁻¹. The adopted relative source abundances of the most abundant isotopes (for particle flux in cm⁻² s⁻¹ (MeV nucleon⁻¹)⁻¹) are ⁴He: 7.199 × 10⁴, ¹²C: 2819, ¹⁴N: 182.8, ¹⁶O: 3822, ²⁰Ne: 312.5, ²²Ne: 100.1, ²³Na: 22.84, ²⁴Mg: 658.1, ²⁵Mg: 82.5, ²⁶Mg: 104.7, ²⁷Al: 76.42, and ²⁸Si: 725.7. These values are relative to the proton normalization N_p for a proton source abundance 1.06 × 10⁶, but this is only formal since the antiprotons, secondary positrons, and gamma rays were computed from an independent fit to proton and He data. N_p is used only to normalize C and O to the data, via the ratios given above (other data like N are not used explicitly).

Special handling is required to treat the solar modulation of the propagated CR spectra, for which we introduce an extra nuisance parameter for each of the data set we consider. The motivation and choice of the Gaussian priors, with mean and standard deviation as given in Table 1, is described in Section 4.3. In addition, we also introduce a set of parameters τ designed to mitigate the possibility that the fit be dominated by undetected systematic errors in the data, as explained in detail in the next section. Overall, we thus fit a total of 16 free parameters, including 7 model parameters, 4 modulation parameters, and 5 observational variance rescaling factors.

For demonstration of the method we use the most accurate CR data sets available preferably taken near solar minimum.¹²

The B/C ratio is well measured by a number of space- and balloon-borne missions. The HEAO-3 data (Engelmann et al. 1990) remain the most accurate to date in the energy range 0.6–35 GeV nucleon⁻¹ and have been recently confirmed by PAMELA (R. Sparvoli 2010, private communication). At higher energies, from 30 GeV nucleon⁻¹–1 TeV nucleon⁻¹, we use ATIC-2 (Panov et al. 2008) and CREAM-1 data (Ahn et al. 2008). At low energies, the Voyager 1 and 2 (Lukasiak et al. 1999), Ulysses (Duvernois et al. 1996), and ACE (de Nolfo et al. 2006) data agree with each other, while the ACE data (50–200 MeV nucleon⁻¹) have the smallest statistical error. Therefore, we use the ACE measurements corresponding to the solar minimum conditions (George et al. 2009).

The ¹⁰Be/⁹Be ratio is most accurately measured (70–145 MeV nucleon⁻¹) by ACE (Yanasak et al. 2001), which we include in our fit. Those measurements are in agreement with Voyager 1 and 2 (Lukasiak et al. 1999), and Ulysses (Connell 1998) data. At higher energies (per nucleon), there are only two data points by ISOMAX (Hams et al. 2004) with very large error bars, which we also include in the fit.

We also use the carbon and oxygen spectra as measured by ACE at the solar minimum (George et al. 2009) and by HEAO-3 (Engelmann et al. 1990).

As mentioned above, a very important issue is the treatment of the heliospheric modulation. We fit to the CR data in the whole energy range from some 10 MeV nucleon⁻¹ to TeV energies. However, a comparison of calculated CR spectra, the elemental and isotopic ratios with low-energy data (below ∼ 20 GeV nucleon⁻¹) measured inside the heliosphere requires care as the calculated spectra depend significantly on the treatment of the heliospheric modulation. As mentioned in Section 2, the modulation can be realistically treated with full three-dimensional models, but application of such models to the current study does not seem feasible since the number of free parameters and the computing requirements would considerably increase. Currently, it is only possible to use a simple force-field approximation (Gleeson & Axford 1968), which is characterized by the value of the modulation potential. However, directly using the modulation potentials from different experiments is problematic because they cannot be interpreted independently from the experiments themselves. The derived values depend on the choices of interstellar spectra used for their analyses, which differ from experiment to experiment (and are sometimes not provided).

To deal with this type of uncertainty, instead of fixing a collection of a priori values to the modulation potential, we allow some flexibility to the fits and include the modulation potentials as free nuisance parameters in our inference, with one free parameter per experiment (i.e., ACE, HEAO-3, ISOMAX, and CREAM-1). To avoid unphysical/implausible values, we adopt Gaussian priors with mean and standard deviation as given in Table 1, which are motivated by the estimated ballpark values of the modulation by the experimentalists. Note that no modulation parameter is given for ATIC-2 as we only use high-energy data for that experiment and modulation is not relevant.

For a given set of the CR model parameters Θ and the modulation potential parameters ϕ (where ϕ = {ϕ₁, ..., ϕ₄}, with a different choice of the modulation potential for each data set employed) we can compute via GALPROP the ensuing CR spectrum, as a function of energy, Φ_X(E, Θ, ϕ) for species X. Assuming Gaussian noise on the observations, we take the likelihood function for each observation of species X at energy E_i to be of the form

$$\begin{equation} P\big(\hat{\Phi }_X^{ij} | \Theta, \phi\big) = \frac{1}{\sqrt{2\pi }\sigma _{ij}} \exp \left(-\frac{1}{2}\frac{\big(\Phi _X(E_i, \Theta, \phi) - \hat{\Phi }_X^{ij} \big)^2}{\sigma _{ij}^2} \right), \end{equation} \tag{ 6 }$$

where Φ_X(E_i, Θ, ϕ) is the prediction from the CR propagation model for species X at energy E_i, $\hat{\Phi }_X^{ij}$ is the measured spectrum, and σ_ij is the reported standard deviation. The sub/superscript i runs through the data points within each of the data sets j. We assume bins to be independent, so that the full likelihood function is given by the product of the terms of the form given above:

$$\begin{equation} P({\bf D}| \Theta,\phi) = \prod _{j=1}^{5}\prod _{i=1}^{N_j} P\big(\hat{\Phi }_X^{ij} | \Theta, \phi\big). \end{equation} \tag{ 7 }$$

However, a careful analysis of a plot of the data points for each CR species reveals that there are fairly strong discrepancies between different data sets. This might point to either an underestimation of the actual experimental error bars or to undetected systematic errors between data sets. If some or all of the reported error bars are significantly underestimated, this would lead to a handful of data points incorrectly dominating the global fit, introducing systematic bias in the reconstructed value of the parameters. To mitigate against undetected systematics, we follow the procedure described in, e.g., Barnes et al. (2003). For each data set we introduce in the likelihood a parameter τ_j (j = 1, ..., 5), whose function is to rescale the variance of the data points in order to account for possible systematic uncertainties. Therefore, Equation (6) is modified:

$$\begin{equation} P\big(\hat{\Phi }_X^{ij} | \Theta, \phi, \tau\big) {=} \frac{\sqrt{\tau _j}}{\sqrt{2\pi }\sigma _{ij}} \exp \left({-}\frac{1}{2}\frac{ \left(\Phi _X(E_i, \Theta, \phi) - \hat{\Phi }_X^{ij} \right)^2}{\sigma _{ij}^2/\tau _j} \right)\!. \end{equation} \tag{ 8 }$$

The role of the set of parameters τ = {τ₁, ..., τ₅}, which we call "error bar rescaling parameters," is to allow for the possibility that the error bars reported by each of the experiments underestimate the true noise. We then add τ to the set of parameters Θ and sample over it, too, thus allowing the data themselves to decide whether there are systematic discrepancies in the reported error bars. A value τ_j < 1 means that the data prefer a systematically larger value for the errors for data set j. Note that τ_j not only appears in the exponential of the Gaussian in Equation (8), but also in the pre-factor, which, being proportional to $\sqrt{\tau _j}$, ensures that τ_j never attains a value arbitrarily close to 0 (implying infinite error bars). Furthermore, the variance scaling parameters τ also take care of all aspects of the model that are not captured by the reported experimental error: this also includes theoretical errors (i.e., the model not being completely correct), errors in the cross section normalizations, etc.

The full posterior distribution for the CR propagation model parameters Θ, the variance rescaling parameters τ, and the modulation parameters ϕ is given by

$$\begin{equation} P(\Theta, \phi, \tau, | {\bf D}) \propto P({\bf D}| \Theta, \phi, \tau) P(\Theta)P(\tau)P(\phi), \end{equation} \tag{ 9 }$$

where the likelihood P(D|Θ, τ, ϕ) is given by Equations (7) and (8).

The priors P(Θ), P(ϕ), and P(τ) in Equation (9) are specified as follows. We take the prior on a set of model parameters, P(Θ), to be uniform on Θ with ranges as given in Table 1. As shown below, the posterior is reasonably well constrained and close to Gaussian for Θ, hence we expect our results to be fairly independent of the prior choice.

This conclusion is strengthened by the inspection of the profile likelihood, which is obtained from our samples by maximizing the value of the likelihood along the hidden dimensions rather than integrating over the posterior. The profile likelihood statistic is independent of the priors, provided the parameter space has been sampled with sufficient resolution, and thus it constitutes a cross-check for the presence of large volume effects coming from the priors. Such volume effects are typically important when the priors play a major role in the inference, while they are usually negligible when the posterior is dominated by the likelihood (in which case the prior influence is minimal; see, e.g., Trotta et al. 2008 for an illustration). We have found the profile likelihood to be in excellent agreement with the posterior pdf presented below, and therefore we do not consider it further in our results below. This means that the prior influence is small and our results can be considered to be robust with respect to reasonable changes in priors.

Regarding the modulation parameters, we adopt a Gaussian prior on each of them, informed by the values reported by each experiment (see Table 1), in order to avoid physically implausible values. The description of the experimental CR data sets can be found in Section 4.3.

The τ_j are scaling parameters in the likelihood, and thus the appropriate prior is given by the Jeffreys' prior, which is uniform on log τ_j (see Barnes et al. 2003 or Jaynes & Bretthorst 2003 for a justification). Therefore, we adopt the proper prior

$$\begin{equation} P(\log \tau _j) = \left\lbrace \begin{array}{@{}l l}2/3 & {\rm for } -3/2 \le \log \tau _j \le 0 \\ \ms 0 &{\rm otherwise} \end{array} \right. \end{equation} \tag{ 10 }$$

that corresponds to a prior on τ_j of the form

$$\begin{equation} P(\tau _j) \propto \tau _j^{-1}. \end{equation} \tag{ 11 }$$

The inclusion in our analysis of the nuisance parameters ϕ and τ (which are then marginalized over; see Equation (17)) has two main effects on our inference about the CR parameters of interest, Θ. First, it increases the robustness of our fit, since the nuisance parameters account for potential systematic effects in the data (τ) and approximately capture the impact of solar modulation on the measurements (ϕ). Second, it makes our CR constraints more conservative, since our marginalized errors on Θ fully account for all possible values of the nuisance parameters compatible with the data.

A powerful and efficient alternative to classical MCMC methods has emerged in the last few years in the form of the so-called nested sampling algorithm, invented by John Skilling (Skilling 2004, 2006; Feroz & Hobson 2008; Feroz et al. 2009). Although the original motivation for nested sampling was to compute the evidence integral of Equation (4), the recent development of the MultiNest algorithm (Feroz & Hobson 2008; Feroz et al. 2009) has delivered an extremely powerful and versatile algorithm that has been demonstrated to be able to deal with extremely complex likelihood surfaces in hundreds of dimensions exhibiting multiple peaks.

As samples from the posterior are generated as a by-product of the evidence computation, nested sampling can also be used to obtain parameter constraints in the same run as computing the Bayesian evidence. In addition, multi-modal nested sampling exhibits an efficiency that is almost independent of the dimensionality of the parameter space being explored, thus beating the "curse of dimensionality."

The essential element of nested sampling is that the multi-dimensional evidence integral is recast into a one-dimensional integral. This is accomplished by defining the prior volume X as dX ≡ P(Θ)dΘ so that

$$\begin{equation} X(\lambda) = \int _{{\mathcal L}(\Theta)>\lambda } P(\Theta) d \Theta, \end{equation} \tag{ 12 }$$

where the integral is over the parameter space enclosed by the iso-likelihood contour ${\mathcal L}(\Theta) = \lambda$. So X(λ) gives the volume of parameter space above a certain level λ of the likelihood. Then the Bayesian evidence, Equation (4), can be written as

$$\begin{equation} P({\bf D}) = \int _0^1 {\mathcal L}(X) d X, \end{equation} \tag{ 13 }$$

where ${\mathcal L}(X)$ is the inverse of Equation (12). Samples from ${\mathcal L}(X)$ can be obtained by uniformly drawing samples from the likelihood volume within the iso-contour surface defined by λ. The one-dimensional integral of Equation (13) can be obtained by simple quadrature, thus

$$\begin{equation} P({\bf D}) \approx \sum _i {\mathcal L}(X_i) W_i, \end{equation} \tag{ 14 }$$

where the weights are $W_i = \frac{1}{2}(X_{i-1} - X_{i+1})$ (see Skilling 2004, 2006; Feroz & Hobson 2008; Feroz et al. 2009; Mukherjee et al. 2006, for details). It has been shown in the context of CMB data analysis in cosmology and in supersymmetry phenomenology studies that this technique reduces the number of likelihood evaluations by over an order of magnitude with respect to conventional MCMC.

In this paper, we adopt the publicly available MultiNest algorithm (Feroz & Hobson 2008), as implemented in the SuperBayeS code (Trotta et al. 2008; Ruiz de Austri et al. 2006), that we have interfaced with the GALPROP code. First, we performed an exploratory scan with MultiNest, adopting 4000 live points in our 16-dimensional parameter space, with the aim of scouting the structure and degeneracy directions of the posterior distribution. An important characteristic of MultiNest, which sets it apart from conventional MCMC methods, is its ability to sample reliably from multi-modal distributions. Therefore, it is highly desirable to employ MultiNest to perform scans of parameter spaces that have not been investigated before, as MultiNest will make it less likely to miss important substructure in the probability distribution when the latter is multi-modal.

Our exploratory MultiNest scan gathered ∼10⁵ samples from the posterior, with an overall efficiency of about 10% and a total computational effort of ≈13 CPU years. This initial scan revealed a well-behaved, unimodal distribution over the prior ranges given in Table 1. We then computed the parameter-set covariance matrix and adopted this as a Gaussian proposal distribution for a conventional Metropolis–Hastings MCMC scan. Since the proposal distribution was well matched to the posterior, our MCMC scan reached an efficiency of ∼15%, and we built 10 parallel chains with 14,000 samples each (after burn-in), for a total of 1.4 × 10⁵ samples. We checked that the Gelman & Rubin mixing criterion (Gelman & Rubin 1992) is satisfied for all of our parameters (i.e., R ≪ 0.1, where R is the inter-chain variance divided by the intra-chain variance).

We verified that the posterior distribution obtained with MCMC was in excellent agreement with the one obtained from MultiNest, which validates our results as the two sampling schemes are completely different. The results presented in this paper are obtained from the MCMC run, which allows for slightly smoother posterior distributions as it contains 40% more samples than the MultiNest scan.¹³

In order to keep the computational cost within reasonable limits, we carry out our MultiNest and MCMC runs assuming a relatively coarse spatial and energy/nucleon grid for the CR propagation. For each point in parameter space, we adopt δz = 0.2 kpc, δR = 1 kpc, and δE = 1.4. With these parameters, one likelihood evaluation takes approximately 15 s on an eight-way 2.4 GHz Opteron CPU machine. We then reprocessed the MCMC samples using importance sampling, decreasing the spacing of the spatial grid by a factor of two in each direction while increasing the energy resolution to δE = 1.1. The computational cost per likelihood evaluation is increased by a factor of ∼20, but allows a more precise computation of the CR spectra. The statistical distribution is adjusted accordingly, thus obtaining a posterior distribution that is close to what would have been obtained by running the scan at the higher resolution initially. The advantage of using posterior sampling in this context is that the resampling of the points can be done in a massively parallel way (using 800 CPUs, the resampling step takes only a few hours). Although the best-fit point shifts somewhat after importance sampling, we have verified that the bulk of our probability distributions remain stable. We therefore conclude that the results presented here are robust with respect to increases in the spatial and energy resolution of the scan.

Once a sequence {Θ⁽⁰⁾, Θ⁽¹⁾, ..., Θ^(M−1)} of samples from the posterior pdf has been gathered, obtaining Monte Carlo estimates of expectations for any function of the parameters is straightforward. For example, the posterior mean is given by (〈 · 〉 denotes the expectation value with respect to the posterior)

$$\begin{equation} \langle \Theta \rangle = \int P(\Theta |{\bf D})\Theta d\Theta \approx \frac{1}{M} \sum _{t=0}^{M-1} \Theta ^{(t)}, \end{equation} \tag{ 15 }$$

where the second equality follows because the samples Θ^(t) are generated from the posterior by construction. In general, the expectation value of any function of the parameters f(Θ) is obtained as

$$\begin{equation} \langle f(\Theta) \rangle \approx \frac{1}{M}\sum _{t=0}^{M-1} f(\Theta ^{(t)}). \end{equation} \tag{ 16 }$$

It is useful to summarize the results of the inference by giving the one-dimensional marginal probability for the jth element of Θ, Θ_j, obtained by integrating out all other parameters from the posterior:

$$\begin{equation} P(\Theta _j|{\bf D}) = \int P(\Theta |{\bf D}) d \Theta _1 \ldots d \Theta _{j-1} d \Theta _{j+1} \dots d \Theta _{n}, \end{equation} \tag{ 17 }$$

where P(Θ_j|D) is the marginal posterior for the parameter Θ_j. From the posterior samples (obtained either by MCMC or MultiNest) it is straightforward to obtain and plot the marginal posterior on the left-hand side of Equation (17): since samples are drawn from the full posterior by construction, their density reflects the value of P(Θ|D). It is then sufficient to divide the range of Θ_j into a series of bins and count the number of samples falling within each bin, ignoring the coordinates values Θ_i (for i ≠ j). The two-dimensional posterior is defined in an analogous fashion. The one-dimensional, two-tail symmetric α% credible region is given by the interval (for the parameter of interest) within which α% of the samples are found, obtained in such a way that a fraction (1 − α)/2 of the samples are outside the interval on either side. In the case of a one-tail upper (lower) limit, we report the value of the quantity below (above) where α% of the samples are found.

In Figures 1 and 2, we show one-dimensional marginalized posterior probabilities for the propagation and nuisance parameters of the model. The red cross represents the best fit, the vertical thin line the posterior mean, and the horizontal bar the 68% and 95% error ranges (yellow/blue, respectively). Two-dimensional constraints on some parameter combinations are presented in Figure 3. The best-fit point and posterior ranges are summarized in Table 2. For reference the galdef parameter definition files with the best-fit parameter values presented in this study are available as a supplementary material to this paper. These give precise definitions of the model used, which can be reproduced as required. The final MCMC chains from which Figures 1–3 were produced are also available.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We see that δ and v_Alf are quite well constrained, with the posterior mean δ = 0.31 ± 0.02 being very close to the canonical value of 1/3 for Kolmogorov diffusion. The Alfvén speed, v_Alf = 38.4 ± 2.1 km s⁻¹, is higher than in earlier studies, but this is dictated by the fit to the ACE data on the B/C ratio at low energies. The posterior intervals on the values of D₀ and z_h are fairly large, D₀ = (8.32 ± 1.46) × 10²⁸ cm² s⁻¹ and z_h = 5.4 ± 1.4 kpc. The typical value of 4 kpc adopted in many studies is at the lower end of the viable range, but still within the 95% interval, z_h∈ [3.2, 8.6] kpc. We can see from the D₀ versus z_h panel in Figure 3 that the diffusion coefficient and the halo size are positively correlated, as expected.

Other parameters exhibit less pronounced correlations. The injection indexes ν₁ and ν₂ are tightly constrained and almost uncorrelated (Figure 3), but this reflects the fact that the position of the injection spectral break is fixed in this analysis, so that the indices are determined by δ and v_Alf with their narrow ranges. The value of the injection index ν₂ = 2.38 ± 0.04 provides a consistency check as the value of the sum ν₂ + δ should be close to the spectral indices of directly measured carbon and oxygen spectra ∼2.70, and indeed we find that ν₂ + δ = 2.69 ± 0.05. Comparison to the value of the injection index ν₁ = 1.92 ± 0.04 shows that the spectral break required is 0.46 ± 0.05. The pdf for the spectral break is plotted in the bottom right panel of Figure 1. While at face value the break appears very statistically significant, it has to be kept in mind that the value found is dependent on the break energy, which was fixed in this analysis. Future analyses will allow more freedom in the form of the spectrum.

In general, there is remarkable agreement between the "by-eye" fitting in the past (e.g., Strong & Moskalenko 1998, 2001; Moskalenko et al. 2002; Ptuskin et al. 2006) and the parameter constraints found using the refined Bayesian inference analysis described in this paper. The posterior mean values of the diffusion coefficient D₀ = (8.32 ± 1.46) × 10²⁸ cm² s⁻¹ at 4 GV and the Alfvén speed v_Alf = 38.4 ± 2.1 km s⁻¹ are also in fair agreement with earlier estimates of 5.73 × 10²⁸ cm² s⁻¹ and 36 km s⁻¹ (Ptuskin et al. 2006), respectively. The posterior mean halo size is 5.4 ± 1.4 kpc, also in agreement with our earlier estimated range z_h = 4–6 kpc (Strong & Moskalenko 2001), although our best-fit value of z_h = 3.9 kpc is somewhat lower, due to the degeneracy between D₀ and z_h. However, the well-defined posterior intervals produced by the present study are significantly more valuable than just the best-fit values themselves as they provide an estimate of associated theoretical uncertainties and may point to a potential inconsistency between different types of data.

We now assess the quality of our best-fit model. Define the χ² as

$$\begin{equation} \chi ^2 \equiv \sum _{j=1}^5 \sum _{i=1}^{N_j} \frac{\big(\Phi _X(E_i, \Theta, \phi) - \hat{\Phi }_X^{ij} \big)^2}{\sigma _{ij}^2/\tau _j}, \end{equation} \tag{ 18 }$$

i.e., we compute the χ² using the rescaled error bars for the data points (note that the χ² ≠ −2log P(D|Θ, ϕ, τ), i.e., the χ² is not minus twice the log-likelihood because of the pre-factor containing τappearing in Equation (8)). There are N = 76 total data points and M = 16 fitted parameters, including both the modulation and the error rescaling parameters. Therefore, the number of degrees of freedom (dof) is 60, and for the best-fit model we find χ² = 69.3, which leads to a reduced chi-squared χ²/dof = 68/60 = 1.15. This is not surprising, since by construction the error bar rescaling parameters, τ, are adjusted dynamically during the global fit to achieve this. A more detailed breakdown of the contribution to the total χ² by data set is given in Table 3.

The predictions for the fitted CR spectra of the best-fit model parameters are shown in Figures 4–6, including an error band delimiting the 68% and 95% probability regions. The species shown are B/C and ¹⁰Be/⁹Be ratios, and the spectra of carbon and oxygen. In each plot, we show the spectrum modulated with the potential corresponding to our best-fit parameters from our global fits for each of the data sets employed. We also show the data sets, each with error bars enlarged by the best-fit value of our scaling parameters, τ, as given in Table 2. The yellow/blue band delimits regions of 68% and 95% probability, and is modulated according to the potential given in the each panel. We emphasize that the power of our statistical technique is such that we can, for the first time, provide not only a best-fit model but also an error band with a well-defined statistical meaning.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In order to better visualize the comparison of our best-fit model to the fitted data, we plot in the bottom part of each panel the best-fit residuals, i.e., the difference between data and best-fit model, divided by the experimental error bar (enlarged by the correct error scaling parameter):

$$\begin{equation} {\mathcal R}_{ij} = \frac{ \hat{\Phi }_X^{ij} - \Phi _X(E_i, \Theta, \phi)}{\sigma _{ij}/\sqrt{\tau _j}}. \end{equation} \tag{ 19 }$$

Visual inspection of the residuals for the B/C and the ¹⁰Be/⁹Be ratios (see Figures 4 and 5) shows that our best-fit model gives an excellent fit to those data, with the distribution of the residuals approximately symmetric around 0. This indicates that there is no systematic bias in our best fit. The contribution to the overall χ² from those data sets is, if anything, smaller than would be expected statistically: Table 3 indicates that each datum contributes about ∼0.3 units to the χ². This could point to a degree of overfitting, or to our error bar rescaling parameters being too small. However, the origin of this slight overfitting becomes clear when one considers the oxygen and carbon spectra, and their residuals (Figures 6). Residuals here are significantly larger, especially at low energies, E < 3 GeV, and the average contribution to the total χ² by each datum is much larger, of order ∼1.4; see Table 3. Therefore, the error bars on carbon and oxygen seem to require enlargement in order for our model to provide a good fit. Note from the shape of the residuals in Figure 6 that there is no indication of systematic bias in our fit, i.e., the enlargement of the error bars does not come about because the model cannot reproduce the data, rather, because the data themselves seem to show an amount of scatter that is incompatible with a smooth spectrum (unless the error bars are enlarged sufficiently).

As a consequence, we can conclude that it is the carbon and oxygen spectra that are driving the value of the error bars rescaling parameters for HEAO-3 and ACE. Therefore, the error bars are correspondingly enlarged in the B/C and ¹⁰Be/⁹Be spectra, since we are using one single error bar rescaling parameter for each experiment. This results in much smaller residuals for the latter spectra. Introducing a larger number of error bar rescaling parameters, one for each experiment and for each CR species, and fitting them independently could resolve this issue. Then, the rescaling will be less important for B/C and ¹⁰Be/⁹Be, leading to tighter constraints from those data sets.

We also use our best-fit model to calculate secondary antiprotons, electrons, positrons, and diffuse γ-rays which were not fitted, but provide a useful consistency check. For this calculation the spectra of CR protons and He were adjusted to the BESS data (Sanuki et al. 2000). Secondary antiprotons were calculated using the same formalism as in Moskalenko et al. (2002). Calculation of secondary electrons, positrons, and diffuse γ-rays is described in Section 3. The spectra of secondaries are shown in Figures 7–9. The primary electron injection spectrum is based on fitting to pre-Fermi electron data (conventional model; Strong et al. 2004; Ptuskin et al. 2006) and is parameterized as a broken power law with indices 1.6/2.5 below/above 4 GeV and a steepening (index 5) above 2 TeV, and normalized to the Fermi data at 25 GeV (Ackermann et al. 2010). The plot for the CR electron and positron spectrum is shown separately in Figure 7 and compared to relevant measurements, while the plot for total leptons (electrons plus positrons) is displayed in the left panel of Figure 8. Even though the total electron and diffuse emission data were not fitted, they agree well with our best-fit model predictions. The positron fraction, shown in the right panel of Figure 8, does not agree with the PAMELA data (Adriani et al. 2009), but this was expected since secondary positron production in the general ISM is not capable of producing an abundance that rises with energy.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Antiprotons, shown in Figure 9, also not fitted, present an interesting example where the intensity at a few GeV is significantly underpredicted by the reacceleration models. As has been already shown (Moskalenko et al. 2002, 2003), the antiproton flux measurements by BESS taken during the last solar minimum, 1995–1997 (Orito et al. 2000), are inconsistent with reacceleration models at the 40% level at about 2 GeV, while the stated measurement uncertainties in this energy range are 20%. The reacceleration models considered are conventional models, based on local CR measurements, Kolmogorov diffusion, and uniform CR source spectra throughout the Galaxy. Models without reacceleration that can reproduce the antiproton flux, however, fall short of explaining the low-energy decrease in the secondary/primary nuclei ratio. To be consistent with both, the introduction of breaks in the diffusion coefficient and the injection spectrum is required, which may suggest new phenomena in particle acceleration and propagation. Inclusion of a local primary component at low energies, possibly associated with the Local Bubble, could reconcile the data (Moskalenko et al. 2003).

Figure 9 shows that the reacceleration model underproduces antiprotons in the GeV range by ∼ 30% also compared to the PAMELA data taken during the current solar minimum (Adriani et al. 2010). On the other hand, high-energy PAMELA data agree well with the predictions, which may suggest that the excess over the model predictions at low energies can be associated with the solar modulation. Since reacceleration is also most important below a few GeV, at present it does not appear possible to distinguish these effects. However, a more systematic analysis which includes evaluation of other propagation models may help and will be carried out in a future work.

Most other analyses have used analytical propagation models, in particular Donato et al. (2002) and Maurin et al. (2002). The most recent reported results from this approach are in Maurin et al. (2010) and Putze et al. (2010), which used χ² and MCMC techniques to analyze a wide range of semi-analytical models. The nearest case to ours is their diffusion/reacceleration model where they found δ = 0.23–0.24, v_Alf = 70–80 km s⁻¹, and z_h = 4–6 kpc. Better fits are found with convection included as well, with a smaller v_Alf, but this requires a very large δ = 0.8–0.9. In contrast, we find a more plausible range of values for the diffusion/reaccleration model (however, we do not consider convection here).

We also note that they assume an injection spectrum that is a single power law in rigidity, which also includes a factor of $\beta ^{\eta _S}$, ∼ $\beta ^{\eta _S}\rho ^{-\nu }$ (Putze et al. 2010). In the case of proton and He injection spectra in the reacceleration model (their Model II) they use η_S ≈ 1. This is equivalent to our break in the injection spectrum provided the diffusion coefficient has the form of Equation (5). For heavier nuclei, C to Fe, their injection spectrum has η_S ≈ 2, which compensates for the large values of the v_Alf in their fits. Therefore, their best-fit parameters (e.g., v_Alf) are not equivalent to ours and are dependent on a particular choice of the injection spectrum. Their conclusion that they do not need a break in the injection spectrum is thus significantly overstated as they need an ad hoc factor $\beta ^{\eta _S}$ with index η_S different for different groups of nuclei.