PAPER • OPEN ACCESS

Universal properties of primary and secondary cosmic ray energy spectra

Marco Smolla1,2,4, Benjamin Schäfer³, Harald Lesch² and Christian Beck³

Published 1 September 2020 • © 2020 The Author(s). Published by IOP Publishing Ltd on behalf of the Institute of Physics and Deutsche Physikalische Gesellschaft

New Journal of Physics, Volume 22, September 2020

Author e-mails: smolla@mpe.mpg.de, b.schaefer@qmul.ac.uk, c.beck@qmul.ac.uk

Author affiliations
¹ Max-Planck-Institut für extraterrestrische Physik, Giessenbachstrasse, 85741 Garching, Germany
² Ludwig-Maximilians-Universität München, Universitäts-Sternwarte München, Scheinerstrasse 1, 81679 München, Germany
³ School of Mathematical Sciences, Queen Mary University of London, London E1 4NS, United Kingdom
⁴ Author to whom any correspondence should be addressed.

Dates
Received 24 April 2020; Revised 21 July 2020; Accepted 28 July 2020

Citation
Smolla et al 2020 New J. Phys. 22 093002

Keywords
cosmic rays, high-energy physics, generalized statistical mechanics

DOI: https://doi.org/10.1088/1367-2630/abaa03

Abstract

Atomic nuclei appearing in cosmic rays (CRs) are typically classified as primary or secondary. However, a better understanding of their origin and propagation properties is still necessary. We analyse the flux of primary (He, C, O) and secondary nuclei (Li, Be, B) detected with rigidity (momentum/charge) between 2 GV and 3 TV by the alpha magnetic spectrometer on the International Space Station. We show that q-exponential distribution functions, as motivated by generalized versions of statistical mechanics with temperature fluctuations, provide excellent fits for the measured flux of all nuclei considered. Primary and secondary fluxes reveal a universal dependence on kinetic energy per nucleon for which the underlying energy distribution functions are solely distinguished by their effective degrees of freedom. All given spectra are characterized by a universal mean temperature parameter ∼200 MeV which agrees with the Hagedorn temperature. Our analysis suggests that QCD scattering processes together with nonequilibrium temperature fluctuations imprint universally onto the measured CR spectra, and produce a similar shape of energy spectra as high energy collider experiments on the Earth.

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

1. Introduction

A fundamental challenge of current cosmic ray (CR) research is to understand the origin of highly energetic CRs, their abundance in terms of different particle types, and to identify the processes at work for acceleration and propagation. Collectively these processes determine the energy dependent flux of CRs, that is their energy spectra. Because charged particles gyrate around the magnetic field lines of the interstellar medium (ISM), the directional information about the source is ultimately lost, leading to a roughly isotropic distribution observed here at Earth. The atomic nuclei among the CRs are classified as primary CRs, usually thought to be expelled by supernovae explosions and accelerated in shock fronts of supernova remnants, and secondary CRs, which result from particle collisions in the ISM. Here, we consider the flux of six different nuclei, namely the primaries He, C, O and the secondaries Li, Be, B as observed with the alpha magnetic spectrometer (AMS) on the International Space Station [1, 2].

It is commonly accepted that the major fraction of He, C, O can be classified as primary CRs whereas Li, Be, B are secondary CRs because their relative abundance exceeds the chemical composition of the ISM by a few orders of magnitude [3]. Some progress has been made in explaining CR acceleration (e.g. at supernova remnant shocks) [4] and propagation (e.g. diffusion confinement) [5] which allows to better investigate the specific processes responsible for the observed distributions. Nevertheless, considering the multitude of physical processes involved, our understanding remains incomplete and theoretical models accounting for the given nuclei spectra contain many unknown parameters and are currently under debate [6].

As measured CR energy spectra decay in good approximation with a power law over many orders of magnitudes, it is reasonable to apply a generalized statistical mechanics formalism (GSM) [7] which generates power laws rather than exponential distributions as the relevant effective canonical distributions. Canonical Boltzmann–Gibbs (BG) statistics is only valid in an equilibrium context for systems with short-range interactions, but it can be generalized to a nonequilibrium context by introducing an entropic index q, where q > 1 accounts for heavy-tailed statistics and q = 1 recovers BG statistics [8–10]. The occurrence of the index q can be naturally understood due to the fact that there are spatio-temporal temperature fluctuations in a general nonequilibrium situation, as addressed by the general concept of superstatistics, a by now standard statistical physics method [11]. Since the flux distribution as a function of energy in CRs evidently does not decay exponentially, it is reasonable not to use BG statistics but rather GSM, which has been successfully applied to CRs before in [12–14] and also applied to particle collisions in LHC experiments [15–17]. Other applications of this superstatistical nonequilibrium approach are Lagrangian [18] and defect turbulence [19], fluctuations in wind velocity and its persistence statistics [20, 21], fluctuations in the power grid frequency [22, 23] and air pollution statistics [24].

