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.