Universal properties of primary and secondary cosmic ray energy spectra

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.

We investigate the CR flux, given as differential intensity with respect to kinetic energy per nucleon, defined as E = (E_total − m)/A, with total energy ${E}_{\text{total}}=\sqrt{{p}^{2}+{m}^{2}}$, momentum $p=\vert \overrightarrow {p}\vert$, 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

$$\begin{equation}{P}_{E}\left(E\right)=C\rho \left(E\right){\mathrm{e}}_{q}^{-bE},\end{equation} \tag{ 1 }$$

which corresponds to the following differential intensity of flux

$$\begin{equation}{J}_{E}\left(E\right)=v\left(E\right){P}_{E},\end{equation} \tag{ 2 }$$

where C, q > 1 and T = b⁻¹ > 0 are free parameters and the q-exponential is defined as ${\mathrm{e}}_{q}^{x}\equiv {\left(1+\left(1-q\right)x\right)}^{\frac{1}{1-q}}$ which implies ${\mathrm{e}}_{q\to 1}^{x}={\mathrm{\text{e}}}^{x}$. ρ(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 $\rho \left(E\right)={p}^{2}\frac{\mathrm{d}p}{\mathrm{d}E}=\left(E+u\right)\sqrt{E\left(E+2u\right)}$ which leads to the flux derived from our superstatistical model

$$\begin{equation}{J}_{E}^{\mathrm{m}\mathrm{o}\mathrm{d}}\left(E\right)=CE\left(E+2u\right){\mathrm{e}}_{q}\left(-bE\right),\end{equation} \tag{ 3 }$$

which we compare to the observed flux ${J}_{E}^{\text{obs}}$.

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

$$\begin{equation}\underset{E\to \infty }{\mathrm{lim}}\enspace {\mathrm{e}}_{q}^{-bE}\propto {E}^{\gamma },\end{equation} \tag{ 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 $\left({J}_{E}^{\mathrm{m}\mathrm{o}\mathrm{d}}-{J}_{E}^{\text{obs}}\right)/\sigma$, 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.

Zoom In Zoom Out Reset image size

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 q_prim = 1.2109 (n = 3.5) and q_sec = 1.1969 (n = 4.2), where n can be interpreted as degrees of freedom of temperature fluctuations as outlined below.

Zoom In Zoom Out Reset image size

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

$$\begin{equation}\beta =\sum _{i=1}^{n}{X}_{i}^{2}\end{equation} \tag{ 5 }$$

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

$$\begin{equation}g\left(\beta \right)=\frac{1}{{\Gamma}\left(n/2\right)}{\left(\frac{n}{2{\beta }_{0}}\right)}^{\frac{n}{2}}{\beta }^{\frac{n}{2}-1}\enspace \mathrm{exp}\left(-\frac{n\beta }{2{\beta }_{0}}\right).\end{equation} \tag{ 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

$$\begin{equation}{p}_{E}\left(E\vert \beta \right)=\frac{1}{Z\left(\beta \right)}\rho \left(E\right){\text{e}}^{-\beta E}.\end{equation} \tag{ 7 }$$

In order to normalize our conditional distribution function we need to integrate over all possible energy states, obtaining the normalization constant $Z\left(\beta \right)={\int }_{0}^{\infty }\rho \left(E\right){\text{e}}^{-\beta E}\mathrm{d}E$. The marginal distribution P_E(E) (the unconditioned distribution of energies) can be computed by integrating the conditional distribution p_E(E|β) over all inverse temperatures β weighted with g(β). In the relativistic limit (neglecting mass terms) this yields

$$\begin{equation}{P}_{E}\left(E\right)={\int }_{0}^{\infty }g\left(\beta \right){p}_{E}\left(E\vert \beta \right)\mathrm{d}\beta \sim \rho \left(E\right){\mathrm{e}}_{q}^{-bE}\end{equation} \tag{ 8 }$$

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

$$\begin{equation}{\beta }_{0}={\int }_{0}^{\infty }\beta g\left(\beta \right)\mathrm{d}\beta {:=}\frac{1}{{T}_{0}}.\end{equation} \tag{ 9 }$$

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

$$\begin{equation}n=\frac{2}{q-1}-6.\end{equation} \tag{ 10 }$$

Considering the physical meaning of the random variable X_i 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 X_i 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.

For the temperature parameter ${T}_{0}={\beta }_{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 T_H. 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 P_E(E) decay asymptotically as

$$\begin{equation}{P}_{E}\left(E\right)\sim {E}^{-1-\frac{n}{2}}.\end{equation} \tag{ 11 }$$

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

$$\begin{equation}\langle E\rangle ={\int }_{0}^{\infty }E{P}_{E}\left(E\right)\mathrm{d}E,\end{equation} \tag{ 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(β)

$$\begin{equation}\langle {\beta }^{-1}\rangle ={\int }_{0}^{\infty }{\beta }^{-1}g\left(\beta \right)\mathrm{d}\beta .\end{equation} \tag{ 13 }$$

The above mean only exists for n > 2, since the integrand behaves as ${\beta }^{\frac{n}{2}-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.

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.

Zoom In Zoom Out Reset image size

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.

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'.

We derive the distribution function (1), that is ${P}_{E}=C\rho \left(E\right){\mathrm{e}}_{q}^{-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 p_E(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\left(\beta \right)={\int }_{0}^{\infty }\rho \left(E\right){\text{e}}^{-\beta E}\mathrm{d}E$. 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:

$$\begin{equation}{P}_{E}\left(E\right)={\int }_{0}^{\infty }g\left(\beta \right){p}_{E}\left(E\vert \beta \right)\mathrm{d}\beta \end{equation} \tag{ A.1 }$$

$$\begin{equation}\sim {\left(\frac{n}{2{\beta }_{0}}\right)}^{\frac{n}{2}}\rho \left(E\right){\int }_{0}^{\infty }{\beta }^{\frac{n}{2}+2}\enspace {\text{e}}^{-\beta \left(E+\frac{n}{2{\beta }_{0}}\right)}\mathrm{d}\beta \end{equation} \tag{ A.2 }$$

$$\begin{equation}\sim \rho \left(E\right){\left(\frac{n}{2{\beta }_{0}}\right)}^{\frac{n}{2}}{\left(E+\frac{n}{2{\beta }_{0}}\right)}^{-3-\frac{n}{2}}.\end{equation} \tag{ 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:

$$\begin{equation}{P}_{E}\left(E\right)\sim \rho \left(E\right){\left(E+\frac{n}{2{\beta }_{0}}\right)}^{-3}{\left(\frac{\frac{n}{2{\beta }_{0}}}{E+\frac{n}{2{\beta }_{0}}}\right)}^{\frac{1}{q-1}-3}\end{equation} \tag{ A.4 }$$

$$\begin{equation}\sim \rho \left(E\right){\left(\frac{n}{2{\beta }_{0}}\right)}^{-3}{\left(\frac{1}{1+E\frac{2{\beta }_{0}}{n}}\right)}^{\frac{1}{q-1}}\end{equation} \tag{ A.5 }$$

$$\begin{equation}\sim \rho \left(E\right){\mathrm{e}}_{q}^{-bE}.\end{equation} \tag{ 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

$$\begin{equation}{P}_{E}\left(E\right)\sim \rho \left(E\right){\int }_{0}^{\infty }\frac{g\left(\beta \right)}{Z\left(\beta \right)}\enspace {\text{e}}^{-\beta E}\mathrm{d}\beta .\end{equation} \tag{ A.7 }$$

We provide a thorough derivation of equation (2), that is J_E = v(E)P_E, 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=\vert \overrightarrow {p}\vert$, [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 J_R(R) → J_E(E) such that J_R(R)dR = J_E(E)dE is conserved. This is a simple transformation of variables and yields

$$\begin{equation}{J}_{E}\left(E\right)=\frac{A}{Ze}\frac{E+u}{\sqrt{E\left(E+2u\right)}}{J}_{R}\left(R\right),\end{equation} \tag{ B.1 }$$

with [J_E(E)] = [m⁻² sr⁻¹ s⁻¹ GeV⁻¹]. For better visibility of the accuracy of our fits, we multiplied the flux with E^2.7, such that the units for the flux in the presented plots are [GeV^1.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

$$\begin{equation}\frac{\mathrm{d}N}{{\mathrm{d}}^{3}x}\sim {P}_{E}\left(E\right)\mathrm{d}E\sim {P}_{p}\left(E\right){\mathrm{d}}^{3}p.\end{equation} \tag{ B.2 }$$

Analogously to ${P}_{E}\left(E\right)\sim \rho \left(E\right){e}_{q}^{-bE}$ in the previous section one can derive that ${P}_{p}\left(E\right)\sim {\mathrm{e}}_{q}^{-bE}$. Thus the density of states ρ(E) can be calculated from the conservation condition (B.2). Using $E=\left(\sqrt{{p}^{2}+{m}^{2}}-m\right)/A$, which implies that the energy depends only on the magnitude of the momentum, simplifies d³p = 4πp²dp, and therefore

$$\begin{equation}{P}_{E}\left(E\right)\sim \rho \left(E\right){\mathrm{e}}_{q}^{-bE}\sim {p}^{2}\frac{\mathrm{d}p}{\mathrm{d}E}\enspace {\mathrm{e}}_{q}^{-bE}.\end{equation} \tag{ B.3 }$$

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

$$\begin{equation}\rho \left(E\right)\sim \left(E+u\right)\sqrt{E\left(E+2u\right)}.\end{equation} \tag{ 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²P_p] = [P_E] = [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 J_E, corresponding to the distribution function P_E, which yields

$$\begin{equation}{J}_{E}^{\mathrm{m}\mathrm{o}\mathrm{d}}\left(E\right)\sim v\left(E\right){P}_{E}\left(E\right)\sim v\left(E\right)\rho \left(E\right){\mathrm{e}}_{q}^{-bE}.\end{equation} \tag{ 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 [vP_E] = [eV m⁻² s⁻¹].

In order to express the velocity in terms of E we use p = γmv with $\gamma =\frac{1}{\sqrt{1-v}}$ (in c = 1 convention), $p=A\sqrt{E\left(E+2u\right)}$ and m = Au to find

$$\begin{equation}v\left(E\right)=\frac{\sqrt{E\left(E+2u\right)}}{\left(E+u\right)}.\end{equation} \tag{ B.6 }$$

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

$$\begin{equation}{J}_{E}^{\mathrm{m}\mathrm{o}\mathrm{d}}\left(E\right)=CE\left(E+2u\right){\mathrm{e}}_{q}^{-bE}\end{equation} \tag{ B.7 }$$

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