Here, we apply GSM and superstatistical methods to the observed CR flux of atomic nuclei to infer the physical parameters of the underlying energy distributions, which turn out to be nearly identical for all primaries and secondaries, respectively. The universal properties of the two CR types can be distinguished by a single parameter, the entropic index q, which we relate to the effective degrees of freedom of temperature fluctuations that are relevant in a GSM description. The average temperature parameter that fits all nuclei spectra turns out to be universal as well and is given by about 200 MeV, coinciding with the Hagedorn temperature. This suggests that QCD scattering processes play a dominant role in shaping the spectrum of observed CRs. The spectra are indeed similar to observed momentum spectra in high energy proton–proton-collider experiments on the Earth, which are known to generate q-exponential power laws [15, 17].

The paper is organized as follows: in section 2 we demonstrate that CR nuclei spectra, as measured by AMS, are well described by q-exponential distribution functions. We show that the spectra exhibit data collapse if they are related to the kinetic energy per nucleon. In section 3 we interpret the observed spectra in terms of temperature fluctuations occurring during the production process of the CRs, based on χ² superstatistics. We relate the power law spectral index to the relevant degrees of freedom contributing to the temperature fluctuations. Finally, a possible physical explanation of the universal properties of the observed spectra is given in section 4.

2. Results

We investigate the CR flux, given as differential intensity with respect to kinetic energy per nucleon, defined as E = (Etotal − m)/A, with total energy Etotal=p2+m2, momentum p=|p→|, rest mass m = Au, mass number A, atomic mass unit u = 0.931 GeV and [m] = [p] = [GeV] in c = 1 convention.

In order to infer physical parameters from the energy distribution fit to the observed CR flux, we employ an established GSM model [13], modified slightly by replacing the total energy by the kinetic energy per nucleon E. This choice of variable is common practice in CR literature because the kinetic energy per nucleon and the charge to mass ratio of a given particle provides the essential properties that decide how the particle’s trajectory will be modified by the presence of magnetic field lines. In the appendix we rigorously derive the distribution function

PE(E)=Cρ(E)eq−bE,(1)

which corresponds to the following differential intensity of flux

JE(E)=v(E)PE,(2)

where C, q > 1 and T = b⁻¹ > 0 are free parameters and the q-exponential is defined as eqx≡1+(1−q)x11−q which implies eq→1x=ex. ρ(E) is a phase space factor which describes the density of states, i.e. how many energy states can be taken on in a given range. For our fits we used ρ(E)=p2dpdE=(E+u)E(E+2u) which leads to the flux derived from our superstatistical model

JEmod(E)=CE(E+2u)eq(−bE),(3)

which we compare to the observed flux JEobs.

The entropic index q determines the high-energy (i.e. the tail) behaviour of the distribution since the q-exponential asymptotically approaches a power law

limE→∞eq−bE∝Eγ,(4)

with spectral index γ = 2 − 1/(q − 1) for q > 1. The parameter T = b⁻¹ represents a temperature in energy units that constrains the low-energy regime of maximum flux. Since our analysis focuses on the spectral shape of the energy distribution we collect all global factors, which do not depend explicitly on the energy, in the amplitude C, which is merely a gauge for the absolute magnitude of the flux.

Figure 1 illustrates that most data points are fitted by our model within a single standard deviation for all six nuclei. We determined the best fit by applying χ² minimisation with (JEmod−JEobs)/σ, meaning deviation of model from data weighted by the respective measurement uncertainty, where the standard deviation σ is the sum of measurement errors for a specific energy bin. For most of the data the error is of the order of a few percent whereas the uncertainty tends to increase with energy up to the largest uncertainty of 89 percent associated with the beryllium flux measured in the highest energy bin.

Figure 1. The particle flux of each CR species was fitted with (3) using three parameters C, T, q. The vertical axis in this log–log plot was multiplied with E2.7 for better visibility. The fit’s accuracy can be quantified by the deviation from modelled Jmod to observed flux Jobs weighted by the respective measurement error σ. Evidently, almost all data points fall within the uncertainty range of ±σ illustrated as grey shaded area. The mean temperature T₀ is defined by (9). The amplitude C has dimensions [C] = [m⁻² sr⁻¹ s⁻¹ GeV⁻³].

