Implications of Lithium to Oxygen AMS-02 spectra on our understanding of cosmic-ray diffusion

We analyze recent AMS-02 comic-ray measurements of Lithium, Beryllium, Boron, Carbon, Nitrogen and Oxygen. The emphasis of the analysis is on systematic uncertainties related to propagation and nuclear cross sections. To investigate the uncertainties in the propagation scenario, we consider five different frameworks, differing with respect to the inclusion of a diffusion break at a few GV, the presence of reacceleration, and the presence of a break in the injection spectra of primaries. For each framework we fit the diffusion equations of cosmic rays and explore the parameter space with Monte Carlo methods. At the same time, the impact of the uncertainties in the nuclear production cross sections of secondaries is explicitly considered and included in the fit through the use of nuisance parameters. We find that all the considered frameworks are able to provide a good fit. In particular, two competing scenarios, one including a break in diffusion but no reacceleration and the other with reacceleration but no break in diffusion are both allowed. The inclusion of cross-section uncertainties is, however, crucial, to this result. Thus, at the moment, these uncertainties represent a fundamental systematic preventing a deeper understanding of the properties of CR propagation. We find, nonetheless, that the slope of diffusion at intermediate rigidities is robustly constrained in the range $\delta\simeq0.45-0.5$ in models without convection, or $\delta\simeq0.4-0.5$ if convection is included in the fit. Furthermore, we find that the use of the AMS-02 Beryllium data provides a lower limit on the vertical size of the Galactic propagation halo of $z_\mathrm{h}\gtrsim3$ kpc.


I. INTRODUCTION
Over the last decade AMS-02 has acquired precise measurements of cosmic-ray electrons and positrons as well as multiple nuclei from protons to iron [1]. These data are provided with unprecedented precision at the level of a few percent in a large energy range between 1 GeV and a few TeV and have triggered numerous studies [2][3][4][5][6][7][8][9][10][11][12][13][14]. The propagation of cosmic rays (CR) in the Galaxy can be modeled as a diffusive process within a halo extending a few kpc above and below the Galactic plane. While electrons and positrons lose energy quickly and can thus be used to study the local Galactic environment (up to few kpc) CR nuclei sample a much larger volume of the diffusion halo up to about 10 kpc. Beside diffusion and energy losses, other ingredients like convective winds or reacceleration might be important. These processes can be constrained by using secondary CRs like Boron (B), Beryllium (Be) and Lithium (Li), i.e., CRs produced by the fragmentation of primary CRs, mainly Carbon (C), Nitrogen (N) and Oxygen (O), on the interstellar medium (ISM) during propagation 1 . Historically, B has been the most investigated secondary nucleus. Nonetheless, AMS-02 has now provided excellent measurements of Li and Be [15], which have been indeed studied in recent analyses. Beryllium is particularly important because its isotope 10 Be has a life time of ∼ 1 Myr, similar to the propagation time scale of CR nuclei. Precise measurements of this isotope would, thus, allow to disentangle the well-known degeneracy between the normalization of the diffusion coefficient and the height of the diffusive halo. Current measurements of this isotope are still not very precise and bound to very low energies. Nonetheless, in the absence of precise 10 Be measurements also the total Be flux might allow to draw some conclusions on the height of the diffusive halo.
Deuterium and Helium 3 are also secondaries sensitive to propagation, with 3 He also recently measured by AMS-02 [16]. Further heavier secondaries are important, in particular sub-Iron species like Scandium (Sc), Titanium (Ti), Vanadium (V), Chromium (Cr) and Manganese (Mn), whose measurements from AMS-02 are not yet available.
The aim of the present analysis is to combine the information of multiple primary CRs C, N and O and secondary CRs Li, Be and B in order to derive constraints on the diffusion properties of CR nuclei and systematically explore the role of the additional components which might impact in the process of CR propagation. We will consider two main complementary propagation frameworks. In the first, CR diffusion is described by a smoothly broken power law in rigidity with a break at about ∼ 5-10 GeV and no reacceleration. The second one, instead, includes reacceleration but it does not assume a break in diffusion, which is instead replaced by a break in the injection spectrum of primaries. Beside these two scenarios, we consider three further frameworks, making it a total of five, that extend the previous analyses in various ways. In the most general case we explore a scenario with both reacceleration, a break in diffusion, and a break in injection.
The other ingredient which we will investigate thoroughly is the systematic uncertainty in the nuclear cross-section production of secondaries. Most of these cross sections are poorly measured and thus the effect of their uncertainties on the prediction of the abundance of secondary elements is significant. This aspect has been stressed in several recent works on CR nuclei [6-11, 14, 17] as well as for CR antiprotons [18][19][20][21]. Recently, there has been some experimental effort to gather new measurements of these cross sections [22,23]. Furthermore, there is some effort to measure nuclear cross sections directly with the AMS-02 detector itself, see [1]. As we will show, taking into account these uncertainties is crucial. In particular, propagation scenarios which seem excluded at first remain viable once these uncertainties are considered. The important element of novelty which we will include in the following will be the simultaneous exploration of the propagation and cross-section uncertainties in a global fit, with the latter uncertainties included via a parametrization of the cross sections themselves. The exploration of the large joint parameter space given by propagation and cross-section parameters is a delicate task which we perform through the use of Monte Carlo scanning techniques. Furthermore, we check the impact of possible correlations in the AMS-02 data on our results.
The manuscript is structured in the following way: In Sec. II we give a more general introduction to CR propagation of nuclei in our Galaxy. Then, in Sec. III we summarize the CR data used in this analysis. Section IV explains in more detail the importance of fragmentation cross sections for this analysis. The methods of the global fits using Monte Carlo scanning techniques are detailed in Sec. V. In this work we will consider 5 different scenarios for CR propagation, which are specified in Sec. VI. Finally, results are presented in Sec. VII. More specific results on the halo half-height z h are given in VIII. Before concluding in Sec. X we discuss and compare our analysis to other works in Sec. IX.

