Solar and nuclear physics uncertainties in cosmic-ray propagation
Recent data released by the AMS experiment on the primary spectra and secondary-to-primary ratios in cosmic rays (CRs) can pose tight constraints to astrophysical models of CR acceleration and transport in the Galaxy, thereby providing a robust baseline of the astrophysical background for dark matter search via antimatter. However, models of CR propagation are affected by other important sources of uncertainties, notably from solar modulation and nuclear fragmentation, that cannot be improved with the sole use of the AMS data. The present work is aimed at assessing these uncertainties and their relevance in the interpretation of the new AMS data on the boron-to-carbon (B/C) ratio. Uncertainties from solar modulation are estimated using improved models of CR transport in the heliosphere constrained against various types of measurements: monthly resolved CR data collected by balloon-born or space missions, interstellar flux data from the Voyager-1 spacecraft, and counting rates from ground-based neutron monitor detectors. Uncertainties from nuclear fragmentation are estimated using semiempirical cross-section formulae constrained by measurements on isotopically resolved and charge-changing reactions. We found that a proper data-driven treatment of solar modulation can guarantee the desired level of precision, in comparison with the improved accuracy of the recent data on the B/C ratio. On the other hand, nuclear uncertainties represent a serious limiting factor over a wide energy range. We therefore stress the need for establishing a dedicated program of cross-section measurements at the (100 GeV) energy scale.
pacs:98.70.Sa, 96.50.S, 96.50.sh, 13.85.Tp, 95.35.+d
An important challenge in astroparticle physics is the determination of the injection and transport parameters of Galactic cosmic rays (CRs), that establish the relation between the measured fluxes near Earth and the properties of their sources. Astrophysical models of CR propagation account for particle acceleration mechanisms in Galactic sources, diffusive transport processes in the turbulent magnetic fields, and production of secondary particles from CR collisions with the gas. The vast majority of Galactic particles, such as protons, He, and C-N-O nuclei are of primary origin, originating by diffusive-shock acceleration in supernova remnants (SNRs). Rarer particles such as H, He or Li-Be-B elements are called ssecondary, as they are predominantly generated by interactions of primary CRs in the interstellar medium (ISM). Along with the spectra of primary CRs, the use of secondary-to-primary nuclear ratios –and most notably the boron-to-carbon (B/C) ratio– is of paramount importance for constraining the key parameters that regulate acceleration and transport properties of these particles (Grenier et al.2015, ). Understanding CR propagation is essential to predict the secondary production of antimatter particles , , or , as it constitutes the astrophysical background for the search of new-physics signals from the annihilation of dark matter particles. Along with the background, the determination of the CR transport parameters is also important for modeling the signal of dark-matter induced antiparticles. The B/C ratio in CRs has been precisely measured in the 0.5-1000 GeV/n energy range by the Alpha Magnetic Spectrometer (AMS) experiment in the International Space Station (Aguilar et al., 2016b). Recent results from AMS on the fluxes of protons and helium (Aguilar et al., 2015a), antiparticles (Aguilar et al., 2016), and preliminary measurements on high-energy Li-Be-B spectra have generated novel ideas (Serpico, 2015; Tomassetti & Donato, 2015; Amato & Blasi, 2017; Blasi, 2017) and advanced analysis efforts (Feng et al., 2016) aimed at the determination of subtle effects in CR propagation.
With these new standards of experimental accuracy, however, we consider critical to make an assessment of those unavoidable sources of model uncertainties that cannot be directly constrained by the AMS data: solar-physics uncertainties, from the modulation effect of CRs in the heliosphere, and nuclear-physics uncertainties, from fragmentation processes of light elements in the ISM.
The goal of this paper is to assess these uncertainties and their relevance in the interpretation of the B/C data. Modeling solar modulation is of crucial importance to study the CR transport processes that reshape the CR spectra and secondary-to-primary ratios in the 0.1-20 GeV/n energy region. This problem is often underestimated in CR astrophysics, with the frequent use of ultrasimplified approaches based on poorly justified physics assumptions. In this respect, a substantial advance can be done thanks to the availability of numerical solvers for CR transport in the heliosphere and to the recent release of very precious sets of time-resolved data. The second issue deals with nuclear physics inputs in CR propagation models. The determination of the transport parameters relies on the calculation of the secondary production rate of Be and B nuclei, for which several energy-dependent cross-section estimates are required. Propagation models make use of semiempirical cross-section formulae that are based on accelerator data at the sub-GeV/n energy scale, and then extrapolated at the relevant energies. The accuracy of the inferred propagation parameters is clearly linked to the quality of the available measurements on nuclear fragmentation. Both sources of uncertainties are expected to affect the interpretation of secondary CR data, as they blur the connection between experimental observations and CR propagation effects, therefore hiding valuable pieces of information that are encoded in the data. A proper assessment of systematic uncertainties in the models is therefore essential, in CR physics analysis, for a full exploitation of the AMS data.
In this section, we outline the various types of measurements that are used in this work. A reference model of CR propagation, which we use to provide local interstellar (IS) particle fluxes and nuclear ratios, is constructed using recent primary CR data on proton and helium from AMS (Aguilar et al., 2015a) in the GeV/n - TeV/n energy region, new measurements from CREAM-III (Yoon et al., 2017) at multi-TeV/n energies, and direct IS flux data from Voyager-1 (Cummings et al., 2016) in the sub-GeV/n energy window. To model the B/C ratio, we make use of data released recently by AMS and Voyager-1 (Aguilar et al., 2016b; Cummings et al., 2016). The AMS measurements provide precise B/C ratio between 0.5 and 1000 GeV/n of kinetic energy under a period of medium-high solar activity (from May 2011 to May 2016) corresponding to 2.310 and 8.310 boron and carbon events, respectively. The Voyager-1 data are available only at 0.1 GeV/n energies. Nonetheless, these data are very precious because they constrain directly the low energy B/C ratio outside the heliosphere.
In order to model the solar modulation effect of CRs in the heliosphere, we made use of a large variety of data collected at Earth, in high-altitude balloons, or in space missions. In particular, we use the data provided from the worldwide network of neutron monitor (NM) detectors. We have retrieved monthly averaged measurements from NM stations in Newark, Oulu, Apatity, and Jungfraujoch (Steigies, 2015). These data consist in energy-integrated counting rates, corrected for detector efficiency and atmospheric pressure, between January 2005 and January 2017. NM rates are used in combination with direct measurements of the CR proton flux performed by space missions or high-altitude balloons projects. We use monthly resolved proton data from the PAMELA experiment (Adriani et al., 2013) collected under solar minimum conditions between July 2006 and January 2010, yearly-resolved data recently released by the EPHIN/SOHO satellite (Kühl et al., 2016), collected over the period 2000-2016, and proton flux measurements from the BESS Polar-I (Polar-II) mission between 13 and 21 December 2004 (between 23 December 2007 and 16 January 2008) (Abe et al., 2012).
Finally, we have made use of a large compilation of accelerator data on nuclear fragmentation reactions. Most of the data consist in energy-dependent cross-sections for the production of B, B, Be, Be, and Be isotopes from the collisions of B-C-N-O elements off hydrogen. These reactions have been measured between 50 to 5000 MeV/n by a number of experimental projects : Read & Viola 1984 (); : Webber et al. 1998 (); : Webber et al. 1990 (); : Olson et al. 1983 (); : Fontes 1977 (); : Korejwo et al. 2000 (); : Korejwo et al. 2001 (); : Radin et al. 1979 (); : Ramaty et al. 1997 (); : Webber et al. 1998 (); : Raisbeck1971 & Yiou 1971 () (symbols are used later). Along with reactions on isotopically separated projectiles and fragments, we also use measurements on charge-changing reactions, from Webber et al. 2003 (), which are available from 0.1 GeV/n to nearly GeV/n energies. Finally, cross-sections involving heavier nuclei, such as Ne, Si, or Fe, are extracted from GALPROP code of CR propagation (Strong & Moskalenko, 1998).
Iii Reference model of CR propagation
In this section, we define a reference model of CR diffusive propagation into an homogeneous cylindrical-shaped Galactic halo. Models of CR propagation in the Galaxy employ fully analytical (Thoudam & Horandel, 2014; Tomassetti, 2012a). semianalytical (Jones et al., 2001; Maurin et al., 2001), or fully numerical calculation frameworks (Strong & Moskalenko, 1998; Evoli et al., 2008; Kissmann, 2014). We utilize semianalytical calculations implemented under Usine, a global toolkit that allows to solve the CR propagation equation for given nuclear and astrophysical inputs such as ISM gas distribution, source distribution and spectra, fragmentation cross-sections, boundary conditions (Maurin, 2015). We implement a Kraichnan-like diffusion model with minimal-reacceleration as a benchmark. The propagation equation for a -type CR particle is written as:
where is the CR density of the species per unit of rigidity . The source term, , includes the primary injection spectra, from SNRs, and the terms describing secondary CR production in the ISM or decays. The rigidity dependence of the primary term is traditionally modeled as a power law, , with a universal injection index for all primary elements. This picture, however, has recently been challenged by the observation of puzzling features in the spectrum of primary CRs (Serpico, 2015; Amato & Blasi, 2017): a difference in spectral index between protons and nuclei, with /He at - GV, and a common spectral hardening of all fluxes at 375 GV (Aguilar et al., 2015a). Hence, we adopt a primary source term of the type:
|Injection, low energy shaping||0|
|Injection, smoothing factor||4|
|Diffusion, normalization||[kpc myr]||0.026|
|Diffusion, scaling index||1/2|
|Diffusion, low energy shaping||-1|
|Reacceleration, Alfvén speed||[km s]||15|
where is particle dependent, and in particular is higher for proton, and describes a universal hardening occurring in all primary spectra at rigidity . The parameter describes the behavior of the spectrum at low energy, with from the simplest linear diffusive-shock acceleration. The composition factors are normalized to the primary CR abundances at reference rigidity .
The secondary production term describes the products of decay and spallation of heavier CR progenitors. For decay of radioactive nuclei with corresponding lifetime , the source term is . The occurrence of fragmentation processes is described by:
where is the differential rate of -type particle production at kinetic energy , from collisions of -type CR particles, with density , with -type targets of the ISM component (with number density ). The ISM is essentially composed by hydrogen H and helium He atoms, with and . Nuclear interactions couple the equation of each -type nucleus to those of all heavier -nuclei. The system is resolved by starting from the heavier nucleus, which is assumed to be purely primary, and proceeding downward in mass. The fragmentation loop is repeated twice. The diffusion coefficient is taken as spatially homogeneous and rigidity dependent:
where sets its normalization at reference rigidity , and the scaling index describes its rigidity dependence. The factor allows for a low rigidity change in the diffusion regime. In particular, negative values for the parameter are used to to effectively account for a faster CR diffusion in the sub-GV rigidity region, which may be expected, e.g., from wave damping (Ptuskin et al., 2006). This effect is reflected by the observed peak in the B/C ratio at 1 GeV/n. The characteristic shape of the B/C ratio in the GeV/n energy region, however, may be also ascribed to other processes such as advection of strong diffusive reacceleration, if not to a change in spallation cross-sections (Maurin et al., 2001; Jones et al., 2001). Reacceleration is described with a diffusion coefficient in rigidity space which is linked to spatial diffusion by the relation:
The parameter describes the average speed of Alfvén waves in the ISM medium. We consider modest reacceleration ( 15 km/s) under a scenario with Iroshikov-Kraichnan diffusion (1/2). Other parameters are listed in Table 1. The equation is solved in steady-state conditions at the boundary of the diffusion region. The IS fluxes entering the solar system are then computed for each species:
where and are the CR charge and mass number, and the fluxes are given in units of kinetic energy per nucleon . The calculated ISs of proton and helium, along with the B/C ratio, are shown in Fig. 1 in comparison with new data from Voyager-1, AMS, and CREAM-III. The connection between model parameters and physics observables is discussed in several works (Maurin et al., 2001; Maurin et al. 2010, ; Grenier et al.2015, ). In the purely diffusive regime, primary spectra behave as , while secondary-to-primary ratios follow , giving direct constraints to the parameters , , , and . In the presence of energy changes, the connection between parameters and physics observables is blurred by other effects, e.g., reacceleration, which significantly reshape the B/C ratio at 0.1-50 GeV/n. To extract information in this energy region, our ability in modeling the unavoidable effect of solar modulation is then essential, along with calculations of fragmentation reactions that lead to the production of boron nuclei. The interpretation of primary spectra and secondary-to-primary ratios may be different in scenarios where CR anomalies are ascribed to source components of Li-Be-B nuclei, nearby SNRs, or changes in diffusion (Tomassetti & Donato, 2015; Blasi, 2017; Tomassetti, 2012a; Maurin et al. 2010, ; Serpico, 2015). In this work, however, these issues are not addressed. The remainder of this paper is devoted to assess solar and nuclear physics uncertainties. For this task, the reference model present here will serve as the input model.
|Location||Newark, Delaware||Oulu, Finland||Apatity, Russia||Jungfraujoch, Switzerland|
|Coordinates||39.68 N 75.75 W||65.05 N, 25.47 E||67.55 N, 33.33 E||46.55 N, 7.98 E|
|Altitude||50 m||15 m||177 m||3570 m|
|Cutoff||2400 MV||810 MV||480 MV||4500 MV|
|a||(697.5 4.5) MV||(680.3 4.8) MV||(709.8 4.6) MV||(689.9 4.5) MV|
|b||0.39 0.01||0.40 0.01||0.42 0.01||0.38 0.01|
Iv Solar physics uncertainties
When entering the heliosphere, CRs travel across the outflowing solar wind embedded in its turbulent magnetic field, where they undergo convection, diffusion, energy losses, and drift motion. To give a complete description of the solar modulation of CRs, all these effects have to be carefully modeled. Our work takes advantage of recent progress in experimental measurements (Bindi et al., 2017), theoretical developments (Potgieter, 2017), and the availability of numerical codes for solar modulation (Kappl et al., 2016).
iv.1 Solar modulation calculations
The goal of solar modulation calculations is to predict the modification of the CR energy spectrum in heliosphere, at a given epoch , via a transformation of the type:
The problem can be described in terms omnidirectional phase-space CR density , where the evolution is governed by the Krymsky-Parker equation:
where the CR number density is given by . The various terms of Eq. 8 represent convection with the solar wind, of speed , which is caused by the expansion, particle drift motion with velocity , caused by gradients and curvature effects, spatial diffusion with tensor representing its symmetric part, and energy losses due to adiabatic deceleration (Potgieter, 2017). For each CR particle species, Eq. 8 is set to zero in order to calculate steady-state solutions , which is a reasonable assumption for long-term modulation studies where the key parameters change gradually. The initial conditions are placed at the boundary of the modulation region, where the reference model IS fluxes are provided from Eq. 6. A simplified approach that is commonly used in the CR astrophysics is the so-called force-field (FF) approximation. It consists in solving Parker’s equation for a radially expanding wind, of speed , a fully isotropic diffusion coefficient with spatial-dependent part , after neglecting drift and loss terms. The FF approximation provides an analytical one-to-one correspondence between arrival fluxes at and IS fluxes at in terms of a lower rigidity shift, , where:
The parameter is called the modulation potential. In terms of CR energy spectra , at a given epoch, kinetic energy per nucleon is shifted by , so that the modulated flux is given by:
In this work, we do not make the FF approximation to compute the CR modulation, but we use quantity as input parameter. We solve Eq. 8 using a 2D numerical model with azimuthally symmetric spherical coordinates: radius and helio-colatitude (Kappl et al., 2016). The solar wind is taken as radially flowing with speed 400 km s. The parallel component of the diffusion tensor is taken as , in units of cm s, and its perpendicular component is . The adimensional scaling factor accounts for the time-dependence of the problem. The regular magnetic field (HMF) follows the usual Parker model, with , where , rad s is the Sun rotation speed, 3.4 nT sets the HMF at 1 AU, and sets the magnetic polarity state of the Sun. The polarity is negative (positive) when the HMF points inward (outward) in the Northern Hemisphere. This model accounts for drift effects, implemented as in (Strauss et al., 2012; Kappl et al., 2016). Drift is important near the Heliospheric current sheet (HCS) where polarity changes from north to south. When CR particles travel close to the HCS, their drift speed is a large fraction of their total speed, and 10 times larger than . Away from the HCS, is of the order of . The signs of its components depend on the product , which is at the origin of the charge-sign dependence of CR modulation (Strauss et al., 2012; Kappl et al., 2016). Because of drift, particles and antiparticles sample different regions of heliosphere and their role interchange with polarity reversal. The angular extension of the HCS, which sets its waviness, is described by the tilt angle . To obtain steady-state solutions of Eq. 8 we employ a differential stochastic approach. In particular we use the SolarProp numerical engine (Kappl et al., 2016). With this method, Eq. 8 is resolved by means of Monte Carlo generation of large samples of pseudoparticle trajectories. In practice, pseudoparticles are backward propagated from Earth to the modulation boundaries, that we set at the Heliopause HP, AU. We also account for the termination shock (TS), placed at distance 85 AU (Langner et al., 2003), by smooth increase (decrease) of the HMF (wind speed) radial profiles of a factor 3 across the shock position, where is the TS compression ratio (Potgieter, 2017). Due to the TS, CRs in the outer heliosphere follow slower diffusion and weaker convection (Langner et al., 2003; Strauss et al., 2012). The IS fluxes , which state the boundary conditions at the HP, are calculated within the reference model tuned against Voyager-1 and AMS (Cummings et al., 2016; Aguilar et al., 2015a). Further investigations on the phenomenology of this model will be presented in a future work.
iv.2 Uncertainties from solar modulation
To assess the uncertainties on the B/C calculations that arise from solar modulation modeling, we perform a global data-driven reconstruction of the time evolution of the CR flux based on NM counting rates and CR proton data. With this reconstruction, we constrain the key input parameters of the solar modulation model and their time dependence, which allows us to estimate their influence on the predicted time evolution of the B/C ratio. This procedure provides a time-dependent uncertainty band for our model prediction, calculated on monthly basis, which can therefore be averaged over a the period of the AMS observation time.
From the implementation described in Sec. IV.1, two time-dependent parameters have to be determined: the normalization of the diffusion coefficient, , and the tilt angle of the HCS, . To model the tilt angle, we adopt a smooth interpolation of the time-series based on the “radial model” reconstruction (Hoeksema, 1995). This reconstruction is provided on 10-day basis by the Wilcox Solar Observatory (WSO, 2017). To determine the time evolution of the diffusion coefficient normalization, we make use of the monthly resolved series of NM rates. For this task, we convert monthly average NM counting rates into a time-series of modulation potential that, as discussed, arise from the FF solution of Eq. 8. An inverse relationship between and is suggested by Eq. 9, from which . Here we simply write , where the coefficients and are free parameters that we determine using CR flux data, and is determined using NM counting rates. We recall that, although we make use of the parameter obtained in the context of the FF approximation, we do not make the FF approximation in our flux calculations.
The procedure to determine the time-series follows early works (Tomassetti, 2017; Ghelfi et al., 2016). For a given NM detector , located altitude and geomagnetic cutoff , the link between the counting rate and top-of-atmosphere CR fluxes (with , He) is expressed by:
where is a transmission function around the cutoff value , parametrized as a smoothed step function (Smart & Shea, 2005; Tomassetti, 2015), and is the -particle dependent detector response function (Maurin et al., 2015). In particular we use a factorized form, , where accounts for time and energy dependencies of the NM response, including the development of hadronic cascades (Cheminet et al., 2013), and the factor set the absolute normalization and its altitude dependence. To compute Eq. 11, Equation 11 is calculated using monthly averaged NM rates, corrected for detector efficiency and pressure, provided between January 2005 and January 2017 from various stations (Steigies, 2015). The NM detectors considered in this works are listed in Table 2. For a given station , the parameter is obtained from the request of agreement between the calculated rate and the observed rate , together with the requirement that where the integration is performed over the observation periods (Tomassetti, 2017). For all detectors we consider 12 years, between January 2005 and January 2017. The time-series of monthly reconstructed resulting from this procedure is shown in Fig. 2a using the NM detectors of Table 2. Uncertainty bands are shown for each station, corresponding to 25 MV (Ghelfi et al., 2016). This is also the level of discrepancies among the various NM responses. The -series calculated from Ref. (Usoskin et al., 2011) is superimposed as dashed lines, for reference. This time-series is in agreement with our reconstruction.
From the reconstructed time-series of modulation potential, the parameters and that specify the time-evolution of the diffusion scaling . are determined using CR proton data from PAMELA, EPHIN/SOHO, and BESS-Polar, for a total of to 3993 CR data points. Using CR flux measurements , collected at given epoch and energy , we build a global estimator:
The proton flux calculations are carried out from numerical solutions of Eq. 8 using input parameters and . Both time-series and have been smoothly interpolated with monthly resolution.
The best fit parameters and are obtained using standard -minimization techniques. The minimization required a large scan over the - parameter space, for a total simulation of 1.9 billion particle trajectories. This procedure provides the time and energy dependence of the proton flux and its associated uncertainties. The time profile of the CR proton flux at GeV is shown in Fig. 2b, in comparison with the data. The procedure was repeated using the four time-series of , giving slightly different –but consistent– best fit parameters. The best fit results are listed in Table 2 for the various NM detectors considered. Nonetheless, the robustness of the predictions is mostly related to the constraints provided by the CR proton flux data, which is plotted as shaded band in Fig. 2b. From the figure, one can also notice a clear anticorrelation between CR flux and modulation potential. The modulated fluxes presented in the figure are calculated on monthly based between 2005 and 2017, therefore covering a 11-year period corresponding to a complete solar cycle. This period also includes magnetic reversal, which occurred in early 2013, where the Sun’s polarity state switched from to , shown in the figure as vertical shaded bar. It can be seen that the calculations agree well with the AMS proton data collected between May 2011 and November 2013, although these data have not been included in the estimator. We also show, in Fig. 3, the inferred diffusion scaling obtained by performing single fits to the CR energy spectra plotted against the modulation potential reconstructed at the same epoch using data from the NEWK station. The solid line represents the relation , for the NEWK station, and its uncertainty band. This figure shows the inverse relationship between that the two quantities. Once the parameters are constrained, the model can be used to predict CR flux of other species including heavier nuclei or antiparticles. In Fig. 4 we plot the time evolution of the energy spectra of boron (a) and carbon (b), along with their ratio (c) calculated between 2005 and 2017. The figure represents the flux (ratio) values as a function of energy and time. In our calculations for the B/C ratio, we have accounted for various sources of model uncertainties, including uncertainty on the input IS fluxes, uncertainty on the input parameters and on the global fitting procedure. The resulting error is mildly time-dependent. The resulting uncertainties on the B/C ratio are presented in Sec. VI after performing a time average over the AMS observation period.
iv.3 On the force-field approximation
In our calculations, we made use of the modulation parameter as input. This parameter is obtained in the context of the FF approximation, although we adopted a numerical approach to solve the Krymsky-Parker equation. On may wonder if a fully FF-based approach may reach an adequate level of accuracy. In Fig. 2b, NM-driven FF modulation potential derived from NEWK is compared with the -values determined from direct fits to CR spectra (CR-driven) of the various data sets. While the overall time profile is well reproduced in both approach, some inconsistencies are apparent, e.g., in the PAMELA data. Symmetrically, in Fig. 2c, numerical calculations (thick line) are shown together the NM-driven FF calculations (thin dashed line) obtained using the parameter of Fig. 2a to FF-modulate the proton IS. Again, some inconsistency can be observed with the PAMELA data. Similar tensions were also noted in other studies (Cholis et al., 2016). The origin of this tension is that, with the integration of the NM rates of Eq. 11, the modulation parameter is determined at the (10 GeV) energy scale. In contrast, with the CR-driven approach, the fit is performed over the whole energy spectra, i.e., down to 80 MeV for PAMELA, which is where the FF approximation breaks down. More in general, the larger discrepancies are expected during periods (when drift is relevant) and especially at low energy or during high level of solar activity (where the modulation effect is stronger). Nonetheless, when the FF is directly applied for fitting CR data (with as free parameter) the flux calculations can be described reasonably well. This is illustrated in Fig. 5. Along with the IS proton flux, solar modulated spectra are shown from numerical calculations, from NM-driven FF calculations (with determined from NEWK data), and from CR-driven FF calculations (where the flux is directly fit to the CR data of the figure). The CR-driven approach leads to a fairly good description of the data, within 10-15 % of precision for the fluxes of Fig. 5. The main weakness of this approach lies in its limited predictive power. For instance, FF-calculations calibrated with proton data cannot be used to describe antiprotons, which require different -values. Similarly, FF-calculations calibrated with NM data under a given polarity state cannot be applied to the opposite polarity. Hence the FF method can provide an adequate description of the CR spectra, but it cannot be used in a predictive way, e.g., for estimating the astrophysical background of CR antiparticles. Having a predictive model is also necessary to describe the time-unresolved collected over long exposures. For instance, the AMS data represent a time average over years of observation time. The CR flux evolves during this period, which is only covered by NMs, on much shorter timescales. To account for these issues of the simple FF method, recent works proposed easy-to-use generalizations of the FF formula through the introduction of a charge- or rigidity-dependent parameter (Cholis et al., 2016).
V Nuclear physics uncertainties
The determination of the CR transport parameters using the B/C ratio relies obviously on calculations of the boron production rate from the disintegration or decay of the heavier nuclei. The accuracy of secondary production calculations therefore depends on the reliability of the production and destruction cross-sections employed.
v.1 Cross-section calculations
The general source term of secondary production of CR nuclei is given in Eq. 13. This expression is further simplified by the straight-ahead approximation, according to which it the secondary nucleus is assumed to be ejected with the same kinetic energy per nucleon of the fragmenting projectile, i.e., . Furthermore, the composition of the interstellar gas is dominated by hydrogen (90 %) and helium (10 %), so that reactions involving heavier target can be safely neglected. Hence the source term for -type particle production becomes:
In CR propagation models, semiempirical formulae are used to predict the cross-sections at any energy and for all relevant projectilefragment (PF) combinations. These algorithms are based on observed systematics in mass yields, charge dispersion, and energy dependence of the fragments produced in nucleus-nucleus collisions. Popular algorithms are YIELDX (Silberberg et al., 1998) and WNEW (Webber et al., 1990b; Webber et al. 2003, ). The GALPROP code makes use of parametric formulae obtained from fits to the data or from nuclear codes CEM2k and LAQGSM, eventually normalized to the data (Strong & Moskalenko, 1998; Moskalenko & Mashnik, 2003). Along with mass-changing and charge-changing reactions, cross-sections for isotopically separated targets or fragments have been measured by several experiments (see Sec.II). In spite of such large collection of PF reactions, the data are available only at 0.1 5 GeV/n energies. In this work, we mostly are interested in the production of B, B, Be, Be, and Be isotopes from interactions of B, C, N and O off hydrogen and helium targets. These reactions account for 90% of the Be and B production. The remaining reactions involve fragmentation of heavier nuclei such as Ne, Si, or Fe, which give a minor contribution to the boron production uncertainty. It is important to stress that the fragmentation network of CRs in the Galaxy may involve several steps where the disintegration or decay of intermediate particles contributes to the abundance of a given final-state CR nucleus. If the intermediate nucleus is stable or long lived (i.e., its lifetime is larger in comparison with the CR propagation timescale) this species has to be properly accounted in the CR propagation network. Examples of this kind are OCB or ONB reaction processes. The multistep nature of CR fragmentation is illustrated in the alluvial diagram of Fig. 6 which links the secondary Be-B isotopes (left blocks) to their purely primary progenitors (right blocks) via two stages of fragmentation calculated at 1 GeV/n of kinetic energy. The stream fields between the blocks represent the fragmentation reaction channels (to be read from right to left) sized accordingly to their contribution to secondary abundances. From this graph it can be seen, for example, that the B isotope is mostly generated from C and O as expected, while a non-negligible contribution is represented by fragmentation of secondary species such as B, N, or C. Thus, the usual view of purely primary CRs fragmenting into secondaries is only an idealized approximation. The actual calculations account for a complex multistep reaction network to determine the abundance of Li-Be-B particles.
Another situation is when the intermediate nucleus is unstable with a short lifetime to the CR propagation timescale. An example is given by B particles that can be generated either by “direct” one-step processes (e.g., CB) or from production of intermediate isotopes (e.g., C) with subsequent decay into B. In all these cases, intermediate short-lived nuclei are treated as “ghost” so that they are accounted using effective cross-section formulae:
where is the branching ratio for the decay channel (Maurin et al., 2001). Ghost nuclei do not experience CR propagation within their lifetime, but they can be detected in laboratory measurements, depending on the time-resolution of the experiment. An instructive example is given in Fig. 7 for the CB reaction channel. In this reaction, measurements and calculations are shown for the “direct” B production (green line) and for the cumulative production (purple line) which includes C ghost nuclei. Experimentally, the latter is measured by the sole mass identification of end-products. This processes is the relevant one for CR propagation and represents the effective cross-section of Eq. 14. All fragmentation cross-sections used in this work have to be meant as effective ones. Finally, to properly compute the abundance of boron at the GeV/n energy scale, it is also important to model the production of beryllium isotopes via iterative calculations. Beryllium and boron share the same progenitors, but the small component of Be contained in it the Be flux decays radioactively into B, contributing to nearly 10 % in the B/C ratio at low energy.
v.2 Uncertainties from nuclear fragmentation
We utilize a large compilation of cross-section data on isotopically resolved reactions. An extensive survey of the literature can be found in other works (Tomassetti, 2015; Moskalenko & Mashnik, 2003; Moskalenko et al., 2013). To determine the cross-section uncertainties using the data, we perform a data-driven renormalization of the GALPROP parametrizations . Our approach follows the procedure of earlier studies on H – He isotopes (Tomassetti, 2012b) and heavy nuclei (Moskalenko et al., 2013; Tomassetti, 2015), although we introduce some improvements. For all PF channels, parametric cross-section formulae have been redetermined via minimization of the following quantity:
where is the experimental error on the -th data point for the considered PF channel.
|Proj Frag||GALPROP||WNEW||YIELDX||Data sets|
|O B||27.34 mb||14.57 mb||14.41 mb||(25.66 1.06) mb||1.34||(34.37 1.43) mb||0.94 0.04||1.01 0.02|
|O B||11.00 mb||9.52 mb||8.81 mb||(11.92 0.52) mb||1.34||(15.97 0.70) mb||1.08 0.04||0.98 0.03|
|N B||26.12 mb||26.12 mb||21.71 mb||(30.63 2.48) mb||1.31||(40.04 3.24) mb||1.17 0.08||0.67 0.51|
|N B||9.56 mb||8.81 mb||7.63 mb||(9.69 0.77) mb||1.31||(12.66 1.01) mb||1.01 0.07||1.54 0.05|
|N B||29.22 mb||29.98 mb||26.66 mb||(29.80 1.08) mb||1.21||(35.94 1.30) mb||1.02 0.03||0.95 0.02|
|N B||10.44 mb||10.64 mb||9.18 mb||(10.15 0.84) mb||1.21||(12.24 1.01) mb||0.97 0.07||0.99 0.05|
|C B||56.88 mb||54.86 mb||52.83 mb||(54.73 2.57) mb||1.29||(70.50 3.32) mb||0.96 0.04||0.90 0.01|
|C B||12.30 mb||16.21 mb||11.59 mb||(12.05 0.58) mb||1.29||(15.52 0.74) mb||0.98 0.04||1.00 0.02|
|O Be||2.14 mb||1.34 mb||2.07 mb||(1.90 0.13) mb||1.47||(2.79 0.19) mb||0.88 0.05||0.90 0.01|
|O Be||3.48 mb||3.35 mb||3.51 mb||(3.40 0.22) mb||1.47||(4.99 0.32) mb||0.97 0.06||0.98 0.03|
|O Be||10.00 mb||8.75 mb||8.92 mb||(8.97 0.29) mb||1.47||(13.16 0.42) mb||0.89 0.03||0.98 0.02|
|N Be||1.75 mb||1.06 mb||1.81 mb||(1.73 0.21) mb||1.43||(2.47 0.30) mb||0.99 0.09||0.95 0.07|
|N Be||10.10 mb||7.46 mb||8.47 mb||(7.90 0.47) mb||1.43||(11.29 0.67) mb||0.78 0.04||1.05 0.08|
|C Be||3.94 mb||2.05 mb||3.41 mb||(3.61 0.27) mb||1.41||(5.07 0.38) mb||0.91 0.06||0.92 0.01|
|C Be||6.76 mb||5.31 mb||4.98 mb||(6.63 0.29) mb||1.41||(9.32 0.41) mb||0.98 0.04||1.00 0.02|
|C Be||9.58 mb||10.32 mb||10.76 mb||(8.88 0.30) mb||1.41||(12.48 0.42) mb||0.93 0.03||1.09 0.19|
|B B||38.91 mb||42.58 mb||38.91 mb||(37.83 9.25) mb||1.29||(48.63 11.89) mb||0.97 0.21||1.10 0.14|
|B Be||12.95 mb||5.90 mb||4.56 mb||(7.39 0.90) mb||1.40||(10.36 1.26) mb||0.57 0.06||0.90 0.17|
|B Be||10.00 mb||15.27 mb||8.01 mb||(7.22 1.08) mb||1.40||(10.13 1.51) mb||0.72 0.09||0.91 0.12|
|B Be||4.48 mb||4.48 mb||3.63 mb||(4.68 0.49) mb||1.40||(6.57 0.69) mb||1.05 0.10||0.90 0.11|
The parameters and represent normalization and energy scale, respectively, and are determined from the available data for each PF channel. For some channels, measurements are available only in narrow energy ranges so that, for these reactions, the parameter is poorly determined. Nonetheless, in all channels, the procedure returns new cross-sections along with uncertainties corresponding to one-sigma confidence intervals. The fitting results are summarized in Table 3, where the renormalized cross-sections are compared with the predictions of other algorithms at 10 GeV/n. The best fit parameters and are also listed. The symbols corresponding to the various data sets are encoded in Sec.II. It can be noted that the renormalized cross-sections are often close to the original values, because the original GALPROP model is built with the help of a large set of data. In contrast, the WNEW and YIELDX models show large discrepancies for some reactions (e.g., O B or B Be). In fact they are normalized to incomplete sets of data. If one switches among the different cross-section models, the resulting discrepancy amounts to about 20 % (Maurin et al. 2010, ). But this discrepancy is not fully representative of the systematic errors in cross-section, because the available data give tighter constraints. For instance, if the fit procedure is repeated using WNEW or YIELDX as the base, the results agree well at the percent level at energies above a few 100 MeV/n. Larger discrepancies persist in the lowest energy region ( 100-200 MeV/n), because at these energies many PF reactions are resonant shaped, but this feature is accounted for only by the GALPROP cross-sections.
A subdominant contribution comes from CR spallation off helium target, which constitutes a 10% fraction of interstellar gas. To account for this component, all the production and destruction cross-sections off He target are obtained using the scaling factor (Ferrando et al., 1988), which ranges from 1.2 to 1.5. This factor is listed in Table 3, along with the cross-section values for CR collisions off the He target. Due to the lack of data, uncertainties cross-sections reactions off helium cannot be directly estimated via a measurement-validated procedure. In this work, we conservatively assume a full correlation with the uncertainties in hydrogen production. Improving measurements on nucleus-helium collisions is clearly helpful for reducing nuclear uncertainties, but not critical. Reactions off He-target give a second order contribution to the total uncertainty on the secondary production terms. Finally, CR propagation models require estimates of destruction processes, for which we have employed dedicated parametrizations (Allison et al., 2016). These reactions are known with better precision in comparison to partial fragmentation reactions. Furthermore, it can be estimated that only a small fraction (10 %) of light CR nuclei (8) is involved in destructive processes, because the diffusion timescale of these particles is always dominant over spallation (Tomassetti, 2015c). This translates into negligible uncertainties for the Be-B equilibrium fluxes.
An important issue that has to be examined is the influence of systematic energy biases in the cross-section reactions. In the existing parametrizations, all cross-sections are assumed to be asymptotically energy independent above the energy of a few GeV/n energies. Models of CR propagation rely heavily on these extrapolations. However, the possible –unaccounted– presence of energy-dependent biases in these formulae may cause corresponding biases on the calculated B/C ratio which, in turn, would lead to a misdetermination of the CR transport parameters. For instance, energy-dependent bias in the boron production cross-section would directly affect the B/C-driven determination of the scaling parameter . Examples of this effect were illustrated in Maurin et al. 2010 (), and it was noticed that the application of a factor to the WNEW parametrization is fairly consistent with the data for 0.05. A slow increase of nucleus-nucleus cross-sections at high energy is also expected from recent developments based on the Glauber-Gribov model, and experimentally measured in proton-nucleus reactions. An example is shown in Fig. 8 for He+C and C+C reactions, where calculations based on the Glauber-Gribov model are compared with older calculations used in the GEANT-4 hadronic generator (Allison et al., 2016). Energy-dependent biases are not accounted for in the fitting method of Eq. 15, so that it returns energy independent errors. On the other hand, the presence of high-energy biases cannot be tested for the isotopically resolved channels of Table, 3, due to the lack of multi-GeV data. To tackle the issue, we make use of data on charge-changing reactions that are available to nearly 100 GeV/n of kinetic energy. Owing the scarcity of these data, we accumulate all charge-changing measurements involving light elements from to . We compute the discrepancy between cross-section data and best fit calculations as function of kinetic energy. Then, we adopt a criterion to determine the size of the allowed bias on statistical basis. From the relation with 0.1 GeV/n, we obtain . Thus, cross-section data are essentially consistent with a constant behavior (), bounded by a one-sigma error of which sets the level of energy-dependent systematic uncertainty. Charge-changing cross-sections are shown in Fig. 8 (c to f) for various reactions involving Be to Si nuclei, with the estimated uncertainty band. Each charge-changing reaction reflects the combination of several isotopic channels, where each channel carries its own errors on normalization and energy scale. Although these measurements can be partially correlated among channels, (due to, e.g., common sources of systematic errors from the experimental setup), this information is not available to us. In contrast, energy-dependent systematics are assumed to be “universal” for all reactions because, with the existing data, we were unable to perform a channel-by-channel determination. Nonetheless, we expect this contribution to be highly correlated among the various channels, featuring common dynamical aspects of nucleus-nucleus collisions (Allison et al., 2016).
Vi Results and discussions
Calculations on CR modulation (Sec. IV.1) and nuclear fragmentation (Sec. V.1) are now utilized to provide final uncertainty estimates for the B/C ratio. From the procedure outlined in Sec. IV.2, time-dependent uncertainties from solar modulation can be obtained for the CR fluxes of any nuclear species. boron and carbon fluxes are calculated by performing a time average over the AMS observation period , from May 2011 to May 2016, and using time bins of 1 month:
A similar procedure was used to obtain a time averaged and energy-dependent error band for the modulated B/C ratio near Earth. To assess the influence of the estimated cross-section errors in the near-Earth fluxes and in the B/C ratio, we proceed in calculating the uncertainties in the source terms of secondary isotopes:
Along with the source term, the estimated uncertainties have to be propagated to the near-Earth equilibrium fluxes after accounting for the effects of Galactic transport and Heliospheric modulation. To first approximation, one has at any energy , although the connection between and is slightly blurred by ionization, adiabatic cooling, or reacceleration effects. Thus we have performed a large number of propagation runs where the production terms of all secondary CRs are randomly varied, using Monte Carlo generation techniques, according to the Gaussian-like probability distribution of width . This provides Gaussian uncertainties for the secondary CR fluxes and in the B/C ratio calculated at different energies. Typical uncertainties are found to be 5-8 % for B production in the 10 GeV/n energy scale. All errors increase with energy because, as discussed in Sec. V.2, we have accounted for the possible presence of a common energy-dependent bias.
The solar-modulated B/C ratio is shown in Fig. 9a along with the two uncertainty bands. In Fig. 9b we plot the residuals between data () and model (), i.e., the quantity . Although the physical modeling of Galactic CR processes may be improved in several directions, we note, from this figure, that the discrepancy between data and model lies fairly well within the systematic uncertainties. At 100 GeV/n, the B/C ratio has a slight tendency to harden which is not recovered by the model and may point to unaccounted effects in CR propagation (Genolini et al., 2017). In the low energy region, the comparison of Fig. 9b with the data from Voyager-1 is made using IS flux calculations. With decreasing energy, below 1 GeV/n, the ratio decreases rapidly due to reacceleration or possible changes in diffusion (Ptuskin et al., 2006). Our calculations break at MeV/n. Below these energies, and down to 8 MeV/n, Voyager-1 measured a nearly energy independent B/C ratio that is at tension with CR propagation models (Cummings et al., 2016). This behavior can be further investigated with complementary data, such as charge and isotopic composition measurements on He, H or Li-Be, that are currently being measured by AMS. It worth noticing, in fact, that the Voyager-1 probe is better optimized for the energy determination of low charge nuclei.
The relative errors on the B/C ratio associated with both contributions are plotted in Fig. 10. Uncertainties from solar modulation (short-dashed lines) are estimated at the level of 10 % in the 100 MeV/n energies, then rapidly decreasing with increasing kinetic energy, below 5 % at 1 GeV/n, and below 1 % at 10 GeV/n. This trend reflects the characteristic of particle transport in the heliosphere, the effects of which become weaker and weaker with the increasing energies of CR particles. Uncertainties arising from nuclear fragmentation (long-dashed lines) are found to follow an opposite behavior. These errors lie at the level of 6-8 % in the low energy region, i.e., where cross-section data are available, and then following a progressive increase with energy that reflects the lack of high-energy measurements. In the figure, these curves are compared with the experimental uncertainties from the AMS and Voyager-1 measurements on the B/C ratio. The AMS data (filled circles) can be regarded as the potential level of precision to which CR injection and transport can be understood. Such a potential has seen a dramatic improvement in precision and energy range, thanks to AMS, turning from 15 % (in pre-AMS data) to the level 3 % at 50 GeV/n of kinetic energy per nucleon. As seen in the figure, it is reassuring to notice that solar physics uncertainties become rapidly subdominant at 1 GeV/n. Also, in the energy window 0.1-1 GeV/n, the experimental errors on the B/C provided by other experiments (George et al., 2009) are larger than those from solar modulation. The experimental errors from Voyager-1 (open squares) are shown in the figure for reference, as the uncertainty line is calculated for AMS. Furthermore, Voyager-1 data represent a direct measurement of the IS ratio B/C. These data are extremely useful, e.g., to rule out CR propagation scenarios without modeling the effects of solar modulation.
The major limiting factor is represented by nuclear fragmentation inputs that dominate the uncertainties up to 500 GeV/n of kinetic energy. The highest energy points of the AMS data are statistically dominated; thus we expect that the experimental errors will be reduced after longer exposure times. This region will also be covered, very soon, from the expected release of TeV/n nuclear data by the space experiments CALET, DAMPE or ISS-CREAM (Serpico, 2015).
A further point of future discussion (not inspected in this work) is the validity of the straight-ahead approximation, according to which a narrow distribution of the fragments energies is assumed. In previous studies on the Li-Be-B production, the change introduced by this approximation was found to be small in comparison with the precision of the data (Kneller et al., 2003). Within the accuracy of the new data, however, the legitimacy of this approximation may deserve a reevaluation.
Clearly, both sources of model uncertainties considered in this work affect the evaluation of the astrophysical background for dark matter searches, i.e., the secondary fluxes of , , or particles that have to be calculated with B/C-driven models of CR propagation. Along with B/C-driven uncertainties, however, astrophysical background calculations are directly affected by uncertainties on production and transport of antiparticles. In particular, modeling solar modulation of antinuclei demands to account for charge-sign dependent effects (and uncertainties). Thus, a numerical model of solar modulation including for drift motion, such as the one utilized here, is essential to get reliable estimates. Similarly, nuclear uncertainties on the production of antiparticles are an important factor, for estimating the background level, that is getting increasing attention (Feng et al., 2016; Kappl et al., 2015).
The investigation of energy spectra and composition of CRs addresses fundamental science questions: origin of CRs: source properties, nature of injection and acceleration processes; CR diffusion properties and interstellar turbulence; spatial extent and properties of the CR confinement region; presence of dark matter signatures in the CR spectrum. In spite of much effort in the development of CR propagation models, our knowledge of the underlying physics processes has been for long time plagued by large uncertainties and unresolved tensions between inferred parameters and theoretical expectations. The need for high-energy data on secondary-to-primary ratios has often been invoked in literature (Grenier et al.2015, ). In this respect, the past years can be considered a golden age for CR measurements. In comparison with previous measurements, the new AMS data on the B/C ratio provide a dramatic improvement corresponding to a factor 5 in the experimental accuracy. To exploit these new data, however, it is very important to precisely model the effect of solar modulation and to assess the corresponding uncertainties. Fortunately, a large wealth of useful data on CR modulation has recently become available (Bindi et al., 2017). New releases of local IS data and monthly resolved data of the CR flux are of invaluable help for modeling solar modulation. Using these observations in combination with the NM counting rates, we have developed a data-driven approach using a numerical model of solar modulation in order to predict the near-Earth B/C ratio and to assess its uncertainties. The input parameter of our model consist in a time-series of measurements on the HCS tilt angle and NM rates (where the latter have been converted into a more practical time-series of modulation potential ), hence our model is highly predictive. With the time-series, one may also employ the simple FF approximation for describing the data, but as discussed, the FF breaks down at low energies or in condition. Direct FF fits to CR data provides a fairly good description, but this approach introduces degeneracies for the CR astrophysics investigation and falls short of predicting antiparticle fluxes. The precision of new CR data demands a reliable modeling of solar modulation. It is also important to stress that the CR data reported by new-generation experiments are collected over ultra-long exposures (e.g., 5 years for the AMS B/C ratio) during which the solar activity evolves appreciably. Rather than using a static picture for the heliosphere with fixed parameters describing mean diffusion properties or tilt angle, we have adopted a quasi-steady state approach where a time-series of input parameters is converted into a time-series of modulated CR fluxes, to be subsequently averaged over the desired observation period. It is therefore essential, for a reliable modeling of the modulation effect, having time-resolved data over the period of interest. In this respect, the present work would strongly benefit from monthly resolved measurements of low energy fluxes over the period of the AMS observations. The AMS experiment is now probing the descending phase of the solar activity cycle toward the next minimum. Measurements in solar minimum conditions can be very precious for CR propagation studies, because the modulation effect during this period can be more reliable modeled. Another important point is that solar-modulation uncertainties are much larger for individual B-C flux calculations, while the effect is partially suppressed when analyzing the ratios of these fluxes. More in general, uncertainties from solar modulation are minimized with the use of secondary-to-primary ratios between species of similar charge or mass, e.g., d/, He/He, B/C, and F/Ne.
The other important source of uncertainty is that related to nuclear fragmentation processes of CRs in the ISM.
Models of CR propagation rely heavily on high-energy extrapolations of parametric formulae for total and partial reactions.
Cross-section formulae are based on low energy nuclear data (up to a few GeV/nucleon) mostly collected in the
90’s. These data are essentially untested in the energy range relevant for the current CR physics investigation.
In this work, we have performed a data-driven reevaluation of isotopic cross-section data for several PF channels involving the production of B-Be nuclei.
Along with normalization and energy scale, we have accounted for uncertainties from energy-dependent biases in cross-section formulae.
For this task we used data on charge-changing reactions.
As we have shown, nuclear-physics uncertainties are found to be a major limiting factor for the interpretation of the new B/C data from AMS.
Hence, in the analysis of CR data, these uncertainties cannot be ignored any longer.
In line with recent studies, all calling for more precise nuclear data for CRs
(Gondolo , 2014; Maurin et al. 2010, ; Moskalenko et al., 2013; Tomassetti, 2015, 2012b),
we therefore stress once more the urgent need for establishing a program
of cross-section measurements at the (100 GeV) energy scale.
To conclude, we would like to quote Chen et al. 1997 () from the Transport Collaboration:
“With the shutdown of the LBL Bevalac and the pending closure of the Saclay Saturne accelerators, opportunities for obtaining
cross-section measurements relevant to the interpretation of CR data are rapidly dwindling worldwide.
Thus, future experiments will rely heavily upon cross-section predictions, and it is important to update our formulae using data (…)
to ensure that the solutions to some astrophysical problems are not dominated by cross-section inaccuracies rather than by CR measurements”.
Nowadays, the solutions to topical problems in CR physics seems actually dominated by cross-section inaccuracy rather than by CR measurements.
I thank David Maurin for sharing the code Usine of CR propagation calculations. This work has also benefit from the public codes GALPROP (for fragmentation cross-sections) and SolarProp (for solar modulation). CR data are retrieved from online databases SSDC-DB at the Italian Space Agency and CRDB at LPSC-Grenoble (Maurin2014, ). NM data are taken from NMDB (Steigies, 2015). The author acknowledges the European Commission for support under the H2020-MSCA-IF-2015 action, grant No. 707543 MAtISSE – Multichannel Investigation of Solar Modulation Effects in Galactic Cosmic Rays.
- (1) I. A. Grenier, J. A. Black, A. W. Strong, Annu. Rev. Astron. Astrophys., 53, 199 (2015); A. W. Strong, I. V. Moskalenko, V. S. Ptuskin, Ann. Rev. of Nucl. & Part. Sci. 57, 285-327 (2007);
- Aguilar et al. (2016b) M. Aguilar, L. Ali Cavasonza, G. Ambrosi, et al., PRL 117, 231102 (2016);
- Aguilar et al. (2015a) M. Aguilar, D. Aisa, B. Alpat, et al., PRL 114, 171103 (2015); M. Aguilar, D. Aisa, B. Alpat, et al., PRL 115, 211101 (2015);
- Aguilar et al. (2016) M. Aguilar, L. Ali Cavasonza, B. Alpat, et al., PRL 117, 091103 (2016); M. Aguilar, D. Aisa, A. Alvino, et al., PRL 113, 121102 (2014);
- Serpico (2015) P. D. Serpico, Proc. 34 ICRC - The Hague, PoS(ICRC2015)009 (2015); P. Maestro, Proc. 34 ICRC - The Hague, PoS(ICRC2015)016 (2015);
- Tomassetti & Donato (2015) N. Tomassetti & F. Donato, ApJ 803, L15 (2015); N. Tomassetti, PRD 92, 081301(R) (2015);
- Amato & Blasi (2017) E. Amato & P. Blasi, ASR, in press (2017); P. Blasi, Astron. Astrophys. Rev. 21, 70 (2013);
- Blasi (2017) P. Blasi, Mon. Not. R. Astron. Soc. 471, 1662 (2017);
- Feng et al. (2016) J. Feng, N. Tomassetti, A. Oliva, PRD 94, 123007 (2016); M. J. Boschini, S. Della Torre, M. Gervasi, et al., ApJ 840, 115 (2017); M. Korsmeier & A. Cuoco, PRD 94, 123019 (2016); G. Jóhannesson, R. Ruiz de Austri, A. C. Vincent, et al., ApJ 824, 16 (2016);
- Yoon et al. (2017) S. Yoon, T. Anderson, A. Barrau, et al., ApJ 839, 1 (2017);
- Cummings et al. (2016) A. C. Cummings, E. C. Stone, B. C. Heikkila, et al., ApJ 831, 18 (2016);
- Steigies (2015) C. Steigies, Proc. 34 ICRC - The Hague, PoS(ICRC2015)225 (2015); H. Mavromichalaki, A. Papaioannou, C. Plainaki, et al., ASR 47, 2210-2222 (2011) [see also http://www.nmdb.eu];
- Adriani et al. (2013) O. Adriani, G. C. Barbarino, G. A. Bazilevskaya, et al., ApJ 765, 91 (2013);
- Kühl et al. (2016) P. Kühl, R. Gómez-Herrero, B. Heber, Sol. Phys. 291, 965-974 (2016);
- Abe et al. (2012) K. Abe, H. Fuke, S. Haino, et al., PRL 108, 051102 (2012); K. Abe, H. Fuke, S. Haino, et al., PLB 670, 103 (2008);
- (16) S. M. Read, & V. E. Viola, At. Data Nucl. Data Tables 31, 359-397 (1984);
- (17) W. R. Webber, J. C. Kish, J. M. Rockstroh, et al., ApJ 508, 949-958 (1998);
- (18) W. R. Webber, J. C. Kish, D. A. Schrier, PRC 41, 547 (1990);
- (19) D. L. Olson, B. L. Berman, D. E. Greiner, et al., PRC 28-4, 1602-1613 (1983);
- (20) P. Fontes, PRC 15-6, 2159-2168 (1977);
- (21) A. Korejwo, T. Dzikowski, M. Giller, et al., J. Phys. G: Nucl. Part. Phys. 26, 1171-1186 (2000);
- (22) A. Korejwo, M. Giller, T. Dzikowski, et al., Proc 27 ICRC - Hamburg, 2, 1371-1374 (2001);
- (23) J. R. Radin, E. Gradsztajn, A. R. Smith, PRC 20, 787 (1979);
- (24) R. Ramaty, B. Kozlovsky, R. E. Lingenfelter, H. Reeves, ApJ 488, 730-748 (1997);
- (25) W. R. Webber, A. Soutoul, J. C. Kish, et al., PRC 58, 3539 (1998);
- (26) G. M. Raisbeck & F. Yiou, PRL 27, 875 (1971);
- (27) W. R. Webber, A. Soutol, J. K. Kish, J. M. Rockstroh, ApJ 144, 153-167 (2003);
- Strong & Moskalenko (1998) A. W. Strong & I. V. Moskalenko, ApJ 509, 212-228 (1998);
- Thoudam & Horandel (2014) S. Thoudam & J. Horandel, A&A 567, A33 (2014);
- Tomassetti (2012a) N. Tomassetti, ApJ 752, L13 (2012); N. Tomassetti & F. Donato, A&A 544, A16 (2012);
- Jones et al. (2001) F. C. Jones, A. Lukasiak, V. Ptuskin, W. R. Webber, ApJ 547, 264 (2001);
- Maurin et al. (2001) D. Maurin, F. Donato, R. Taillet, P. Salati, 2001, ApJ, 555, 585-596
- Evoli et al. (2008) C. Evoli, D. Gaggero, D. Grasso, L. Maccione, JCAP, 10, 18 (2008);
- Kissmann (2014) R. Kissmann, Astropart. Phys. 55, 37-50 (2014);
- Maurin (2015) D. Maurin, Proc.34 ICRC, The Hague, PoS(ICRC2015)484 (2015);
- Ptuskin et al. (2006) V. S. Ptuskin, I. V. Moskalenko, F. C. Jones, et al., 2006, ApJ, 642, 902-916;
- (37) D. Maurin, A. Putze, L. Derome, A&A 516, A67 (2010); Y. Genolini, A. Putze, P. Salati, P. D. Serpico, A&A 580, A9 (2015);
- Bindi et al. (2017) V. Bindi, C. Corti, C. Consolandi, et al., ASR 60, 865-878 (2017);
- Potgieter (2017) M. S. Potgieter, ASR 60, 848-864 (2017); M. S. Potgieter, Living Rev. Sol. Phys., 10, 3 (2013);
- Kappl et al. (2016) R. Kappl, Comp. Phys. Comm. 207, 386-399 (2016);
- Strauss et al. (2012) R. D. Strauss, M. S. Potgieter, I. Büshing, A. Kopp, Astrophys. Space Sci. 339, 223 (2012); Bobik, P., G. Boella, M. J. Boschini, et al., ApJ 745, 132 (2012); M. S. Potgieter, ASR 53, 1415â1425 (2014); R. A. Burger & M. S. Potgieter, ApJ 339, 501-511 (1989);
- Langner et al. (2003) U. W. Langner, M. S. Potgieter, W. R. Webber, J. Geophys. Res. 108, A10, 8039 (2003); U. W. Langner & M. S. Potgieter, J. Geophys. Res. 109, A01103 (2004);
- Hoeksema (1995) J. T. Hoeksema, Space Sci. Rev. 72, 137-148 (1995);
- WSO (2017) http://wso.stanford.edu
- Tomassetti (2017) N. Tomassetti, ASR 60, 815-825 (2017);
- Ghelfi et al. (2016) A. Ghelfi, D. Maurin, A. Cheminet, et al., ASR 60, 833-847 (2017);
- Smart & Shea (2005) D. F. Smart & M. “P.” Shea, ASR 36, 2012-2020 (2005);
- Tomassetti (2015) N. Tomassetti, PRC 92, 045808 (2015);
- Maurin et al. (2015) D. Maurin, A. Cheminet, L. Derome, et al., ASR 55, 363 (2015);
- Cheminet et al. (2013) A. Cheminet, G. Hubert, V. Lacoste, et al., J. Geophys. Res. 118, 7488-7496 (2013);
- Usoskin et al. (2011) I. G. Usoskin, G. A. Bazilevskaya, G. A. Kovaltsov, J. Geophys. Res. 116, A2 (2011);
- Cholis et al. (2016) Cholis, I., Hooper, D., & Linden, T., PRD 93, 043016 (2016); C. Corti, V. Bindi, C. Consolandi, K. Whitman, ApJ 829, 8 (2016); J. Gieseler, Proc. 34 ICRC (The Hague) PoS(ICRC2015)118 (2015);
- Silberberg et al. (1998) R. Silberberg, C. H. Tsao, A. F. Barghouty, 1998, ApJ, 501, 911–919
- Webber et al. (1990b) W. R. Webber, J. C. Kish, D. A. Schrier, PRC 41, 533-546 (1990); W. R. Webber, J. C. Kish, D. A. Schrier, PRC 41, 566-571 (1990);
- Moskalenko & Mashnik (2003) I. V. Moskalenko & S. G. Mashnik, Proc 28 ICRC, Tsukuba, 2, 1969–1972 (2003);
- Moskalenko et al. (2013) I. V. Moskalenko, A. E. Vladimirov, T. Porter, A. W. Strong, Proc. 33 ICRC, Rio de Janeiro (2013); Moskalenko, I. V., Vladimirov, A. E., Porter, T. A., Proc. 32 ICRC, Beijing, 238 (2011);
- Tomassetti (2012b) N. Tomassetti, Astrophys. Space Sci. 342, 131 (2012);
- Ferrando et al. (1988) P. Ferrando, W. R. Webber, P. Goret, et al., PRC 37–4, 1490-1502 (1988);
- Allison et al. (2016) J. Allison, K. Amako, J. Apostolakis, et al., Nucl. Inst. Meth. A 835, 186-225 (2016); R. K. Tripathi, M. G. Cucinotta, J. W. Wilson, Nucl. Inst. Meth. B 117, 347-349 (1996);
- Tomassetti (2015c) N. Tomassetti, PRD 92, 063001 (2015); ApJ 815, L1 (2015);
- Genolini et al. (2017) Y. Genolini, P. Serpico, M. Boudaud, et al., arXiv:1706.09812 (2017);
- George et al. (2009) J. S. George, K. A. Lave, M. E. Wiedenbeck, et al., ApJ 698, 1666-1681 (2009);
- Kneller et al. (2003) J. P. Kneller, J. R. Phillips, T. P. Walker, ApJ 589, 217-224 (2003);
- Kappl et al. (2015) R. Kappl, A. Reinert, M. W. Winkler, ; JCAP 10, 34 (2015); Kachelrieß, M., Moskalenko, I. V, Ostapchenko, S., ApJ 803, 54 (2015); F. Donato, M. Korsmeier, M. di Mauro, PRD 96, 043007 (2017);
- Gondolo (2014) P. Gondolo, Nuclear Data Sheets 120, 175-179 (2014);
- (66) C. X. Chen, S. Albergo, Z. Caccia, et al., 1997, ApJ 479, 504;
- () V. Di Felice, C. Pizzolotto, D. D’Urso, et al., Proc. 35 ICRC - Bexco, PoS 1073 (2017) [SSDC-DB: https://tools.asdc.asi.it/CosmicRays]; D. Maurin, F. Melot, R. Taillet, A&A, 569, A32 (2014) [LPSC-DB: https://lpsc.in2p3.fr/cosmic-rays-db];