Figure 2 reveals the universal properties of the primary (He, C, O) and secondary (Li, Be, B) CR fluxes when rescaling each nuclei flux with a suitable global factor such that all data points collapse to a single line in the low energy range. Fixing the global amplitude parameter to C = 1 and T = 0.240 GeV, which is the average value for the temperatures inferred from the individual best fits in figure 1, allows us to do a best fit with q as the only free parameter for the collapsed data of primaries and secondaries. This yields qprim = 1.2109 (n = 3.5) and qsec = 1.1969 (n = 4.2), where n can be interpreted as degrees of freedom of temperature fluctuations as outlined below.

Figure 2. Each particle flux was rescaled with a suitable factor such that the data points (roughly) collapse to a single line at the low energy end and the universal properties of primary and secondary CR nuclei spectra become visible. For larger energies the spectrum splits into primaries and secondaries which can be distinguished by a single parameter, the entropic index q which can be interpreted by the underlying effective degrees of freedom.

3. Interpretation in terms of temperature fluctuations

We consider the observed CR spectra to be the result of many different high energy scattering processes, each having a different local temperature β⁻¹ in the local scattering volume. This idea was previously worked out in detail for collider experiments using LHC data, e.g. in [15]. There are strong fluctuations of temperature in each scattering event, which can be described by superstatistics [11], a standard method in the theory of complex systems. For CRs, we need to generate asymptotic power laws and this can be achieved by so-called χ² superstatistics. As is generally known (see, e.g. [13]) the probability density function for a fluctuating β of the form

β=∑i=1nXi2(5)

with independent and identically distributed Gaussian random variables Xi is a χ² distribution given by

g(β)=1Γ(n/2)n2β0n2βn2−1exp−nβ2β0.(6)

It is well-known in the formalism of superstatistics that superimposing various subsystems with different temperature weighted with g(β) leads to q-exponential statistics. For each scattering event we apply ordinary statistical mechanics locally, i.e. the conditional probability density of a kinetic energy state E in a given scattering event for a given temperature is

pE(E|β)=1Z(β)ρ(E)e−βE.(7)

In order to normalize our conditional distribution function we need to integrate over all possible energy states, obtaining the normalization constant Z(β)=∫0∞ρ(E)e−βEdE. The marginal distribution PE(E) (the unconditioned distribution of energies) can be computed by integrating the conditional distribution pE(E|β) over all inverse temperatures β weighted with g(β). In the relativistic limit (neglecting mass terms) this yields

PE(E)=∫0∞g(β)pE(E|β)dβ∼ρ(E)eq−bE(8)

with b = β₀/(4 − 3q) and mean inverse temperature

β0=∫0∞βg(β)dβ≔1T0.(9)

The effective degrees of freedom n are related to the entropic index q via

n=2q−1−6.(10)

Considering the physical meaning of the random variable Xi defined in equation (5) there are different possible interpretations, see [13, 18, 25, 26] for details. The main conclusions in our given analysis are independent of the particular interpretation chosen but it is worth emphasizing that it is physically plausible to understand Xi as a measure for the fluctuating effective energy dissipation, with n = 3 representing the three spatial degrees of freedom as minimum value, which is increased to n = 4 when including variations in time.

In the appendix we provide a more detailed derivation of the above results, which show that q-exponential statistics follows naturally from summing up ordinary Boltzmann distributions with χ²-distributed inverse temperatures.

4. Physical interpretation and possible reason for universality

For the temperature parameter T0=β0−1, defined in (9), we get the value T₀ ∼ 600 MeV for each of the six CR species in our fits. Hence, the average effective temperature per quark is of the order T₀/3 ∼ 200 MeV, i.e. we recover approximately the observed value of the Hagedorn temperature which is roughly known to be in the range 140 to 200 MeV [27–29] and represents a universal critical temperature for the quark gluon plasma and for high energy QCD scattering processes. Remarkably, the fitted value of T₀/3 in the fits is observed to be the same for all six nuclei, i.e. for both primary and secondary CRs within a range of about one tenth of its absolute value.

