A Markov Chain Monte Carlo for Galactic Cosmic Ray physics: I. Method and results for the Leaky Box Model

0.88 MB·English

by  A. Putze

Astronomy & Astrophysics manuscript no. mcmcI c ESO 2009
(cid:13)
January 21, 2009
A Markov Chain Monte Carlo for Galactic cosmic-ray physics
I. Method and results for the Leaky-Box Model
A. Putze1, L. Derome1, D. Maurin2, L. Perotto1,3, and R. Taillet4,5
1 Laboratoire de PhysiqueSubatomiqueet de Cosmologie lpsc, 53 avenuedes Martyrs, Grenoble, 38026, France
9
2 Laboratoire de PhysiqueNucl´eaire et desHautes Energies lpnhe,Tour 33, Jussieu, Paris, 75005, France
0
3 Laboratoiredel’acc´el´erateurlin´eairelal,Universit´eParis-Sud11,Bˆatiment200,B.P.34,91898OrsayCedex,France
0
4 Laboratoire de PhysiqueTh´eorique lapth,Chemin de Bellevue BP 110, 74941 Annecy-le-Vieux,France
2
5 Universit´ede Savoie, Chamb´ery,73011, France

n
a Received / Accepted
J
ABSTRACT
1
2
Context. Propagation of charged cosmic-rays in the Galaxy depends on the transport parameters, whose number can

