---
abstract: 'This article aims at establishing new benchmark scenarios for Galactic cosmic-ray propagation in the GV-TV rigidity range, based on fits to the AMS-02 B/C data with the propagation code. We employ a new fitting procedure, cautiously taking into account data systematic error correlations in different rigidity bins and considering Solar modulation potential and leading nuclear cross section as nuisance parameters. We delineate specific low, intermediate, and high-rigidity ranges that can be related to both features in the data and peculiar microphysics mechanisms resulting in spectral breaks. We single out a scenario which yields excellent fits to the data and includes all the presumably relevant complexity, the model. This model has two limiting regimes: (i) the  model, a minimal diffusion-only setup, and (ii) the  model, a convection-reacceleration model where transport is tuned by non-relativistic effects. All models lead to robust predictions in the high-energy regime ($\gtrsim10$ GV), i.e. independent of the propagation scenario: at $1\sigma$, the diffusion slope $\delta$ is $[0.43-0.53]$, whereas $K_{10}$, the diffusion coefficient at 10 GV, is $[0.26-0.36]$ kpc$^2$ Myr$^{-1}$; we confirm the robustness of the high-energy break, with a typical value $\Delta_h\sim 0.2$. We also find a hint for a similar (reversed) feature at low rigidity around the B/C peak ($\sim 4$ GV) which might be related to some effective damping scale in the magnetic turbulence.'
author:
- 'Y. Génolini'
- 'M. Boudaud'
- 'P.-I. Batista'
- 'S. Caroff'
- 'L. Derome'
- 'J. Lavalle'
- 'A. Marcowith'
- 'D. Maurin'
- 'V. Poireau'
- 'V. Poulin'
- 'S. Rosier'
- 'P. Salati'
- 'P. D. Serpico'
- 'M. Vecchi'
bibliography:
- 'BCprd.bib'
title: 'Cosmic-ray transport from AMS-02 B/C data: benchmark models and interpretation'
---

