Universal properties of primary and secondary cosmic ray energy spectra

Atomic nuclei appearing in cosmic rays 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 (AMS) 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 $\sim$ 200 MeV which agrees with the Hagedorn temperature. Our analysis suggests that QCD scattering processes together with nonequilibrium temperature fluctuations imprint heavily and universally onto the measured cosmic ray spectra, and produce a similar shape of energy spectra as high energy collider experiments on the Earth.

Atomic nuclei appearing in cosmic rays 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 (AMS) 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 heavily and universally onto the measured cosmic ray spectra, and produce a similar shape of energy spectra as high energy collider experiments on the Earth.

I. 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. Collectiveley 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 * smolla@usm.lmu.de † b.schaefer@qmul.ac.uk ‡ c.beck@qmul.ac.uk spectra contain many unknown parameters and are currently under debate [6]. As measured cosmic ray 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][9][10]. The occurrence of the index q can be naturally understood due to the fact that there are spatiotemporal 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 cosmic rays before in [12][13][14] and also applied to particle collisions in LHC experiments [15][16][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 tempera-ture 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 cosmic rays. 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 II we demonstrate that cosmic ray 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 III we interpret the observed spectra in terms of temperature fluctuations occurring during the production process of the cosmic rays, based on χ 2 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 IV.

II. RESULTS
We investigate the cosmic ray flux, given as differential intensity with respect to kinetic energy per nucleon, defined as E = (E total − m)/A, with total energy E total = p 2 + m 2 , 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 cosmic ray flux, we employ an established GSM model [13], modified slightly by replacing the total energy by the kinetic energy per nucleon E. In the Appendix we rigorously derive the distribution function which corresponds to the following differential intensity of flux where C, q > 1 and T = b −1 > 0 are free parameters and the q-exponential is defined as 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) = p 2 dp dE = (E + u) E(E + 2u) which leads to the flux derived from our superstatistical model which we compare to the observed flux J obs E . The entropic index q determines the high-energy (i.e. the tail) behavior of the distribution since the qexponential asymptotically approaches a power law with spectral index γ = 2 − 1/(q − 1) for q > 1. The parameter T = b −1 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. Fig. 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 χ 2 minimisation with (J mod E − J obs E )/σ, 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. Fig. 2 reveals the universal properties of the primary (He, C, O) and secondary (Li, Be, B) cosmic ray 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 Fig. 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.

III. INTERPRETATION IN TERMS OF TEMPERATURE FLUCTUATIONS
We consider the observed cosmic ray spectra to be the result of many different high energy scattering processes, each having a different local temperature β −1 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 cosmic rays, we need to generate asymptotic power laws and this can be achieved by socalled χ 2 superstatistics. As is generally known (see, e.g. [13]) the probability density function for a fluctuating β of the form with independent and identically distributed Gaussian random variables X i is a χ 2 distribution given by  (3) using three parameters C, T, q. The vertical axis in this log-log plot was multiplied with E 2.7 for better visibility. The fit's accuracy can be quantified by the deviation from modelled J mod to observed flux J obs 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 T0 is defined by (9). The amplitude C has 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 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 −βE dE. 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 with b = β 0 /(4 − 3q) and mean inverse temperature The effective degrees of freedom n are related to the en-  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 cosmic ray 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.
tropic index q via 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 χ 2 -distributed inverse temperatures.