be large depending on the propagation model under scrutiny. A standard approach for determining these parameters
]
h is a manual scan, leading to an inefficient and incomplete coverage of theparameter space.
p Aims. In analyzing the data from forthcoming experiments, a more sophisticated strategy is required. An automated
- statistical toolisused,whichenablesafullcoverageoftheparameterspaceandprovidesasounddetermination ofthe
o transport and source parameters. The uncertainties in these parameters are also derived.
r Methods.WeimplementaMarkovChainMonteCarlo(MCMC),whichiswellsuitedtomulti-parameterdetermination.
t
s Itsspecificities (burn-inlength, acceptance, and correlation length) are discussed in thecontext of cosmic-ray physics.
a Itscapabilities and performances are explored in thephenomenologically well-understood Leaky-Box Model.
[ Results. From a technical point of view, a trial function based on binary-space partitioning is found to be extremely
2 efficient, allowing a simultaneous determination of up to nine parameters, including transport and source parameters,
v such as slope and abundances. Our best-fit model includes both a low energy cut-off and reacceleration, whose values
7 are consistent with those found in diffusion models. A Kolmogorov spectrum for the diffusion slope (δ = 1/3) is ex-
cluded.The marginalised probability-densityfunction for δ and α (theslope of thesourcespectra) are δ≈0.55−0.60
3
and α≈2.14−2.17, dependingon thedataset used and thenumberof free parameters in thefit. Allsource-spectrum
4
parameters (slope and abundances) are positively correlated among themselves and with the reacceleration strength,
2
but are negatively correlated with the otherpropagation parameters.
.
8 Conclusions. The MCMC is a practical and powerful tool for cosmic-ray physicanalyses. It can be used to confirm hy-
0 pothesesconcerningsourcespectra(e.g.,whetherαi 6=αj)and/ordeterminewhetherdifferentdatasetsarecompatible.
8 A forthcoming study will extendour analysis to more physical diffusion models.
0
: Key words.Methods: statistical – ISM: cosmic-rays
v
i
X
r 1. Introduction sentedtheB/Cratioat0.5 50GeV/n(Panovetal.2007),
a and for H to Fe fluxes at 1−00 GeV 100 TeV (Panov et al.
One issue of cosmic-ray (CR) physics is the determination −
2006). At higher energy, two long-duration balloon flights
of the transport parameters in the Galaxy. This determi-
willsoonprovidespectrafor Z=1-30nuclei.The TRACER
nationisbasedontheanalysisofthesecondary-to-primary
collaboration has published spectra for oxygen up to iron
ratio (e.g., B/C, sub-Fe/Fe), for which the dependence on in the GeV/n-TeV/n range (Boyle et al. 2007; Ave et al.
the source spectra is negligible, and the ratio remains in-
2008). A second long-duration flight took place in summer
stead mainly sensitive to the propagation processes (e.g.,
2006, during which the instrument was designed to have a
Maurin et al. 2001 and references therein). For almost 20 wider dynamic-range capability and to measure lighter B,
years, the determination of these parameters relied mostly
C, and N elements. The CREAM experiment (Seo et al.
on the most constraining data, namely the HEAO-3 data,
2004) flew a cumulative duration of 70 days in December
taken in 1979, which covered the 1 35 GeV/n range 2004 and December 2005 (Seo et al. 2006, and preliminary
∼ −
(Engelmann et al. 1990).
results in Marrocchesi et al. 2006 and Wakely et al. 2006),
For the first time since HEAO-3, several satellite or and again in December 2007. A fourth flight was sched-
balloon-borne experiments (see ICRC 2007 reporter’s talk uledforDecember20081.Excitingdatawillarrivefromthe
Blasi 2008) have acquired higher quality data in the same PAMELAsatellite(Picozzaetal.2007),whichwassuccess-
energyrangeorcoveredascarcelyexploredrange(interms fully launched in June 2006 (Casolino et al. 2008).
of energy, 1 TeV/n PeV/n, or in terms of nucleus): from With this wealth of new data, it is relevant to ques-
−
the balloon-borne side, the ATIC collaboration has pre- tion the method used to extract the propagation parame-
Send offprint requests to: Antje Putze, [email protected] 1 http://cosmicray.umd.edu/cream/cream.html
2 A.Putze et al.: A Markov Chain Monte Carlo for Galactic cosmic-ray physics
ters. The value of these parameters is important to many wherethemodelfitsthedata)onlyfillatinyvolume.Thisis
theoretical and astrophysical questions, because they are where the MCMC approach, based on the Bayesian statis-
linked, amongst others, to the transport in turbulent mag- tics,issuperiortoagridapproach.Asinthegridapproach,
netic fields, sources of CRs, and γ-raydiffuse emission (see oneend-productoftheanalysisistheχ2surface,butwitha
Strong et al. 2007 for a recent review and references). It more efficient sampling of the region of interest. Moreover,
also proves to be crucialfor indirect dark-matter detection asopposedtoclassicalstatistics,whichisbasedonthecon-
studies(e.g.,Donatoetal.2004,andDelahayeetal.2008). struction of estimators of the parameters, Bayesian statis-
The usage in the past has been based mostly on a manual tics assumes the unknown parameters to be random vari-
orsemi-automated—hencepartial—coverageoftheparam- ables. As such, their full distribution—the so-called condi-
eter space (e.g., Webber et al. 1992, Strong & Moskalenko tional probability-density function (PDF)—given some ex-
1998,andJonesetal.2001).Morecompletescanswereper- perimental data (and some prior density for these parame-
formed in Maurin et al. (2001, 2002), and Lionetto et al. ters, see below) can be generated.
(2005), although in an inefficient manner: the addition of Tosummarise,the MCMCalgorithmprovidesthePDF
a single new free parameter (as completed for example in of the model parameters, based on selected experimental
Maurin et al. 2002 compared to Maurin et al. 2001) re- data (e.g., B/C). The mean value and uncertainty in these
mains prohibitive in terms of computing time. To remedy parameters are by-products of the PDF. The MCMC en-
these shortcomings, we propose to use the Markov Chain ables the enlargementof the parameterspace ata minimal
Monte Carlo (MCMC) algorithm, which is widely used in computingtimecost(althoughtheMCMCandMetropolis-
cosmological parameter estimates (e.g., Christensen et al. Hastings algorithms used here are not the most efficient
2001, Lewis & Bridle 2002, and Dunkley et al. 2005). One one).Thetechnicalities ofthe MCMC arebriefly described
goal of the paper is to confirm whether the MCMC algo- below. The reader is referred to Neal (1993) and MacKay
rithm can provide similar benefits in CR physics. (2003) for a more substantial coverage of the subject.
The analysis is performed in the framework of the
Leaky-Box Model (LBM), a simple and widely used prop-
agation model. This model contains most of the CR phe- 3. Markov Chain Monte Carlo (MCMC)
nomenology and is well adapted to a first implementation
Considering a model depending on m parameters
of the MCMC tool. In Sect. 2, we highlight the appropri-
ateness of the MCMC compared to other algorithms used θ θ(1), θ(2), ..., θ(m) , (1)
in the field. In Sect. 3, the MCMC algorithm is presented. ≡{ }
In Sect. 4, this algorithm is implemented in the LBM. In
weaimtodeterminetheconditionalPDFoftheparameters
Sect.5,wediscusstheMCMCadvantagesandeffectiveness given the data, P(θ data). This so-called posterior proba-
inthefieldofCRphysics,andpresentresultsfortheLBM. |
bility quantifies the change in the degree of belief one can
We present our conclusions in Sect. 6. Application of the
have in the m parameters of the model in the light of the
MCMC technique to a more up-to-date modelling, such as
data.Appliedtotheparameterinference,Bayestheoremis
diffusion models, is left to a forthcoming paper.
P(dataθ) P(θ)
P(θ data)= | · , (2)
2. Link between the MCMC, the CR data, and the | P(data)
model parameters
where P(data) is the data probability (the latter does not
Various models describe the propagation of CRs in the in- depend on the parametersand hence, can be consideredto
terstellarmedium(Webberetal.1992;Bloemenetal.1993; beanormalisationfactor).Thistheoremlinkstheposterior
Strong & Moskalenko 1998; Maurin et al. 2001; Berezhko probability to the likelihood of the data (θ) P(dataθ)
et al. 2003; Shibata et al. 2006; Evoli et al. 2008). Each andtheso-calledpriorprobability,P(θ),Lindica≡tingthed|e-
model is based on his own specific geometry and has its gree of belief one has before observing the data. To extract
own set of parameters, characterising the Galaxy proper- information about a single parameter, θ(α), the posterior
ties. The MCMC approach aims to study quantitatively density is integrated over all other parameters θ(k6=α) in
how the existing (or future) CR measurements can con- a procedure called marginalisation. Finally, by integrating
strain these models or, equivalently, how in such models the individual posterior PDF further, we are able to deter-
the set of parameters (and their uncertainties) can be in- mine the expectation value, confidence level, or higher or-
ferred from the data. dermodeoftheparameterθ(α).Thisillustratesthe techni-
In practice, a given set of parameters in a propagation caldifficulty ofBayesianparameterestimates: determining
model implies, e.g., a given B/C ratio. The model param- the individual posterior PDF requires a high-dimensional
eters are constrained such as to reproduce the measured integration of the overall posterior density. Thus, an effi-
ratio. The standard practice used to be an eye inspection cient sampling method for the posterior PDF is manda-
ofthe goodnessoffit to the data.This wasreplacedby the tory. For models of more than a few parameters, regular
χ2 analysis in recent papers: assuming the χ2 statistics is grid-samplingapproachesarenotapplicableandstatistical
applicable to the problem at stake, confidence intervals in techniques are required (Cowan 1997).
these parameters can be extracted (see App. A). Among these techniques, MCMC algorithms have been
The main drawback of this approach is the computing fully tried and tested for Bayesian parameter inference
time required to extend the calculation of the χ2 surface (MacKay 2003; Neal 1993). MCMC methods explore any
to a wider parameter space. This is known as the curse of targetdistributiongivenbyavectorofparametersp(θ),by
dimensionality, due to the exponential increase in volume generating a sequence of n points (hereafter a chain)
associated with adding extra dimensions to the parame-
ter space, while the good regions of this space (for instance θ θ , θ , ..., θ . (3)
i i=1,...,n 1 2 n
{ } ≡{ }
A.Putze et al.: A Markov Chain Monte Carlo for Galactic cosmic-ray physics 3
Eachθ isavectorofmcomponents[asdefinedinEq.(1)]. Burn-inlength Theburn-indescribesthepracticeofremov-
i
In addition, the chain is Markovian in the sense that the ing some iterations at the beginning of the chain to elim-
distribution of θ is influenced entirely by the value of inate the transient time needed to reach the equilibrium
n+1
θ .MCMCalgorithmsaredevelopedsothatthetimespent or stationary distribution, i.e., to forget the starting point.
n
by the Markov chain in a region of the parameter space is The burn-in length b is defined to be the number of first
proportionaltothetargetPDFvalueinthis region.Hence, samples θ of the chain that must be discarded.
i i=1,...,b
{ }
fromsuchachain,onecanobtainanindependentsampling The stationary distribution is reached when the chain en-
of the PDF. The target PDF as well as all marginalised ters the most probable parameter region corresponding to
PDF are estimated by counting the number of samples the regionwherethetargetfunctioniscloseto itsmaximal
within the related region of parameter space. value. To estimate b, the following criterionis used: we de-
Below, we provide a brief introduction to an MCMC finep tobethemedianofthetargetfunctiondistribution
1/2
usingtheMetropolis-Hastingsalgorithm(seeNeal1993and obtained from the entire chain of N samples. The burn-
MacKay2003,chapter29forfurtherdetailsandreferences). in length b corresponds to the first sample θ , for which
b
p(θ )>p (see App. C for an illustration).
b 1/2
3.1. The algorithm
Correlation length By construction [see Eq. (5)], each step
The prescription that we use to generate the Markov
of the chain depends on the previous one, which ensures
chains from the unknown target distribution is the so-
that the steps of the chain are correlated. We can obtain
called Metropolis-Hastings algorithm. The Markov chain
independent samples by thinning the chain, i.e. by select-
increases by jumping from the current point in the pa-
rameter space θ to the following θ . As said before, the ing only a fraction of the steps with a periodicity chosen
i i+1 toderiveuncorrelatedsamples.Thisperiodisestimatedby
PDF of the new point only depends on the current point,
i.e. (θ θ ,...,θ )= (θ θ ). This quantity defines computing the autocorrelation functions for each parame-
theTtranis+it1i|on1 probabiilityTforis+t1a|teiθ from the state θ . ter. For a parameter θ(α) (α = 1,...,m), the autocorrela-
i+1 i tion function is given by
The Metropolis-Hastings algorithm specifies to ensure
T
thatthestationarydistributionofthechainasymptotically 2
tends to the target PDF one wishes to sample from. E θ(α)θ(α) E θ(α)
Ateachstepi(correspondingtoastateθi),atrialstate c(jα) = h i j+ii−(cid:16) 2h i i(cid:17) , (7)
θ is generatedfrom a proposal density q(θ θ ). This E θ(α)
trial trial| i i
proposaldensityischosensothatsamplescanbeeasilygen- (cid:20)(cid:16) (cid:17) (cid:21)
erated(e.g.,aGaussiandistributioncentredonthecurrent
which we calculate with the Fast Fourier Transformation
state).Thestateθtrial isacceptedorrejecteddependingon (FFT). The correlation length l(α) for the α-th parameter
the following criterion. By forming the quantity
is defined as the smallest j for which c(α) < 1/2, i.e. the
p(θ )q(θ θ ) j
a(θ θ )=min 1, trial i| trial , (4) values θ(α) and θ(α) of the chain that are considered to be
trial| i p(θ ) q(θ θ ) i i+j
(cid:18) i trial| i (cid:19) uncorrelated. The correlationlength l for the chain, for all
thetrialstateisacceptedasanewstatewithaprobabilitya parameters, is defined to be
(rejectedwithprobability1 a).Thetransitionprobability
is then − l max l(α), (8)
(θ θ )=a(θ θ )q(θ θ ). (5) ≡α=1,...,m
i+1 i trial i trial i
T | | |
Ifaccepted,θ =θ ,whereasifrejected,thenewstate which is used as the periodof the thinning (see App. C for
i+1 trial
is equivalent to the current state, θ =θ . This criterion an illustration).
i+1 i
ensures that once at its equilibrium, the chain samples the
targetdistribution p(θ). If the proposaldensity q(θ θ )
trial| i Independentsamplesandacceptance Theindependentsam-
is chosen to be symmetric, it cancels out in the expression ples ofthe chainarechosento be θ ,where k is an
i i=b+lk
of the acceptance probability, which becomes: { }
integer.ThenumberofindependentsamplesN isdefined
ind
p(θ ) to be the fraction of steps remaining after discarding the
trial
a=min 1, p(θ ) . (6) burn-in steps and thinning the chain,
(cid:18) i (cid:19)
We note that the process requires only evaluations of Ntot b
N = − . (9)
ratios of the target PDF. This is a major virtue of this al- ind l
gorithm, in particular for Bayesian applications, in which
thenormalisationfactorinEq.(2),P(data)= P(dataθ) Theindependentacceptancefind istheratioofthenumber
P(θ)dθ is often extremely difficult to compute. Hence,|the· ofindependentsamplesNind tothetotalstepnumberNtot,
R
ratio of the target PDF, i.e. the posterior of the param-
N
ind
eter for our problem, can be calculated directly from the find = . (10)
N
likelihood of the data and the priors. tot
3.3. Choice of the target and trial functions
3.2. Chain analysis
3.3.1. Target function
The chain analysis refers to the study of several proper-
ties ofthe chains.The followingquantities areinspected in As already said, we wish to sample the target func-
order to convert the chains in PDFs. tion p(θ) = P(θ data). Using Eq. (2) and the fact that
|
4 A.Putze et al.: A Markov Chain Monte Carlo for Galactic cosmic-ray physics
the algorithm is insensitive to the normalisation factor, probability and are rejected, leading to a low acceptance
this amounts to sampling the product P(dataθ) P(θ). and a long correlation length. Conversely, for too small a
Assuming a flat prior P(θ) = cst, the target d|istri·bution width, the chain will take a longer time to reach the inter-
reduces to esting regions. Eventually, even if the chain reaches these
p(θ)=P(dataθ) (θ), (11) regions of high acceptance, only a partial coverage of the
| ≡L
PDF supportwillbe sampled(also leading to a longcorre-
and here, the likelihood function is taken to be
lation length).
χ2(θ) In practice, we first define σα (α = 1,...,m) equal to
(θ)=exp . (12) the expected range of the parameter. In a subsequent it-
L − 2
(cid:18) (cid:19) eration, σ is set to be 2√2ln2 2.3 times σcalc, i.e. the
α ≈ α
The χ2(θ) function for n data is FWHM of the PDF obtained with the first iteration. The
data
result is actually insensitive to the numerical factor used.
ndata (yexp ytheo(θ))2
χ2(θ)= k − k , (13)
σ2 Covariance matrix The proposal density is taken to be an
k=1 k
X N-dimensional Gaussian of covariance matrix V
where yexp is the measuredvalue, ytheo is the hypothesised
k k 1
value for both a certain model and the parameters θ, and q(θ ,θ ) exp (θ θ )TV−1(θ θ ) . (15)
trial i trial i trial i
σk istheknownvarianceofthemeasurement.Forexample, ∝ (cid:18)−2 − − (cid:19)
yexp and ytheo represent the measured and calculated B/C
k k The covariance matrix V is symmetric and diagonalisable
ratios.
(D isadiagonalmatrixofeigenvaluesandP representsthe
The link betweenthe targetfunction, i.e., the posterior
change in the coordinate matrix),
PDFoftheparameters,andtheexperimentaldataisestab-
lishedwiththehelpofEqs(11)to(13).Thislinkguarantees V =PTDP,
the proper sampling of the parameter space using Markov
and where again Eq. (6) holds. The parameters θ are
chains, which spend more time in more relevant regions of trial
hence found to be
parameter space, as described above.
θ =θ +PTDx,
trial i
3.3.2. Trial function where x is a vector of m random numbers following a
Despite the effectiveness of the Metropolis-Hastings algo- Gaussian distribution centred on zero and with unit vari-
rithm, to optimise the efficiency of the MCMC and min- ance.
imise the number of chains to be processed, trial functions The covariance matrix V is estimated, e.g., from a pre-
should be as close as possible to the true distributions. We vious iteration using the Gaussian step. The advantage of
useasequenceofthreetrialfunctionstoexploretheparam- this trial function with respect to the previous one is that
eter space. The first step is a coarse determination of the it takes accountof the possible correlationsbetween the m
parameter PDF. This allows us to calculate the covariance parameters of the model.
matrixleadingtoabettercoverageofparameterspace,pro-
vided that the target PDF is sufficiently close to being an
BinarySpacePartitioning(BSP) Athirdmethodwasdevel-
N-dimensional Gaussian. The last step takes advantage of
oped to define a proposal density for which the results of
a binary-space partitioning (BSP) algorithm.
the Gaussian step or the covariance matrix iterations are
used to subdivide the parameter space into boxes, in each
Gaussian step For the first iteration, the proposal density of which a given probability is affected.
q(θ ,θ ),requiredtoobtainthetrialvalueθ fromθ The partitioning of the parameter space can be organ-
trial i trial i
is written as ised using a binary-tree data structure known as a binary-
spacepartitioningtree(deBergetal.2000).Therootnode
2
θ(α) θ(α) of the tree is the m dimensionalbox corresponding to the
q(θ ,θ ) exp 1 trial− i . (14) entire parameter sp−ace. The binary-space partitioning is
trial i ∝α=1,...,m −2(cid:16) σα2 (cid:17)  then performed by dividing each box recursively into two
Y   childboxesifthepartitioningsatisfiesthefollowingrequire-
 
ment: a box is divided only if the number of independent
TheserepresentmindependentGaussiandistributionscen-
tred on θ . The distribution is symmetric, so that the ac- samplescontainedinthisboxishigherthanacertainnum-
i
ceptance probability a follows Eq. (6). The varianceσ2 for ber (here we used a maximum of between 3% and 0.1%
α
(α) of the total number of independent samples). When a box
each parameter α is to be specified. Each parameter θ
trial has to be divided, the division is made along the longer
is hence calculated to be
side of the box (the box-side lengths aredefined relative to
θ(α) =θ(α)+σ x, the root-box sides). For each end node (i.e. node without
trial i α· any children), a probability, defined as the fraction of the
where x is a random number obeying a Gaussian distribu- number of independent samples in the box to their total
tion centred on zero with unit variance. number, is assigned. For empty boxes, a minimum proba-
Itisimportanttochooseanoptimalwidthσ tosample bility is assigned and all the probabilities are renormalised
α
properly the posterior (target) distribution. If the width is so that the sum of all end-node probabilities equals 1.
toolarge,assoonasthechainreachesaregionofhighprob- The proposal density q(θ ) is then defined, in each
trial
ability,mostofthetrialparametersfallintoaregionoflow end-node box, as a uniform function equal to the assigned
A.Putze et al.: A Markov Chain Monte Carlo for Galactic cosmic-ray physics 5
marise, the initial values of the propagation parameters
θ are chosen randomly in their expected range to crank
0
up each Markov chain. The interstellar (IS) CR fluxes are
then calculated for this set of parameters (see, e.g., Fig. 1
in Maurin et al. 2002for further details of the propagation
steps).TheISfluxismodulatedwiththeforce-fieldapprox-
imation and the resulting top-of-atmosphere (TOA) spec-
trum is compared with the data, which allows us to calcu-
late the χ2 value [Eq.(13)],hence the likelihood[Eq.(12)].
This likelihood (in practice the log-likelihood) is used to
compute the acceptance probability [Eq. (4)] of the trial
vectorofparametersθ (asgeneratedbyoneofthethree
trial
trial functions described in Sect. 3.3.2). Whether the trial
vector is accepted or rejected implies whether θ = θ
1 trial
or θ = θ . This procedure is repeated for the N steps
1 0
of the chain. Obviously, when θ = θ , the propagation
i+1 i
step does not need to be repeated. Because of the nature
of the MCMC algorithm, several chains can be executed
in parallel. Once completed, these chains are analysed (see
Sect. 3.2)—discarding the first step belonging to the burn-
ing length and thinned according to the correlation length
l [Eq. (8)]—and combined to recover the desired posterior
PDF P(θ data).
|
In this procedure, the user must decide i) the data to
be used, ii) the observable to retain in calculating the like-
Faivge.c1t.orFl[oEwq.c(h1a)r]toofftthhee iαmp=le1m,e.n.t.e,dmMfrCeeMpCaraalgmoertitehrsmo:fθtihies lviehcotoodr,θa)nfdoriiwi)hitchhewneumsebeekrthofefproeestepraiorarmPeDteFr.s m (of the
model,evaluatedateachstepi,andp(θi)isthetargetfunction
given by Eq. (11). See text for details.
4.2. Leaky-Box Model (LBM)
The LBM assumes that all CR species are confined within
probability. The sampling of this proposal density is sim-
the Galaxy with an escape rate that equals N/τ , where
esc
ple and efficient: an end node is chosen with the assigned
the escape time τ is rigidity-dependent, and is written
esc
probability and the trial parameters are chosen uniformly
as τ (R). This escape time has two origins. First, CRs
esc
in the corresponding box. In comparison to the other two
canleakoutthe confinementvolumeandleavethe Galaxy.
proposaldensities, this proposaldensity basedon a BSP is
Second, they can be destructed by spallation on interstel-
asymmetric, because it is only dependent on the proposal
lar matter nuclei.This latter effect is parameterisedby the
state q(θtrial). Hence, Eq. (4) must be used. grammage x (usually expressed in g cm−2), defined as the
columndensityofinterstellarmatterencounteredbyapath
followedby aCR. The CRs thatreachEarthhavefollowed
4. Implementation in the propagation model
different paths, and can therefore be described by a gram-
TheMCMCwiththethreeabovemethodsareimplemented mage distribution N(x) dN/dx. The LBM assumes that
≡
intheUSINEpackage2,whichcomputesthepropagationof
Galactic CR nuclei and anti-nuclei for several propagation N(x) exp−λesc(R)x , (16)
∝
models (LBM, 1D and 2D diffusion models). The reader is
referred to Maurin et al. (2001) for a detailed description where the mean grammage λesc(R) = x is related to the
h i
for the nuclear parameters (fragmentation and absorption mass m, velocity v and escape time τesc(R) by means of
cross-sections),energylosses(ionisationandCoulomb),and λesc(R)=m¯nvτesc(R).
solar modulation (force-field) used. The function λesc(R) determines the amount of spalla-
We briefly describe how the MCMC algorithm is im- tionsexperiencedbyaprimaryspecies,andthusdetermines
plemented in the propagation part (Sect. 4.1), using a thesecondary-to-primaryratios,forinstanceB/C.Froman
LBM—the procedure would be similar for any other experimentalistpointofview,λesc(R)isaquantitythatcan
model. The LBM and its parameters are briefly discussed be inferred from measurements of nuclei abundance ratios.
(Sect. 4.2) as well as the input spectrum parameters Thegrammageλesc(R)isknowntoprovideaneffectivede-
(Sect.4.3).Additionalinformationaboutthedataaregath- scriptionofdiffusionmodels(Berezinskiietal.1990):itcan
ered in App. B. be related to the efficiency of confinement (which is deter-
mined by the diffusion coefficient and to both the size and
geometry of the diffusion volume), spallative destruction
4.1. Flow chart (which tends to shorten the average lifetime of a CR and
thus lower λ ), and a mixture of other processes (such as
A flow chart of the Metropolis-Hastings MCMC algorithm esc
convection, energy gain, and losses).
used in the context of GCRs is given in Fig. 1. To sum-
In this paper, we compute the fluxes in the framework
2 A public version will be released soon (Maurin, in prepara- of the LBMwith minimal reaccelerationby the interstellar
tion). turbulence, as described in Osborne & Ptuskin (1988) and
6 A.Putze et al.: A Markov Chain Monte Carlo for Galactic cosmic-ray physics
Seo & Ptuskin (1994). The grammage λ (R) is parame- The pattern of the source abundances observed in the
esc
terised as cosmic radiationdiffers fromthat ofthe solarsystem.This
is due to a segregation mechanism during the acceleration
λ (R)= λ0βR0−(δ−δ0)R−δ0 when R<R0, (17) stage.Twohypothesesaredisputedintheliterature:oneis
esc (λ0βR−δ otherwise; based on the CR composition controlled by volatility and
mass-to-chargeratio(Meyeretal.1997;Ellisonetal.1997),
where we allow for a break,i.e. a different slope below and and the other one is based on the first ionisation potential
above a critical rigidityR . The standardform used in the (FIP)ofnuclei(e.g.,Cass´e&Goret1973).Inthiswork,for
0
literature is recovered by setting δ = 0. For the entire set each configuration, the source abundances are initialised
0
of n nuclei, a series of n equations (see Maurin et al. 2001 to the product of the solar system abundances (Lodders
for more details) for the differential densities Nj=1,...,n are 2003), and the value of the FIP taken from Binns et al.
solvedata givenkinetic energyper nucleonE (E is the (1989). The final fluxes are obtained by an iterative calcu-
k/n
total energy), i.e. lationofthepropagatedfluxes,rescalingtheelementabun-
dances—keepingfixedtherelativeisotopicabundances—to
d dNj match experimental data at each step until convergence is
AjNj(E )+ BjNj Cj =Sj(E ). (18)
k/n dE − dE k/n reached(seeFig.1inMaurinetal.2002forfurtherdetails).
(cid:18) (cid:19)
The result is thus insensitive to the input values (more de-
Inthisequation,ther.h.s.termisthesourcetermthattakes tails about the procedure are given in App. B.1).
into account the primary contribution (see Sect. 4.3), the The measurement of all propagated isotopic fluxes
spallative secondary contribution from all nuclei k heavier should characterise all source spectra parameters com-
than j, and the β-decay of radioactive nuclei into j. The pletely, i.e. the q and α parameters should be free.
j j
first energy-dependent factor Aj is given by However,onlyelementfluxesareavailable,whichmotivates
the above rescaling approach. In Sect. 5.3, a few calcula-
1 1
Aj = + n vjσj+ISM+ . tions are undertaken to determine self-consistently, along
τesc ISM=H,He ISM inel τβj with the propagation parameters, i) α and η, and ii) the
X
source abundances for the primary species C, O, and the
The two other terms correspondto energy losses and first-
mixedNelements(themaincontributorstotheboronflux).
order reaccelerationfor Bj and to second-orderreaccelera-
tion for Cj. Following Osborne & Ptuskin (1988) and Seo
& Ptuskin (1994),
5. Results
dE
B = +(1+β2)β2EK and C =β4E2K , We first examine the relative merits of four different pa-
dt ion,coul. pp pp
rameterisations of the LBM, and determine the statistical
where(cid:10) (cid:11) significance of adding more parameters. These models cor-
Kpp = 43Va2δ(4 δτ2es)c(4 δ). (19) respond to {θα}α=1,...,m≤5 with
− − – Model I = λ , R ,δ , i.e. no reacceleration ( = 0)
The strength of the reacceleration is mediated by the { 0 0 } Va
and no break in the spectral index (δ =0).
pkmseusd−o1 Akplfcv−´e1n.icThspiseeids rVelaatoedf tthoeastcrautetesrperesedingivuennitsinoaf – ModelII ={λ0, δ, Va},i.e. no critical0rigidity (R0 =0)
and no break in the spectral index (δ =0).
diffusion model with a thin disk h and a diffusive halo L – ModelIII= λ , R , δ, ,i.e.nobr0eakinthespectral
bAyssmumeainngs toyfpVicaal=vaVluae×s o(hfLh)−=1/02.1(Skepoc&anPdtLusk=in1019k9p4c),. – Mindoedxel(IδV0 ==0{)λ.0, R0, δ V, aδ,} .
the value of a can be directly transposed and compared { 0 0 0 Va}
V
to a true speed V , as obtained in diffusion models.
a VarioussubsetsofB/Cdataareusedtoinvestigatewhether
To summarise, our LBM with reacceleration may in-
old data are useful or just add confusion to the PDF de-
volve up to five free parameters, i.e. the normalisation λ ,
0 termination. We note in Sect. 5.2 that no useful constraint
the slopes δ and δ below or above the cut-off rigidity R ,
0 0 can be drawn from p data alone.
and a pseudo-Alfv´en velocity related to the reaccelera-
Va We also consider additional free parameters (Sect. 5.3)
tion strength.
related to the source spectra, for a self-consistent determi-
nation of the propagation and source properties. Since we
4.3. Source spectra show that a break in the slope (Model IV) is not required
by currentdata, we focus onModel III (for the description
We assume that the primary source spectrum Q (E) for
j of the propagationparameters), defining:
each nuclear species j is given by (β =v/c)
Qj(E)≡dQj/dE =qjβηjR−αj, (20) – sMloopdeelαIIisI+a1fr=ee{pλa0r,aRm0e,tδe,r.Va}+{α},where the source
where q is the source abundance, α is the slope of the – ModelIII+2= λ , R , δ, + α,η ,whereboththe
j j 0 0 a
species j, and the term βηj manifests our ignorance about source slope α a{nd the expVon}ent{η [o}f β, see Eq. (20)]
the low-energy spectral shape. We further assume that are free parameters.
α α for all j, and unless stated otherwise, η η = 1 – ModelIII+4= λ , R , δ, + α, q , q , q ,where
j j 0 0 a C N O
in o≡rder to recover dQ/dp p−α, as obtained fr≡om acc−el- theabundances{q ofthemoVst}sig{nificantlycont}ributing
i
∝
eration models (e.g., Jones 1994). The constraints existing elements are also free parameters.
on η are explored in Sect. 5.3. – Model III+5 = λ , R , δ, + α, η, q , q , q .
0 0 a C N O
{ V } { }
A.Putze et al.: A Markov Chain Monte Carlo for Galactic cosmic-ray physics 7
This allows us to investigate the correlations between pa- 0.2 l 0 [g/cm2] Binary space partitioning step
rameters further and into potential biases in the propaga-
tion parameter determination. Chain analysis
More details about the practical use of the trial func- 0.1 Burn-In length: 1.9
tions can be found in App. C. In particular, the sequential Correlation Length: 1.5
use of the three sampling methods (Gaussian step, covari- find: 149975/200000 = 74.987%
c2 /dof: 1.43
ance matrix step, and then binary-space partitioning) is 0 min
20 30 40
found to be the most efficient: all results presented here- d
after are based on this sequence. 0.6 20
0.55
5.1. Fitting the B/C ratio 10
0.5
5.1.1. HEAO-3 data alone
0
20 30 40 0.5 0.55 0.6
WefirstconstrainthemodelparameterswithHEAO-3data 0.06 V [km/s/kpc]
a
only(Engelmannetal.1990).Thesedataarethemostpre- 100 100
cise data available at the present day for the stable nuclei 0.04
ratio B/C of energy between 0.62 to 35GeV/n. 80 80
The results for the models I, II, and III are presented 0.02
60 60
in Figs. C.2, 2 top, and 2 bottom. The inner and outer
contoursaretakentoberegionscontaining68%and95%of 20 30 40 0.5 0.55 0.6 0 60 80 100
the PDF respectively (see App. A.1). The first observation
one can make for the LBM without reacceleration (Model 0.1 l 0 [g/cm2] Binary space partitioning step
I, Fig. C.2), is that the marginal distributions of the three Chain analysis
LBM parameters are mostly Gaussian. The tail for small 0.05 Burn-In length: 4.0
Correlation Length: 6.9
lvoawlu-eesneorfgRy0daistadu(e<t1oGtheiVs/pna)r:amtheetreerabreeinngocHonEsAtrOai-n3eddabtya 020 30 40 cfm2iinnd/:d o3f3:7 511./32000000 = 16.875%
atlowenergy,soallR0 valuesbelow3GVareequiprobable 4 R0 [GV]
(this remains true for Model III).
0.5
As seeninFig.2,a morecomplicatedshape for the dif-
2
ferent parameters is found for Model II (top panel), and
even more so for Model III (bottom panel). This induces
a longer correlation length (1.5 and 6.9 steps instead of 1 020 30 40 00 2 4
d
step) and hence reduces the efficiency of the MCMC (75% 10
formodelIIand17%formodelIII).Physically,thecorrela- 0.6 0.6
tionbetween the parameters,as seenmostclearlyin Fig.2 5
(bottom), is understoodasfollows.First,λ ,R ,andδ are 0.5 0.5
0 0
positively correlated. This originates in the low-energy re- 20 30 40 0 2 4 0 0.5 0.6
lation λesc ∝ λ0R0−δ, which should remain approximately 100 100 100 0.04 Va [km/s/kpc]
constant to reproduce the bulk of the data at GeV/n en-
80 80 80
ergy. Hence, if R or δ is increased, λ also increases to 0.02
0 0
60 60 60
balance the product. On the other hand, is negatively
a
correlatedwith δ (and hence with all the pVarameters):this 40 40 40 0
20 30 40 0 2 4 0.5 0.6 40 60 80 100
is the standard result that to reach smaller δ (for instance
Fig.2.Posterior distributionsforModelII(top)andModelIII
to reach a Kolmogorov spectrum), more reacceleration is
(bottom) using HEAO-3 data only. For more details, refer to
required.ThiscanalsobeseenfromEq.(19),whereatcon-
caption of Fig. C.2.
stantτ ,K 2/f(δ),wheref isadecreasingfunction
esc pp ∝Va
of δ: hence, if δ decreases, f(δ) increases, and then has
a
V
to increase to retain the balance. The reaccelerationmechanism was invoked in the liter-
The values for the maximumof the PDF for the propa- ature to decrease the spectral index δ toward its preferred
gation parameters along with their 68% confidence inter- valueof1/3givenbyaKolmogorovspectrumofturbulence.
vals (see App. A) are listed in Table 1. The values ob- InTable1,theestimatedpropagationparametervaluesfor
tained for our Model I are in fair agreement with those the models II and III are indeed slightly smaller than for
derived by Webber et al. (1998), who found λ , R , δ = ModelI,buttheKolmogorovspectralindexisexcludedfor
0 0
38.27,3.6, 0.7 . The difference for λ cou{ld be rela}ted all of these three cases (using HEAO-3 data only). This
0
{to the fact th}at Webber et al. (1998) rely on a mere result agrees with the findings of Maurin et al. (2001), in
eye inspection to extract the best-fit solution or/and use whichamorerealistictwo-dimensionaldiffusionmodelwith
a different set of data. For example, comparing Model I reacceleration and convection was used. We note that the
with a combination of HEAO-3 and low-energy data values for a 80 km s−1 kpc−1, should lead to a true
V ∼
(ACE+Voyager1&2+IMP7-8, see Sect. 5.1.2) leads to speedV = √hL 80kms−1 in adiffusion modelfor
a a
V × ∼
λ , R , δ = 52, 5.3, 0.69 , slightly changing the values which the thin disk half-height is h=0.1 kpc and the halo
0 0
{ } { }
of the Model I preferred parameters (compared to the first size is L = 10 kpc: this is consistent with values found in
line of Table 1). Maurin et al. (2002).
8 A.Putze et al.: A Markov Chain Monte Carlo for Galactic cosmic-ray physics
Model λ0 R0 δ Va χ2min/dof Model λ0 R0 δ Va χ2min/dof
g cm−2 GV km s−1kpc−1 Dataset g cm−2 GV km s−1kpc−1
I 54+2 4.2+0.3 0.70+0.01 - 3.35 III-A 30+5 2.8+0.6 0.58+0.01 75+10 1.30
−2 −0.9 −0.01 −4 −0.8 −0.06 −13
II 26+2 - 0.52+0.02 88+6 1.43 III-B 28+2 2.6+0.4 0.53+0.02 85+9 1.09
−2 −0.02 −11 −3 −0.7 −0.03 −8
III 30+5 2.8+0.6 0.58+0.01 75+10 1.30 III-C 27+2 2.6+0.4 0.53+0.02 86+9 1.06
−4 −0.8 −0.06 −13 −2 −0.7 −0.03 −5
Table 1. Most probable values of the propagation parameters
(aftermarginalising overtheotherparameters)formodelsI,II, III-D 26+−22 3.0+−00..45 0.52+−00..0022 95+−76 4.15
andIIIusingHEAO-3alone(14datapoints)andtheB/C con- III-E 30+2 3.7+0.2 0.57+0.01 88+3 6.08
−2 −0.3 −0.02 −6
straint. The uncertainty in the parameters correspond to 68%
CL of the marginalised PDF (see App. A). The last column Table 2. Same as in Table 1, but testing different data sets
shows the minimum χ2/dof obtained for each model (the asso- withModelIII:A=HEAO-3(14datapoints),B=HEAO-3+
ciated best-fit parameters are gathered in Table 3). ACE(20datapoints),C=HEAO-3+ACE+Voyager1&2+
IMP7-8 (22 data points), D = HEAO-3 + all low-energy data
(30 data points), E = all B/C data (69 data points).
o
ati 0.34 TOA
C r (F =250 MV) HEAO-3
B/ 0.32
e.g., the HEAO-3 data, and adding low-energydata allows
0.3 us to more reliably constrain R . We note that by fitting
0
only the low-energy data, only the grammage crossed in
0.28
very narrow energy domain would be constrained.
0.26 In a first step, we add only the ACE (CRIS) data (de
Nolfo et al. 2006), which covers the energy range from
0.24
8 10−2GeV/n to 2 10−1GeV/n, and which is later
0.22 Model I (c m2in/dof=3.35) r∼efer·redtoasdataset∼B(th·edatasetAbeingHEAO-3data
0.2 Model II (c 2 /dof=1.43) alone). The resulting posterior distributions are similar for
min the datasets B and A (B is not shown, but A is given in
0.18 Model III (c 2 /dof=1.30) Fig. 2, bottom). Results for datasets A and B are com-
min
pletely consistent(firstandsecondline ofTable 2),but for
0.16
the latter, propagation parameters are more tightly con-
10-1 1 Ek/n10 (GeV/n) strainedandthefitisimproved(χ2min/dof=1.09).TheACE
(CRIS)dataarecompatiblewithR =0,butthe preferred
0
Fig.3.Best-fitratio forModelI(bluedotted),II(reddashed), critical rigidity is 2.47GV.
andModelIII(blacksolid) usingtheHEAO-3dataonly(green
All other low-energy data (ISEE-3, Ulysses, IMP7-8,
symbols). The curves are modulated with Φ = 250GV. The
Voyager1&2, ACE) are then included (dataset D). The
corresponding best-fitparameters are gathered in Table 3.
resulting values of the propagation parameters are left un-
changed. However, a major difference lies in the higher
χ2 /dof of 4.15, which reflects an inconsistency between
The final column in Table 1 indicates, for each model, min
the best χ2 value per degreeof freedom,χ2 /dof. This al- the different low-energy data chosen for the MCMC. If
min the data point from the Ulysses experiment is excluded,
lows us to compare the relative merit of the models. LB
χ2 /dof decreases to a value of 2.26, and by excluding
models with reacceleration reproduce the HEAO-3 data min
more accurately (with χ2/dof of 1.43 and 1.30 for the also the ISEE-3 data points (dataset C) it decreases fur-
Models II and III respectively compared to χ2/dof = 4.35 ther to 1.06 (see Table 2). Since the set of low-energy
data have different modulation parameters, the difference
for Model I). The best-fit model B/C fluxes are shown
in the results for the various data subsets becomes clearer
with the B/C HEAO-3 data modulated at Φ = 250 MV
after the data have been demodulated. The force-field ap-
in Fig. 3. Physically, the origin of a cutoff R in λ at
0 esc
proximation provides a simple analytical one-to-one corre-
lowenergycanbe relatedto convectionindiffusionmodels
spondence between the modulated top-of-the atmosphere
(Jones 1979). Hence, it is a distinct process as reacceler-
(TOA)andthedemodulatedinterstellar(IS)fluxes.Foran
ation. The fact that Model III performs more successfully
isotopex,the ISandTOAenergiespernucleonarerelated
than Model II implies that both processes are significant,
byEIS =ETOA+Φ(Φ=Z/A φisthemodulationparam-
as found in Maurin et al. (2001). k k ×
eter), and the fluxes by (p is the momentum per nucleon
In the following, we no longer consider Model I and II, x
of x)
and inspect instead, the parameter dependence of Model
III on the dataset selected. pIS 2
ψIS EIS = x ψTOA ETOA+Z/A ψ . (21)
x k pTOA x k ×
5.1.2. Additional constraints from low-energy data (cid:18) x (cid:19)
(cid:0) (cid:1) (cid:0) (cid:1)
The B/C ratio results from a combination of various iso-
The actual data sets for the B/C ratio (see, e.g., Fig. 5)
topes, and assumingthe same Z/Afor allisotopes, we find
showaseparationintotwoenergydomains:the low-energy
that
range extends from 10−2GeV/n to 1GeV/n and the
high-energy range go∼es from 1GeV/∼n to 102GeV/n. B IS B TOA
∼ ∼ EIS = ETOA+Z/A φ . (22)
The spectral index δ is constrained by high-energy data, C k C k ×
(cid:18) (cid:19) (cid:18) (cid:19)
(cid:0) (cid:1) (cid:0) (cid:1)
A.Putze et al.: A Markov Chain Monte Carlo for Galactic cosmic-ray physics 9
o o
ati HEAO-3 (F =250 MV) ati 0.34 HEAO-3
B/C r0.34 VAoCyEa g(CerR1I&S2) ((FF ==222255 MMVV)) B/C r 0.32 ACE (CRIS)
0.32 IMP7-8 (F =250 MV) || Voyager1&2
ISEE-3 (F =370 MV) 0.3 F =250 MV IMP7-8
TOA
0.3 Ulysses (F =420 MV) IS 0.28 F =225 MV (HEAO)
(ACE+Voyager+IMP)
0.28 0.26
0.24
0.26
0.22
Model III-A (c 2 /dof=1.30)
0.24 min
0.2
Model III-C (c 2 /dof=1.06)
min
0.22 0.18
10-1 1 Ek/n (GeV/n) 10-1 1 10
Ek/n (GeV/n)
Fig.4.HEAO-3andACE(CRIS)modulated(TOA,solid line) Fig.5. Best-fit B/C flux (Model III) for datasets A (thin red
and demodulated (IS, dashed line) data points have been con- curves) and C (thick black curves). Above 300 MeV/n, B/C is
nected to guide the eye. Filled symbols (modulated and de- modulatedtoΦ=250MV(solidlines)appropriateforHEAO-3
modulated) correspond to HEAO-3, ACE (CRIS), IMP7-8 and data, whereas below, it is modulated to Φ = 225 MV (dashed
Voyager1&2. On the TOA curve, the empty red stars and the lines). Model III-B, not shown, overlaps with III-C. The corre-
blueuppertriangle correspond to ISEE-3 and Ulysses. sponding propagation valuesare gathered in Table 3.
Model λbest Rbest δbest δbest Vbest χ2/dof
0 0 0 a
Data g cm−2 GV kms−1kpc−1 Voyager1&2 constraints slightly shifts all of the parame-
ters to lower values.
I-A 54.7 4.21 - 0.702 - 3.35 In a final try, we take into account all available data
II-A 25.8 - - 0.514 88.8 1.43 (dataset E,final line of Table 2). Many data are clearly in-
III-A 31.7 2.73 - 0.564 73.0 1.30 consistent with each other (see Fig. 8), but as for the low-
energy case, although the χ2 /dof is worsened, the pre-
III-C 26.9 2.45 - 0.527 88.5 1.06 min
ferredvaluesofthepropagationparametersarenotchanged
IV-C 32.7 2.38 -0.97 0.572 70.5 0.86
drastically(comparewithdatasetsB,C,andDinTable2).
We await forthcoming data from CREAM, TRACER, and
Table 3. Best-fit values (corresponding to χ2 )
min PAMELA to be able to confirm and refine the results for
for B/C data (A=HEAO-3 data alone, C=HEAO-
HEAO-3 data.
3+Voyager1&2+ACE+IMP7-8).
5.1.3. Model IV: break in the spectral index
The modulatedanddemodulatedlow-energyB/Cdataare Wehavealreadymentionedthattherigiditycut-offmaybe
shown in Fig. 4 (see caption for details). The ISEE-3 and associatedwiththe existenceofagalacticwindindiffusion
Ulysses data points, as just underlined, are clearly incon- models. By allowing a break to occur in the spectral index
sistent with other data. To be consistent, Φ =200 MV for of λesc [see Eq. (17)], we search for a deviations from a
ISEE-3 and Φ = 200 MV for Ulysses would be required. single power law (δ0 =δ) or from the cut-off case (δ0 =0).
Significant uncertainties ∆Φ 25 50 GV are quoted in Adding a new parameter δ0 (Model IV) increases the
general, so that it is difficult∼to c−onclude whether there correlationlengthoftheMCMC,sinceR0 andδ0 arecorre-
are systematics in the measurement or if the modulation lated [see Eq. (17)]. The acceptance find [Eq. (9)] is hence
quoted in the papers is inappropriate. Some experiments extremely low. For Model IV-C (i.e. using dataset C, see
have also accumulated the signal for several years, peri- Table 2), we find find = 2%. The PDF for δ0 is shown
ods during which the modulationchanges.It is beyondthe in the left panel of Fig. 6. The most probable values and
scope of this paper to discuss this issue further. Below, we 68% confidence intervals obtained are λ0,R0,δ0,δ,Va =
{ }
discard both ISEE-3 and Ulysses data in selecting an ho- 30+2,2.2+0.4, 0.6+0.2,0.55+0.04,76+9 , which are con-
{ −2 −0.6 − −1.3 −0.02 −11}
mogeneous low-energy data set, which includes the most sistent with values found for other models, as given in
recent ACE (CRIS) data. Tables1and2:addingalow-energyspectralbreakonlyal-
The resulting best-fit models, when taking low-energy lowsustobetteradjustlow-energydata(figurenotshown).
dataintoaccount,aredisplayedinFig.5.TheB/Cbest-fit The best-fit parameters, for which χ2 = 0.86, are re-
min
ratioisdisplayedforModelIII andthe datasetA(redthin ported in Table 3. The small value of χ2 (smaller than
min
lines) and C (black thick lines). Model III-B (not shown) 1)mayindicateanover-adjustment,whichwoulddisfavour
yieldssimilarresultstoModelIII-C.Solidanddashedlines the model.
correspond to the two modulations Φ = 250 MV (HEAO- It is also interesting to compel δ to be positive,
0
3 and IMP7-8) and Φ = 225 MV respectively (ACE and to check whether δ = 0 (equivalent to Model III),
0
Voyager1&2). Although the fit from HEAO-3 alone pro- δ = δ (equivalent to Model II), or any value in-between
0
videsagoodmatchatlowenergy,addingACE(CRIS)and that is preferred. We find the most probable values to
10 A.Putze et al.: A Markov Chain Monte Carlo for Galactic cosmic-ray physics
0.6
d d > 0
Model IV 0 0 c 2/dof
2 3
0.4
2.5
1 2 Model III-C (c 2 /dof=1.06)
0.2 min
c 2 /dof + 0.273 (cid:222) 68% CL
1.5
min
c 2 /dof + 0.513 (cid:222) 95% CL
min
0 1
-4 -2 0 0 0.5 1
Fig.6. Marginalised PDF for the low-energy spectral index δ
0 0.5
in Model IV-C. The parameter δ is either free to span both
0
positiveandnegativevalues(leftpanel)orconstrainedtoδ >0
0
0
(right panel). 1 1.2 1.4 1.6 1.8 2
Fig.7.χ2/dofnormaliseddistributionforModelIII-C.The68%
be λ ,R ,δ ,δ,V = 23+1,1+2,0+0.6,0.49+0.01,102+4 . and 95% CL of the distribution are shown respectively as the
The{c0orre0spo0ndinga}PDF{fo−r1δ −is1shown in th−e0.r0i1ght p−an5}el red and black area.
0
of Fig. 6. The maximum occurs for δ = 0, which is also
0
found to be the best-fit value; we checked that the best-fit
o
pCcoa.rrrAaemsspeeotcenordnsidnmagrayttcopheeMaskotdhaeoplspeIeIag.rivsTehanteiδna0sTs≈oacb0ila.e5te,3dsfuoχcr2hMasofdoδer0l≈tIhIIiδs- B/C rati 00.3.35 TOA || F =250 MV BBBBaaaalllllllloooooooonnnn ####1234
min IMP7-8
configuration is worse than that obtained with δ0 = 0, in ISEE-3
agreement with the conclusion that Model III provides a 0.25 Spacelab-2
closer description of the data than Model II. F =225 MV HUElyAssOes-3
0.2 Voyager1&2
ACE (CRIS)
5.1.4. Summary and confidence levels for the B/C ratio 0.15 ATIC-2
Model III-C
Inthepreviousparagraphs,wehavestudiedseveralmodels
and B/C datasets. The two main conclusions that can be 0.1 Best fit
drawnarei)thebest-fitmodelisModelIII,whichincludes 68% CL
0.05
reaccelerationand a cut-off rigidity, and ii) the most likely 95% CL
valuesofthepropagationparametersarenottoodependent
0
on the data set used, although when data are inconsistent 10-1 1 10 102 103
with each other the statistical interpretation of the good- Ek/n (GeV/n)
ness of fit of a model is altered (all best-fit parameters are
Fig.8. Confidence regions of the B/C ratio for Model III-C as
gatheredinTable3).Thevaluesofthederivedpropagation calculatedfromallpropagationparameterssatisfyingEq.(A.2).
parameters are close to the values found in similar studies The blue-dashed line is the best-fit solution, red-solid line is
and the correlation between the LB transport parameters 68% CL and black-solid line 95% CL. Two modulation param-
are well understood. eters are used: Φ = 225 MV below 0.4 GeV/n (adapted for
Taking advantage of the knowledge of the χ2 distribu- ACE+Voyager1&2+IMP7-8 data) and Φ = 250 MV above
tion, we can extract a list of configurations, i.e. a list of (adapted for HEAO-3data).
parameter sets, based on CLs of the χ2 PDF (as explained
in App. A.2). The χ2 distribution is shown for our best
model(e.g.,diffusionmodel),forwhichthesituationmight
model, i.e. Model III, in Fig. 7. The red and black areas
not be so simple.
correspondtothe68%and95%confidenceintervals,which
From the same lists, we can also derive the range al-
are used to generate two configuration lists, from which
lowed for the source abundances of elements (we did not
68% and 95% CLs on, e.g., fluxes, can be derived3.
try to fit isotopic abundances here, although this can be
The B/C best-fit curve (dashed blue), the 68% (red
achieved, e.g., as in Simpson & Connell 2001 and refer-
solid), and 95% (black solid) CL envelopes are shown in
ences therein). The element abundances are gathered in
Fig. 8. For the specific case of the LBM, this demonstrates
Table 4, for elements from C to Si (heavier elements were
thatcurrentdataarealreadyabletoconstrainstronglythe
not used in this study). They can be compared with those
B/Cflux(asreflectedbythegoodvalueχ2 =1.06),even
min found in Engelmann et al. (1990) (see also those derived
at high energy. This provides encouraging support in the
from Ulysses data, Duvernois et al. 1996). For some ele-
discriminating power of forthcoming data. However, this
ments, the agreement is striking (F, Mg), and is otherwise
conclusionmustbe confirmedbyanalysisofa morerefined
fair. The difference for the main progenitors of boron, i.e.
3 Forinstance,itcanbeusedtopredictthepord¯background C, N, and O, is a bit puzzling, and is probably related
to a difference in the input-source spectral shape. This is
fluxtolook foradark-matterannihilatingcontribution,e.g., as
inDonatoetal.(2004).Notehoweverthatthestatisticalproce- discussedfurtherinSect.5.3,wherewealsodetermineself-
duredescribed hereis farmorerobust thanthecrudeapproach consistently the propagationparameters along with the C,
used by Maurin et al. (2001); Donato et al. (2004). N, and O abundances.