Let us provide some arguments on why we consider T₀/3 as a relevant temperature parameter. In general, there are two alternative formulations of superstatistics, defined as type A and type B in [11], which yield the same form of distribution functions but differ in their definitions for T₀ because type A uses the unnormalized whereas type B uses the normalized Boltzmann factor for deriving the generalized canonical distribution, as we present in detail for type B in the appendix . We consider type B superstatistics as physically more plausible because we understand the generalized canonical distribution as originating from a superposition of many CR ensembles associated with a normalized canonical distribution respectively. Since the Hagedorn temperature is associated with the kinetic energy of particles interacting via the strong nuclear force, we divide the average temperature T₀ by the number of quarks, namely three for atomic nuclei. For a quark–gluon plasma one can either define a temperature for single quarks, or—after hadronization—for mesonic or baryonic states. At the critical Hagedorn phase transition point, where both states exist, this is mainly a question of definition [30].

The emergence of the Hagedorn temperature (at least as an order of magnitude) in our fits suggests that CR energy spectra might originate from high energy scattering processes taking place at the Hagedorn temperature TH. Very young neutron stars, initially formed in a supernova explosion, indeed have a temperature ∼10¹² K to 100 MeV of comparable order of magnitude as the Hagedorn temperature [31]. As the Hagedorn temperature is universal, so is the average kinetic energy per quark of the CRs nuclei, assuming they are produced in a Hagedorn fireball, either during the original supernova explosion, or later in collision processes of highly energetic CR particles with the ISM.

Our observation that the kinetic energy per nucleon (or per quark) yields universal behaviour of the spectra is indeed pointing towards QCD processes as the dominant contribution that shapes the spectra (see also [14]): were there mainly electromagnetic processes underlying the spectra, one would expect invariance under rescaling with Z, but we observe invariance (universality) under rescaling with A. At the LHC one observes similar q-exponentials for the measured transverse momentum spectra, as generated by QCD scattering processes, with a temperature parameter b⁻¹ that is of similar order of magnitude (150 MeV) as in our fits for the CRs, see table 4 in [17]. That paper also shows that hard parton QCD scattering leads to power law spectral behaviour.

Note that while the entire energy spectrum in figure 1 is well fitted by a q-exponential, the residuals tend to oscillate. A similar oscillatory behaviour of the residuals (logarithmically depending on the energy) has been observed in the transverse momentum distribution for high energy pp collision experiments at the LHC [17, 32]. The similarity of these log-periodic oscillations for our CR data and for collider experiments on the Earth is indeed striking, and once again supports our point that both phenomena could have similar roots based on high energy scattering processes.

After having analysed the average temperature, let us now concentrate on the fluctuations of temperature in the individual scattering events, described by the parameter n, which determines the entropic index q and thereby the tail behaviour. One readily notices that in our GSM model the marginal distributions PE(E) decay asymptotically as

PE(E)∼E−1−n2.(11)

In order to calculate the expectation of the fluctuating energy E,

⟨E⟩=∫0∞EPE(E)dE,(12)

one notices that the integrand decays as EP(E) ∼ E−n/2. Thus the expectation value is only well defined if n > 2.

In the absence of further effects, like an energy dependent cross section, we could explain n exclusively by the underlying statistics and thus associate n with the number of Gaussian random variables contributing to the fluctuating β in equation (5). Since for CR propagation energy dependent processes affect the spectral shape the derived value for n will represent both the statistical properties and the energy dependent processes. For this reason effective non-integer values for n are possible. Because the above argument about the existence of the expectation value should apply more generally, and thus even in the absence of additional spectral modifications, we conclude that n = 3 is the minimum value for the degrees of freedom.

A similar argument applies if one looks at the existence of the mean of the temperature β⁻¹ as formed with the probability density g(β)

⟨β−1⟩=∫0∞β−1g(β)dβ.(13)

The above mean only exists for n > 2, since the integrand behaves as βn2−2 for β → 0. We are thus naturally led to the minimum value n = 3 as the strongest fluctuation state of the Hagedorn fireball for which a mean energy ⟨E⟩ and a mean temperature ⟨β⁻¹⟩ is well-defined.

For secondary CRs, there is an additional degree of freedom as an additional collision process at a later time is needed to produce secondary CRs. Thus it is plausible that for secondary CRs n is larger than the minimum value n = 3. The next higher value of n, which can be regarded as an excited state of temperature fluctuations, n = 4, corresponds to q = 1.2000. Indeed, based on our fits (see figure 1), secondary CRs are well approximated by this q and therefore n = 4.