Introduction {#sec:intro}
============

The last decade in direct cosmic-ray (CR) detection experiments has been characterized by a major improvement in the precision of the data available, and by an extension of the covered dynamical range [@AhnEtAl2009; @PanovEtAl2009; @AdrianiEtAl2014a; @Ting2013; @AdrianiEtAl2017; @AmbrosiEtAl2017]. In particular, with the AMS-02 data the community has to deal for the first time with percent level precision and a welcomed redundancy in the measurements.

But, as well known, [*great responsibility inseparably follows from great power*]{} [^1]: since theoretical predictions are very far from attaining that level of precision, both due to ignorance of the detailed underlying microphysics (CR acceleration and transport) and because of irreducible limitations (e.g. due to the intrinsic stochasticity of the sources [@Mertsch2011; @GenoliniEtAl2017a]), a preliminary question that should be addressed is that of the best strategy to take advantage of such a wealth of data. In this paper, we primarily focus on the AMS-02 B/C data [@AguilarEtAl2016a] and investigate how much they can constrain CR transport, aiming at defining new benchmark models.

An ambitious approach would be to proceed with global fits of all available data, attempting an overall and simultaneous understanding of CR sources (all species) and propagation. However, this approach is prone to mixing uncertainties of different nature, with the risk of devaluing the actual strength of the data by introducing poorly controlled parameters (for an illustration, see [@CosteEtAl2012]). Since our current understanding of CR measurements has more firm elements in the propagation part than in the source one, factorizing out propagation effects from source effects, while inspecting their physical plausibility a posteriori, seems justified. In this article, we proceed through partial tests of key aspects of the current propagation paradigm, with the goal of validating it or highlighting its breakdown. An important and relatively new issue in this area is that systematic errors are often dominant over statistical ones. This requires a change of perspective in well-established practices of analyzing the data, as well as new standards of rigor. It calls for establishing a satisfactory protocol for analyzing the data on a relatively simple and homogeneous data sample.

This article represents an important pillar in our CR data analysis based on the overall philosophy sketched above. As long established [@Owens1976a; @GinzburgEtAl1980; @ProtheroeEtAl1981; @StrongEtAl1998; @MaurinEtAl2001], a flux ratio of elements present but in traces in the solar system material and interstellar medium (ISM), such as Lithium, Beryllium, Boron (“secondaries”), to abundant species like Carbon or Oxygen (“primaries”) is extremely sensitive to propagation parameters. It was shown to be also almost insensitive to the energy spectrum of the injected primary species, notably if those are described by a common power-law in rigidity [@MaurinEtAl2002a; @PutzeEtAl2011; @GenoliniEtAl2015], a rather generic prediction of studies of CR acceleration at sources [@MalkovEtAl2001; @CaprioliEtAl2014; @Amato2014a; @MarcowithEtAl2016]. In particular, at high rigidities we expect the B/C ratio (currently the most precisely measured) to be dominantly affected by diffusive propagation and nuclear cross sections.

In [@GenoliniEtAl2017], we performed an analysis of the high-rigidity range of the AMS-02 B/C ratio [@AguilarEtAl2016a], finding evidence for a diffusive origin of the observed spectral break, at the same rigidity scale inferred from a similar feature in the proton and helium CR fluxes [@AguilarEtAl2015], i.e. $\sim 300$ GV. Actually, recent years have been characterized by the observational establishment of “spectral anomalies” (for reviews, see [@Serpico2015; @Serpico2018]), in particular of spectral breaks in primary species [@PanovEtAl2009; @AhnEtAl2010; @AdrianiEtAl2011a; @AguilarEtAl2015; @AguilarEtAl2015a; @AguilarEtAl2017]. In turn, there has been growing evidence in favor of their interpretation in terms of a high-rigidity break in the diffusion coefficient [@GenoliniEtAl2017; @ReinertEtAl2018; @XueEtAl2019], notably after the first AMS-02 publications of nuclear CR fluxes [@AguilarEtAl2017; @AguilarEtAl2018]. In this article, we move several steps beyond our previous analysis [@GenoliniEtAl2017], presenting a [*complete*]{} analysis aiming at constraining CR propagation and at proposing new benchmark setups: First, we rely on an improved analysis of the B/C data by the AMS-02 collaboration [@AguilarEtAl2018]. We further benefit from additional data on the primary species to constrain the break independently from the B/C ratio—using the C and O fluxes [@AguilarEtAl2017] which are most contributing species to B production [@GenoliniEtAl2018], but were not available to Ref. [@GenoliniEtAl2017]. Second, we follow the new methodology proposed in [@DeromeEtAl2019] to analyze the AMS-02 data, carefully accounting for a number of subtle (but highly important) effects which are usually ignored, notably the (partial) correlations in systematic errors. This approach enables a straightforward and sound statistical interpretation of the models (best-fit models have $\chi^2/{\rm dof}\sim 1$), also allowing for their inter-comparison. Third, we propose a new generic propagation model (dubbed  in the following), with a number of parameters that should be sufficient to describe all key features currently present in the data. In addition to a high-rigidity break, a modification of the diffusion coefficient at low rigidity is enabled ($\lesssim 5$ GV), with two limiting cases (dubbed  and ): this allows us to assess the relative discriminating strength of the data in this energy range and to shed new light on propagation in the low-rigidity regime, where a second diffusion break might be present.

The paper is organized as follows: In , the one-dimensional (1D) propagation model and the essential physical effects involved in CR propagation are presented, before introducing our three benchmark scenarios (, , and ). In , we describe the specific iterative procedure used for the B/C analysis, and this procedure is validated and checked in two appendices: assesses the robustness of the derivation of the high-energy break, by taking advantage or not of the C and O fluxes; further discusses the dependence of the fit parameters upon the lower rigidity cut, to better illustrate and give meaning to the terms “low-rigidity” and “high-rigidity”parameters. Best-fit results for our three scenarios are presented in , where these scenarios are also tentatively interpreted in terms of the underlying microphysics. In , we report our conclusions and mention natural follow-up works. Note that all results are obtained for a 1D model of our Galaxy, with the size of the diffusive halo fixed. In order to allow for a broader usage of our results, as they may have some consequence in predicting the fluxes of other secondary species, provides a scaling of the high-rigidity parameters with $L$, whereas reminds the reader of the effectiveness of this description, and provides a “dictionary” to interpret the results in terms of a two-dimensional (2D) model with different halo sizes.

Transport models {#sec:def_bench}
================

In this section, we introduce the generic propagation equation that we further solve semi-analytically in the framework of the  code [@Maurin2018] – for fully numerical frameworks complementary to ours, we refer the reader to Refs. [@StrongEtAl1998; @EvoliEtAl2008; @EvoliEtAl2017; @Kissmann2014; @KissmannEtAl2015]. We also set a generic CR transport configuration motivated by theoretical arguments on the microphysics of CRs, with a focus on possibly important low-energy processes. This generic setup will itself be used as a benchmark configuration, of which we shall explore two limiting regimes. These three cases characterize new benchmark models (dubbed , , and —see ) that are aimed at capturing different theoretical assumptions, while still being data driven. These configurations will be shown to provide excellent fits to the current B/C data assuming simple power-law primary CR spectra. Not only may this stimulate further microphysical interpretations, but it also offers a basis for a description of other CR data, like for instance the positron [@DelahayeEtAl2009; @DelahayeEtAl2010; @BoudaudEtAl2017a] and antiproton fluxes [@DonatoEtAl2001; @GiesenEtAl2015; @ReinertEtAl2018; @BoudaudEtAl2019]. Eventually, these benchmarks will be instrumental in characterizing and hopefully reducing theoretical uncertainties entering searches for exotic phenomena ( [@DelahayeEtAl2008; @BoudaudEtAl2015a; @BoudaudEtAl2015; @BoudaudEtAl2017; @BoudaudEtAl2018; @BoudaudEtAl2018b; @ReinertEtAl2018]).

The [*min*]{}, [*med*]{}, and [*max*]{} benchmark values proposed in Ref. [@DonatoEtAl2004], all based on the same model, were defined to roughly bracket the theoretical uncertainties on dark matter-induced antiproton flux predictions assuming the best-fitting propagation parameters of Ref. [@MaurinEtAl2001]. However, these values have been challenged by a series of complementary constraints [@LavalleEtAl2014; @BoudaudEtAl2017a; @ReinertEtAl2018], and are anyway no longer consistent with the B/C data [@GenoliniEtAl2017]. The revised reference models we propose here rely on different assumptions on the microphysics of CR transport instead of different values of parameters within the same configuration. This change in philosophy stems from the fact that with increasing precision in the (multiwavelength and multimessenger) observational data and improvements on the theory side, we expect to arrive soon to a much better understanding and description of the CR microphysics itself than it was possible two decades ago [@MarcowithEtAl2016; @AmatoEtAl2018].

Transport description {#ssec:transport}
---------------------

### Transport equation {#sssec:prop_eq}

The general formalism that provides a powerful description of the transport of CRs in the Milky Way derives from the seminal textbook by Ginzburg and Syrovatskii [@GinzburgEtAl1964] (see also [@BerezinskiiEtAl1990; @Schlickeiser2002]), and relies on the following diffusion-advection equation for a CR species of index $\alpha$, here in the steady-state approximation and in energy space (rather than rigidity or momentum space):

\[eq:prop\] & - \_[**x**]{} { K(E)\_[**x**]{}\_- \_[c]{} \_} + { b\_[tot]{}(E)\_- \^2 K\_[pp]{} } + \_v\_ n\_[ism]{} \_+ \_\_\
& = q\_+\_ { \_v\_n\_[ism]{}+\_ }\_.

This equation describes the spatial and energy evolution of the differential interstellar CR density per unit energy $\psi_\alpha\equiv dn_\alpha/dE$, assuming a net primary injection rate of $q_\alpha$, and a secondary injection rate arising from inelastic processes converting heavier species of index $\beta$ into $\alpha$ species (with a production rate $ 
\sigma_{\beta\to \alpha}v_\beta n_{\rm ism}$ on the ISM density $n_{\rm ism}$, or a decay rate $\Gamma_{\beta\to\alpha}$). This source term is balanced by several other terms, among which the decay rate $\Gamma_{\alpha}$ (if relevant). The central piece of the propagation equation is the spatial diffusion coefficient $K$, that we discuss in more detail in . The other processes are mostly relevant at low rigidity, but may still affect the determination of higher-energy parameters: convection is featured by a velocity $\vec{V}_{\rm c}$, diffusive reacceleration is parameterized by the energy-dependent coefficient $K_{pp}$, and the inelastic destruction rate is given by $\sigma_{\alpha}v_\alpha n_{\rm ism}\, \psi_\alpha$, with the $\sigma$’s being energy-dependent nuclear cross sections; energy losses are characterized by the rate $b_{\rm tot}\equiv dE/dt$, which includes ionization and Coulomb processes as prescribed in [@MannheimEtAl1994; @StrongEtAl1998], as well as adiabatic losses induced by convection and reacceleration, see [@MaurinEtAl2002a; @PutzeEtAl2010].

Finally, we can switch from the interstellar (IS) CR flux predictions to the top-of-atmosphere (TOA) ones by means of the force-field approximation [@GleesonEtAl1968a; @Fisk1971], for which we only indicate the Fisk potential $\phi_{\rm F}$. The latter is constrained from Ref. [@GhelfiEtAl2016] for the AMS-02 data taking period.

### Geometry and cross sections {#sssec:geom}

We assume a 1D propagation model, as introduced in  [@BulanovEtAl1974; @Ptuskin1974; @Jones1979; @PtuskinEtAl1990; @JonesEtAl2001; @MaurinEtAl2001], where the magnetic halo confining the CRs is an infinite slab in the radial direction and of half-height $L$. Indeed, the radial boundary has only a minor quantitative impact on other transport parameters when the diffusion coefficient is taken as a scalar function (see  [@MaurinEtAl2002a; @PutzeEtAl2010]), and neglecting the radial dependence allows us to more efficiently probe the entire available parameter space without significant loss of generality. Therefore, we consider the vertical coordinate $z$ to be the only relevant spatial coordinate. See however for some considerations on the correspondence between 1D and 2D models.

The sources of CRs and the ISM gas which they scatter off are taken homogeneous in an infinitely thin disk at $z=0$, with an effective half-height $h=100$ pc. Energy losses are also considered to be localized in the disk, $b_{\rm tot}\propto 2\,h\,\delta(z)$. The ISM density is set to $n_{\rm ism}=2\,h\delta(z)\,n_0$, where $n_0=1$ cm$^{-3}$, corresponding to a surface density of $\Sigma_{\rm ism}= 2\,h\,n_0\simeq 6\times 10^{20}$ cm$^{-2}$ consistent with observations [@Ferriere2001]. We assign 0.9 and 0.1 of this budget to hydrogen and helium (in number), respectively. We do not indulge here in a more detailed discussion of the determination of these parameters and of their error from independent observations, since they are largely degenerate with the normalization of the diffusion coefficient (see, [*e.g.*]{}, [@MaurinEtAl2010]).

For the nuclear production and spallation cross sections, we use as reference the sets of tables from the [](https://galprop.stanford.edu/) package, while allowing for normalization, energy shift and low-energy slope changes according to the NSS method described in [@DeromeEtAl2019].

### Diffusion in real space and momentum space, and convection {#sssec:kdiff}

An important physical ingredient to all propagation models is the diffusion coefficient, which describes the scattering of CRs off magnetic turbulence. We assume that it can be taken as a scalar function, homogeneous and isotropic all over the magnetic slab. This is likely a good approximation in the context of B/C analyses because the CR flux is locally isotropic and the magnetic field configuration of the Milky Way exhibits relative fluctuations $\delta B/B\sim 1$ over all relevant spatial scales $\lambda$ (or wavenumbers $k=2\pi/\lambda$) [@CasseEtAl2002; @SunEtAl2008a; @JanssonEtAl2012a]. Since at the energies of interest, the CR flux is expected to be contributed to by many sources mostly located many “magnetic domains” away, an isotropic diffusion should provide at least a reasonable effective description of the data [@BerezinskiiEtAl1990]. For discussions on anisotropic models, see  Refs. [@DeMarcoEtAl2007; @EvoliEtAl2012a].

On the theory side, it is expected that the magnetic turbulence responsible for CR diffusion has different scaling behaviors in $k$-space, as a consequence of various possible phenomena. For instance, the turbulence power spectrum can be dominated by different sources depending on the dynamical range, with the resulting “two-zone” models known to provide good fits to the data [@Tomassetti2012; @GuoEtAl2018] (see also [@SeoEtAl1994; @JonesEtAl2001]). A very appealing scenario is proposed in Refs. [@BlasiEtAl2012a; @AloisioEtAl2013; @AloisioEtAl2015; @EvoliEtAl2018], relying on streaming instability [@Wentzel1969; @Holmes1975; @Skilling1975], where at rigidities beyond a few hundreds of GV, CRs diffuse on the turbulence injected on large scales by supernova bubbles. This turbulence cascades down until crossing the rigidity scale where the turbulence induced by CRs themselves takes over. This naturally generates a break like the one observed in the CR spectra and discussed just above. On the other hand, it is known that the CR spectra observed at low rigidity by the Voyager I spacecraft [@CummingsEtAl2016; @WebberEtAl2017] have a spectral slope rather different from the slope at intermediate rigidities. Due to the CR-wave coupling, any phenomenon with a low-rigidity characteristic scale, affecting either propagation or injection, may thus be at the origin of correlated changes in the CR spectra and the diffusion coefficient. For propagation, such a scale might arise due to the decrease of the CR pressure as CRs get closer and closer to the nonrelativisitic regime, and/or be related to some dissipation of the turbulence power spectrum [@YanEtAl2004; @PtuskinEtAl2005; @PtuskinEtAl2006a; @ShalchiEtAl2010; @EvoliEtAl2014; @XuEtAl2016]. In the following, while remaining agnostic on these specifics, and in contrast to previous B/C studies performed in the context of semi-analytical models ( [@MaurinEtAl2001; @MaurinEtAl2002a; @PutzeEtAl2010; @PutzeEtAl2011]), we want to capture the possibility that the diffusion coefficient departs from a single power law. This is justified by both theoretical arguments and observational evidence, as recalled above.

Starting from general considerations arising in the quasi-linear theory ( [@Jokipii1966; @BerezinskiiEtAl1990; @Schlickeiser2002; @Shalchi2009]), the diffusion coefficient is expected to be linked to the magnetic turbulence spectrum $|\delta B/B|$ through \[eq:k\_linear\_th\] K(E) = , where $v$ is the CR speed, $l_{\rm mfp}$ is the mean free path length, $r_{\rm L}\propto R/B$ is the Larmor radius defined from the rigidity $R=p/Z$, and $k_{\rm L}\propto 1/r_{\rm L}$ is the turbulence mode in resonance with the CR Larmor radius. Consequently, we propose a general form for the diffusion coefficient that can account for breaks in both the high-rigidity range and the low-rigidity range (hidden in the factor $r_{\rm L}/|\delta B/B|_{k_{\rm L}}^2$ above), which reads

\[eq:def\_K\] K(R) = \_ K\_[10]{} \_ \_ \_.

In the above equation, $\beta=v/c$ is the dimensionless CR speed, and $R_{\rm l/h}$ is the location of the low/high-rigidity break, while $R_{\rm 10}$ is an intermediate rigidity (here taken at 10 GV on purpose) such that $R_{\rm l}<R_{10}<R_{\rm h}$ ($R_{\rm l}\ll R_{\rm h}$). We then get the scaling $K(R)\propto \beta^\eta R^{\delta_{\rm l}}$ in the limit $R\ll R_{\rm l}$, and the scaling $K(R)\propto R^{\delta_{\rm h}}$ in the limit $R\gg R_{\rm h}$. Therefore, $\delta_{\rm l}$, $\delta$, and $\delta_{\rm h}$ simply describe the diffusion spectral indices in the low-, intermediate-, and high-rigidity regime, respectively. The parameter $s_{\rm l}$ ($s_{\rm h}$) characterizes how fast the spectral change proceeds around $R_{\rm l}$ ($R_{\rm h}$), and is inspired by the need to describe the very smooth hardening of the B/C data showing up at high rigidity. Indeed, we recall that the previous high-rigidity analysis performed in Ref. [@GenoliniEtAl2017] provided support to a softening of the diffusion coefficient to explain this feature, such that we can already anticipate that $\delta_{\rm h}<\delta$. The normalization of the diffusion coefficient $K_{10}$ (which carries the physical units) is another free parameter. Mind the difference with the convention used in most past analyses, where the normalization was instead $K_0$ and was defined at a rigidity $R_0=1\,{\rm GV}$. Note also that $K_{10}\simeq K(10\,{\rm GV})$, not a strict equality, because of the influence of the other terms.

We further introduce the spectral-change parameters \[eq:Delta\] \_[l]{} &=& - \_[l]{} ,\
\_[h]{} &=& - \_[h]{}, As already mentioned above, $\Delta_{\rm h}$ is expected to be positive. A positive $\Delta_{\rm l}$ is also expected from damping arguments and from the flattening of the primary CR spectra observed by Voyager I, as CRs may diffuse mostly on self-generated turbulence—see the discussion above. Notice that in the low-rigidity regime, additional non-relativistic processes might further be considered in an effective way by raising the velocity $\beta$ to the power $\eta$, an effective index which—it has been argued—might take negative values in some regimes [@PtuskinEtAl2005; @PtuskinEtAl2006a].

A comment on $\eta$ is in order: since the rigidity range of CR data analyzed in this article is always relativistic, sizable departures from $\eta=1$ (the natural value from quasi-linear theory, see ) and/or large values of $V_{\rm A}$ (which allows for energy redistribution) may be needed to affect appreciably CRs whose $v\simeq c$. In a certain sense, $\eta$ is thus not a very valuable [*effective*]{} parameter for the problem at hand. Nonetheless, we keep the $\eta$ parameter in the discussion for historical reasons, since in combination with strong reacceleration it used to be an important ingredient in past studies of B/C data, notably at low energies [@SeoEtAl1994; @StrongEtAl1998; @JonesEtAl2001; @MaurinEtAl2001; @DonatoEtAl2004]. Sufficiently large negative values of $\eta<-\delta_{\rm l}$ (or, similarly, of $\delta_{\rm l}<-1$ if $\eta=1$) can also imply superluminal diffusion [@DunkelEtAl2007; @AloisioEtAl2009] in the non-relativistic regime. In this sense, we caution the reader never to extrapolate a-critically the functional forms obtained here too far from the rigidity range over which the fits have been obtained.

Let us now be more specific about reacceleration. It turns out that spatial diffusion can rather generically be linked to diffusion in momentum space (aka reacceleration) in most (but not all) cases [@BerezinskiiEtAl1990]. We include diffusion in momentum space through an additional diffusion coefficient $K_{pp}$—see . We follow the reacceleration model proposed in Refs. [@OsborneEtAl1988; @SeoEtAl1994; @JonesEtAl2001], which is implemented in  such that $K_{pp}(R,\vec{x})=2\,h\,\delta(z)\,K_{pp}(R)$, and \[eq:Kpp\] K(R)K\_[pp]{}(R) = V\_[A]{}\^2 , where $V_{\rm A}$ is an effective Alfvénic speed characterizing the magnetic turbulence—$\delta$ is the diffusion spectral index in the intermediate inertial regime. Since it appears explicitly only as a normalization factor, we stick to this formula even when spatial diffusion exhibits several spectral regimes. The fact that reacceleration is [*effectively*]{} localized in the disk allows us to partly solve analytically, which significantly speeds up the numerical exploration of the parameter space [@Maurin2018]. While this “pinching” is a fair approximation for ionization and Coulomb processes, it is only a convenient approximation for adiabatic losses induced by convection and reacceleration. Hence, care should be taken when comparing inferred values of the parameters $V_{\rm c} \equiv | \vec{V}_{\rm c}|$ and $V_{\rm A}$ with theoretical expectations. Loosely speaking, one can expect to recover the phenomenology of a more extended reacceleration zone by a rescaling of $V_{\rm A}^2$ by a factor $h/z_{\rm A}$ [@MaurinEtAl2002], where $z_{\rm A}$ is the half-height over which reacceleration would spread in the magnetic slab [@JonesEtAl2001]. So, for $h/z_{\rm A}\simeq {\cal O}(h/L)$, our fitted value of $V_{\rm A}$ should be scaled by a factor $\sqrt{L/h}$ before any comparison against theoretical or observational constraints [@ThornburyEtAl2014; @DruryEtAl2017].

Finally, convection also arises quite naturally in the framework discussed above. We include convection in the standard way by means of the convection velocity \[eq:vc\] \_[c]{}(z) = V\_[c]{}\_z, where $z$ is the vertical coordinate and $\vec{e}_z$ the unit vector along the vertical axis crossing the magnetic slab of extension $[-L,L]$ along that axis.

Benchmark models {#ssec:benchmarks}
----------------

In the most general case, the free parameters featuring the propagation modeling that we have introduced above are the following: $L$ for the magnetic halo size; $K_{10}$, $\delta$, $\eta$, $R_{\rm l}$, $\delta_{\rm l}$ (equivalently $\Delta_{\rm l}$), $s_{\rm l}$, $R_{\rm h}$, $\Delta_{\rm h}$ (equivalently $\delta_{\rm h}$), and $s_{\rm h}$ for the diffusion coefficient; $V_{\rm A}$ for reacceleration; $V_{\rm c}$ for convection. This is a 12-parameter space, hence a huge configuration volume to explore.

Based on previous studies, we can further fix $L$ which is highly correlated with $K_{10}$, see discussion in . Unless specified otherwise, we will set $L$ to 10 kpc in the following. Moreover, as anticipated in , the determination of (an interval for) the three parameters describing the high-rigidity break benefits from fits including primary species, see . Finally, without loss of generality, we fix the smoothing low-rigidity break parameter $s_{\rm l}=0.05$, which amounts to consider a fast transition. This is however not critical to the fit. Hence, we are left with 7 free parameters.

From these 7 parameters, we design three different benchmark propagation models which may be related to quite different limiting regimes of the underlying microphysics. The first, most generic, model includes the whole setup introduced above: let us name it the  model. The second one is much simpler as it is free of convection and reacceleration, hence with much less free parameters, while providing fits to the data comparable to the previous one (see ); let us call it the  model. The third and last one includes both reacceleration and convection, but relates the possible change in the propagation at low rigidities to a change originating specifically in the non-relativistic regime ($\eta$), instead of a more generic low-rigidity break in the diffusion coefficient. This scenario provides a slightly worse fit to the data compared to the previous ones, at the expense of a large reacceleration $V_{\rm A}$. However, it allows us to connect the current analysis to the strong reacceleration models that were popular in the past; let us dub it the  model. Both the  and  models are actually particular cases of the  model, but put the emphasis on different physical processes at low rigidity. In the following, we provide the details of these three configurations.

### : the paradigmatic model {#sssec:big}

The  model includes a double-break diffusion coefficient, as well as convection and reacceleration. Its minimal version fixes the non-relativistic parameter $\eta=1$, while a non-minimal configuration may allow $\eta$ to vary. The latter case will actually help justify the former one independently from theoretical arguments. Therefore, the  model stands for the most general configuration describing the propagation equation, , and which allows us to probe the low-rigidity processes with the largest flexibility and complexity. This model has a total of 6 (7) parameters in the minimal (non-minimal) configuration, which are recalled in .

### : the minimal (double-break diffusion) model {#sssec:slim}

The  model is a subpart of , which discards convection and reacceleration as major players at low rigidity ($V_{\rm A}=V_{\rm c}=0$ km/s), but instead insists on relating low-rigidity features to changes in the magnetic turbulence properties. It also assumes a standard scaling in the non-relativistic regime, with $\eta=1$. This model, though very minimal, will be shown to provide an excellent fit to the data. It has 4 free parameters which are summarized in . Note that an important advantage of this model is that it comes with a fully analytical solution to the transport equation. This is particularly attractive in the context of dark matter predictions [@MaurinEtAl2002; @LavalleEtAl2012].

### : the ‘old-fashion’ strong reacceleration model {#sssec:quaint}

Our last benchmark model is the  model, which is also a subpart of the  model, and which aims at describing the low-rigidity features mostly in terms of reacceleration and convection. This model is actually the direct descendant of the [*min*]{}-[*med*]{}-[*max*]{} models [@MaurinEtAl2001; @DonatoEtAl2004] as it relies on almost the same configuration space, except for the high-rigidity break in the diffusion coefficient (which was not observed at the time of its ancestors and will anyway be treated as a nuisance parameter in the statistical analysis). Large $V_{\rm A}$, in combination with a non-trivial value of $\eta\lesssim 0$ is needed to provide decent fits to the data. A large $V_{\rm A}$ in turn couples low-rigidity and high-rigidity features, maximizing parameter correlations. The model has 5 free parameters, made explicit in . In practice, the diffusion coefficient associated with the  model is that of without the low-rigidity term.

   Free parameters / Models   [**BIG**]{}   [**SLIM**]{}   [**QUAINT**]{}
  -------------------------- ------------- -------------- ----------------
           $K_{10}$                                       
           $\delta$                                       
            $\eta$               1 or            1        
       $\delta_{\rm l}$                                         N/A
         $s_{\rm l}$             0.05           0.05            N/A
         $R_{\rm l}$                                            N/A
         $V_{\rm A}$                            N/A       
         $V_{\rm c}$                            N/A       

  : Free parameters of the three benchmark models , , and . The first block of parameters is associated with the diffusion coefficient in the intermediate regime, and is common to all models. The second block is related to a potential low-rigidity break in the diffusion coefficient, or to purely non-relativistic effects. The last block is related to reacceleration and convection.[]{data-label="tab:free_params"}

Fitting strategy {#sec:fit_bc}
================

In this section, we explain the fitting strategy used to extract the benchmark propagation parameters for the models presented above (, , and ). Fits are performed with the  package [@JamesEtAl1975] interfaced with the  code [@Maurin2018], and in particular, asymmetric error bars on the parameters rely on the  algorithm. For more technical details and subtleties on the setup and the analysis, we refer the reader to Ref. [@DeromeEtAl2019].

### Modeling uncertainties

For each run, the fluxes of the elements from Beryllium (Be) to Silicon (Si) are computed assuming that $^{10}$B, $^{11}$B (and $^{10}$Be, decaying into $^{10}$B) are pure secondary species and that all the heavier elements contain a secondary and a primary component. We assume the primary injection to follow a universal power law in rigidity with index $\alpha$. The secondary component is computed by a full spallation network using the  cross-section parameterization (see appendices of [@GenoliniEtAl2018]). It has been shown in [@DeromeEtAl2019] that this parameterization provides the best agreement with the data, and that uncertainties on spallation cross sections are satisfactorily taken into account using only the $^{12}{\rm C}+{\rm H}\rightarrow {}^{11}{\rm B}$ production cross section as nuisance parameter with the “normalization, slope and shape” (NSS) strategy. For each run, the [*initial*]{} default procedure is to fix the normalization of the primary components of all elements to the 10.6 GeV/nuc data point of HEAO-3 [@Maurin2018], except for the CNO elements which affect more directly the B/C ratio: The latter ones are normalized to the C, N, O data of AMS-02 at a rigidity of 50 GV. The power-law index $\alpha$ is first set to 2.3, and fixed later via the iterative procedure explained below. The solar modulation of CRs is described in the force-field approximation, for which the Fisk potential $\phi_{\rm F}$ is averaged over the AMS-02 B/C data taking period. Based on [@GhelfiEtAl2016], we set $\phi_{\rm F}$ as a nuisance parameter of mean value 730 GV and dispersion $\sigma_{\phi_{\rm F}}=$100 MV.

### Data errors

The AMS-02 collaboration does not provide users with the covariance error matrix of the data. In this case, it is common practice to estimate the *total* errors by summing systematics and statistics in quadrature. This procedure is however inappropriate when systematics dominate, i.e. below $\sim$ 100 GV for the AMS-02 data, and for which correlations in energy are expected to be important. A major novelty of the present analysis is to perform fits taking into account these correlations with a parametric form of the covariance matrix. The matrix was built thanks to the information provided in the Supplemental Material of the AMS-02 B/C analysis [@AguilarEtAl2018]. In particular, the different systematics, which are associated with different physics processes in the detector, have different correlation lengths, and the covariance matrix built reflects this complexity. For more precision we refer the reader to [@DeromeEtAl2019].

### C and O primaries

It has been noted already that recent data show an indication for a high-rigidity break in the diffusion coefficient. However, these data are at present still far for providing us with the precise characteristics of this feature. In fact, in order to start gaining statistical confidence in the very existence of this break, the typical strategy until now has been to combine the B/C data with independent indications for the break. For instance, in [@GenoliniEtAl2017] we used the AMS-02 $p$ and He data to that purpose. As anticipated in , we do not any longer focus our B/C analysis on the high-rigidity regime. Instead, we want to provide reference values for the parameters controlling the low- and intermediate-rigidity regimes. Consequently, it is a natural choice to use the high-rigidity break parameters as [*nuisance*]{} parameters. However, in order to establish the plausible range over which to vary them, it is recommended to resort to complementary and “independent” input. To minimize possible biases due for instance to possibly different origins of the different species, we choose to limit ourselves to the C and O fluxes because:

They are by far the main progenitors of the B and C fluxes entering the B/C ratio.

Fitting them allows us to determine a plausible value of the common spectral index of nuclei $\alpha$ as well as to check their [*consistency*]{} with the parameters obtained with the B/C analysis. Indeed, although we focus here on the B/C observable to determine the propagation parameters, we still want to make sure that our results are consistent with the observed primary fluxes. In , we make the important sanity check of neglecting this external input, relying solely on the B/C data to determine both the high-rigidity spectral break (and the other propagation parameters). We show that the obtained results are perfectly consistent with the “factorizing” procedure sketched above, at the obvious price of a worse determination of the propagation parameters. This is indeed not surprising, since the B flux, which dominates the B/C statistical error at high rigidity, is more than one order of magnitude scarcer than the C and O ones.

### The fitting procedure

The technical implementation of the fits proceeds by iteration. After fixing the (low- and intermediate-rigidity) propagation parameters with a first fit of the B/C ratio (as described above, i.e. with $\phi_{\rm F}$ and the $^{12}{\rm C}+{\rm H}\rightarrow {}^{11}\rm B$ production cross section taken as nuisance parameters), we perform a combined fit of the AMS-02 C and O fluxes keeping the following parameters as free parameters: source-term normalizations, power-law dependence in rigidity $\alpha$, and break parameters ($R_{\rm h}$, $s_{\rm h}$, and $\Delta_{\rm h}$). We then use the best-fit values of the break parameters and associated covariance matrix as nuisance parameters in a new B/C fit, keeping also $\alpha$ fixed to its best-fit value. In practice, only a couple of iterations are needed to get the parameters compatible between two consecutive iterations. The results discussed below are the outcome of this procedure.

![Combined best fits of C (top) and O (bottom) fluxes in the three benchmark models , , and  defined by B/C fits (see ). A panel below each plot reports the $Z$-score, corresponding to the residuals normalized by the [*total*]{} errors $\sigma_{\rm tot}$.[]{data-label="fig:best_prim_modelB"}](local_fluxes_AND_score_C_R_nuis_Galp_final.pdf "fig:"){width="\columnwidth"} ![Combined best fits of C (top) and O (bottom) fluxes in the three benchmark models , , and  defined by B/C fits (see ). A panel below each plot reports the $Z$-score, corresponding to the residuals normalized by the [*total*]{} errors $\sigma_{\rm tot}$.[]{data-label="fig:best_prim_modelB"}](local_fluxes_AND_score_O_R_nuis_Galp_final.pdf "fig:"){width="\columnwidth"}

For the fits of the C and O fluxes, a simpler yet sufficient approximation is to assume uncorrelated *total* errors $\sigma_{\rm tot}$, i.e. statistical and systematic errors summed in quadrature: on the one hand, only statistical uncertainties dominate around the high-energy break position, so that this is a reasonable approximation. On the other hand, this fit only enters the B/C analysis via the treatment of the high-rigidity parameters as nuisance.

In we report the fits of the source and high-rigidity break parameters to the C (top panel) and O (bottom panel) fluxes for our three benchmarks, the , , and  models. It is clear even by visual inspection that the fits with a simple, common power-law index $\alpha$ are excellent: The fits fall within one $\sigma_{\rm tot}$ and never beyond two $\sigma_{\rm tot}$’s from all intermediate and high-rigidity points, showing that our consistency check is successful. Some minor discrepancy at low rigidity is noticeable, but not worrisome for our purposes. In fact, should one aim at describing C and O primary fluxes in detail down to low rigidities, a more accurate fitting procedure treating cross-section parameters as nuisance and accounting for bin-to-bin correlations of the systematic errors (as done for B/C) would certainly reduce these minor disagreements. This is beyond our goals here, but will be of interest for future more global analyses.

Results {#sec:results}
=======

In this section, we illustrate our results for the (low- and intermediate-rigidity) propagation parameters and discuss their implications. Initialization files used for the analysis, along with the resulting best-fit values and covariance matrix of best-fit parameters will be provided with the forthcoming new release .

Best-fit values and $1\sigma$ uncertainties
-------------------------------------------

The best-fit values and errors on the three model parameters (, , and ) are reported in . In the first block, we report the diffusion parameters $\delta$ and $K_{10}$ common to all models, which control the *intermediate-rigidity* regime. We then report the *low-rigidity* parameters, which are different (both in nature and number) between , on one side, and,  and , on the other. The high-rigidity break parameters, fixed following the nuisance procedure, are reported at the bottom of the Table. The range over which we scan for them will be discussed in , since their determination is affected by the inclusion of external data (in our case, C and O absolute fluxes).

In all these fits, nuisance parameters vary within reasonable pre-assigned intervals. The solar modulation parameter $\phi_{\rm F}$ attains a value of 731, 734 and 725 MV in the best-fit model , , and , respectively. Concerning the nuisance of the spallation cross section $\rm ^{12}C+H \rightarrow {}^{11}B $, its best *normalization* is found to be 12%, 13% and 11% above the reference GP17 value in the best-fit model , , and , respectively. [The preferred *slope* encoding the low energy shape is of 0.12, 0 and 0.16, for the same models. The induced spectral distorsions in  and correspond to a slight decrease of the cross section at low energy.]{}

[l c c c]{} Parameters & & &\
$\chi^2/{\rm dof}$ & $61.7/61\!=\!1.01$ & $61.8/63\!=\!0.98$ & $62.1/62\!=\!1.00$\
\
[$K_{10}$]{}\[kpc$^2$Myr$^{-1}$\] & $0.30_{-0.04}^{+0.03}$ & $0.28_{-0.02}^{+0.02}$ & $0.33^{+0.03}_{-0.06} $\
$\delta$ & $0.48_{-0.03}^{+0.04}$ & $0.51_{-0.02}^{+0.02}$ & $0.45_{-0.02}^{+0.05} $\
\
$V_{\rm c}$ \[kms$^{-1}$\] & $0^{+7.4}$ & N/A & $0.0^{+8}$\
$V_{\rm A}$ \[kms$^{-1}$\] & $67^{+24}_{-67}$ & N/A & $101_{-15}^{+14}$\
$\eta$ & 1 (fixed) & 1 (fixed) & $-0.09_{-0.57}^{+0.35}$\
$\delta_{\rm l}$ & $-0.69^{+0.61}_{-1.26}$ & $-0.87^{+0.33}_{-0.31}$ & N/A\
$R_{\rm l} $ \[GV\] & $3.4_{-0.9}^{+1.1}$ & $4.4_{-0.2}^{+0.2}$ & N/A\
\
\
$\Delta_{\rm h}$ & $0.18$ & $0.19$ & $0.17$\
$R_{\rm h} $ \[GV\] & $247$ & $237$ & $270 $\
$s_{\rm h}$ & $0.04$ & $0.04$ & $ 0.04$\

![Best fit B/C curve for models , , and . Results for the best fit parameter values are given in . The bottom panel shows the $Z$-score.[]{data-label="fig:best_fit_model"}](local_fluxes_AND_score_BC_R_mod_all_nuis_Galp_final.pdf){width="\columnwidth"}

Our best-fit curves are reported in for the three models. Note that all models lead to analogous curves and fit quality, only differing in the fine features of the spectral shape at low rigidity. The inset displays the $Z$-score, i.e. the residuals normalized to the *total* errors $\sigma_{\rm tot}$. Note that this has only a qualitative purpose, since technically the $\chi^2$ is computed accounting for correlations in the systematics of B/C data, a major novelty of this analysis. The similar fit quality of the  and  models indicates that the additional free parameters present in the former are actually unnecessary to describe the data: If the fit allows for a low-rigidity break, there is but a minor and currently unnecessary role played by $V_{\rm c}$ and $V_{\rm A}$. We note a tiny and statistically insignificant preference for model  (and a fortiori ) with respect to , which is only worth noticing since  has one free parameter more than in . In fact, we stress that if we had fixed $\eta=1$ in the  model, its fit quality would have degraded, and it would have been rejected at $>2\,\sigma$ with respect to the  and  models. [Finally, we note that, compare to , the benchmark  and  have respectively a weaker and no break at low rigidity, althought the latter is partly mimicked by the spectral distorsions of the cross section in nuisance. This tends to provide additional support to the possible presence of a low-rigidity break in the diffusion coefficient.]{}

![[]{data-label="fig:diffusion"}](K_diffusion.pdf){width="\columnwidth"}

Also, it is important to notice that the parameters common to the three models are found with values compatible within $\sim 1\sigma$. This suggests that the diffusive properties at intermediate rigidities are constrained rather robustly by the data [(see  for an illustration of this)]{}, independently of the specific scenario within which the low-rigidity behavior is interpreted and fitted. This conclusion is rather encouraging when one considers interpretations of the high-rigidity spectral break.

The  model comes out with a few apparently surprising features: at face value, the best fit for $V_{\rm A}$ is rather large, perhaps even more so in the light of the value found for $V_{\rm c}$ which is compatible with zero. A too large value for $V_{\rm A}$ would lead to the surprising conclusion that the power in diffusive reacceleration is comparable to the total CR luminosity. This would for instance imply that the bulk of CRs energy rather originates from ISM turbulence than from shocks in SNRs, as customarily assumed (see e.g.[@DruryEtAl2017]). This apparent conundrum is alleviated once accounting for the effective nature of our parameterization of the reacceleration term, which for technical reasons is artificially pinched to the thin disk, rather than being present in the whole propagation halo. The actual power in the turbulence in the whole Galaxy is thus reduced by the ratio $h/z_{\rm A}\simeq {\cal O}(h/L)$ (see ), hence resulting roughly consistent with expectations, and also more in line with the allowed range for $V_{\rm c}$. Another perhaps surprising outcome is the value of $\eta$, whose best fit is [*negative*]{}, and anyway rather away from typically considered values $\simeq 1$. This conclusion is also qualitatively valid in models  and : there is a slight preference for the diffusion coefficient below about 4.5 GV to [*increase with lower rigidity*]{}. Note that, within the allowed range for $\delta_l$, extrapolation in the non-relativistic regime could lead to nonphysical results, as soon as $\delta_l+1=\delta-\Delta_{\rm l}+1<0$. Be that as it may, neither this caveat nor the previous one for $\eta$ in the  model should be overstated, since they arise at best at the $1\sigma$ level. All model fits are consistent with a perhaps more physically acceptable flat behavior, or a rigidity-independent diffusion coefficient at low-$R$. Furthermore, it is worth noting how the  model is closer to the relatively unproblematic regime $\delta_{\rm l}>-1$ than the  one. The introduction of some reacceleration and convection (both physically expected) tends to yield more reasonable values for the low-rigidity slope.

Possible interpretation and microphysics
----------------------------------------

In any case, the most obvious interpretation of these results is that there are less and less waves onto which CRs can scatter at low rigidity. One possible reason is that turbulence dissipation effects lead to a parallel diffusion coefficient which decreases with increasing rigidity, with turning point at $\sim 3\;$GV for plausible choices for the parameters [@ShalchiEtAl2010]. Another possibility arises in models where the CRs scatter onto self-generated turbulence below some rigidity (see e.g. [@BlasiEtAl2012a]). The energy density (and the pressure) carried by CRs peaks at the few GV scale; above this rigidity, the induced diffusion coefficient increases with rigidity as customarily assumed, because of the relatively steep CR power-law spectrum. Below this rigidity, however, the lower the rigidity (or Larmor radius), the smaller the turbulence with respect to extrapolations, simply because there are less and less CRs that can generate it by streaming instability due to their spectral inflection. The order of magnitude of the break in the low-energy CR spectrum seems to be in the right ballpark, but these qualitative arguments deserve a more detailed investigation, which we postpone to future work.

Robustness of low-, intermediate-, and high-rigidity parameters
---------------------------------------------------------------

A very encouraging finding is that, within uncertainties, the diffusive properties at intermediate rigidities do not depend on the specific scenario considered at low-rigidity. The value found for $\delta$ appears closer to a Kraichnan turbulence spectrum ($\delta\simeq 0.5$) than to a Kolmogorov one ($\delta\simeq 1/3$), although this conclusion should not be overstated since the model involves an [*effective*]{} isotropic diffusion coefficient. An indirect implication of this robustness is to increase the credibility in any deviation found at high rigidity, of course.

Concerning the low-rigidity regime, however, there are several important caveats, which suggest some prudence to avoid over-interpreting the values found. First of all, while there is a clear indication for a different regime of propagation at low rigidity, the “hardest” parameters to interpret ($\eta$ and $\delta_{\rm l}$) are actually heavily influenced by the one or two lowest-rigidity points. This is illustrated in more detail in , where one can compare the behavior of $R_{\rm l}$ vs. $R_{\rm min}$ with respect to $\delta_{\rm l}$ vs. $R_{\rm min}$, $R_{\rm min}$ being the rigidity above which the fit is performed. There is simply not enough of a baseline at low rigidity in the AMS-02 data to unambiguously measure the slope in this range. Another point to keep in mind is that the low-rigidity range is quite influenced by the uncertainties in the nuclear cross sections and the treatment of solar modulation. Indeed, including the nuisance parameters for the production cross section increases the $1\sigma$ uncertainties on $\eta$ () by 50%, and on $R_l$ and $\delta_l$ () by 90%. In our fits, including solar modulation is a second order effect, since it increases the low-energy parameters uncertainties by order 5%.

The only model-independent conclusion that we can safely make on the low-rigidity range is that multiple models can account for the observations, with rather different physical interpretations possible. So, statements such as “the reacceleration/convection velocity determined from the B/C data is …” should be taken with a grain of salt, since they appear [*very*]{} model dependent, if compared, for instance, with the determination of $\delta$. The fitted values should [*only*]{} be used as references in the same model used to fit them, and extrapolations at lower rigidities (below the range covered by the data) are not guaranteed to be physical.

Summary, conclusions, and perspectives {#sec:conclusion}
======================================

This article has set the stage for the propagation scenarios that we want to test, challenge and refine with further AMS-02 data, defining benchmark models and ranges of parameters. We have validated the first step of this program with a statistically more sound analysis of the AMS-02 B/C data, going beyond state-of-the-art in the modern literature, and checking different theoretical frameworks differing in the treatment of transport at low rigidities with a major (model ) or a negligible (model ) role played by reacceleration. Both models are limiting cases of a more general model (). We have made sure that issues like numerical stability, the effects of cross sections uncertainties, the bin-to-bin correlation of systematic errors are handled sufficiently well not to bias significantly the conclusions.

For the time being, either model can describe with comparable performances the low-rigidity regime, with a statistically insignificant preference for model . The parameters describing intermediate rigidities are consistently determined in either case. At low rigidity, degeneracies with nuisance parameters impact both the best fit and uncertainties, in particular the ones controlling the energy shape of cross sections and solar modulation. This means that qualitatively different models offer almost equally good description of the data, so that inferring the physics of the propagation at low rigidity is challenging, and we must content ourselves with one or another “effective” description. This lesson on the shaky discrimination power among models with mild differences at low rigidities is likely to apply more generally, even to alternative models not tested here, because it partially relies on the effect of the nuisance parameters. Obviously, any discrimination between the two sets of models must be based on complementary data or arguments, such as the (astro)physical plausibility of the parameters found, an issue which we also briefly discussed. However, finding a break in the diffusion coefficient at low rigidity should not come so much as a surprise, since this feature, possibly related to some damping in the turbulence spectrum and the subsequent increase of the CR mean free path, is expected from theoretical grounds [@YanEtAl2004; @PtuskinEtAl2006a; @ShalchiEtAl2010; @EvoliEtAl2014]. In this respect, a careful study of low energy data complimentary to AMS-02 ones (e.g. from ACE-CRIS [@LaveEtAl2013] and Voyager I [@CummingsEtAl2016; @WebberEtAl2017]), together with a more realistic account of the systematic error correlations (based on further information provided by the experimental collaborations), could certainly help in drawing more robust conclusions on the properties and the nature of this break.

Besides extracting reference propagation parameters and uncertainty ranges from B/C data, which are intended for references for further studies, we have also performed a first test of the consistency of the obtained results with simple CR source spectra (power laws). We also confirmed and strengthened our conclusions in [@GenoliniEtAl2017], that the high-rigidity data can be consistently interpreted as a consequence of a break in the diffusive coefficient, in agreement with AMS-02 high-energy primaries spectra. Indeed, we show that this preference does persist in a generalized analysis extending to the whole rigidity range, and for alternative propagation setups, notably with/without reacceleration.

Following this study, the most pressing issue is of course to test the reference models provided here against other secondary data (e.g., Li, Be, pbar, positrons). In particular, our forthcoming publication will focus on the antiproton channel [@BoudaudEtAl2019]. There has been a recent interest in the possibility that these data hide a signal of dark matter annihilation, see e.g. [@CuocoEtAl2019; @CholisEtAl2019], and it is interesting and important to re-examine those claims within our analysis framework. Finally, it is also known that putative dark matter signals are sensitive to the diffusive halo size, hence an important and motivated follow-up project analysis will involve other secondaries, including isotopes such as the radioactive species (e.g. ${}^{10}$Be). To that purpose, we provide the reader in (1D models) and in (2D models) with the scaling relations that allow to extrapolate our benchmark models (derived assuming $L=10$ kpc) to a range of $L$ between 4 and 18 kpc.

[*Note added:*]{} As we were completing this study, we became aware of [@VittinoEtAl2019], where the authors find support for the presence of multiple breaks in the diffusion coefficient, based on CR electron and positron data. This result is complementary to (and consistent with) ours, while interestingly based on independent datasets concerning other CR species.

This work has been supported by the “Investissements d’avenir, Labex ENIGMASS", by the French ANR, Project DMAstro-LHC, ANR-12-BS05-0006, and by Univ. de Savoie, AAP “DISE”. It has also benefited from the support of the ANR project GaDaMa (ANR-18-CE31-0006), the OCEVU Labex (ANR-11-LABX-0060), the CNRS IN2P3-Theory/INSU-PNHE-PNCG project “Galactic Dark Matter”, and European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreements N$^\circ$ 690575 and No 674896. The work of Y.G. is supported by the IISN, the FNRS-FRS and a ULB ARC. The work of M.B. is supported by the European Research Council ([ERC]{}) under the EU Seventh Framework Program (FP7/2007-2013) / [ERC]{} Starting Grant (Agreement No 278234 – [“NewDark”]{} project).

On the high-rigidity break from C, O and the fitting procedure {#app:consistency}
==============================================================

The fitting procedure described in makes use of the C and O fluxes:

1.  As a sanity check for the actual diffusion parameters inferred (see ).

2.  To determine $\alpha$, the common spectral index for all nuclei, although its value is irrelevant for the B/C calculation (and for the transport parameter determination).

3.  Above all, for cornering a plausible window for the nuisance of the the high-energy break parameters.

Below, we provide some consistency checks as well as some comments on these ancillary results.\

Consistency check
-----------------

In , we report the results of the reference fits of C, O fluxes corresponding to the B/C model fits discussed in .

       Parameters                 [BIG]{}                        [SLIM]{}                           [QUAINT]{}
  --------------------- ---------------------------- --------------------------------- ------------------------------------
   $\chi^2/{\rm dof}$         $75.7/129= 0.59$                $73.2/129=0.57$                    $80.3/129= 0.62$
        $\alpha$                    2.35                           2.33                                2.36
    $\Delta_{\rm h}$       $0.18_{-0.05}^{+0.13}$         $0.18^{+0.11}_{-0.04}$              $0.18_{-0.01}^{+0.18}$
   $R_{\rm h} $ \[GV\]       $244_{-52}^{+198}$             $236_{-51}^{+152}$             $282_{-89}^{\textit{+349}} $
       $s_{\rm h}$       $0.04_{\it -0.04}^{+0.11}$   $0.03_{\textit{-0.03}}^{+0.09}$   ${ 0.04}^{+0.15}_{\textit{-0.04}}$

  : The $\chi^2/{\rm dof}$, the best-fit value for $\alpha$ as well as error range (used as nuisance parameters and ranges in the B/C analysis) for the *high-rigidity* parameters, coming from the combined fit to absolute C, O fluxes in the iterative procedure described in . Values in italics means that the fit reached the border of the interval. []{data-label="keyresH"}

All propagation models inferred from B/C appear to provide excellent fits to the C, O fluxes as well. The $\chi^2$ cannot be used at face value as a quantitative estimator of the quality of the fit, since [*total*]{} errors have been used in the C, O fits: Hence, we likely underestimate the contribution to the $\chi^2$, notably those of the intermediate- and low-rigidity data mostly influenced by systematic errors and their correlations (see also the companion paper [@DeromeEtAl2019]). Nonetheless, a relative preference seems to emerge for the [BIG]{} and [SLIM]{} models, compared to the [QUAINT]{} model, which is interesting as the same trend is also present from the more rigorous B/C analysis.

Concerning $\alpha$, the values found are intriguingly similar to the ones found in the fit of the He flux, which is performed in [@BoudaudEtAl2019], another reassuring consistency test of our procedure. We are thus consistent with the current universality of the spectra of nuclei (while the proton flux seems to be somewhat steeper). This is an interesting observable to keep an eye on in the future, of course. Note that the fit yields a nominal error on the parameter $\alpha$ at the sub-percent level, since $\delta$ is kept fixed in the iteration. Realistic uncertainties on $\alpha$ are however comparable to the ones of $\delta$ reported in .

In all cases, the indication for a high-rigidity break $\Delta_{\rm h}$ is rather significant ($\gtrsim 4\,\sigma$), again consistently with AMS-02 results, but here referring to the underlying diffusive coefficient (i.e. a break in the model space, not in the flux spectral index). Also, the values found are consistent within the errors with those found from $p$, He analyses (e.g. [@GenoliniEtAl2017]), although a bit higher, i.e. indicating a slightly more pronounced break. It will be interesting to follow-up on this in the light of further analyses of both light and intermediate/heavy nuclei, to see if the situation will relax towards a more common value or point to some discrepant hardening.\
Let us briefly develop further on the significance of the high-rigidity break in the light of B/C data only.

          Parameters                       [BIG]{}                           [SLIM]{}                       [QUAINT]{}
  -------------------------- ------------------------------------ ------------------------------ ---------------------------------
           $\delta$                 $0.55_{-0.04}^{+0.20}$           $0.55 ^{+0.09}_{-0.03}$            ${\it 0.9}_{-0.23}$
   $K_{10}$ \[kpc$^2$/Myr\]         $0.26_{-0.2}^{+0.05}$            $0.26 _{-0.01}^{+0.07}$          $0.10^{+0.07}_{-0.01}$
     $V_{\rm A}$ \[km/s\]                 $0^{+64}$                             NA                        $71_{-7}^{+20}$
     $V_{\rm c}$ \[km/s\]                 $0^{+16}$                             NA                        $19_{-5}^{+3}$
           $\eta_t$                       1 (fixed)                         1 (fixed)                 $-0.30_{-0.75}^{+0.54}$
       $\delta_{\rm l}$            $-0.84^{+0.32}_{-0.36}$           $-0.87^{+0.35}_{-0.33}$                    NA
     $R_{\rm l} $ \[GV\]             $4.4_{-2.1}^{+0.46}$              $4.4^{+0.2}_{-0.2}$                      NA
       $\Delta_{\rm h}$             $0.27_{-0.12}^{+0.22}$         $0.27_{{\it -0.12}}^{+0.21}$       $0.56_{-0.24}^{+0.09}$
     $R_{\rm h} $ \[GV\]           $158_{{\it -58}}^{+235}$          $159^{+240}_{{\it -59}}$           $ {\it 100}^{+96}$
         $s_{\rm h}$          $0.10_{{\it -0.10}}^{{\it +0.20}}$   $0.11_{-0.1}^{{\it +0.19}}$    $0.26_{{\it-0.26}}^{\it +0.04}$
      $\chi^2/{\rm dof}$                $58.6/58=1.01$                    $58.7/60=0.98$                  $59.7/59=1.01$

  : [****]{}Best fit parameters for models [BIG]{}, [SLIM]{}, and [QUAINT]{}, if fitting the high-energy break of the diffusion coefficient as well on the B/C data only. Values in italics means that the fit reached the border of the interval.

\[tab:results\_break\_free\]

If we were to use solely B/C data to fit [*also*]{} the high-rigidity parameters (i.e. [*without*]{} relying on the C, O flux data), we would obtain the results listed in Table IV. The low and intermediate rigidity propagation parameters are consistent with our reference one (see ), with larger error bars, as expected since we now are determining more parameters from a more restricted set of data. Similar considerations apply to high-rigidity parameters, compare with . The largest departures are seen in the [QUAINT]{} model, where one suffers from a partial degeneracy of the (large) $V_{\rm A}$ parameter with the others, including $\delta$. Also, in this case parameters tend to drift towards the borders of the “plausible” interval fixed beforehand, which puts into question how physically meaningful this model results really are. Still, in all cases there is an evidence for a high-rigidity break (at $\gtrsim 2\,\sigma$ level, naively speaking) from the B/C alone, which [*a posteriori*]{} is a justification for our choice of the parameterization of the diffusion coefficient, .

Break vs no-break
-----------------

In Tab. V we report the best fit propagation parameters without high-rigidity break in the diffusion coefficient. Note how the values of $\delta$ would be biased (at the $1\div 2\,\sigma$ level), resulting in a harder diffusion coefficient. In the same spirit as [@GenoliniEtAl2017], we compute the $\Delta\chi^2$ with respect to our results in  (break parameters fit to C and O fluxes). In the [QUAINT]{} model, as intuitively expected, the presence of a large $V_{\rm A}$ can partially mimic the break, but not completely, and the “no break” case is still disfavored (at $\sim 2\,\sigma$ level). For the [BIG]{} and [SLIM]{} models, which are refinements of the purely diffusive intermediate/high-rigidity model considered in [@GenoliniEtAl2017], we find $\Delta\chi^2>10$, confirming (and thus reinforcing the robustness of) the results presented in [@GenoliniEtAl2017].

          Parameters                  [BIG]{}                  [SLIM]{}                 [QUAINT]{}
  -------------------------- ------------------------- ------------------------- ------------------------
           $\delta$           $0.48^{+0.02}_{-0.02}$    $0.48^{+0.02}_{-0.02}$    $0.42^{+0.03}_{-0.02}$
   $K_{10}$ \[kpc$^2$/Myr\]   $0.29^{+0.02}_{-0.02}$    $0.29^{+0.02}_{-0.02}$    $0.36^{+0.02}_{-0.04}$
     $V_{\rm A}$ \[km/s\]           $0^{+115}$                    NA                $113^{+7}_{-15} $
     $V_{\rm c}$ \[km/s\]            $0^{+12}$                    NA                   $0 ^{+4.1}$
           $\eta_t$                  1 (fixed)                 1 (fixed)           $0.6^{+0.3}_{-0.5} $
       $\delta_{\rm l}$       $-0.88_{-0.30}^{+0.31}$   $-0.88_{-0.30}^{+0.32}$             NA
     $R_{\rm l} $ \[GV\]       $4.4^{+0.23}_{-2.4}$      $4.4^{+0.24}_{-0.21}$              NA
      $\chi^2/{\rm dof}$          $72.8/61=1.19$            $72.8/63=1.16$            $67.1/62=1.08$
        $\Delta\chi^2$                 11.1                      11.0                      5.0

  : Best fit parameters for models [BIG]{}, [SLIM]{}, and [QUAINT]{}, with no high-rigidity break in the diffusion coefficient.

\[tab:results\_without\_break\]

Fit parameters dependence upon low-rigidity cutoff {#app:low_rig}
==================================================

In , we present the evolution of the best-fit parameters as a function of a low-rigidity cut $R_{\rm min}$ above which the fit is performed, for the three models considered. We note that the value of the parameters $\delta$ and $K_{10}$ remains essentially unchanged whatever $R_{\rm min}$. In contrast, $\delta_{\rm l}$, $R_{\rm l}$ (for  and ), $V_a$, and $\eta$ (for ) depend crucially on the first data points, notably those below $\sim$4 GV. This explains their denomination of *low-rigidity* parameters. For the  model, note how the error on $\delta_{\rm l}$ crucially depends on the first couple of AMS-02 points, and the evidence for a change of slope (a determination of $R_{\rm l}$) is stronger than the actual value of the slope at low rigidity. Finally, it is worth commenting on $V_{\rm A}$: this parameter is (anti)correlating with low-energy ones (in particular $\eta$ for ) and, to a minor extent, also with $\delta$. This is not very surprising since large values of $V_{\rm A}$ imply “cross-talk” among energy bins. \[app:primaries\]

![image](scan_Rmin_model_Model_AB.pdf){width="0.64\columnwidth"} ![image](scan_Rmin_model_Model_B.pdf){width="0.64\columnwidth"} ![image](scan_Rmin_model_Model_A.pdf){width="0.64\columnwidth"}

Scaling of propagation parameters with $L$ in 1D model {#app:Ldiff}
======================================================

The benchmark parameters for , , and  where derived assuming the Galactic magnetic halo shapes as a 1D-slab of half-thickness $L=10$ kpc. By fitting the B/C ratio in these models, it is well known that the normalization of the diffusion coefficient $K_{10}$ and the halo thickness $L$ are degenerated so that the ratio $K_{10}/L$ is constant. We have checked that this was still the case given the higher sensitivity of AMS-02 data, and found the following scaling relations for values of $L$ within \[4,18\] kpc:

\[eq:sc\_relation\] $$\begin{aligned}
&&\text{\BIG:  }\quad\frac{K_{10}}{L}=0.030^{+0.003}_{-0.004}\;\rm kpc/Myr\;.\\
&&\text{\SLIM:  }\quad\frac{K_{10}}{L}=0.028^{+0.002}_{-0.002}\;\rm kpc/Myr\;,\\
&&\text{\QUAINT:  }\quad\frac{K_{10}}{L}=0.033^{+0.003}_{-0.006}\;\rm kpc/Myr\;.\end{aligned}$$

Dictionary to use 1D propagation parameters in 2D models. {#app:dict}
=========================================================

A 1D-slab geometry for the magnetic halo does not allow one to account for CRs escaping radially from the Galaxy (see, e.g. refs [@Taillet:2002ub; @Maurin:2002uc]). This choice could be thought as an over-simplification. A more realistic 2D geometry commonly used is to consider the Galactic halo as a cylindrical box of radius 20 kpc and half-thickness $L=10$ kpc, where CR sources lie uniformly in the disk and the Earth is set at 8.5 kpc from its center (for an illustration see e.g. Fig.10 of [@GenoliniEtAl2015]). However, using this geometry we have found that variations of the best fit values () for all parameters, except $K_{10}$, are within their respective uncertainties. In fact, in this 2D (uniform disk) case, a degeneracy between $K_{10}$ and $L$ is still present, but is no longer described by the relation ; the escape from the radial boundaries increases with increasing $L$. Starting from , for each benchmark, we summarize below our empirical prescription to go from 1D to 2D, only for the parameters which drift with $L$. Note that the preferred value for $V_{\rm A}$ in the  model is now zero, although the uncertainty on this parameter is quite large.

$$\begin{aligned}
&& \text{\BIG: }\,\frac{K^{\rm 2D}_{10}}{\tanh{(L^{1.1}/10.1)}}=0.25^{+0.04}_{-0.02}\,{\rm kpc^2/Myr}, \quad
\delta^{2D}\times \tanh{(L^{0.4}/0.77)}=0.50^{+0.02}_{-0.04}\quad\text{and} \quad \begin{cases} V_{\rm A}=0^{+80} \,{\rm km.s^{-1}} \\ R_{\rm l}=4.4_{-0.2}^{+0.2} \,{\rm GV}  \\ \delta_{\rm l}=-0.83_{-0.3}^{+0.3} \end{cases}\,,\\
&& \text{\SLIM:  }\,\frac{K^{\rm 2D}_{10}}{\tanh{(L^{1.1}/10.1)}}=0.25^{+0.02}_{-0.02}\,{\rm kpc^2/Myr}\quad\text{and} \quad \delta^{2D}\times \tanh{(L^{0.4}/0.77)}=0.50^{+0.02}_{-0.02}\,,\\
&& \text{\QUAINT:  }\,\frac{K^{\rm 2D}_{10}}{\tanh{(L^{1.1}/10.1)}}=0.30^{+0.02}_{-0.07}\,{\rm kpc^2/Myr}\quad\text{and}\quad\delta^{2D}\times \tanh{(L^{0.4}/0.77)}=0.44^{+0.06}_{-0.02}\,.\end{aligned}$$

[^1]: Or, in the original 1793 French version: “[*Une grande responsabilité est la suite inséparable d’un grand pouvoir*]{}”. French Revolution Parliamentary Archives, “Tome 64 : Du 2 au 16 mai 1793, Séance du mardi 7 mai 1793, page 287, available e.g. at https://frda.stanford.edu/fr/catalog/wx067jz0783\_00\_0293.