II. COSMIC-RAY PROPAGATION
Galactic CRs are well described by diffusion within a cylindrical halo extending a few kpc above and below the Galactic plane. In general, we distinguish primary and secondary CRs. While primaries are produced and injected into the diffusion halo by astrophysical sources like supernova remnants or pulsars, the latter are produced during the process of CR propagation, mostly by fragmentation of heavier CR nuclei on the ISM. The ISM mostly consists of Hydrogen (∼ 90%) and Helium (∼ 10%) gas located in the Galactic disc. The whole process of CR propagation can be modeled by a diffusion equation in phase-space densities [24]: In addition to diffusion several further processes are incorporated in this equation. We will briefly mention and define each process in the following. The term q i (x, p) denotes the source term of each CR species i. For primaries, we factorize the source term into a space-and a rigidity-dependent part. Assuming that our Galaxy is cylindrically symmetric we obtain q i (x, p) = q i (r, z, R) = q 0,i q r,z (r, z) q i,R (R) , where r and z are cylindrical coordinates with respect to the Galactic center and R = p/Z is the rigidity. The rigidity-dependent part, q i,R (R), is modeled by a smoothly broken power-law with a break positioned at R inj,0 . The spectral indices above and below the break are denoted γ 1,i and γ 2,i , respectively: . ( The parameter s allows for a smoothing of the power law around the break: the larger s the smoother is the transition from γ 1,i to γ 2,i . For completeness, we state also the adopted spatial form of the source term, which is the same for all primaries: with parameters α = 0.5, β = 2.2, r s = 8.5 kpc, and z 0 = 0.2 kpc. 2 On the other hand, the source term of secondary CRs is given by: Here σ k+j→i denotes the fragmentation cross section of species k into i on the ISM component j. Furthermore, φ k is the primary CR flux, and n ISM,j is the density of the ISM gas component j. The rigidities of the secondary CR i and primary CR k are related by the assumption that the kinetic-energy-per nucleon (T /A) is conserved during the fragmentation, which is a good approximation [27]. In the above expression for simplicity the spatial dependence of the various quantities has been suppressed but it is taken into account in the calculations.
At the heart of Eq. (1) we find the spatial diffusion term. The diffusion coefficient D xx is modeled by a smoothly broken power law as function of rigidity R with up to two breaks: Here δ l , δ , δ h are the three power law indices below the first break at R D,0 , between the two breaks, and above the second break at R D,1 . Furthermore, s D,0 and s D,1 describe the amount of smoothing of the two breaks, respectively, and β = v/c is the velocity of the CRs. The diffusion coefficient is normalized to D 0 at a reference rigidity of 4 GV, D xx (R = 4 GV) = D 0 . The first break, if included in the model, is typically in the range of 1-10 GV. The existence of a second break at higher rigidities is suggested by the recent observations of secondary lithium, beryllium, and boron [15]. We will thus include a break around R D,1 ∼ 300 GV with a negligible amount of smoothing, s D,1 → 0. Next to diffusion, reacceleration, convection and energy losses play an important role. The amount of diffusive reacceleration is parametrized in terms of the velocity v A of Alfvèn magnetic waves as [28,29]: The parameter δ is the power-law index between the breaks from Eq. (6). We assume that the velocity of convective winds, V (x), is constant and orthogonal to the Galactic plane, V (x) = sign(z) v 0,c e z . Finally, we consider energy losses. These might be continuous, adiabatic or catastrophic. Continuous energy losses are encoded in the term dp/dt and include processes like ionization and Coulomb losses, while adiabatic energy losses originate from a nonzero gradient of convective winds. Catastrophic energy losses are encoded in the lifetime for fragmentation or decay as described by the parameters τ f and τ r . We note that the catastrophic loss terms give rise to the source term of secondaries (see Eq. (5)).
All in all, Eq. (1) provides a chain of coupled differential equations, where the flux of heavier nuclei (partly) define the source term of lighter secondary CR species. It can be solved by starting with heavier species gradually evolving to lighter nuclei. We use the Galprop code 3 [30,31] in order to solve the diffusion equation numerically. More specifically, we use Galprop version 56.0.2870 combined with Galtoollibs 855 4 as basis for our analysis and implement some custom modifications. We employ a grid in kinetic energy per nucleon and in the two spatial dimensions r and z, namely we assume cylindrical symmetry. 5 It is assumed that CRs are in a steady state. Our diffusion halo has a maximal extension of r = 20 kpc and z = ±z h . Given the well-known degeneracy between z h and D 0 we will keep z h fixed to a benchmark value of 4 kpc for most of the analysis. In a specific paragraph we will vary also z h and discuss its constraints.
Beside propagation in the Galaxy, CRs also propagate in the Solar system. This affects CR spectra only at relatively low rigidities, below ∼ 50 GV, through the effect of solar winds and the solar magnetic field. This is commonly referred to as solar modulation and varies in strength in a 22-year cycle. A propagation equation similar to the one of interstellar CR propagation, but adjusted to the heliospheric conditions, can be used to describe Solar modulation. The equation can be solved numerically [37][38][39][40]. There are also approaches to solve approximated versions with semi-anlytical methods [41]. Progress in the understanding of CR propagation in the heliosphere has been made [3,42] particularly thanks to the data of the Voyager I probe which left the heliosphere a few years ago and is thus sampling CRs fluxes before they are affected by the solar modulation. Beside this, time-dependent measurements of CR fluxes both by PAMELA and AMS-02 [43,44] are also very useful to constrain the properties of heliospheric propagation. Nonetheless, a detailed understanding is still missing at the moment. For this reason we resort to the commonly employed force-field approximation [45], where the CR flux near Earth is given by where e is the elementary charge, and Z i , m i , and ϕ SM,i are the atomic number, mass, and solar modulation potential of species i, respectively. E LIS,i is the CR energy of the local interstellar spectrum, while E ⊕,i is the energy effectively measured at the top of the Earth's atmosphere after propagation in the Solar system. In the force-field approximation ϕ SM,i is species-independent. Nonetheless, it might be necessary to consider different potentials if the measurements for different species are not simultaneous. Fortunately, the measurements provided by AMS-02 of Li, Be, B, C, N, O which we will use in the analysis all refer to the same period of time, thus we will use a single modulation potential for all of the species. N AMS-02 66 [47] O AMS-02 67 [46] B/C AMS-02 67 [15] Li/C AMS-02 67 [15] Be/B AMS-02 66 [15] TABLE II. Summary of cross section related nuisance parameters included in the CR fits. The first column contains the name of the effective fit parameter name. The second column lists the physical parameters related to the effective one and effectively varied in the fit (notation corresponding to Eq. 10). The third column states the sampling procedure, i.e., whether the parameter is sampled by MultiNest or on-the-fly. The last column reports the range of values scanned in the fit. Each nuisance parameter is only included if the product species is included in the CR fit. We assume that the fragmentation cross sections X + H → Y and X + He → Y are proportional to each other. Thus, we omit the '+H' and '+He' from the notation.
fit parameter nuisance parameters sampling prior δXS→B δ 16 In the following we will use the recent AMS-02 measurements of lithium, beryllium and boron [15], carbon and oxygen [46], and nitrogen [47]. All the data refer to the same period of time, namely the first five years of operation. This, as mentioned above, allows us to simplify the treatment of solar modulation and assume a single modulation potential for all the species.
For Li, Be and B absolute, fluxes are available. Nonetheless, we will perform fits on the ratios B/C, Li/C and Be/B since in the ratio some measurement systematics cancel out and the overall error is smaller. Furthermore, while the absolute secondary spectrum depends on the spectrum of the primary and thus on the injection parameters, the secondary over primary ratios are very weakly dependent on them, and this helps to speed-up the convergence of the fit. In Table I we report a list of the employed data. In the following also VOYAGER data on carbon, oxygen and boron-over-carbon ratio [48] are shown, but it is only used for plotting purposes and not actually included in the analysis.

IV. NUCLEAR CROSS-SECTION
The precision of fragmentation cross sections to produce secondary CRs is in many cases significantly worse compared to the precision of recent CR measurements provided by the AMS-02 experiment [17]. Uncertainties are very often at the level of 20-30% or even more in the cases of cross sections for which data are very scarce. Thus, we allow some flexibility in the default fragmentation cross sections in order to take into account the related uncertainties. More specifically, we allow for a freedom in the overall normalization and in the low-energy slope of each fragmentation cross section: Here A k+j→i is an overall renormalization factor and δ k+j→i adjusts the slope of the cross section below a reference kinetic energy-per-nucleon chosen to be T ref /A = 5 GeV/nuc. This choice for T ref /A is justified by the fact that all the cross-section models predict a break around this energy, and a flat behavior above it. As a default model for the cross section we use the Galprop parametrization. 6 The network of reactions involved in the production of the secondary species is extremely large [17]. In order to keep the analysis feasible we thus allow additional freedom only in the cross section of the main reactions contributing to the secondaries, namely the ones involving fragmentation of carbon and oxygen [17], which are the most abundant primary species. Below we briefly list the reactions considered for each species. We also provide approximate numbers for the relative isotopic abundances of each species, which result from a typical run of Galprop using the default network of cross sections.
Lithium has two stable isotopes 6 3 Li and 7 3 Li, roughly equally abundant in CRs. We thus consider the four reactions producing these two isotopes from the fragmentation of 16 8 O and 12 6 C. There is one subtlety, in fact one should consider eight fragmentation reactions: four on H and four on He. The uncertainty on the fragmentation cross sections on He are even larger than the ones involving H. On the other hand, the abundance of He in the ISM is a factor of 10 smaller than the one of H. We follow the typical assumption, as implemented in Galprop, that the two reactions on H and He are related by a rescaling factor. 7 Thus, we only introduce one common nuisance parameter for the spallation on H and He.
Beryllium has 3 isotopes 7 4 Be, 9 4 Be and 10 4 Be. We notice that 7 4 Be decays by electron capture and is very short lived in atomic form. However, the nuclear form as in CRs is stable, apart from a possible capture of electrons in the ISM. The electron capture process is implemented in Galprop and we tested that it has a negligible impact on the Be spectrum. 7 4 Be and 9 4 Be are the most abundant CR isotope making between 30% and 50% of total beryllium each. The less abundant is 10 4 Be composes only about 10% at low energies and 20% above 100 GV. Thus, we consider nuisance parameters for the production cross sections only for the two most abundant species, 7 4 Be and 9 4 Be. Boron has two stable isotopes 10 5 B and 11 5 B contributing about 1/3 and 2/3, respectively. We consider the four related production channels.
CR C is mainly made of the isotope 12 6 C with a subdominant part of up to 10% of 13 6 C, and it is mainly primary. Nonetheless, it has a small secondary component, at the level of ∼ 20 %, coming from fragmentation of oxygen. The primary part is almost entirely made of 12 6 C since 13 6 C is produced in star synthesis or primordial nucleosynthesis only at the percent level. The secondary part, though, has roughly equal abundances of the two isotopes. We thus consider the two reactions producing the two carbon isotopes from oxygen.
Nitrogen also has two stable isotopes, 14 7 N and 15 7 N and has both a primary and secondary component roughly equal in abundance below 50 GV. At higher energies 14 7 N is dominant. Similarly to what happens for C, the primary component is almost entirely made of 14 7 N since 15 7 N is only found in nature at the sub-percent level. The secondary component coming from the O fragmentation though has contributions to both isotopes. We thus consider the two fragmentation reactions.
Finally, O is almost purely the isotope 16 8 O and can be considered almost entirely as primary, with very small traces coming from the fragmentation of heavier nuclei (at the level of less than percent). Thus, we do not consider any freedom in the reactions producing oxygen, although in the network we include all the reactions producing it which involve nuclei up to silicon. In this respect we stress that all the reactions listed above are the ones in which we allow freedom in the cross sections, and which give the main contribution to the secondaries, typically at the level of 70-80%. The remaining reactions in the network (which are up to several hundreds, but give overall a subdominant contribution) are all included, although the related cross sections are fixed to the original model. Table II lists the nuclear reactions which we will vary in the fit and the related normalization and slope parameters. We note that even though we have reduced the number of reactions to a few crucial ones, the number of related parameters is still too large for the analysis to be feasible. We thus employ a further simplification and we force the parameters relative to the same element to be all equal. In this way we have in total five normalizations and five slopes, i.e. one for each element which is secondary or has a significant secondary component. This simplification is justified by the fact that, at the moment, we only have measurements of the global abundance of each element, but not yet separate spectra for each isotope. Varying independently the parameters for each isotope would thus give rise to strong degeneracies in the parameter space with no added physical value. These parameters will be sampled in the fit simultaneously with the CR propagation parameters. Just for three of them, listed in the last column of Tab. II a simplified treatment is possible. In fact, the cross-section normalizations of Li, Be and B are basically equivalent to the overall normalization of the spectra itself. Thus, they can be marginalized on-the-fly during the fit. This procedure is explained in more detail in the next section.

V. METHODS
The methods used to perform CR fits in this work are very similar to those introduced in [2] and [49]. While [2] and several followup applications [49][50][51][52][53] have focused on fitting the light CR nuclei, i.e. protons, helium and antiprotons, here we apply this methodology for the first time to heavier nuclei. In particular, we will focus on the nuclei from lithium to oxygen which have been measured precisely by the AMS-02 experiment. We will briefly remind the main ingredients of our method as detailed in [49] and then describe the necessary extensions required for this analysis.
As log-likelihood we use a simple χ 2 summed over the various CR species: where φ s,i is the experimentally measured flux of the CR species s at the rigidity R i and φ (m) s is the corresponding model prediction. The covariance matrix V (s) contains the uncertainty of the flux measurement. In the our fiducial setup we assume uncorrelated uncertainties, i.e., V where σ s,i is the error measurement reported by AMS-02 for species s in the rigidity bin i. More precisely, we sum in quadrature the statistical and the systematic error. However, we will also explore, in a few selected cases, the impact of correlated systematics in the experimental uncertainties. To derive the full correlation matrix for these cases we will follow the approach of Refs. [53]. 8 A similar approach was also used in [54]. We will not summarize the procedure here and refer the reader to Refs. [53,54] for further details. The χ 2 and φ (m) have a dependence on the model parameters that for simplicity we have suppressed in the above formula. There are two types of parameters. The first set pertains to CR propagation and is given by a total of 17 parameters. They are partly described in Sec. II, but, for convenience, we list them again below. Four parameters are used to describe the injection spectrum of primaries, i.e., the slopes below and above the rigidity break, γ 1 , γ 2 , the rigidity break R 0 and a smoothing parameter s. Nine more parameters describe propagation, i.e., the normalization D 0 and the slopes δ l , δ, and δ h of the diffusion coefficient, the break positions and their smoothing R D,0 , R D,1 , and s D,0 , 9 the velocity of Alfvèn magnetic waves, v A , and the convection velocity, v 0c . Furthermore, there are 3 parameters to determine the isotopic abundance of primary C, N, and O, measured relative to the hydrogen abundance, whose value is arbitrarily fixed to 1.06 × 10 6 . The last parameter is the solar modulation potential ϕ AMS-02 . During several preliminary fits we left these parameters completely free to vary and we found that the preferred value is about 600 MV. For better stability of the subsequent fits, we then decided to apply a Gaussian prior on ϕ AMS-02 , i.e., we add to the main likelihood the term −2 log(L SM ) = (ϕ AMS-02 − 600 MV) 2 /σ 2 ϕ where σ ϕ = 30 MV. Not all of these parameters are present in each fit. Depending on the propagation scenario under consideration (see below) only a subset might be present. These parameters, together with the prior ranges explored in the fit are reported in Tab. III.
The second set of parameters is related to the nuclear fragmentation cross sections. As described in Sec. IV we consider a total of 10 cross-section parameters. Also in this case, not all of the parameters are necessarily present in each fit. A fit with BCNO for example only will use 6 parameters, while the fit with all of the species will include all 10 parameters. As mentioned in Sec. IV the 3 normalization parameters of Li, Be and B are treated as on-the-fly marginalization parameters, as explained further below.
To scan this large combined propagation-cross-section parameter space (up to 27-dimensional in the largest case) we use MultiNest [55]. For the MultiNest setup we use 400 live points, an enlargement factor efr=0.7, and a stopping criterion of tol=0.1. As mentioned above some parameters are treated in a simplified way, namely the modulation potential and three cross-section normalizations. We profile over these four parameters together on-thefly at each MultiNest likelihood evaluation following [56]. More precisely, for each evaluation in the fit within the parameter space scanned by MultiNest the likelihood is maximized over the four remaining simplified parameters. This maximization is performed with Minuit [57]. We remark that the parameter scans performed in this analysis require a large computational effort. A typical MultiNest scan requires 0.5 to 1 M likelihood evaluations. However, depending on the complexity of the fit this number can increase significantly. For example, when the half-height of the diffusion halo is further included as free parameter (cf. Sec. VIII), the fit has to explore the additional D 0 − z h degeneracy. Thus, in some cases the number of required likelihood evaluations can exceed 3 M. A single likelihood evaluation takes between 150 and 200 sec on a single core.
MultiNest formally explores the Bayesian posterior and thus naturally provides sample from it. Nonetheless, these samples can also be used to perform a frequentist statistical analysis of the results. As default choice to interpret the scan result we use a frequentist framework, and we build one and two-dimensional profile likelihoods in the different parameters, from which we derive contours which are shown in the various figures in the following. Nonetheless, we also check the analysis results given by a Bayesian approach. We find that two interpretations give almost indistinguishable constraints indicating that the data are constraining enough to make the results basically prior-independent, and statical interpretation independent. An explicit example comparison is given in the Appendix.

VI. FIVE FRAMEWORKS TO DESCRIBE CR PROPAGATION
In Sec. II we have described a very general approach to model CR propagation in this work. However, one of the key motivations of this analysis is to check whether the new precise AMS-02 measurements of several CR secondary spectra can point to a more specific framework of CR propagation. Therefore, we define, select and study five propagation scenarios, four of which are sub-cases of the most general one. The parameters of the various scenarios together with the range of variation explored in the fits are listed in Tab. III.
The most simple framework used in this analysis, i.e., the one with the smaller number of parameters, is called BASE. In this framework there is no break in the injection spectrum of primary CRs which thus follows a simple power law in rigidity with spectral index γ 2 . On the other hand, the diffusion coefficient is modeled by a double-broken power law. The first break is at a few GVs. A break in this range is motivated by various theoretical studies [58,59].
While a second break around few hundreds of GVs is favoured by the recent measurements of AMS-02 of Li, Be and B [15]. Finally, in this scenario we also allow for convection, while reacceleration is forced to be zero. In total this model has 9 free propagation parameters. To them we then need to add C, N, and O abundances, the solar modulation potential and the cross-section uncertainty parameters.
The BASE scenario is gradually complicated by adding further ingredients. We consider the BASE+v A scenario which is identical to the BASE one but further includes reacceleration. It thus have one additional parameter. The BASE+inj scenario, instead, it goes beyond the BASE scenario including a break in the injection spectrum of primaries. It thus has three parameters more, namely γ 1 , R inj,0 and s. The most general model is given by an extension of the BASE framework by both reacceleration and a break in the injection spectrum which we dub BASE+inj+v A .It has four more parameters with respect to BASE. Finally, the last scenario we consider includes reacceleration and a break in the injection spectrum but we remove the break in the diffusion coefficient at rigidities of a few GV. We refer to this framework as BASE+inj+v A −diff.brk. This last scenario corresponds to the one explored in our previous analyses [2,[49][50][51][52][53] as well as in other works [10,60,61]. Regarding the dataset employed, every fit includes AMS-02 data on C, N and O. Fits can then be grouped into two categories: In the first one we include the B/C ratio dubbing this configuration BCNO, while in the second group we include also the Li/C ratio and the Be/B ratio, i.e., all the nuclei considered in this work, dubbing this configuration as LiBeBCNO. All in all, this gives 10 CR fits corresponding to 5 different frameworks of CR propagation and 2 CR data set configurations. Finally, only for the BASE and BASE+inj+v A −diff.brk frameworks with BCNO and LiBeBCNO data configurations we consider four more fits, to study the effect of correlated systematic uncertainties. In total we will thus perform 14 fits.

VII. RESULTS
The results of the fits are summarized in Tab. IV, which, for each fit, reports the total χ 2 and the contribution to the χ 2 for each single species, the number of degrees of freedom, the best-fit value for each parameter and the 1 σ error. We can draw several conclusions: • We can see that the fits are all very good with typical χ 2 /dof of about 0.3 for the BCNO fits and 0.45 for LiBeBCNO fits. For better readability the same results are reported in graphical form in Fig. 1, while Fig. 2 and Fig. 3, instead, show the data for the spectra of each species together with the best-fit model predictions for each framework, with the bottom panels of each figure showing the residuals. The quality of the fits can also be seen from the residuals which are all compatible with the size of the error bars themselves and thus at the level of 3-4%. Actually, the χ 2 /dof is a bit small, a result also seen in previous analyses [2,49] and related to the fact that the systematic uncertainty reported by AMS-02 is probably slightly overestimated or correlated. While this typically indicates that the derived constraints should be conservative, one can also try to be more aggressive and model more explicitly the systematic uncertainty and its energy correlation, as we did for the four fits also reported in Tab. IV and labeled "(corr)". One can see that indeed in these cases the χ 2 /dof increases to values a bit above one while the constraints on the parameters are slightly stronger. A similar behavior is indeed also observed in Refs. [49,53]. For better clarity we do not display the residuals for the fits of correlated systematic uncertainties in Fig. 2 and 3. As can be seen in the main panels, the best-fit lines of the fits with correlated uncertainties are often offset with respect to the data points. However, this is expected and due to the fact that correlations in the data allow for an offset in normalization and/or small tilts. Therefore, in these cases is more difficult or even misleading to judge the fit from the residuals and one should rather rely on the χ 2 values themselves reported in the tables and Fig. 1.
• We are able to obtain good fits with very similar best-fit χ 2 for all the five framework explored. This means that the frameworks with fewer parameters should be preferred. For example, we see that the BASE framework in the BCNO case has a best-fit χ 2 of 72.4. Adding one parameter more in the BASE+v A case leaves the χ 2 basically unchanged. Adding three parameters more in the BASE+inj case improves the χ 2 only by ∼ 4.5, while with four parameters more in the BASE+inj+v A case the χ 2 improves only by ∼ 5.0. Using a likelihood ratio test, which in this case is equivalent to calculating the ∆χ 2 , shows that all the improvements are at the level of about only 1 σ, and thus not significant. The BASE framework and the BASE+inj+v A −diff.brk cannot be compared directly since they are not nested models. We can, however, compare the BASE+inj+v A −diff.brk and BASE+inj+v A cases, which are nested. In this case, BASE+inj+v A has three parameters more than BASE+inj+v A −diff.brk and the χ 2 is better by ∆χ 2 = 7. Based on the likelihood ratio and a χ 2 distribution with 3 degrees of freedom, this is at the level of less than 2 σ, and thus again not significant. All in all, the two simplest models compatible with the data which emerge from this discussion are the BASE one and BASE+inj+v A −diff.brk. Checking the χ 2 values of the LiBeBCNO cases one can see that similar conclusions are reached. Also assuming correlated systematic uncertainties does not change these conclusions. Indeed, as can be seen from Fig. 5 which collects all the 2D contours related to the propagation parameters as well as the 1D profile χ 2 , the parameters of the BASE and BASE+inj+v A −diff.brk framework are generally well-constrained, while the more general framework BASE+inj+v A , where all the parameters are considered, only appears to produce several degeneracies among the parameters and loose constraints, indicating, indeed, that the model has too many parameters with respect to the constraining power of the data. These two models, as we already mentioned, are quite complementary and physically very different, since in the BASE case a break in diffusion is present at few GV which is not present in the BASE+inj+v A −diff.brk framework. This break is the crucial feature which in the BASE case allows us to fit the flattening and turning down observed at few GVs in the spectra of secondary-over-primary ratios as can be seen in Figs. 2-3. At the same time, contrary to the BASE case, the BASE+inj+v A −diff.brk framework includes reacceleration, which is in this case the important feature allowing to fit the low-energy secondary-over-primary ratios. It is quite remarkable that these two very different models can both provide good fits and are thus both valid models. Still, there is actually an important catch to this result regarding cross-section uncertainties, as discussed further below.   Comparison of AMS-02 measurements with the best-fit spectra for the BCNO data set. The solid lines correspond to best-fit model for each of the five CR propagation frameworks. The dotted lines are the best-fit models for the BASE and the BASE+inj+vA−diff.brk scenario with correlated systematic uncertainties. Dashed lines represent the interstellar CR spectra, i.e., the spectra before solar modulation is applied. Residuals are only shown for the cases with uncorrelated AMS-02 uncertainties. For comparison we display Voyager data at low energies.

A. Results on CR propagation parameters
• From Fig. 1 we can see that all the propagation parameters are compatible in all the 14 fits, which is quite reassuring regarding the robustness of the result. There is one notable exception, regarding the D 0 − δ plane which we will discuss in more detail below. Apart from this, it means that the three secondary nuclei Li, Be and B point to compatible propagation parameters and can be combined in a single global fit, a result which is certainly not trivial.
• As mentioned above, the BASE and BASE+inj+v A −diff.brk models point to compatible parameter values for all parameters they have in common. At first glance, though, there seems to be a slight offset in the D 0 − δ contours, as shown in the triangle plot in Fig. 5. However, this is mostly related to the normalization rigidity scale of the diffusion coefficient chosen at 4 GV, which is very close to the rigidity of the break in the BASE framework. In Fig. 4 we show the diffusion coefficient as function of rigidity for the BASE, BASE+inj+v A ,  and BASE+inj+v A −diff.brk models. Above 10 GV the diffusion coefficients of all models agree well within the uncertainty. Furthermore, the diffusion coefficients derived from the BCNO and LiBeBCNO data sets match well, underlining again the compatibly of the three secondary CR species considered in this work. We note that Fig. 4 is produced for a fixed z h of 4 kpc. If instead z h were allowed to vary in the fits, the uncertainty band around the best-fit diffusion coefficient would be enlarged significantly due to the well-known degeneracy between D 0 and z h .
• We can see that for all of the fits convection is compatible with zero. This is actually due to a strong degeneracy between v 0,c , D 0 , δ and, partly, s D,0 , at least for the models with a break in diffusion, as can be seen in Fig. 5. For this reason, values of v 0,c up to 50 km/s are still possible at the price of larger ranges allowed for D 0 and δ. The framework BASE+inj+v A −diff.brk, instead, prefers no convection allowing values only up to ∼ 10 km/s. At the moment, it is thus more economical to consider a model with no convection, since it fits the data equally well. In this case lower values of D 0 are selected together with higher values of δ (since the degeneracy between v 0,c and δ has a negative slope). We note, nonetheless, that we are exploring only a simple model of convection with constant winds, but more complex models with wind gradients as function of z h are possible and might give different results. A systematic study of convection goes beyond the aims of the current analysis.
• Due to the degeneracy discussed above for the BASE case, δ is constrained in the fairly large range 0.4-0.5 when convection is included, while assuming zero convection would give a preference for δ 0.5 with a quite small error of the order of 0.01. For the BASE+inj+v A −diff.brk framework this degeneracy is weaker and the error on δ is thus smaller, but smaller values of δ are preferred, namely δ 0.45 and δ 0.41 for the BCNO and LiBeBCNO fits, respectively. We thus see that without convection the framework with reacceleration tends to prefer smaller values of δ. These values of δ are in line with the results of [2,49] where the BASE+inj+v A −diff.brk model is analyzed, but with different CR data, namely proton, helium and antiproton. Also in this case, it is remarkable and not trivial that light nuclei and heavier nuclei point to a compatible propagation scenario. Other recent analyses [5] of AMS-02 data, which consider a break in diffusion, also tend to find higher values of δ with respect to a scenario with reacceleration and no break in diffusion. From the point of view of diffusion theory this means that the latest AMS-02 data point to Kraichnan (δ ∼ 0.5) model of turbulence while Kolmogorov turbulence (δ = 1/3) seems disfavored.
• The models with a break in diffusion find a position of the break at around 4 GV, compatible with the predictions

FIG. 5.
Triangle plot with the fit results for the CR propagation frameworks BASE (blue), BASE+inj+vA (green), and BASE+inj+vA−diff.brk (red), and using the BCNO data set. The triangle plot shows only the subset of fit parameters related to CR propagation. The 2D plots show the 1σ and 2σ contours for each combination of two parameters, which are derived from the 2D χ 2 -profiles. The diagonal shows the 1D χ 2 -profiles (the y-axis is always the ∆χ 2 and ranges from 0 to 10).
of theoretical studies [58,59]. The slope δ l below the break is not well constrained. There is a preference for a negative value, but, conservatively only an upper limit of about δ l < ∼ 0 can be derived. A negative value of δ l is present, for example, in models of turbulence with damping of Alfven waves at low energies [58]. All of the frameworks include a second break in diffusion at higher energies and this is robustly detected in all of the fits, with a break around 200 GV and amount of breaking of about ∆δ ∼ 0.15. This confirms the findings of the AMS-02 analysis [15] where this result is derived modeling directly the observed spectra rather than from a diffusion model. It is also in agreement with the results of [5]. In absolute value the slope of diffusion at high energy becomes δ h ∼ 0.33, which seems to indicate a transition from Kraichnan turbulence to Kolmogorov   6. Same as Fig. 5, but for the subset of cross section nuisance parameters.
turbulence at few hundreds GVs.
• Reacceleration is a common ingredient in CR propagation models. However, typically, the amount of reacceleration required to fit the CR data is fairly large and might be in disagreement with energy considerations, as discussed recently in [62]. It is thus important to understand whether CR spectra do indeed require reacceleration for a good fit. In the BASE+inj+v A −diff.brk framework reacceleration is quite tightly constrained with a value of about 20 km/s. The BASE framework does not include reacceleration but is nonetheless able to fit the data well. The more general framework BASE+inj+v A where both reacceleration and break in diffusion are included also allows values of v A up to 20 km/s, but due to the presence of degeneracies, v A is also compatible with zero. To understand if a certain amount of reacceleration is indeed physically present would thus require us to demonstrate the model BASE+inj+v A −diff.brk is preferred over the ones without reacceleration. However, at the moment, both models provide a good fit to the data, thus a definitive conclusion is not possible yet.
• The abundances of primaries C, N, O are all well-constrained in the fit. There is, nonetheless, an expected degeneracy between the C and O normalization. No clear correlation between the abundances of the primaries at the source and CR propagation parameters is observed, although it would be expected at some level. This is probably due to complicated degeneracies in the high-dimensional parameter space, which blur single 2D correlations. However, we note that there is an anti-correlation with the B production cross-section A XS → B.
The values of C and O abundances preferred by the fit are about 20-30% larger than typical values used in Galprop [63]. While the N abundance is larger by a factor of ∼ 2. The result reflects the fact that the C and O flux normalizations measured by AMS-02 are about 20-30% larger than pre-AMS-02 determinations (e.g. [64], [46] and references therein). These abundances can be also compared with measurements from the Sun or meteorites [65,66]. In this respect the O/N∼ 10 ratio found here is very similar to the value observed in the Solar system given in [65,66].

B. The role of cross-section uncertainties
• A crucial role is played by the cross-section uncertainties that we have included in the fits. The triangle plots containing the 2D contours and 1D χ 2 profiles of the relevant parameters are shown in Fig. 6 (and also for the LiBeBCNO in the appendix: Fig. 9). One can see that in all of the frameworks containing a break in diffusion, and thus in particular in the BASE framework, all of the cross-section parameters are compatible with having no deviations from the benchmark cross-section model (with a couple of notable exceptions to be discussed better below), i.e., the normalizations A are compatible with 1 and the slopes δ XS are compatible with 0. Thus, this means that the CR spectra can be fit while only assuming the benchmark cross-section model, which is in this case the Galprop one, without any modification. The important point to notice is that this is not valid for the reacceleration model BASE+inj+v A −diff.brk. In this case significant deviations are seen. For example, the slope for B production needs to be modified to a value δ XS → B 0.15 and it is incompatible with zero. Results are similar for the other slopes: δ XS → N 0.1, δ XS → Be 0.25, and δ XS → Li 0.2. The normalization parameters A, on the other hand, are compatible among all of the models. The conclusion is that it would not have been possible to achieve a good fit for the BASE+inj+v A −diff.brk framework if cross-section uncertainties were not included in the analysis. However, we cannot conclude that this favors the BASE model, since these deviations are only at the level of 20% and are within the systematic uncertainties present in the cross sections. Thus, the conclusion is rather the fact that the cross-section uncertainties are at the moment significantly limiting the inference of ISM propagation properties from CRs. If the uncertainties were much smaller, at the level of few percent, i.e., at the same level as the AMS-02 uncertainties, we would be able to select a specific propagation scenario over the others. As also discussed in other works [6,7,9,11,12,[17][18][19], this further confirms that a community effort is needed in order to carry out specific campaigns of measurements to reduce nuclear cross-section uncertainties, in such a way to have then the possibility to exploit the full potential of CR data to study the properties of the ISM.
• As mentioned above, some deviation of the cross-section parameters from the benchmark model is present in all of the frameworks explored. The most evident case is the normalization of the cross-section for the secondary production of carbon, A XS → C, which always saturates the lower bound of our prior, i.e., a value of 0.5. Although uncertainties in the cross sections are typically quite large, a downward deviation at the level of 50% or more of the cross sections involving the production of carbon, in particular the ones coming from fragmentation of oxygen, seems unlikely. This deviation thus seems to indicate some underlying problem in fitting carbon. A possible solution is to allow different injections for carbon and oxygen. We verified with a specific fit that in this case A XS → C becomes indeed again compatible with 1, and the fit gives slightly different injections for C and O. A different injection for C and O may be justified in the light of the fact that also H and He are already observed to give indications of different injections. We stress, nonetheless, that in both cases this has a negligible impact on the propagation parameters, which is the main focus of this analysis.
Another parameter which presents a systematic deviation from the benchmark in all the fits is the normalization of the Li production cross section, A XS → Li, with a preferred value of 1.3 ± 0.15. In this case, however, the deviation is within the expected uncertainties and thus acceptable. Similarly, although less pronounced, a systematic deviation is present in A XS → N with a value of 1.1 ± 0.05, again perfectly acceptable.

VIII. CONSTRAINTS ON THE HALO HALF-HEIGHT z h
In the previous section we have assumed a fixed halo half-height of z h = 4 kpc, which is justified because of the well-known degeneracy between z h and D 0 . As discussed in the introduction, this degeneracy can be broken with the use of radioactive clocks (most notably 10 4 Be). Measurements of 10 4 Be from AMS-02 are not yet available, however, the total Be spectrum measured by AMS-02 has small error bars at the level of few percent. The contribution of the 10 4 Be isotope to the total Be spectrum is at a level between 10% and 20%. So, given the small uncertainties of AMS-02 measurement, the total Be spectrum might be able to break the D 0 -z h degeneracy. In this section we thus explore the possible constraints on z h . To this end we perform a total of 8 additional fits using the same setup as in the previous fits but with a further free parameter, namely z h , left free to vary in the range between 2 kpc and 10 kpc. For these fits we only consider the BASE and the BASE+inj+v A −diff.brk frameworks with the two data sets BCNO and LiBeBCNO. Furthermore, we consider both the case with uncorrelated uncertainties and the case with correlated ones. In total this gives eight different fits to explore.
The χ 2 -profiles as function of z h for all fits are shown in Fig. 7. As expected, the BCNO data set is not able to provide constraints on the halo height, and the χ 2 -profiles have a flat behavior and do not rise above the 2σ level both for the BASE and BASE+inj+v A −diff.brk scenario. The picture changes with the LiBeBCNO dataset, i.e., when the Be spectrum is included also in the fits. In the BASE scenario (left panel of Fig. 7) a robust lower limit on z h > ∼ 4 kpc at 2σ (i.e., ∆χ 2 = 4) can be derived, which is valid both when correlations in the data are considered and when they are not included. No upper bounds on the halo height can be placed. The situation is slightly different for the BASE+inj+v A −diff.brk scenario. With uncorrelated uncertainties a lower bound on z h > ∼ 3 kpc at 2σ is still present, which is a bit less stringent than the BASE case. With correlated uncertainties, however, the lower bound appears to further weaken going down to 2 kpc. A weakening of the constraints when correlated uncertainties are considered is in principle possible since correlations introduce some freedom in the normalization and long rigidity-scale tilt in the data which can be somehow degenerate with the effect of varying some parameters, such as z h . Indeed, this effect seems to be present also in the BASE case, although in that case it seems to be less prominent. Finally, an upper bound of z h < ∼ 8 kpc at 2σ is given in the correlation case, which also does not appear to be very robust since it is not present in the uncorrelated case. Overall, we conclude that the constraints on z h appear to be both propagationscenario dependent and dependent on the assumptions regarding the correlation properties of data systematics, and thus they should be taken with some care. Nonetheless, from the above analysis we can infer a strong indication, although not a conclusive one, for a lower limit of about 3 kpc on z h . The situation might likely change and become more constraining once AMS-02 data on the 10 4 Be isotope becomes available.

IX. DISCUSSION AND COMPARISON WITH OTHER WORKS
In this section we discuss and compare our results with previous works. The analysis in Refs. [6,7], similarly to this work, takes into account nuclear cross-section uncertainties. A main difference is the fact that we include primaries in the fit to constrain the diffusion framework while Refs. [6,7] consider only secondaries for this purpose. Another difference is the use of a simplified semi-analytical framework for propagation as opposed to our fully numerical treatment through the Galprop code. Nonetheless, overall our results are in reasonable agreement. We both find a value of δ ∼ 0.5 in models without reacceleration and with a break in diffusion, and a tendency to get a smaller value of δ in the complementary case with reacceleration but without break in diffusion. A further common conclusion is that the inclusion of cross-section uncertainties is crucial to achieve a good fit and that sizable deviations from the benchmark models are observed. Ref. [7] constrains the z h to 5 ± 3 kpc at 1σ. The lower limit is in agreement with our findings.
In Refs. [8,9,14] a simplified analytical approach similar to that of Refs. [6,7] is used. They perform the fit only above 20 GV to exclude the rigidity region most strongly affected by Solar modulation. They do not include nuclear cross section uncertainties, although they employ an updated parametrization from a refit of the relevant processes. They find a constrain on z h in the range 3-6 kpc when using only the AMS-02 statistical error bars. When systematic uncertainties are included, the constraint is softened to approximately z > ∼ 2 kpc. They find a best fit for δ = 0.54 (in Ref. [9] fitting to B/C, Be/C, B/O, and Be/O data) and δ = 0.63 (in Ref. [8] fitting to B/C, C, N, and O data) which is in principle incompatible with our findings, although they do not report any error so a proper comparison is not possible.
In Refs. [10,11] the authors use a propagation framework equivalent to the BASE+inj+v A −diff.brk one together with a fully numerical treatment both for Galactic propagation using the Galprop code and for Solar System propagation using the Helmod code. Nuclear cross-section uncertainties are not included. The main finding is the need to consider primary Li in order to achieve a reasonable fit for Li CR data. This is possibly related to the choice of not allowing freedom in the cross-section normalization. As we have shown, we are able to fit Li well but with some substantial deviation in the Li production cross section from the benchmark model. A similar deviation for the Li cross-section is indeed found also in Refs. [6,7] and Refs. [12,13]. Finally, they find a value of δ = 0.415 ± 0.025 which is compatible with the value δ = 0.41 ± 0.01 we find for our LiBeBCNO fit in the BASE+inj+v A −diff.brk framework.
In Refs. [12,13], the authors study the CR spectra of Li, Be, and B also taking cross section uncertainties into account. They use the Dragon2 code to model CR propagation in a framework similar to BASE+inj+v A −diff.brk. To model cross section uncertainties they compare the two default cross section options of Dragon2 and Galprop and, furthermore, they iteratively adjust the normalization of the secondary production cross sections within uncertainties to obtain a better fit to the AMS-02 data. They confirm that different models of cross sections have a significant impact on the final predicted CR spectra. In agreement with our results, they also report a preferred normalization of the Li production cross section around ∼1.3 when using the Galprop model. As in Refs. [6,7] they use only secondary species to constrain propagation. Regarding the slope of diffusion they find a value of ∼ 0.4 in agreement with what we find for the BASE+inj+v A −diff.brk case.

X. SUMMARY AND CONCLUSIONS
In the present work, we have performed global fits on the recently published AMS-02 measurements of secondary nuclei Li, Be, B, primaries C and O and mixed secondary/primary N, in order to study the propagation properties of CRs in the ISM. We consider five different propagation frameworks differing with respect to the inclusion of a diffusion break at few GVs, the presence of reacceleration, and the presence of a break in the injection spectra of primaries. The first one, which we dub BASE, has the injection spectrum of primaries modeled as a simple power law, no reacceleration and a break in the diffusion coefficient at about few GVs in rigidity. The BASE+v A one, as the name suggests, enlarges the first one including reacceleration. Similarly, BASE+inj includes a break in the injection of primaries, while BASE+inj+v A includes both of them and it is the most general model. Finally, the model BASE+inj+v A −diff.brk has reacceleration, a break in injection also at a few GVs, but no break in the diffusion coefficient. In particular, among these, BASE and BASE+inj+v A −diff.brk are the "minimal" two, in the sense that they have the minimal number of parameters. Furthermore they are also the ones which differ the most from the physical point of view. In the rest of the conclusions we will thus focus on discussing these two models.
It is well known, however, that the production cross section of secondary CRs are affected by sizable uncertainties and this can lead to a bias in the inferred propagation properties. To properly take into account this uncertainty we consider nuisance parameters in the cross sections in the form of overall normalizations and tilts in the low energy regime below few GeVs. We thus extend our formalism in order to perform a global fit where both propagation parameters and cross-section nuisances are included together so that to properly consider all these uncertainties and their correlations. We use Monte Carlo scanning techniques in order to handle the large (> 20) dimensionality of the parameter space. Finally, a further uncertainty we explore is the presence of correlation in the systematic error bars of AMS-02, a possibility which has been discussed recently.
The main result of our analysis is that all the five propagation scenarios we examine are able to provide good fits to the AMS-02 data, including both the BASE and BASE+inj+v A −diff.brk scenarios, which is remarkable, since, as stressed above, they describe quite different physical scenarios. Crucial to this result, however, is the fact that we have taken into account the cross section uncertainties. While in the BASE scenario the preferred values of the cross-section nuisance parameters are in line with the default cross-section model we use, namely the Galprop model, in the BASE+inj+v A −diff.brk scenario sizable deviations are required, especially in the tilt at low energy, at the level up to 20-30%. However, these larger cross sections are currently within the level of uncertainties of these cross sections, and thus acceptable. Thus, the conclusion is that, at the moment, the nuclear cross section uncertainties are significantly limiting our capacity to study the properties of the ISM and the propagation of CRs. In fact, since the AMS-02 data are in principle stringent enough to discriminate among different scenarios, the cross section systematics are able to introduce enough uncertainty to make them degenerate. It is thus desirable to have experimental efforts in order to reduce these uncertainties with proper measurement campaigns in the laboratory.
Nonetheless, several conclusions can still be drawn which are robust to the inclusion of the cross-section uncertainties. For example, the slope of the diffusion coefficient at intermediate energies is well constrained in the range δ 0.45−0.5 in models without convection, or δ 0.4 − 0.5 if convection is included in the fit, pointing to a Kraichnan model of turbulence of the ISM. We also find that convection is not required in any of the five considered CR propagation scenarios, although we did not perform a systematic exploration of different convection models but we only considered a simple model with constant convective winds. The BASE+inj+v A −diff.brk requires reacceleration with Alfven waves of about v A 20 km/s. However, since the BASE model does not require reacceleration and fits the data equally well, at the moment the data are inconclusive about the presence of reacceleration. This is an example of the physics limitation introduced by the nuclear cross-section uncertainties. Finally, since Be data are sensitive to the vertical size of the propagation halo z h , due to the contribution of the radioactive isotope 10 4 Be, in the fits including Be data we are able to derive a lower limit on z h of about z h > ∼ 3 kpc at 2σ. However, while all of the previous results are robust with respect to the inclusion of correlation in the AMS-02 data systematics this last result gets weakened to z h > ∼ 2 kpc when correlations are included, especially in the case of BASE+inj+v A −diff.brk. On the other hand, the BASE scenario provides a more robust and stringent lower limit of z h > ∼ 4 kpc. Given the dependence of these constraints on the propagation scenario and the correlation properties of the data, the Be spectrum does not seem optimal to study z h . Thus, precise measurements of the 10 4 Be isotope are awaited, while a better understanding of the correlation properties of the AMS-02 uncertainties is desirable.

Appendix: Supplements
We list in this appendix further plots related to the main analysis described in the main text. The first two, Figs. 8 and 9 are the equivalents of Figs. 5 and 6, respectively, but for fits with the LiBeBCNO data sets. The results and conclusions are in agreement with those from the fits on the BCNO data set, as discussed in the main text. We note that Figs. 9 contains more cross section nuisance parameters compared to Figs. 6, namely, those related to the Be and Li spectra. As expected, we can observe that the slopes δ XS as well as the normalizations A XS of Li, Be, and B are correlated. Correlation is particularly pronounced between B and Li. This might possibly be due to the on-the-fly sampling strategy applied to these two parameters.       FIG. 9. Same as Fig. 5, but for the fits with the LiBeBCNO data set, and for the subset of cross section nuisance parameters. Figure 10 contains additional results related to the fits which include the halo half-height z h as free parameters (c.f. Sec. VIII). We show the 2D uncertainty contours for all other CR propagation parameters with z h for the BASE and BASE+inj+v A −diff.brk setup. Results are shown both for the BCNO and LiBeBCNO data set. We note that the well-known degeneracy between the normalization of the diffusion coefficient D 0 and halo half-height z h is less pronounced for the BASE framework. The reason for this is that in the BASE scenario the diffusion coefficient contains a break at around 4 GeV. Therefore, D 0 , which is normalized at a rigidity of 4 GV, is much more degenerate also with other parameters like s D,0 and δ l . Consequently the correlation between D 0 and z h is smeared out in the BASE framework, while while it remains well visible for BASE+inj+v A −diff.brk one. v0,c [k Finally, we show an example in which we compare the CR propagation parameters obtained in the Bayesian or frequentist statistical interpretation, namely the BASE framework of CR propagation and the fit to the BCNO data set. In Fig. 11 we compare the two statistical interpretation. For each couple of parameters the frequentist 1 and 2 σ contours (blue) are derived from the 2D χ 2 -profiles. On the other hand, the Bayesian 1 and 2 σ contours (amber) are derived from the marginalized posterior function and correspond to the 68.3% and 95.5% credible level, respectively. On the diagonal we show the marginalized 1-dimensional posterior function for the Bayesian interpretation which is compared to the 1-dimensional profiled likelihood of the frequentist interpretation. As can be seen, the 1D profiled likelihood and 1D marginalized posterior on the diagonal of Fig. 11 match extremely well. We thus conclude that the constraints derived in this work do not depend on the statistical interpretation. We remark, however, that 2D Bayesian contours are slightly larger compared to the 2D frequentist contours meaning that the Bayesian interpretation is a bit more conservative. Triangle plot with the fit results for the CR propagation frameworks BASE and using the BCNO data set. We compare the frequentist (blue) and Bayesian (amber) interpretation. The triangle plot shows only the subset of fit parameters related to CR propagation. The 2D plots show the 1σ and 2σ contours for each combination of two parameters and the diagonal shows the profiled likelihood and marginalized posterior for the frequentist and Bayesian interpretation, respectively.