In the experimental data detected by AMS, it is to be expected that we will not observe the exact values of q = 1.2222 and q = 1.2000 since the spectra are modified by diffusion processes in the galaxy, by energy dependent escape processes from the shock front of the accelerating supernova remnant, and by radiative losses from acceleration. All these effects can alter the spectrum and lead to small changes in the optimum fitting parameters q and T. We think this is the reason why the best fits of the observed spectra correspond to n = 3.5 rather than n = 3, and n = 4.2 rather than n = 4, equivalent to minor negative corrections for the spectral index γ of the order Δγ ≈ −0.1. Also, the effective temperature T may be increased by diffusion processes in the galaxy, which will broaden the distributions. However, it seems these effects are only small perturbations that slightly modify the universal parameters set by the QCD scattering processes.

While the connection of QCD and GSM was emphasized in [17], the model that we implement in our paper is mainly based on a nonequilibrium statistical mechanics approach, as originally introduced in [13]. This approach is based on temperature fluctuations in each (small) interaction volume, where the scattering event takes place. These local temperature fluctuations are at the root of the observed q-exponentials and the associated temperature scale turns out to coincide approximately with the Hagedorn temperature known from QCD. Other authors [33] have emphasized the fractal and hierarchical structure of scattering events and hadronization cascades, or the complexity of long-range interactions in the hadronization process [30, 34] coming to similar conclusions.

5. Discussion on relevance of solar wind modulation

The AMS measurements were taken at about 400 km above Earth’s surface and are thus subject to solar wind modulation which yields a suppressed flux compared to outside the heliosphere, in particular for charged particles with kinetic energies per nucleon below ≲10 GeV [38]. Thus for our given AMS data with kinetic energies per nucleon in the range of 0.4 GeV ≲ E ≲ 1.2 TeV we would like to quantify the effect of solar wind modulation on our given spectra. Using CR propagation models allows to infer the unmodulated flux before CRs are entering the heliosphere, that is the LIS flux.

This was recently done by [35–37] who combined the two CR propagation models HelMod and Galprop and published the calculated flux for all our given atomic nuclei and for the entire energy range covered by AMS. We use their data and interpolate it to match the AMS energy bins definition. We find a maximum deviation of LIS to flux inside heliosphere (AMS) for the lowest energetic particles with LIS/AMS ≲ 4 for E = 0.4 GeV. The two spectra converge for larger energies quickly and we find LIS/AMS ≲ 1.2 for E = 10 GeV such that in fact only the lowest energy range of our spectra is significantly affected. Since the propagation model provides the flux without giving any uncertainty, we assign each estimated flux the same relative error as in the AMS data set. This makes the comparison between the flux inside and outside the heliosphere consistent and allows to put appropriate weight on measurements with smaller uncertainties for our least-square optimization.

Analogously to the steps performed for the given AMS data, we apply our GSM methodology to the LIS data and present the resulting fits and parameters in figure 3. The average temperatures T₀/3 for the different nuclei are about 50 to 80 MeV lower than for the AMS data, namely in the range 129 to 152 MeV. Still, these temperatures are all about the scale of the Hagedorn temperature and in fact coincide with the temperature range 130 to 160 MeV inferred by GSM methods applied to LHC experiments found by [17, 32]. The effective degrees of freedom remain approximately the same. Since the reliability of our methodology ultimately depends on having a large energy range measured for all the different nuclei, the AMS data is the best currently available experimental data set. In contrast, measurements acquired by Voyager outside the heliosphere only cover energies from about 3 MeV to a few hundred MeV [39, 40]. Hence, we apply our analysis to the large range of AMS-measured data and estimate the modulation by the solar wind, rather than using theoretically derived data for unmodulated spectra.

Figure 3. Superstatistical results are robust when correcting for heliospheric impact. Using CR propagation models HelMod and Galprop allows [35–37] to estimate flux outside the heliosphere, that is unmodulated by solar wind representing the local interstellar spectra, in short local interstellar (LIS) flux. Here, we use the data published in [35–37] which we investigate for the given AMS energy bins. We apply JEmod(E)=CE(E+2u)eq(−bE) with eqx≡1+(1−q)x11−q in order to derive the best fit global amplitude C, temperature T and entropic index q. From the best fit parameter T we derive the average temperature per quark as T₀/3. The entropic index q can be translated into effective degrees of freedom n and into the spectral index γ representing the asymptotic power law behaviour limE→∞eq−bE∝Eγ.

6. Conclusion