IV. PHYSICAL INTERPRETATION AND POSSIBLE REASON FOR UNIVERSALITY
For the temperature parameter T 0 = β −1 0 , defined in (9), we get the value T 0 ∼ 600 MeV for each of the six CR species in our fits. Hence, the average effective temperature per quark is of the order T 0 /3 ∼ 200 MeV, i.e. we recover approximately the numerical value of the Hagedorn temperature T H ∼ 180 MeV, which is well known to be a universal parameter for the quark gluon plasma and for high energy QCD scattering processes. Remarkably, the fitted value of T 0 /3 in the fits is observed to be the same for all six nuclei, i.e. for both primary and secondary cosmic rays, within the statistical error.
The emergence of the Hagedorn temperature (at least as an order of magnitude) in our fits suggests that cosmic ray 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 12 K ∼ 100 MeV of comparable order of magnitude as the Hagedorn temperature [25]. As the Hagedorn temperature is universal, so is the average kinetic energy per quark of the cosmic rays nuclei, assuming they are produced in a Hagedorn fireball, either during during the original supernova explosion, or later during collision processes of highly energetic CR particles with the ISM.
Our observation that the kinetic energy per nucleon (or per quark) yields universal behavior 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. In fact, the shape of the observed cosmic ray spectra is very similar to those observed in high energy proton-proton collisions, see e.g. [17]. 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 −1 that is of similar order of magnitude (150 MeV) as in our fits for the cosmic rays, see table IV in [17]. That paper also shows that hard parton QCD scattering leads to power law spectral behavior.
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 behavior. One readily notices that in our GSM model the marginal distributions P E (E) decay asymptotically as In order to calculate the expectation of the fluctuating energy E, one notices that the integrand decays as EP (E) ∼ E −n/2 . Thus the expectation value is only well defined if n > 2. Since n is an integer (given by the number of Gaussian random variables contributing to the fluctuating β in eq.(5)), we thus can only have the values n = 3, 4, 5, . . ., and the smallest n is given by n min = 3.
A similar argument applies if one looks at the existence of the mean of the temperature β −1 as formed with the probability density g(β) The above mean only exists for n > 2, since the integrand behaves as β 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 β −1 is well-defined. This value of n implies q = 11/9 = 1.2222, very close to the fitted value for primary cosmic rays.
For secondary cosmic rays, there is an additional degree of freedom as an additional collision process at a later time is needed to produce secondary cosmic rays. Thus it is plausible that for secondary cosmic rays 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 Fig. 1), secondary cosmic rays 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.

V. CONCLUSION
We provide excellent fits for the measured AMS spectra of primary (He, C, O) and secondary cosmic rays (Li, Be, B) using a simple superstatistical model. The observed q-exponential spectra are interpreted in terms of temperature fluctuations occuring in the Hagedorn fireball during the production process of cosmic rays 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 P E = Cρ(E)e −bE q , using the framework of superstatistics by which we can interpret the best fit parameters with a temperature T = b −1 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(β) = ∞ 0 ρ(E)e −βE dE. In agreement with [13] we apply the ultra-relativistic approximation for the density of states ρ(E) ∼ E 2 in order to calculate Z(β) ∼ β −3 . Given the χ 2 -distributed β, defined by (6), we calculate the generalized canonical distribution as follows: Introducing q = 1+2/(n+6) (equivalent to n/2 = 1/(q − 1)−3) and b = β 0 /(4−3q), allows us to express the result as: (A2) Thus we have derived the distribution function (1), which we used for our fits, building on the framework of generalized statistical mechanics and superstatistics.
Note that the above equations are only valid for the particular case ρ(E) ∼ E 2 and g(β) being a χ 2 distribution. More generally, one has Appendix B: Applying theory to observation 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 = | p|, [R] = [V] and the corresponding flux measured in units [J(R)] = [m -2 sr -1 s -1 GV -1 ]. 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 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 Analogously to P E (E) ∼ ρ(E)e −bE q in the previous section one can derive that P p (E) ∼ e −bE q . Thus the density of states ρ(E) can be calculated from the conservation condition (B2). Using E = ( p 2 + m 2 − m)/A, which implies that the energy depends only on the magnitude of the momentum, simplifies d 3 p = 4πp 2 dp, and therefore Calculating the derivative and using p 2 = A 2 E(E + 2u) we find ρ(E) ∼ (E + u) E(E + 2u).
Note that we generally neglect constant global factors in our equations because we are focusing on the shape of the spectrum rather than its absolute magnitude. Evidently, . 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 [26] provides a detailed overview about the different ways to count particles including this relation. Evidently, it yields the desired physical dimensions since In order to express the velocity in terms of E we use p = γmv with γ = 1 Plugging everything into (B5) reveals the relation between q-exponential distribution function and the observed differential intensity