We provide excellent fits for the measured AMS spectra of primary (He, C, O) and secondary CRs (Li, Be, B) using a simple superstatistical model. The observed q-exponential spectra are interpreted in terms of temperature fluctuations occurring in the Hagedorn fireball during the production process of CRs in their individual scattering events. We provide evidence that the observed spectra of CR nuclei share universal properties: the spectra collapse if the kinetic energy per nucleon is taken as the relevant variable. Primary and secondary CRs can be uniquely distinguished by their respective entropic index q, corresponding to different degrees of freedom associated with the temperature fluctuations. They share the same average temperature parameter, whose order of magnitude coincides with the Hagedorn temperature.

Acknowledgments

This project has received funding from the European Union’s Horizon 2020 Research and Innovation Programme under the Marie Sklodowska-Curie Grant Agreement No. 840825. In addition, we acknowledge support from the DFG Cluster of Excellence ‘ORIGINS’.

Appendix A. Deriving the superstatistical distribution function

We derive the distribution function (1), that is PE=Cρ(E)eq−bE, using the framework of superstatistics by which we can interpret the best fit parameters with a temperature T = b⁻¹ and effective degrees of freedom n.

Superstatistics [11] is a generalization of Boltzmann statistics in the sense that the distribution function can be derived by integrating the conditional probability distribution pE(E|β) = ρ(E)e−βE/Z(β) for all given values of inverse temperature β. The normalization is calculated by summing over all possible energy states, yielding Z(β)=∫0∞ρ(E)e−βEdE. In agreement with [13] we apply the ultra-relativistic approximation for the density of states ρ(E) ∼ E² in order to calculate Z(β) ∼ β⁻³. Given the χ²-distributed β, defined by (6), we calculate the generalized canonical distribution as follows:

PE(E)=∫0∞g(β)pE(E|β)dβ(A.1)

∼n2β0n2ρ(E)∫0∞βn2+2e−βE+n2β0dβ(A.2)

∼ρ(E)n2β0n2E+n2β0−3−n2.(A.3)

Introducing q = 1 + 2/(n + 6) (equivalent to n/2 = 1/(q − 1) − 3) and b = β₀/(4 − 3q), allows us to express the result as:

PE(E)∼ρ(E)E+n2β0−3n2β0E+n2β01q−1−3(A.4)

∼ρ(E)n2β0−311+E2β0n1q−1(A.5)

∼ρ(E)eq−bE.(A.6)

Thus we have derived the distribution function (1), which we used for our fits, building on the framework of GSM and superstatistics.

Note that the above equations are only valid for the particular case ρ(E) ∼ E² and g(β) being a χ² distribution. More generally, one has

PE(E)∼ρ(E)∫0∞g(β)Z(β)e−βEdβ.(A.7)

Appendix B. Applying theory to observation

We provide a thorough derivation of equation (2), that is JE = v(E)PE, which relates the distribution function from our superstatistical model with the observed differential flux intensity measured by AMS.

The AMS data [1, 2] was published in bins of rigidity R = pc/Ze with atomic number Z, electric charge e, momentum p=|p→|, [R] = [V] and the corresponding flux measured in units [J(R)] = [m⁻² sr⁻¹ s⁻¹ GV⁻¹]. Instead of rigidity we have chosen to investigate the spectrum in respect to kinetic energy per nucleon. To convert the flux dependence from rigidity R to kinetic energy per nucleon E, we need to transform the flux JR(R) → JE(E) such that JR(R)dR = JE(E)dE is conserved. This is a simple transformation of variables and yields

JE(E)=AZeE+uE(E+2u)JR(R),(B.1)

with [JE(E)] = [m⁻² sr⁻¹ s⁻¹ GeV⁻¹]. For better visibility of the accuracy of our fits, we multiplied the flux with E2.7, such that the units for the flux in the presented plots are [GeV1.7 m⁻² sr⁻¹ s⁻¹]. For the atomic number A we refer to AMS [1, 2] who inferred the following average abundance of isotopes ⁴He, ¹²C, ¹⁶O, 6.5Li, ⁸Be and 10.7Be among the detected nuclei. The measured flux J represents a differential intensity. Thus it counts the number of particles with energy E (or rigidity R) coming from a unit solid angle that pass through a unit surface per unit of time.

Our superstatistical model builds on a distribution function, denoted as P, which counts the spatial density of particles within a given momentum/energy range as

dNd3x∼PE(E)dE∼Pp(E)d3p.(B.2)

Analogously to PE(E)∼ρ(E)eq−bE in the previous section one can derive that Pp(E)∼eq−bE. Thus the density of states ρ(E) can be calculated from the conservation condition (B.2). Using E=(p2+m2−m)/A, which implies that the energy depends only on the magnitude of the momentum, simplifies d³p = 4πp²dp, and therefore

PE(E)∼ρ(E)eq−bE∼p2dpdEeq−bE.(B.3)

Calculating the derivative and using p² = A²E(E + 2u) we find

ρ(E)∼(E+u)E(E+2u).(B.4)

Note that we generally neglect constant global factors in our equations because we are focussing on the shape of the spectrum rather than its absolute magnitude. Evidently, [p²Pp] = [PE] = [eV⁻¹ m⁻³] does not have the same dimension as the detected flux, given as differential intensity J with [J] = [eV⁻¹ m⁻² s⁻¹ sr⁻¹]. This reminds us that in order to derive the associated differential intensity from a distribution function we have to account for the rate at which particles go through the detector. That is we multiply with the particle’s velocity to obtain the flux JE, corresponding to the distribution function PE, which yields

JEmod(E)∼v(E)PE(E)∼v(E)ρ(E)eq−bE.(B.5)

[41] provides a detailed overview about the different ways to count particles including this relation. Evidently, it yields the desired physical dimensions since [vPE] = [eV m⁻² s⁻¹].

In order to express the velocity in terms of E we use p = γmv with γ=11−v (in c = 1 convention), p=AE(E+2u) and m = Au to find

v(E)=E(E+2u)(E+u).(B.6)

Plugging everything into (B.5) reveals the relation between q-exponential distribution function and the observed differential intensity

JEmod(E)=CE(E+2u)eq−bE(B.7)

which recovers the function we fitted to the data (3).

References

[1]

Aguilar M et al 2017 Observation of the identical rigidity dependence of He, C, and O cosmic rays at high rigidities by the alpha magnetic spectrometer on the international space station Phys. Rev. Lett. 119 251101
Crossref

[2]

Aguilar M et al 2018 Observation of new properties of secondary cosmic rays lithium, beryllium, and boron by the alpha magnetic spectrometer on the International Space Station Phys. Rev. Lett. 120 021101
Crossref

[3]

Gaisser T K, Engel R and Resconi E 2016 Cosmic Rays and Particle Physics (Cambridge: Cambridge University Press)

[4]

Tatischeff V and Gabici S 2018 Particle acceleration by supernova shocks and spallogenic nucleosynthesis of light elements Annu. Rev. Nucl. Part. Sci. 68 377–404
Crossref

[5]

Strong A W, Moskalenko I V and Ptuskin V S 2007 Cosmic-ray propagation and interactions in the galaxy Annu. Rev. Nucl. Part. Sci. 57 285–327
Crossref

[6]

Gabici S, Evoli C, Gaggero D, Lipari P, Mertsch P, Orlando E, Strong A and Vittino A 2019 The origin of galactic cosmic rays: challenges to the standard paradigm Int. J. Mod. Phys. D 28 1930022
Crossref

[7]

Wilk G and Włodarczyk Z 2010 Nonextensive thermal sources of cosmic rays Cent. Eur. J. Phys. 8 726–36
Crossref

[8]

Tsallis C 2009 Introduction to Nonextensive Statistical Mechanics (Berlin: Springer)

[9]

Beck C 2009 Generalized information and entropy measures in physics Contemp. Phys. 50 495
Crossref

[10]

Jizba P and Korbel J 2019 Maximum entropy principle in statistical inference: case for non-Shannonian entropies Phys. Rev. Lett. 122 120601
Crossref

[11]

Beck C and Cohen E G D 2003 Superstatistics Physica A 322 267–75
Crossref

[12]

Tsallis C, Anjos J C and Borges E P 2003 Fluxes of cosmic rays: a delicately balanced stationary state Phys. Lett. A 310 372
Crossref

[13]

Beck C 2004 Generalized statistical mechanics of cosmic rays Physica A 331 173–81
Crossref

[14]

Cigdem Yalcin G and Beck C 2018 Generalized statistical mechanics of cosmic rays: application to positron-electron spectral indices Sci. Rep. 8 1764
Crossref

[15]

Beck C 2009 Superstatistics in high-energy physics Eur. Phys. J A 40 267
Crossref

[16]

Marques L, Andrade E and Deppman A 2013 Nonextensivity of hadronic systems Phys. Rev. D 87 114022
Crossref

[17]

Wong C-Y, Wilk G, Cirto L J L and Tsallis C 2015 From qcd-based hard-scattering to nonextensive statistical mechanical descriptions of transverse momentum spectra in high-energy pp andpp¯ collisions Phys. Rev. D 91 114027
Crossref

[18]

Beck C 2007 Statistics of three-dimensional lagrangian turbulence Phys. Rev. Lett. 98 064502
Crossref

[19]

Daniels K E, Beck C and Bodenschatz E 2004 Defect turbulence and generalized statistical mechanics Physica D 193 208–17
Crossref

[20]

Rizzo S and Rapisarda A 2004 Environmental atmospheric turbulence at florence airport AIP Conf. Proc. 742 176–81
Crossref

[21]

Weber J, Reyers M, Beck C, Timme M, Pinto J G, Witthaut D and Schäfer B 2019 Wind power persistence is characterized by superstatistics Sci. Rep. 9 19971
Crossref

[22]

Schäfer B, Beck C, Aihara K, Witthaut D and Timme M 2018 Non-Gaussian power grid frequency fluctuations characterized by lévy-stable laws and superstatistics Nat. Energy 3 119
Crossref

[23]

Gorjão L R, Anvari M, Kantz H, Beck C, Witthaut D, Timme M and Schäfer B 2020 Data-driven model of the power-grid frequency dynamics IEEE Access 8 43082–97
Crossref

[24]

Williams G, Schäfer B and Beck C 2020 Superstatistical approach to air pollution statistics Phys. Rev. Res. 2 013019
Crossref

[25]

Beck C 2001 Dynamical foundations of nonextensive statistical mechanics Phys. Rev. Lett 87 180601
Crossref

[26]

Gravanis E, Akylas E and Livadiotis G 2020 Physical meaning of temperature in superstatistics Europhys. Lett. 130 30005
Crossref

[27]

Hagedorn R 1965 Statistical thermodynamics of strong interactions at high energies Suppl. Nuovo Cimento 3 147–86

[28]

Rafelski J and Ericson T 2016 The Tale of the Hagedorn Temperature BT–Melting Hadrons Boiling Quarks–from Hagedorn Temperature to Ultra-relativistic Heavy-Ion Collisions at CERN: With a Tribute to Rolf Hagedorn (Berlin: Springer) pp 41–8

[29]

Broniowski W, Florkowski W and Leonid Y 2004 Update of the Hagedorn mass spectrum Phys. Rev. D 70 117503
Crossref

[30]

Beck C 2000 Non-extensive statistical mechanics and particle spectra in elementary interactions Physica A 286 164–80
Crossref

[31]

Lattimer J M 2015 Introduction to neutron stars AIP Conf. Proc. 1645 61–78
Crossref

[32]

Wilk G and Włodarczyk Z 2017 Oscillations in multiparticle production processes Entropy 19 670
Crossref

[33]

Deppman A 2016 Thermodynamics with fractal structure, Tsallis statistics, and hadrons Phys. Rev. D 93 54001
Crossref

[34]

Bediaga I, Curado E M F and de Miranda J M 2000 A nonextensive thermodynamical equilibrium approach in e⁺e⁻ → hadrons Physica A 286 156–63
Crossref

[35]

Boschini M J et al 2017 Solution of heliospheric propagation: unveiling the local interstellar spectra of cosmic-ray species Astrophys. J. 840 115
Crossref

[36]

Boschini M J et al 2018 Deciphering the local interstellar spectra of primary cosmic-ray species with HelMod Astrophys. J. 858 61
Crossref

[37]

Boschini M J et al 2020 Deciphering the local interstellar spectra of secondary nuclei with the Galprop/Helmod framework and a hint for primary lithium in cosmic rays Astrophys. J 889 167
Crossref

[38]

Moraal H 2014 Cosmic rays in the heliosphere: observations Astropart. Phys. 53 175–85
Crossref

[39]

Cummings A C, Stone E C, Heikkila B C, Lal N, Webber W R, Jóhannesson G, Moskalenko I V, Orlando E and Porter T A 2016 Galactic cosmic rays in the local interstellar medium: Voyager 1 observations and model results Astrophys. J. 831 18
Crossref

[40]

Stone E C, Cummings A C, Heikkila B C and Lal. N 2019 Cosmic ray measurements from Voyager 2 as it crossed into interstellar space Nat. Astron. 3 1013–8
Crossref

[41]

Moraal H 2013 Cosmic-ray modulation equations Space Sci. Rev. 176 299–319
Crossref