Partition function and thermodynamic parameters of the all-particle cosmic-ray flux

In Section 3, we study the phase-space probability distribution defining the spectral number and flux densities. We obtain the phase-space representation of the thermodynamic variables and derive the entropy functional of the stationary non-equilibrium densities introduced in Section 2. The spectral densities discussed in Sections 2 and 3 are classical. In Section 4, we quantize these densities in Fermi and Bose statistics depending on the spin of the nuclei. At low fugacity, the quantum distributions become virtually undistinguishable from the classical densities. In Section 6.2, we give estimates of the fugacity of cosmic-ray nuclei, which turns out to be very small, so that the classical flux densities are uniformly applicable over the complete spectral range.

In Section 5, we perform spectral fits to the proton and helium fluxes in the $1 GeV - 10^{4} GeV$ interval as defined by AMS-02 and CREAM flux data. We use the analytic factorizing flux densities introduced in Section 2, which also admit extrapolation to lower energies. As the cosmic-ray flux in this interval consists almost entirely of hydrogen and helium nuclei, we obtain an approximation of the low-energy all-particle flux by adding the proton and helium flux densities. In the ultra-relativistic $10^{5} GeV - 10^{11} GeV$ interval, where nuclear masses become negligible, we perform a fit to the all-particle spectrum covered by data sets from Tibet-III, IceTop-73, KASCADE-Grande, Auger, and Telescope Array. The spectrum in this interval shows several spectral breaks such as knee and ankle and two less pronounced breaks in their interconnecting slope as well as an ultra-high energy cutoff above 10¹⁰ GeV. The low-energy flux density obtained from the proton and helium spectra in the $1 GeV - 10^{4} GeV$ interval is extended to higher energies by iteratively adding factors to the low-energy spectral kernel which are defined by successive spectral breaks, as outlined in Appendix A. In this way, we obtain the flux density of the wideband spectrum ranging from 1 GeV to 10¹¹ GeV.

In Section 6, we show that the all-particle spectrum admits an additive decomposition into four exponentially decaying spectral peaks with crossovers at the spectral breaks. The partial flux densities generating the peaks are iteratively extracted from the factorizing total flux density. We calculate the thermodynamic parameters of the all-particle flux and the partial fluxes, namely the specific densities of number count, internal energy and entropy as well as the pressure. In Section 7, we present our conclusions. In Appendix A, we explain wideband spectral fitting with factorizing flux densities.

2. Flux densities with factorizing spectral kernel

The spectral number density of a relativistic gas in stationary non-equilibrium reads [10] and [11]

equation2.1

d ρ (p) = \frac{4 π s}{{(2 π)}^{3}} e^{- β \sqrt{p^{2} + m^{2}} - α} h (p) p^{2} d p,

Turn MathJax on

where we have put ℏ=c=1

ℏ = c = 1

. β=1/(k_BT)

β = 1 / (k_{B} T)

denotes the temperature parameter, p is the momentum variable and s the spin degeneracy of the gas particles. The fugacity parameter α determines the normalization of the density, the fugacity being z=e^−α

z = e^{- α}

. The all-particle cosmic-ray flux will be assembled with partial densities of type (2.1), with varying spin degeneracy and nuclear mass m, see Section 5.2. The spectral kernel h(p) defines a perturbation of the equilibrium phase-space probability, see Section 3, to be determined from spectral fits; the relativistic Maxwell equilibrium density is recovered by putting h(p)=1

h (p) = 1

in (2.1). The differential flux density per steradian based on number density (2.1) is

equation2.2

\frac{d Φ}{d p} [{(GeV / c)}^{- 1} s r^{- 1} m^{- 2} s^{- 1}] = \frac{υ}{4 π} \frac{d ρ}{d p},

Turn MathJax on

with momentum variable p in units of GeV/c and group velocity υ=dE/dp=p/E

υ = d E / d p = p / E

. The dispersion relation reads

E = \sqrt{p^{2} + m^{2}} = m γ_{L}

with Lorentz factor

γ_{L} = 1 / \sqrt{1 - υ^{2}}

. In Section 5, we will employ the rescaled flux density

equation2.3

F [{(GeV / c)}^{κ - 1} s r^{- 1} m^{- 2} s^{- 1}] = p^{κ} \frac{d Φ}{d p} = A_{Φ} \frac{p^{3 + κ} h (p)}{\sqrt{p^{2} + m^{2}}} e^{- β \sqrt{p^{2} + m^{2}}},

Turn MathJax on

where A_Φ

A_{Φ}

is the flux amplitude [12]

equation2.4

\begin{matrix} A_{Φ} [{(GeV / c)}^{- 3} s r^{- 1} m^{- 2} s^{- 1}] & = & \frac{s e^{- α}}{{(2 π)}^{3}} \frac{1}{c^{2} [m / s] ℏ^{3} [GeV s]} \\ \approx & 1.573 \times 10^{53} s e^{- α} \end{matrix}

Turn MathJax on

after restoring units. The spin degeneracy factor is s=2

s = 2

for protons and s=1

s = 1

for scalar bosons such as alpha particles. The scaling exponent κ in ( 2.3) can be freely chosen, as it does not affect least-squares fits, since the scale factor p^κ drops out of the χ² functional once the data points are rescaled with the same factor.

In order to minimize the χ² functional, a good initial guess of the flux density is needed, uniformly valid over the complete spectral range. In a wideband fit extending over several energy decades, it is crucial to first locate the spectral breaks. The starting point, therefore, is to rescale the flux data so that the breaks become clearly visible in log-log plots. In the case of the all-particle cosmic-ray flux, this is achieved with a scaling exponent κ=2.7 $κ = 2.7$ in (2.3). The rescaling facilitates the systematic construction of an initial guess for the flux density, from spectral break to spectral break, covering a spectral range stretching over eleven orders in energy, as explained in Section 5.1 and Appendix A.

The spectral number and flux densities in (2.1) and (2.3) are related by

equation2.5

d ρ [m^{- 3}] = \frac{4 π}{c} A_{Φ} h (p) e^{- β \sqrt{p^{2} + m^{2}}} p^{2} d p = \frac{4 π}{c} F (p) \frac{\sqrt{p^{2} + m^{2}}}{p^{1 + κ}} d p,

Turn MathJax on

with c[m/s] and p[GeV/c]. The deviation of flux density ( 2.3) from thermal equilibrium is specified by a factorizing spectral kernel in the form of a multiply broken power law with smooth transitions,

equation2.6

h (p) = \prod_{k = 1}^{n} {(1 + {(p / c_{k})}^{γ_{k} / | δ_{k} |})}^{δ_{k}},

Turn MathJax on

where the amplitudes c_k<c_k+1

c_{k} < c_{k + 1}

define the location of the spectral breaks, inferred from a spectral fit like the exponents γ_k > 0 and δ_k. Each factor in ( 2.6) represents a spectral break, and the factors are ordered by increasing amplitude. Wideband fits can iteratively be assembled by successively adding factors to the product (2.6), provided that the spectral breaks are sufficiently separated, c_k<<c_k+1

c_{k} < < c_{k + 1}

, see Section 5.1. The number density (2.1) converges to the Maxwell equilibrium distribution in the non-relativistic limit, since h(p→0)=1

h (p \to 0) = 1

.

To obtain integral representations of the thermodynamic variables, it is convenient to use a slightly more general distribution function than stated in (2.1),

equation2.7

\begin{matrix} d ρ (p) & = & \frac{4 π s}{{(2 π)}^{3}} \exp [- α - β \sqrt{p^{2} + m^{2}} \\ + γ \log h (p) - δ \frac{1}{3} \frac{p^{2}}{\sqrt{p^{2} + m^{2}}}] p^{2} d p, \end{matrix}

Turn MathJax on

which defines the partition function

\log Z = V \int_{0}^{\infty} d ρ (p)

studied in Section 3; V is a volume factor. The parameters γ and δ drop out in variables obtained as derivatives of log Z at γ=1

γ = 1

, δ=0

δ = 0

. The number count N, internal energy U and pressure P read

equation2.8

N = - {(\log Z)}_{, α} = V \int_{0}^{\infty} d ρ (p) = \log Z,

Turn MathJax on

equation2.9

U = - {(\log Z)}_{, β} = V \int_{0}^{\infty} \sqrt{p^{2} + m^{2}} d ρ (p),

Turn MathJax on

equation2.10

V P = - {(\log Z)}_{, δ} = \frac{1}{3} V \int_{0}^{\infty} \frac{p^{2}}{\sqrt{p^{2} + m^{2}}} d ρ (p) .

Turn MathJax on

The integration measure in (2.8)–(2.10) is the number density (2.1), coinciding with density (2.7) at γ=1 $γ = 1$ , δ=0 $δ = 0$ . The integral kernel of the pressure variable (2.10) is the product of momentum p , velocity $p / \sqrt{p^{2} + m^{2}}$ and a geometric factor $\int_{0}^{1} \cos^{2} θ dcos θ = 1 / 3$ taking account of the incidence angles of the gas particles (incident upon an infinitesimal surface element). The entropy rescaled with the Boltzmann constant is obtained as

equation2.11

S/k_B=logZ−α(logZ)_,α−β(logZ)_,β−(logZ)_,γ,

S / k_{B} = \log Z - α {(\log Z)}_{, α} - β {(\log Z)}_{, β} - {(\log Z)}_{, γ},

Turn MathJax on

where log Z is taken at γ=1

γ = 1

, δ=0

δ = 0

. The fugacity parameter α can be extracted from the flux amplitude ( 2.4) and the temperature parameter β from the exponential cutoff in flux density ( 2.3). In Section 3, we will derive the entropy variable (2.11) from the phase-space probability density. The fourth term in (2.11) is the excess entropy

equation2.12

H / k_{B} = - {(\log Z)}_{, γ} = - V \int_{0}^{\infty} \log h (p) d ρ (p),

Turn MathJax on

defining the deviation of S from an equilibrium entropy,

equation2.13

S=k_B(logZ+αN+βU)+H.

S = k_{B} (\log Z + α N + β U) + H .

Turn MathJax on

For the classical densities (2.1), this can be simplified by substituting the classical relation logZ=N $\log Z = N$ according to (2.8). The entropy of the quantum distributions discussed in Section 4 is defined by (2.13), with number density dρ(p ) in the above integral representations replaced by the quantized counterpart. In equilibrium, h(p)=1 $h (p) = 1$ , the excess entropy H vanishes and we can identify βPV=logZ(δ=0) $β P V = \log Z (δ = 0)$ via integration by parts.

3. Partition function and entropy

To derive the classical partition function defined by number density (2.1), we start with the n-particle phase-space measure

equation3.1

d^{3 n} (p, q) = \frac{1}{{(2 π)}^{3 n}} d^{3} p_{1} d^{3} q_{1} d^{3} p_{2} d^{3} q_{2} \dots d^{3} p_{n} d^{3} q_{n},

Turn MathJax on

where the spatial d³q_k integrations range over a volume V and the d³p_k integrations over momentum space. (The normalization with (2π)³ⁿ is obtained by performing the classical limit in the quantized number densities derived in Section 4.) The n-particle probability density leading to number density ( 2.7) reads [13]

equation3.2

P_{n} (p) = \frac{1}{Z} b_{n} (p),

Turn MathJax on

equation3.3

\begin{matrix} b_{n} (p) & = & \exp (- α n - β \sum_{k = 1}^{n} \sqrt{p_{k}^{2} + m^{2}} \\ + γ \sum_{k = 1}^{n} \log h (p_{k}) - δ \frac{1}{3} \sum_{k = 1}^{n} \frac{p_{k}^{2}}{\sqrt{p_{k}^{2} + m^{2}}}), \end{matrix}

Turn MathJax on

where Z is a normalization constant and b_n(p)

b_{n} (p)

is the phase-space analogue to the exponential in (2.7). The spectral function h(p) defined in ( 2.6) quantifies the perturbation of the equilibrium phase-space probability. Density P_n(p)

P_{n} (p)

is independent of the space coordinates as there is no interaction potential included. We will ultimately put γ=1

γ = 1

, δ=0

δ = 0

as in Section 2. The normalization condition is

equation3.4

\sum_{n = 0}^{\infty} \frac{s^{n}}{n!} \int P_{n} (p) d^{3 n} (p, q) = 1,

Turn MathJax on

where the summation is over the fluctuating particle number n. The factorial accounts for indistinguishability, and s is the spin degeneracy. In this way, we find the normalization constant Z in ( 3.2) as

equation3.5

Z = \sum_{n = 0}^{\infty} \frac{s^{n}}{n!} \int b_{n} (p) d^{3 n} (p, q) .

Turn MathJax on

This is the grand partition function in phase-space representation, which can be simplified as stated after (2.7), as the integrals factorize.

Differentiation of Z with respect to the fugacity parameter α in ( 3.3) at γ=1 $γ = 1$ , δ=0 $δ = 0$ gives

equation3.6

- \frac{Z_{, α}}{Z} = \sum_{n = 0}^{\infty} \frac{s^{n}}{n!} \int n P_{n} (p) d^{3 n} (p, q) = 〈 n 〉 .

Turn MathJax on

The expectation value 〈n 〉 coincides with the number count N=logZ $N = \log Z$ in (2.8), after evaluation of the factorizing integrals. The phase-space representation of internal energy, pressure and excess entropy in (2.9), (2.10) and (2.12) is obtained in like manner via differentiation of Z in ( 3.5),

equation3.7

U = - {(\log Z)}_{, β} = \sum_{n = 0}^{\infty} \frac{s^{n}}{n!} \int u_{n} P_{n} (p) d^{3 n} (p, q) = 〈 u_{n} 〉,

Turn MathJax on

equation3.8

V P = - {(\log Z)}_{, δ} = \sum_{n = 0}^{\infty} \frac{s^{n}}{n!} \int v_{n} P_{n} (p) d^{3 n} (p, q) = 〈 v_{n} 〉,

Turn MathJax on

equation3.9

H / k_{B} = - {(\log Z)}_{, γ} = \sum_{n = 0}^{\infty} \frac{s^{n}}{n!} \int h_{n} P_{n} (p) d^{3 n} (p, q) = 〈 h_{n} 〉,

Turn MathJax on

where

equation3.10

u_{n} = \sum_{k = 1}^{n} \sqrt{p_{k}^{2} + m^{2}}, v_{n} = \frac{1}{3} \sum_{k = 1}^{n} \frac{p_{k}^{2}}{\sqrt{p_{k}^{2} + m^{2}}}, h_{n} = - \sum_{k = 1}^{n} \log h (p_{k})

Turn MathJax on

are the random variables of energy, pressure and excess entropy in phase space. All derivatives are taken at γ=1

γ = 1

, δ=0

δ = 0

. The entropy variable is assembled as

equation3.11

S / k_{B} = \sum_{n = 0}^{\infty} \frac{s^{n}}{n!} \int s_{n} P_{n} (p) d^{3 n} (p, q) = 〈 s_{n} 〉,

Turn MathJax on

equation3.12

s_n=−logP_n(p)=logZ+αn+βu_n+h_n,

s_{n} = - \log P_{n} (p) = \log Z + α n + β u_{n} + h_{n},

Turn MathJax on

where log Z and the probability density P_n(p)

P_{n} (p)

are taken at γ=1

γ = 1

, δ=0

δ = 0

. Substitution of the above averages into (3.11) leads to S in ( 2.13). The quantum analogue to the −k_BP_nlogP_n

- k_{B} P_{n} \log P_{n}

phase-space representation of entropy is briefly discussed after (4.7).

4. Quantization of stationary non-equilibrium densities

The quantized version of the classical partition function Z in ( 3.5) is obtained by a trace calculation in fermionic (F) or bosonic (B) occupation number representation [14] and [15],

equation4.1

Z_{F, B} = Tr [\exp (- α \hat{N} - β \hat{E} - γ \hat{H} - δ \hat{P})] .

Turn MathJax on

The particle number operator is $\hat{N} = \sum_{k} {\hat{N}}_{k}$ , ${\hat{N}}_{k} = a_{k}^{†} a_{k}$ , where the $a_{k}^{(†)}$ are fermionic or bosonic creation/annihilation operators. The energy operator reads $\hat{E} = \sum_{k} \sqrt{k^{2} + m^{2}} {\hat{N}}_{k}$ , the excess entropy is quantized as $\hat{H} = - \sum_{k} \log h (k) {\hat{N}}_{k}$ , and the pressure operator is $\hat{P} = (1 / 3) \sum_{k} k^{2} {(k^{2} + m^{2})}^{- 1 / 2} {\hat{N}}_{k}$ . The Hermitian number operators are diagonal, ${\hat{N}}_{k} | n 〉 = n_{k} | n 〉$ . The fermionic occupation numbers n_k attached to a wave vector k can take the values zero and one, whereas for bosons, the n_k are non-negative integers. We employ box quantization, discretizing the wave vectors as k=2πn/L $k = 2 π n / L$ , so that k summations are taken over integer lattice points n, corresponding to periodic boundary conditions on a box of volume V=L³ $V = L^{3}$ .

With these prerequisites, trace (4.1) can be evaluated in Fermi statistics as

equation4.2

Z_{F} = \sum_{{n_{k}} = 0, 1} \exp (- \sum_{k} G (k) n_{k}) = \prod_{k} (1 + e^{- G (k)}),

Turn MathJax on

where we defined the shortcut

equation4.3

G (k) = α + β \sqrt{k^{2} + m^{2}} - γ \log h (k) + δ \frac{1}{3} \frac{k^{2}}{\sqrt{k^{2} + m^{2}}} .

Turn MathJax on

In Bose statistics, we find

equation4.4

Z_{B} = \sum_{{n_{k}} = 0}^{\infty} \exp (- \sum_{k} G (k) n_{k}) = \prod_{k} \sum_{n_{k} = 0}^{\infty} e^{- G (k) n_{k}} = \prod_{k} \frac{1}{1 - e^{- G (k)}} .

Turn MathJax on

Accordingly,

equation4.5

\log Z_{F, B} = \pm \sum_{k} \log (1 \pm e^{- G (k)}) .

Turn MathJax on

Performing the thermodynamic continuum limit L → ∞ and taking into account the spin degeneracy, we can replace the summation in ( 4.5) by the integration sV∫d³k/(2π)³,

equation4.6

\log Z_{F, B} = \pm \frac{4 π s V}{{(2 π)}^{3}} \int_{0}^{\infty} \log (1 \pm e^{- G (k; α, β, γ, δ)}) k^{2} d k .

Turn MathJax on

The quantized spectral number densities dρ_{F, B}(k) can be identified by means of the particle count in ( 2.8),

equation4.7

\begin{matrix} N = - {(\log Z_{F, B})}_{, α} = V \int_{0}^{\infty} d ρ_{F, B} (k), \\ d ρ_{F, B} (k) = \frac{4 π s}{{(2 π)}^{3}} \frac{k^{2} d k}{e^{G (k; α, β, γ, δ)} \pm 1}, \end{matrix}

Turn MathJax on

taken at γ=1

γ = 1

, δ=0

δ = 0

. For instance, density dρ_F/dk with s=2

s = 2

applies to protons and dρ_B/dk with s=1

s = 1

to ⁴He nuclei. The internal energy U=−(logZ_F,B)_,β

U = - {(\log Z_{F, B})}_{, β}

, the excess entropy H/k_B=−(logZ_F,B)_,γ

H / k_{B} = - {(\log Z_{F, B})}_{, γ}

and the pressure VP=−(logZ_F,B)_,δ

V P = - {(\log Z_{F, B})}_{, δ}

are obtained from the quantized partition function (4.6) analogously to the classical case in (2.9), (2.10) and (2.12), and the total entropy S is assembled as in ( 2.13).

The quantized version of the phase-space representation of entropy in (3.11) is $S = - k_{B} Tr [\hat{ρ} \log \hat{ρ}]$ ; the counterpart to the phase-space probability density P_n(p) $P_{n} (p)$ in (3.2) is the density operator $\hat{ρ} = \exp (- \hat{s})$ , and the entropy operator $\hat{s} = \log Z_{F, B} + α \hat{N} + β \hat{E} + \hat{H}$ is the counterpart to the entropy random variable s_n in ( 3.12). The quantized partition function Z_F,B(γ=1,δ=0) $Z_{F, B} (γ = 1, δ = 0)$ is stated in (4.1), and the operators $\hat{N}$ , $\hat{E}$ and $\hat{H}$ are defined after (4.1). The excess entropy $H / k_{B} = 〈 \hat{H} 〉 = Tr [\hat{H} \hat{ρ}]$ is determined by the spectral function h(p) in ( 2.6) which can be obtained from a spectral fit, see Section 5. $\hat{H}$ specifies the deviation of $\hat{ρ}$ from an equilibrium density operator. The counterpart to the phase-space normalization (3.4) is $Tr [\hat{ρ}] = 1$ .

The fermionic/bosonic flux densities dΦ_{F, B}/dk read as in ( 2.2), with the classical density dρ/dp replaced by dρ_{F, B}/dk . (k=p $k = p$ since we have put ℏ=c=1 $ℏ = c = 1$ .) To recover the classical limit (2.1) of the number densities dρ_{F, B}/dk, we only need to replace the factor 1/(e^G(p) ± 1) in (4.7) by $e^{- G (p; α, β, γ = 1, δ = 0)} = h (p) e^{- β \sqrt{p^{2} + m^{2}} - α}$ . The condition for the classical limit is thus e^−G(p)<<1 $e^{- G (p)} < < 1$ , which holds for cosmic-ray nuclei over the complete spectral range, as the fugacity z=e^−α $z = e^{- α}$ of the partial fluxes is very small, see the caption to Table 2. Therefore, we can use the classical number and flux densities in (2.1) and (2.2) to estimate the thermodynamic parameters of cosmic-ray fluxes, see Section 6.2.

5. Proton and helium fluxes and the all-particle energy spectrum

5.1. Product representation of spectral flux densities over multiple spectral breaks

The spectral fits are performed with the rescaled flux density (2.3), F=pκdΦ/dp $F = p^{κ} d Φ / d p$ ,

equation5.1

F [{(GeV / c)}^{κ - 1} s r^{- 1} m^{- 2} s^{- 1}] = \frac{A_{Φ} p^{3 + κ} h (p)}{\sqrt{p^{2} + m^{2}}} e^{- β \sqrt{p^{2} + m^{2}}},

Turn MathJax on

where

h (p) = \prod_{k = 1}^{n} {(1 + {(p / c_{k})}^{γ_{k} / | δ_{k} |})}^{δ_{k}}

is the factorizing spectral kernel introduced in (2.6), a multiply broken power-law with smooth transitions at the spectral breaks. We use the scaling exponent κ=2.7

κ = 2.7

to make spectral breaks visible in log-log plots of the all-particle spectrum covering the complete

1 GeV - 10^{11} GeV

interval, see also the remarks after (2.4). Usually (but not always) the exponents δ_k of the spectral kernel have alternating signs, assuming that the factors of h(p ) are ordered by increasing amplitude, c_k<c_k+1

c_{k} < c_{k + 1}

. The c_k define the location of the spectral breaks, and the exponents γ_k together with the sign of δ_k determine the slopes of the spectral curve between the breaks. A moderate or large |δ_k| leads to a smoothing of the spectral break at p ≈ c_k, that is, to a gradual change of the slope in the vicinity of p ≈ c_k. A |δ_k| close to zero results in a sudden break in the spectral curve, appearing as edge joining two nearly straight (approximately power-law) slopes in log-log plots. This is illustrated by the spectral breaks visible in the figures and the numerical value of the exponents δ_k attached to the respective break points c_k listed in Table 1. The product representation of h(p ) is efficient if the spectral breaks are widely separated, c_k≪c_k+1

c_{k} ≪ c_{k + 1}

. In this case, one can perform an iterative fit over an extended spectral range as outlined in Appendix A. Finally, the nuclear mass in (5.1) becomes negligible if m²/p² ≪ 1. That is, at energies above 10⁴ GeV, in the ultra-relativistic regime, the mass squares in (5.1) can be dropped even for iron nuclei.

Table 1.

Fitting parameters of the proton and helium spectra in Figs. 1 and 2 and of the all-particle spectrum in Fig. 3, Fig. 4, Fig. 5, Fig. 6, Fig. 7 and Fig. 8. The spectral fits are performed with the proton and helium flux densities F_{H, low} and F_{He, low} in (5.3) and the all-particle flux density F_all in (5.4). The amplitudes c_k in the power-law factors of F_all define the location of the spectral breaks, and the exponents γ_k and δ_k determine the slopes and curvature of the spectral curves depicted in the figures, as pointed out after ( 5.1). (c₄ is the knee, c₆ the second knee and c₇ the ankle, see Fig. 8.) The temperature parameter k_BT is extracted from the exponential spectral cutoff of the all-particle flux density ( 5.4); k_B denotes the Boltzmann constant. A_Φ,H $A_{Φ, H}$ and A_Φ,He $A_{Φ, He}$ are the amplitudes of the low-energy proton and helium flux densities (5.3).

k	c_k[GeV/c]	γ_k	δ_k
1H	1.283 ± 0.022	4.910 ± 0.050	−(4.125±0.22)
2H	(5.128 ± 2.3) × 102	0.3709 ± 0.14	0.3698 ± 0.33
1He	2.694 ± 0.084	4.775 ± 0.0094	−(3.848±0.12)
2He	(6.955 ± 1.1) × 102	0.2045 ± 0.024	(6.409±3.2)×10−2
3	(6.925 ± 0.15) × 104	(8.896±0.094)×10−2	−(3.106±0.33)×10−2
4	(3.935 ± 0.15) × 106	0.5220 ± 0.032	−(8.835±1.6)×10−2
5	(1.227 ± 0.14) × 107	0.2080 ± 0.039	2.312×10−4
6	(1.091 ± 0.20) × 108	0.2735 ± 0.026	−3.317×10−3
7	(4.866 ± 0.24) × 109	1.041 ± 0.081	0.1678 ± 0.069
	AΦ,H	AΦ,He	kBT[GeV]
	[(GeV/c)−3sr−1m−2s−1]	[(GeV/c)−3sr−1m−2s−1]
	(8.428 ± 1.6) × 103	(1.095 ± 0.20) × 102	(2.458 ± 0.19) × 1010

Full-size table

Table options

In view of the phase-space reconstruction of the thermodynamic variables in Section 3 and the quantization of the partition function in Section 4, the momentum parametrization of flux density (5.1) is the natural choice. We briefly point out the conversion of the momentum parametrization to a parametrization by kinetic energy per nucleon as well as to a rigidity parametrization [16], both often used in observational papers. We write the energy of the nucleus as E=nE_k+m=mγ_L $E = n E_{k} + m = m γ_{L}$ , with Lorentz factor $γ_{L} = 1 / \sqrt{1 - υ^{2}}$ , where n is the number of nucleons, E_k $E_{k}$ the kinetic energy per nucleon and m the mass of the nucleus. The dispersion relation is $p = \sqrt{E^{2} - m^{2}}$ . We use GeV units for energy, writing p[GeV/c ] (momentum per nucleus) in the figures. The conversion of momentum plots to plots parametrized by kinetic energy per nucleon E_k[GeV/n] $E_{k} [GeV / n]$ is thus performed with $p = \sqrt{n E_{k} (n E_{k} + 2 m)}$ and

equation5.2

E_{k}^{κ} \frac{d Φ}{d E_{k}} = \frac{n (n E_{k} + m) E_{k}^{κ}}{{[n E_{k} (n E_{k} + 2 m)]}^{(κ + 1) / 2}} F (\sqrt{n E_{k} (n E_{k} + 2 m)}),

Turn MathJax on

where F(p)=pκdΦ/dp

F (p) = p^{κ} d Φ / d p

as in (5.1). We will put n=1

n = 1

for hydrogen and n=4

n = 4

for helium, assuming ¹H and ⁴He isotopes [8] and [17]. At high energy, we can neglect the rest mass (using the 0th order of an m/(nE_k)

m / (n E_{k})

expansion of (5.2)), and approximate p∼nE_k

p \sim n E_{k}

,

E_{k}^{κ} d Φ / d E_{k} \sim n^{1 - κ} F (n E_{k})

. The conversion of the momentum parametrization in (5.1) to a parametrization by rigidity R=cp/(Ze)

R = c p / (Z e)

(where e is the electron charge, taken positive, and Z the charge number of the nuclei) is performed with p[GeV/c]=ZR[GV]

p [GeV / c] = Z R [GV]

and RκdΦ/dR=Z1−κF(ZR)

R^{κ} d Φ / d R = Z^{1 - κ} F (Z R)

. For protons, momentum plots in p[GeV/c] thus coincide with rigidity R[GV] plots.

Fig. 1 shows spectral fits to the AMS-02 hydrogen and helium fluxes, with CREAM data points added, in the low-energy interval from 1 GeV to 10⁴ GeV. For protons, the fit is performed with flux density

equation5.3

\begin{matrix} F_{H, low} (p) & = & \frac{A_{Φ, H} p^{3 + κ}}{\sqrt{p^{2} + m_{H}^{2}}} h_{H, low} (p), \\ h_{H, low} (p) & = & \frac{{(1 + {(p / c_{2 H})}^{γ_{2 H} / δ_{2 H}})}^{δ_{2 H}}}{{(1 + {(p / c_{1 H})}^{γ_{1 H} / | δ_{1 H} |})}^{| δ_{1 H} |}}, \end{matrix}

Turn MathJax on

Fig. 1.

Proton and helium flux densities in the low-energy regime $1 GeV - 10^{4} GeV$ . Data points from AMS-02 [7] and [8] and CREAM [9]. The hydrogen and helium spectra are fitted with the flux densities F_{H, low}(p) and F_{He, low}(p) in ( 5.3) (solid curves) and parameters listed in Table 1. The helium CREAM points are rescaled by a factor of 0.86 to obtain a smooth continuation of the AMS-02 data set. The error bands defined by the dotted curves indicate the 1σ confidence limits of the least-squares fits, χ²(H)/dof ≈ 15.1/68, χ²(He)/dof ≈ 12.2/66; the residuals are depicted in the lower panels. By adding the spectral curves, F_H,low+F_He,low $F_{H, low} + F_{He, low}$ , we obtain an approximation of the all-particle spectrum, as heavier nuclei contribute only about one percent to the total flux at energies below 10⁴ GeV. The specific number densities of protons and helium nuclei below 10⁴ GeV are 1.07×¹⁰⁻⁴/m³ $1.07 \times 10^{- 4} / m^{3}$ and 1.6×¹⁰⁻⁵/m³ $1.6 \times 10^{- 5} / m^{3}$ , respectively, obtained by truncating the integral representation of the number density in (6.8) at 10⁴ GeV and substituting the low-energy flux limits (5.3). The spectral curve of the proton flux F_{H, low} remains unchanged in rigidity parametrization, if the coordinate labels are replaced by R [GV] and R2.7dΦ/dR[GV^1.7sr⁻¹m⁻²s⁻¹] $R^{2.7} d Φ / d R [G V^{1.7} s r^{- 1} m^{- 2} s^{- 1}]$ , see after (5.2); this is not the case for the helium flux, as the charge number differs from one.

Figure options

where $m_{H} \approx 0.938 GeV$ is the proton mass and κ=2.7 $κ = 2.7$ the rescaling exponent, see (2.3) and after (5.1). As there is no spectral cutoff visible in the low-energy data sets, we drop the exponential in (5.1), assuming a sufficiently high temperature so that $e^{- β_{H} \sqrt{p^{2} + m_{H}^{2}}} \sim 1$ below 10⁴ GeV. The fitting parameters, labeled by a subscript H in (5.3), are listed in Table 1. The same density (5.3) with subscript H replaced by He is employed in the spectral fit of the low-energy helium flux F_{He, low}. We use $m_{He} \approx 3.727 GeV$ as helium mass assuming alpha particles. The parameters A_Φ,He $A_{Φ, He}$ , c_1He $c_{1 He}$ , etc. of the helium fit in Fig. 1 are also recorded in Table 1. We obtain an approximation of the low-energy all-particle flux density in the $1 GeV - 10^{4} GeV$ interval by adding the proton and helium fits, F_all,low∼F_H,low+F_He,low $F_{all, low} \sim F_{H, low} + F_{He, low}$ , see Fig. 1, since heavier nuclei contribute only about one percent to the total flux in this energy range. In Fig. 2, the momentum parametrization of the proton and helium fluxes is converted to a parametrization by kinetic energy per nucleon by means of (5.2).

Fig. 2.

Spectral densities of the low-energy hydrogen and helium fluxes, parametrized by kinetic energy per nucleon. Data points from AMS-02 [7] and [8] and CREAM [9]. The spectral curves and flux points in Fig. 1 are converted from momentum representation to E_k (kinetic energy per nucleon) representation by way of transformation (5.2), assuming ¹H and ⁴He isotopes. As we consider a free relativistic gas of cosmic-ray nuclei rather than nucleons, a flux parametrization by momentum per nucleus has been adopted in this paper except for this figure; also see the end of Section 5.1 for flux parametrization in the ultra-relativistic regime above 10⁴ GeV.

Figure options

In Fig. 3, we perform a spectral fit to the all-particle flux above 10⁵ GeV. To this end, we extend the low-energy flux F_{all, low} by adding power-law factors defined by successive spectral breaks as explained after (5.1):

equation5.4

\begin{matrix} F_{all} (p) & = & \frac{F_{H, low} (p) + F_{He, low} (p)}{{(1 + {(p / c_{3})}^{γ_{3} / | δ_{3} |})}^{| δ_{3} |}} \frac{{(1 + {(p / c_{5})}^{γ_{5} / δ_{5}})}^{δ_{5}}}{{(1 + {(p / c_{4})}^{γ_{4} / | δ_{4} |})}^{| δ_{4} |}} \\ \times \frac{{(1 + {(p / c_{7})}^{γ_{7} / δ_{7}})}^{δ_{7}}}{{(1 + {(p / c_{6})}^{γ_{6} / | δ_{6} |})}^{| δ_{6} |}} e^{- β p} . \end{matrix}

Turn MathJax on

Fig. 3.

Wideband spectral map of the all-particle flux. Low-energy data points from AMS-02 [7] and [8] and CREAM [9]. High-energy data from Tibet-III (QGSJET+PD) [1], IceTop-73 (rescaled by a factor of 0.8) [2] and [3], KASCADE-Grande [4], Auger (rescaled by 1.1) [5] and Telescope Array (SD and FD Mono, rescaled by 0.9) [6]. The rescaling is required to make the data sets compatible in overlapping intervals so that they define a continuous spectral curve. The solid curve depicts the all-particle flux F_all defined by flux density (5.4) with parameters in Table 1. Below 10⁴ GeV, F_all can be approximated by adding the low-energy proton and helium fluxes in Fig. 1, F_all∼F_H,low+F_He,low $F_{all} \sim F_{H, low} + F_{He, low}$ , as pointed out after (5.3). Above 10⁵ GeV, F_all is inferred from a least-squares fit to the high-energy flux points as explained after (5.4), with χ²/dof ≈ 223/101; also see the close-ups in Fig. 5, Fig. 6 and Fig. 7. The dashed and dotted curves depict extensions of the low-energy proton and helium fluxes in Fig. 1 defined by the partial fluxes F_{H, ext} and F_{He, ext} in (5.5), both of which contain admixtures of heavier nuclei and add up to the total all-particle flux density F_all=F_H,ext+F_He,ext $F_{all} = F_{H, ext} + F_{He, ext}$ of the wideband spectrum ranging from 1 GeV to 10¹¹ GeV. These analytic flux densities also extrapolate the spectrum to energies below 1 GeV, as required in the integral representations (6.8)–(6.11) of the thermodynamic variables.

Figure options

The parameters of the high-energy fit in the $10^{5} GeV - 10^{11} GeV$ interval are (c_k, γ_k, δ_k ), k=3,…,7 $k = 3, \dots, 7$ and the cutoff temperature k_BT defining the decay exponent β=1/k_BT $β = 1 / k_{B} T$ , see Table 1. The low-energy proton and helium fluxes F_{H, low} and F_{He, low} are considered input obtained from the fits in Fig. 1. The factors in (5.4) are ordered by increasing spectral break, ¹⁰⁴<<c_k<<c_k+1<<k_BT $10^{4} < < c_{k} < < c_{k + 1} < < k_{B} T$ , see Table 1, so that F_all approximates F_{all, low} below 10⁴ GeV. In the χ² functional of F_all(p), the parameters defining the densities F_{H, low} and F_{He, low} (labeled with an H or He subscript in (5.3) and Table 1) are kept constant at the values inferred from the above described fits of the hydrogen and helium spectra. The high-energy flux points included in χ² stem from Tibet-III [1], IceTop-73 [2] and [3], KASCADE-Grande [4], Auger [5] and Telescope Array [6], see Fig. 3. The sixteen parameters varied in the least-squares fit of F_all(p) are listed in Table 1 (without an H or He subscript). To find the minimum of this χ², it is essential to have a good initial guess of the flux density; an iterative method to find this is sketched in Appendix A. The flux density F_all defined in (5.3) and (5.4) with parameters in Table 1 reproduces the all-particle energy spectrum in the $1 GeV - 10^{11} GeV$ range as depicted in Fig. 3. This analytic density also provides an extrapolation of the spectrum to energies below 1 GeV as well as an interpolation in the $10^{4} GeV - 10^{5} GeV$ range not yet covered by data points, see the close-up in Fig. 4.

Fig. 4.

Crossover from the low- to the high-energy regime of the all-particle spectrum. This is a close-up of Fig. 3. The partial fluxes F_{H, ext} and F_{He, ext} (dashed and dotted curves) are extensions of the low-energy spectral fits F_{H, low} and F_{He, low} of the proton and helium spectra in Fig. 1, see (5.5) and the caption to Fig. 3. In the $10^{4} GeV - 10^{5} GeV$ decade, the partial fluxes F_{H, ext} and F_{He, ext} start to deviate from the depicted CREAM flux points due to admixtures of heavier nuclei, see also the remarks after (5.6). (Therefore, in the χ² functionals of the low-energy fits in Fig. 1, we only included the first three hydrogen CREAM points and the first five helium CREAM points, respectively.) The all-particle flux density F_all=F_H,ext+F_He,ext $F_{all} = F_{H, ext} + F_{He, ext}$ (solid curve) is the continuation of the low-energy flux F_H,low+F_He,low $F_{H, low} + F_{He, low}$ in Fig. 1 to ultra-relativistic energies above 10⁴ GeV; the analytic representation thereof is given in (5.4) and (5.5). A soft spectral break emerges at c₃ ≈ 6.9 × 10⁴ GeV, see Table 1.

Figure options

At energies below 10 GeV, the proton and helium spectra become time dependent due to solar modulations as exemplified by the BESS-Polar spectra [18]. This time dependence can be accommodated in an adiabatic time variation of the fitting parameters (c_1H, γ_1H, δ_1H) and (c_1He, γ_1He, δ_1He) (which define the first spectral break of the proton and helium fluxes causing the broad peaks in Fig. 1, see (5.3) and Table 1), resulting in a quasi-stationary flux density (5.4). In this paper, we use time-averaged spectra and stationary spectral densities to model the proton and helium fluxes.

The fine-structure of the all-particle spectrum in the $10^{5} GeV - 10^{11} GeV$ interval is illustrated by the close-ups in Fig. 5, Fig. 6 and Fig. 7, the latter two extending over four energy decades each. Between the knee c₄ $c_{4}$ and the ankle c₇ $c_{7}$ (listed in Table 1), there are two soft spectral breaks at $c_{5} \approx 1 . 2 \times 10^{7} GeV$ and $c_{6} \approx 1 . 1 \times 10^{8} GeV$ , which are clearly visible in Fig. 6 and accounted for by the respective factors in flux density (5.4). A close-up of the exponential spectral cutoff generated by the Boltzmann factor in the flux density is shown in Fig. 7. The Telescope Array flux points above 2 × 10¹⁰ GeV have not been included in the fit, as they substantially deviate from the Auger flux, which decays exponentially and allows us to extract the temperature parameter. The decay indicated by the Telescope Array points is less pronounced and may even turn out to be subexponential, of Weibull type for instance, with the Boltzmann factor in (5.4) replaced by exp(−ρ(βp)) $\exp (- {(β p)}^{ρ})$ , the exponent 0 < ρ < 1 being an additional fitting parameter.

Fig. 5.

High-energy spectral fit to the all-particle flux in the $10^{5} GeV - 10^{11} GeV$ interval, a close-up of the wideband spectrum in Fig. 3. Data points as listed in the caption of Fig. 3. The least-squares fit F_all (solid curve, χ²/dof ≈ 223/101) is performed with flux density (5.4), the fitting parameters are recorded in Table 1. In this ultra-relativistic regime, the coordinate labels can be replaced by E [GeV] (energy per nucleus) and E2.7dΦ/dE[GeV^1.7sr⁻¹m⁻²s⁻¹] $E^{2.7} d Φ / d E [Ge V^{1.7} s r^{- 1} m^{- 2} s^{- 1}]$ , as the nuclear rest mass drops out of the dispersion relation. The dotted curves indicate the 1σ error band (better visible in the close-ups of Figs. 6 and 7). The last ten Telescope Array points depicted in the figure have not been included in the fit, see the remarks at the end of Section 5.1 and the caption to Fig. 7. The residuals are shown in the second and third panel. The latter includes the residuals of all Telescope Array points except the last one (with ΔF_all/F_all ≈ 63.5), which is off the scale due to the rapid exponential decay of F_all. The second panel is a close-up of the third.

Figure options

Fig. 6.

Close-up of the high-energy all-particle spectrum in Fig. 5. Data points as in Fig. 3. The solid curve F_all depicts flux density (5.4) with parameters in Table 1, the dotted curves define the 1σ error band. Residuals are shown in the lower panel. The knee is defined by the spectral break at $c_{4} \approx 3.9 \times 10^{6} GeV$ , see also (5.4) and Table 1, and is followed by two less pronounced but still clearly visible breaks at $c_{5} \approx 1 . 2 \times 10^{7} GeV$ and $c_{6} \approx 1 . 1 \times 10^{8} GeV$ in the slope joining knee and ankle, the latter is shown in Fig. 7.

Figure options

Fig. 7.

Decaying spectral tail of the all-particle spectrum. Close-up of Fig. 5, data points as in Fig. 3; the depicted ten highest energy points of the Telescope Array have not been included in the χ² fit, as they are not compatible with the Auger data. The latter define an exponential spectral cutoff, whereas the Telescope Array points suggest subexponential Weibull spectral decay [31], [32], [33] and [34]. The caption to Fig. 5 applies regarding the flux density F_all (solid curve). The 1σ error band is indicated by the dotted curves. The ankle is the spectral break at $c_{7} \approx 4 . 9 \times 10^{9} GeV$ , see Table 1. The exponential spectral decay is caused by the ultra-relativistic Boltzmann factor e^−p/(k_BT) $e^{- p / (k_{B} T)}$ of flux density (5.4) with cutoff temperature k_BT ≈ 2.5 × 10¹⁰ GeV. The residuals are plotted in the two lower panels at different scales.

Figure options

Above 10⁴ GeV, we can identify E ∼ p, as the nuclear mass drops out of the dispersion relation, see after ( 2.2). Therefore, the plots in Fig. 5, Fig. 6 and Fig. 7 do not change in energy parametrization, that is, the coordinate labels can be replaced by E [GeV] and E2.7dΦ/dE[GeV^1.7sr⁻¹m⁻²s⁻¹] $E^{2.7} d Φ / d E [Ge V^{1.7} s r^{- 1} m^{- 2} s^{- 1}]$ in this ultra-relativistic regime. Here, E is the energy per nucleus, not to be confused with E_k labeling Fig. 2, which denotes energy (with nuclear rest mass subtracted) per nucleon, see (5.2).

5.2. Decomposition of the all-particle spectrum into partial fluxes

The all-particle flux density F_all in (5.4) is the superposition of two partial fluxes of type (5.1), F_all=F_H,ext+F_He,ext $F_{all} = F_{H, ext} + F_{He, ext}$ , which are extensions of the low-energy proton and helium fluxes F_{H, low} and F_{He, low} in (5.3):

equation5.5

\begin{matrix} F_{H, ext} (p) & = & \frac{A_{Φ, H} p^{3 + κ}}{\sqrt{p^{2} + m_{H}^{2}}} \frac{{(1 + {(p / c_{2 H})}^{γ_{2 H} / δ_{2 H}})}^{δ_{2 H}}}{{(1 + {(p / c_{1 H})}^{γ_{1 H} / | δ_{1 H} |})}^{| δ_{1 H} |}} \\ \times \frac{1}{{(1 + {(p / c_{3})}^{γ_{3} / | δ_{3} |})}^{| δ_{3} |}} \\ \times \frac{{(1 + {(p / c_{5})}^{γ_{5} / δ_{5}})}^{δ_{5}}}{{(1 + {(p / c_{4})}^{γ_{4} / | δ_{4} |})}^{| δ_{4} |}} \frac{{(1 + {(p / c_{7})}^{γ_{7} / δ_{7}})}^{δ_{7}}}{{(1 + {(p / c_{6})}^{γ_{6} / | δ_{6} |})}^{| δ_{6} |}} e^{- β p}, \end{matrix}

Turn MathJax on

and the helium extension F_{He, ext}(p) is obtained by replacing the subscript H by He. These partial fluxes are depicted in Figs. 3 and 4 as dashed and dotted curves. The factors in (5.5) defined by the break energies c_1H

c_{1 H}

, c_2H

c_{2 H}

and c_k

c_{k}

, k=3,⋯,7

k = 3, \dots, 7

, constitute the spectral kernel h_{H, ext} of F_{H, ext}, see (5.1),

equation5.6

h_{H, ext} (p) = \prod_{j = 1}^{2} {(1 + {(p / c_{j H})}^{γ_{j H} / | δ_{j H} |})}^{δ_{j H}} \prod_{k = 3}^{7} {(1 + {(p / c_{k})}^{γ_{k} / | δ_{k} |})}^{δ_{k}},

Turn MathJax on

and analogously for h_{He, ext}. The numeration of the factors is by increasing spectral break, c_1H<<c_2H<<c_k<<c_k+1

c_{1 H} < < c_{2 H} < < c_{k} < < c_{k + 1}

. The exponential in (5.5) corresponds to

e^{- β \sqrt{p^{2} + m^{2}}}

, β=1/k_BT

β = 1 / k_{B} T

, in (5.1); since the cutoff temperature k_BT is very high, see Table 1, we have dropped the mass square in the exponent. In the first ratio in (5.5), the nuclear mass is essential at low p, in the mildly relativistic regime.

The partial fluxes F_{H, ext} and F_{He, ext} fit the hydrogen and helium spectra in the low energy range $1 GeV - 10^{4} GeV$ , where they approximate the flux limits (5.3). At higher energies, in the $10^{4} GeV - 10^{11} GeV$ range, there are admixtures of heavier nuclei to both components, as their sum F_H,ext+F_He,ext $F_{H, ext} + F_{He, ext}$ constitutes the flux F_all in (5.4) obtained from a fit to the all-particle spectrum. These admixtures of intermediate and heavy nuclei can even dwarf the hydrogen and helium contributions to the all particle flux F_H,ext+F_He,ext $F_{H, ext} + F_{He, ext}$ over sections of the high-energy range, see the remarks after (6.7) and Refs. [19], [20], [21], [22] and [23]. The CREAM flux points define the proton and helium spectra in the $10^{4} GeV - 10^{5} GeV$ interval, although with rather large error bars, see Fig. 4. Comparing these points to the flux extensions F_{H, ext} and F_{He, ext}, there is already a small contribution from heavier nuclei noticeable in this interval, especially in the proton extension F_{H, ext}.

We have split the all-particle flux density (5.4) into two partial fluxes, one depending on the hydrogen mass and one on the helium mass, see (5.5). The masses of heavier nuclei contributing to the flux densities F_{H, ext} and F_{He, ext} do not enter, as the all-particle spectrum in the low-energy range essentially consists of protons and helium and, at high energies above 10⁴ GeV, the rest masses are negligible, since the squared mass-energy ratio m²/p² is very small even for iron nuclei. Above 10⁴ GeV, the momentum variable can be identified with energy for the same reason. In Section 6, we will use the partial fluxes F_{H, ext} and F_{He, ext} to decompose the all-particle flux F_all into a series of spectral peaks and to estimate their thermodynamic parameters.

6. All-particle spectrum as additive superposition of spectral peaks

6.1. Partial flux densities defined by spectral breaks

We start by splitting the factorizing flux density F_{H, ext} in (5.5) into partial fluxes F_{H, ext,i} generating four consecutive, exponentially decaying spectral peaks, $F_{H, ext} = \sum_{i = 1}^{4} F_{H, ext, i}$ . The decay exponent of F_{H, ext,i} is denoted by β_iH, and the peaks are ordered by decreasing decay exponent, β_1H ≥ β_2H ≥ β_3H ≥ β_4H. These exponents are free parameters, except for β_4H which coincides with the temperature parameter of the all-particle flux, β_4H=β=1/(k_BT) $β_{4 H} = β = 1 / (k_{B} T)$ , inferred from the wideband fit. The exponents will be specified in (6.7) by identifying them with spectral breaks.

The partial fluxes constituting the hydrogen extension F_{H, ext} in (5.5) are obtained iteratively. The flux density F_{H, ext,1} generating the first peak is found by dropping all numerators of F_{H, ext} defined by spectral breaks, namely the factors

equation6.1

\begin{matrix} f_{2} & = & {(1 + {(p / c_{2 H})}^{γ_{2 H} / δ_{2 H}})}^{δ_{2 H}}, f_{3} = {(1 + {(p / c_{5})}^{γ_{5} / δ_{5}})}^{δ_{5}}, \\ f_{4} & = & {(1 + {(p / c_{7})}^{γ_{7} / δ_{7}})}^{δ_{7}}, \end{matrix}

Turn MathJax on

and by replacing β → β_1H in the exponential in (5.5):

equation6.2

F_{H, ext, 1} = f_{0} f_{1} e^{- β_{1 H} p}, f_{0} = \frac{A_{Φ, H} p^{3 + κ}}{\sqrt{p^{2} + m_{H}^{2}}},

Turn MathJax on

equation6.3

\begin{matrix} f_{1} & = & \frac{1}{{(1 + {(p / c_{1 H})}^{γ_{1 H} / | δ_{1 H} |})}^{| δ_{1 H} |}} \frac{1}{{(1 + {(p / c_{3})}^{γ_{3} / | δ_{3} |})}^{| δ_{3} |}} \\ \times \frac{1}{{(1 + {(p / c_{4})}^{γ_{4} / | δ_{4} |})}^{| δ_{4} |}} \frac{1}{{(1 + {(p / c_{6})}^{γ_{6} / | δ_{6} |})}^{| δ_{6} |}} . \end{matrix}

Turn MathJax on

The partial fluxes defining the second, third and fourth peak are obtained by successively restoring the numerators (6.1), replacing the decay exponents β_iH→β_(i+1)H $β_{i H} \to β_{(i + 1) H}$ and subtracting the preceding peaks:

equation6.4

F_H,ext,2=f₀f₁f₂e^−β_2Hp−F_H,ext,1=f₀f₁(f₂−e^{−(β_1H−β_2H)p})e^−β_2Hp,

F_{H, ext, 2} = f_{0} f_{1} f_{2} e^{- β_{2 H} p} - F_{H, ext, 1} = f_{0} f_{1} (f_{2} - e^{- (β_{1 H} - β_{2 H}) p}) e^{- β_{2 H} p},

Turn MathJax on

equation6.5

\begin{matrix} F_{H, ext, 3} & = & f_{0} f_{1} f_{2} f_{3} e^{- β_{3 H} p} - \sum_{i = 1}^{2} F_{H, ext, i} \\ = & f_{0} f_{1} f_{2} (f_{3} - e^{- (β_{2 H} - β_{3 H}) p}) e^{- β_{3 H} p}, \end{matrix}

Turn MathJax on

equation6.6

\begin{matrix} F_{H, ext, 4} & = & f_{0} f_{1} f_{2} f_{3} f_{4} e^{- β_{4 H} p} - \sum_{i = 1}^{3} F_{H, ext, i} \\ = & f_{0} f_{1} f_{2} f_{3} (f_{4} - e^{- (β_{3 H} - β_{4 H}) p}) e^{- β_{4 H} p} . \end{matrix}

Turn MathJax on

The subtraction of the preceding peaks does not violate positivity, as the second identities in (6.4)–(6.6) show, since the exponents defining the factors f_i=2,3,4 $f_{i = 2, 3, 4}$ in (6.1) are positive and the decay exponents satisfy β_iH≥β_(i+1)H $β_{i H} \geq β_{(i + 1) H}$ . By adding the four fluxes F_{H, ext,i} defined in (6.2) and (6.4)–(6.6), we recover flux density F_{H, ext} in (5.5), provided that the smallest decay exponent β_4H is identified with the temperature parameter β of F_{H, ext}; the decay exponents of the first three peaks drop out in $F_{H, ext} = \sum_{i = 1}^{4} F_{H, ext, i}$ .

The decomposition of the helium flux extension into spectral peaks is performed analogously, $F_{He, ext} = \sum_{i = 1}^{4} F_{He, ext, i}$ ; we just need to replace the subscript H by He in the flux densities (6.2) and (6.4)–(6.6) and in the factors f_{1, 2} defined in (6.1) and (6.3). The four spectral peaks F_all,i in Fig. 8 constituting the all-particle flux $F_{all} = \sum_{i = 1}^{4} F_{all, i}$ are obtained by adding the corresponding peaks of the extended hydrogen and helium components, F_all,i=F_H,ext,i+F_He,ext,i $F_{all, i} = F_{H, ext, i} + F_{He, ext, i}$ , see the beginning of Section 5.2.

Fig. 8.

All-particle energy spectrum as additive superposition of spectral peaks. Data points as in Fig. 3. The solid curve F_all is the χ² fit of the all-particle flux density (5.4), see also the caption to Fig. 3. The iterative decomposition of this flux density into spectral peaks is explained in Section 6.1. The four strongly overlapping peaks F_all,i (dotted/dashed curves) add up to the all-particle flux $F_{all} = \sum_{i = 1}^{4} F_{all, i}$ . The spectral breaks c_k listed in Table 1 are indicated on the abscissa. All four peaks terminate in exponential decay, see (6.7): F_all,2 is exponentially cut at the knee c₄ and F_all,3 at the second knee c₆; the cutoff of the low-energy peak F_all,1 is defined by the spectral breaks c_2H $c_{2 H}$ and c_2He $c_{2 He}$ in the hydrogen and helium components of the all-particle flux. The double peak F_all,4 coalesces with the exponentially decaying all-particle flux above the ankle c₇. The thermodynamic parameters of the total flux F_all and the partial fluxes F_all,i generating the peaks are recorded in Table 2.

Figure options

We identify the decay exponents β_iH and β_iHe (defining the exponential cutoff of the peaks F_all,i=1,2,3 $F_{all, i = 1, 2, 3}$ ) with spectral breaks listed in Table 1,

equation6.7

β_{1 H} = \frac{1}{c_{2 H}}, β_{1 He} = \frac{1}{c_{2 He}}, β_{2 H} = β_{2 He} = \frac{1}{c_{4}}, β_{3 H} = β_{3 He} = \frac{1}{c_{6}} .

Turn MathJax on

The exponential decay of the fourth peak F_all,4, which is actually a double peak centered at the ankle c₇, is determined by the cutoff temperature of the all-particle flux, β_4H=β_4He=1/(k_BT) $β_{4 H} = β_{4 He} = 1 / (k_{B} T)$ , see Table 1. The knee c₄ and the second knee c₆ emerge at the intersections of overlapping fluxes, F_all,2, F_all,3 and F_all,3, F_all,4, respectively, see Fig. 8.

The first peak F_all,1 in Fig. 8 located below 10⁴ GeV has a light mass composition, as the all-particle spectrum below 10⁴ GeV consists predominantly of proton and helium nuclei. The second peak F_all,2 in Fig. 8 (which is exponentially cut at the knee c₄) has also a light composition below 10⁴ GeV for the same reason. Above 10⁴ GeV, the all-particle spectrum increasingly deviates from the combined CREAM proton and helium spectra, due to admixtures of heavier nuclei, see Section 5.2 and Fig. 4. Composition studies with Tibet-III indicate that the light component accounts for less than 30% of the all-particle flux in the knee region [19] and [20]. This suggests that the second peak F_all,2 has an increasingly heavy composition above 10⁴ GeV, containing at least intermediate-mass nitrogen-like nuclei.

Above the knee, the IceCube Collaboration derives a mean log mass which increases in the interval between 5 × 10⁶ GeV and 3 × 10⁸ GeV where it reaches a plateau and starts to decrease [3]. Composition studies with the Yakutsk array in the $10^{6} GeV - 10^{9} GeV$ range also indicate a maximum of the mean log mass at around 10⁸ GeV followed by a substantial decrease in the subsequent decade [21]. A heavy iron-like composition at the second knee around 10⁸ GeV has also been found by the KASCADE-Grande Collaboration [22] and [23]. Above 10⁹ GeV, a light composition is inferred with the Telescope Array [6], whereas a mixed composition, mainly light to intermediate, is suggested by the Auger Collaboration [24]. Accordingly, the fourth peak F_all,4 in Fig. 8 has a light or light-to-intermediate composition, whereas the third peak F_all,3, which covers the knee c₄ and is exponentially cut at the second knee c₆, has presumably a heavy composition.

Models of diffusive shock acceleration in supernova remnants (SNRs) suggest that protons in SNR shocks can reach energies up to the knee at c₄ ∼ 4 × 10⁶ GeV [25] and [26]. As the maximum energy obtainable by shock acceleration scales with the atomic charge number, the iron component can thus extend to energies of around 10⁸ GeV, reaching the second knee c₆. This is consistent with the above quoted composition studies, which suggest a heavy component extending to the second knee and declining above [21]. This heavy component is manifested as spectral peak F_all,3 in Fig. 8. The existence of high-energy ( ∼ 10⁵ GeV) protons in SNRs is also supported by the TeV γ-ray spectra of several SNRs (enumerated in Refs. [27] and [28]) which could stem from hadronic γ-ray production by pion decay, the pions being generated in collisions of high-energy protons with nuclei. The first three peaks in Fig. 8 could thus be produced by shock acceleration in SNRs. Another possible source mechanism is one-shot acceleration in the induced electric fields of pulsars [29].

The fourth peak F_all,4 in Fig. 8 constituting the ultra-high energy component of the all-particle spectrum is a double peak which, above the ankle c₇, coincides with the all-particle flux. Assuming an extragalactic origin of this peak (possible sources could be shock acceleration in active galactic nuclei, also γ-ray bursts or decaying super-massive particles [25]) and a proton-dominated composition, the exponential cutoff of this peak can be explained by scattering of protons by cosmic microwave background photons (energy loss via pion production, GZK cutoff [25]). In the case of a light-to-intermediate mass composition, the cutoff could be caused by acceleration limits of the sources. The double-peak structure of F_all,4 could be carved out by attenuation effects such as electron-positron pair production, again by scattering of protons by microwave photons, or photo-disintegration of heavier nuclei, also caused by interaction with microwave photons or the extragalactic background light [30].

The cutoff of the third peak F_all,3 also affects the shape of the fourth peak because of the subtraction in (6.6). If the peak F_all,3 is cut off at the ankle, β_3H=β_3He=1/c₇ $β_{3 H} = β_{3 He} = 1 / c_{7}$ , instead of at the second knee c₆ as in (6.7), then the double peak F_all,4 in Fig. 8 shrinks to a single peak with its maximum located above the ankle. The latter then emerges at the intersection of a broad Galactic peak F_all,3 and a small extragalactic peak F_all,4. However, a heavy mass component F_all,3 extending to the ankle is at odds with the mentioned composition studies (e.g. Ref. [21]), which indicate a cutoff of the heavy component at the second knee as shown in Fig. 8.

The multiplicative flux representation (5.4) is more efficient than an additive superposition of spectral peaks when performing wideband fits, especially if the source mechanisms are uncertain and the peaks are substantially overlapping as in Fig. 8. The iterative decomposition of the wideband spectrum into spectral peaks is subsequently carried out without invoking approximations.

6.2. Thermodynamic parameters of cosmic-ray fluxes

We will calculate the specific number and energy densities N and U, the excess entropy density H and the pressure P, based on the integral representations ( 2.8)–(2.10) and (2.12). (The specific densities are obtained by dropping the volume factors there.) By making use of the flux representation of the spectral number density in (2.5), we find

equation6.8

N [m^{- 3}] = \frac{4 π}{c} \int_{0}^{\infty} F (p) \frac{\sqrt{p^{2} + m^{2}}}{p^{1 + κ}} d p,

Turn MathJax on

equation6.9

U [GeV / m^{3}] = \frac{4 π}{c} \int_{0}^{\infty} F (p) \frac{p^{2} + m^{2}}{p^{1 + κ}} d p,

Turn MathJax on

equation6.10

\frac{H}{k_{B}} [m^{- 3}] = - \frac{4 π}{c} \int_{0}^{\infty} F (p) \frac{\sqrt{p^{2} + m^{2}}}{p^{1 + κ}} \log h (p) d p,

Turn MathJax on

equation6.11

P [GeV / m^{3}] = \frac{1}{3} \frac{4 π}{c} \int_{0}^{\infty} F (p) p^{1 - κ} d p,

Turn MathJax on

with c [m/s]. The specific density of the total entropy is assembled as S/k_B=(α+1)N+βU+H/k_B

S / k_{B} = (α + 1) N + β U + H / k_{B}

, see (2.13), with fugacity and temperature parameters α=log(1.573×¹⁰⁵³s/A_Φ)

α = \log (1.573 \times 10^{53} s / A_{Φ})

and β=1/k_BT[GeV]

β = 1 / k_{B} T [GeV]

, where A_Φ

A_{Φ}

is the flux amplitude and s the spin degeneracy of the gas particles as in ( 2.4).

We denote the specific densities and the pressure of the extended hydrogen flux F_{H, ext} in (5.5) by N_{H, ext}, U_{H, ext}, H_{H, ext} and P_{H, ext}, and calculate these quantities with F_{H, ext} and h_{H, ext} in (5.6) substituted into the integral representations (6.8)–(6.11). The total entropy density of this partial flux is S_H,ext/k_B=(α_H+1)N_H,ext+βU_H,ext+H_H,ext/k_B $S_{H, ext} / k_{B} = (α_{H} + 1) N_{H, ext} + β U_{H, ext} + H_{H, ext} / k_{B}$ , with fugacity parameter α_H=log(1.573×¹⁰⁵³s_H/A_Φ,H) $α_{H} = \log (1.573 \times 10^{53} s_{H} / A_{Φ, H})$ , spin degeneracy s_H=2 $s_{H} = 2$ and flux amplitude A_Φ,H $A_{Φ, H}$ listed in Table 1.

The thermodynamic parameters of the extended helium flux F_{He, ext} are obtained analogously and denoted by N_{He, ext}, U_{He, ext}, H_{He, ext}, P_{He, ext} and S_{He, ext}; the fugacity parameter α_He $α_{He}$ is calculated with spin degeneracy s_He=1 $s_{He} = 1$ and amplitude A_Φ,He $A_{Φ, He}$ in Table 1. The specific densities N_all $N_{all}$ , U_all $U_{all}$ , H_all $H_{all}$ , S_all and the pressure P_all $P_{all}$ of the all-particle flux F_all in (5.4) are listed in the first row of Table 2, obtained by adding the respective densities and pressure variables of the partial fluxes F_{H, ext} and F_{He, ext}, e.g. N_all=N_H,ext+N_He,ext $N_{all} = N_{H, ext} + N_{He, ext}$ .

Table 2.

Thermodynamic parameters of the all-particle flux F_all and the partial fluxes F_all,i in Fig. 8. N denotes the specific number density, U the specific internal energy density, H/k_B and S/k_B are the specific densities of excess entropy and total entropy (rescaled with the Boltzmann constant) and P is the pressure, see Section 6.2. In the first row, these quantities are listed for the all-particle flux F_all; their calculation is explained after (6.11). The total entropy density S/k_B of the all-particle flux is assembled with the extended hydrogen and helium fluxes F_{H, ext} and F_{He, ext} in (5.5) and their fugacity parameters α_H=114.1 $α_{H} = 114.1$ and α_He=117.8 $α_{He} = 117.8$ . As the fugacity is very small, z_H,He=exp(−α_H,He) $z_{H, He} = \exp (- α_{H, He})$ , the classical flux densities used in the spectral fits are applicable over the complete spectral range, see the end of Section 4. Also recorded are the thermodynamic parameters of the partial fluxes F_all,i generating the four spectral peaks in Fig. 8; these parameters are assembled as outlined after (6.12) and add up to the corresponding quantities of the all-particle flux F_all.

	N [m⁻³] $[m^{- 3}]$	U $[GeV m^{- 3}]$	H/k_B[m⁻³] $[m^{- 3}]$	S/k_B[m⁻³] $[m^{- 3}]$	P $[GeV m^{- 3}]$
All-part., Fall	1.238×10−4	5.935×10−4	5.937×10−4	1.490×10−2	1.705×10−4
Peak 1, Fall, 1	1.225×10−4	5.544×10−4	5.825×10−4	1.475×10−2	1.576×10−4
Peak 2, Fall, 2	1.217×10−6	3.892×10−5	1.114×10−5	1.509×10−4	1.283×10−5
Peak 3, Fall, 3	1.332×10−10	1.158×10−7	1.353×10−9	1.667×10−8	3.860×10−8
Peak 4, Fall, 4	4.963×10−12	8.425×10−9	5.045×10−11	6.212×10−10	2.808×10−9

Full-size table

Table options

Also recorded in Table 2 are the specific densities N_all,i $N_{all, i}$ , U_all,i $U_{all, i}$ , H_all,i, S_all,i and the pressure P_all,i $P_{all, i}$ of the partial fluxes generating the spectral peaks F_all,i, i=1,⋯,4 $i = 1, \dots, 4$ , depicted in Fig. 8. To assemble these variables, we first express the partial fluxes F_{H, ext,i} in (6.2) and (6.4)–(6.6) by flux densities which are structured as in (2.3) and (2.6): F_iH=f₀h_iHe^−β_iHp $F_{i H} = f_{0} h_{i H} e^{- β_{i H} p}$ , with spectral kernels

equation6.12

h_{1 H} = f_{1}, h_{2 H} = f_{1} f_{2}, h_{3 H} = f_{1} f_{2} f_{3}, h_{4 H} = f_{1} f_{2} f_{3} f_{4} .

Turn MathJax on

The factors f_i=0,…,4 $f_{i = 0, \dots, 4}$ are defined in (6.1)–(6.3) and the exponents β_iH in (6.7). The four partial fluxes F_H,ext,i=F_iH−F_(i−1)H $F_{H, ext, i} = F_{i H} - F_{(i - 1) H}$ (with F_0H=0 $F_{0 H} = 0$ ) add up to the total flux F_{H, ext} in (5.5), see also after (6.6).

We write N(F_iH), U(F_iH), H(F_iH, h_iH) and P(F_iH) for the specific densities (6.8)–(6.11), with F_iH and h_iH substituted for the flux F and the spectral kernel h in the integrands. The specific number density N_{H, ext,i} of the partial flux F_{H, ext,i} is calculated as N_H,ext,i=N(F_iH)−N(F_(i−1)H) $N_{H, ext, i} = N (F_{i H}) - N (F_{(i - 1) H})$ , and analogously U_{H, ext,i}, H_{H, ext,i} and P_{H, ext,i}. The specific total entropy density of F_{H, ext,i} reads S_H,ext,i=S(F_iH)−S(F_(i−1)H) $S_{H, ext, i} = S (F_{i H}) - S (F_{(i - 1) H})$ with

equation6.13

S(F_iH)/k_B=(α_H+1)N(F_iH)+β_iHU(F_iH)+H(F_iH,h_iH)/k_B.

S (F_{i H}) / k_{B} = (α_{H} + 1) N (F_{i H}) + β_{i H} U (F_{i H}) + H (F_{i H}, h_{i H}) / k_{B} .

Turn MathJax on

The fugacity parameter α_H is calculated as explained after (6.11), also see the caption to Table 2. The specific number densities of the partial fluxes F_{H, ext,i} add up to the number density N_{H, ext} of the total flux F_{H, ext} already calculated after (6.11), $N_{H, ext} = \sum_{i = 1}^{4} N_{H, ext, i}$ , and analogously for energy, excess entropy, total entropy and pressure.

The foregoing derivations also apply to the extended helium flux F_{He, ext}, with subscript H replaced by He. The specific densities N_all,i $N_{all, i}$ , U_all,i $U_{all, i}$ , H_all,i $H_{all, i}$ , S_all,i and the pressure P_all,i $P_{all, i}$ of the partial fluxes F_all,i defining the spectral peaks in Fig. 8 are found by adding the corresponding densities of the extended hydrogen and helium components, N_all,i=N_H,ext,i+N_He,ext,i $N_{all, i} = N_{H, ext, i} + N_{He, ext, i}$ , etc. The all-particle flux F_all is the sum of the partial fluxes F_all,i, as pointed out after (6.6). Accordingly, the specific number density N_all of the all-particle flux calculated after (6.11) can also be obtained by adding the number densities N_all,i of the four spectral peaks F_all,i in Fig. 8, see Table 2, and this holds as well for the specific energy and entropy densities and the pressure.

7. Conclusion

Cosmic-ray nuclei constitute a relativistic multi-component gas in stationary non-equilibrium. We have found the partition function of this gas, defined by a spectral flux density which accurately reproduces the spectral range currently accessible, from 1 GeV to 10¹¹ GeV. Spectral maps of the all-particle energy spectrum in log-log representation exhibit, besides knee and ankle, several weaker, widely spaced spectral breaks. In the least-squares fits, we used a multiplicative spectral representation of the flux density, each factor being determined by a spectral break. The wideband fit can iteratively be performed over an extended energy range by adding factors to the spectral kernel of the flux density, which is a multiply broken power-law function defined by successive spectral breaks, see (2.6) and Section 5.1. The factorizing all-particle flux density admits an additive decomposition into four partial fluxes, each of them generating a spectral peak which is exponentially cut off at a spectral break.

We have chosen to model cosmic-ray nuclei as stationary non-equilibrium system and to extract the number density and partition function thereof from available spectra. In this way, we can avoid making specific assumptions on source distributions and interaction mechanisms producing the spectral peaks, which otherwise would have to be quantified in models based on kinetic equations. Instead, the spectral peaks constituting the wideband spectrum are iteratively reconstructed from the fitted flux density. The qualitative mass composition of the peaks and possible Galactic and cosmological source mechanisms are briefly discussed at the end of Section 6.1.

The spectral number density (2.1) employed in the spectral fits can be traced back to the phase-space probability distribution studied in Section 3. Once this distribution is known, one can assemble the partition function of the relativistic nuclei and calculate the thermodynamic variables, as done in Sections 3 and 6.2. We quantized the classical spectral number density in Fermi and Bose statistics, depending on the nuclear spin, see (4.7). In either case, the classical density (2.1) approximates the quantum density in the limit of vanishing fugacity. As the fugacity of cosmic-ray fluxes is extremely small, see the caption to Table 2, the classical formalism of Section 2 applies over the complete spectral range.

Figs. 1 and 2 depict spectral fits to the AMS-02 [7] and [8] and CREAM [9] proton and helium spectra. Both of these low-energy spectra consist of a peak followed by a trough and a slightly increasing slope. The peak and trough are treated as spectral breaks defining the factorizing kernel of the flux density employed in the fits, see (5.3) and Table 1. By adding the proton and helium flux densities, we obtained an approximation to the all-particle flux in the $1 GeV - 10^{4} GeV$ interval, since heavier nuclei contribute only a negligible fraction, roughly one percent of the all-particle flux, in this interval.

The low-energy flux density is then iteratively extended to higher energies by fitting the all-particle spectrum defined by the Tibet-III [1], IceTop-73 [2] and [3], KASCADE-Grande [4], Auger [5] and Telescope Array [6] data sets, see Fig. 3. The high-energy fit in the $10^{5} GeV - 10^{11} GeV$ interval is performed with the analytic flux density defined in (5.4) and Table 1, which is a multiply broken power-law function with smooth transitions at the spectral breaks and an exponential ultra-high energy cutoff by the Boltzmann factor. This flux density also provides an interpolation of the all-particle spectrum in the $10^{4} GeV - 10^{5} GeV$ interval. The crossover from the low- to the high-energy spectrum is shown in Fig. 4. The spectral breaks of the high-energy component are visible in Fig. 5, Fig. 6 and Fig. 7 (which are close-ups of the wideband spectrum in Fig. 3), in particular the two soft breaks between knee and ankle resulting in a slight shift of the spectral slope, see Fig. 6. The exponentially decaying spectral tail is shown in Fig. 7; the temperature parameter in the Boltzmann factor of this stationary non-equilibrium gas of relativistic nuclei is inferred from the spectral fit. The all-particle flux density of the $1 GeV - 10^{11} GeV$ wideband spectrum is an additive superposition of four strongly overlapping spectral peaks with crossovers at the spectral breaks as depicted in Fig. 8. The iterative flux decomposition into exponentially decaying spectral peaks is explained in Section 6.1. The thermodynamic parameters of the total flux and of the partial fluxes constituting the spectral peaks are calculated in Section 6.2 and recorded in Table 2.

Appendix A. Iterative wideband fitting

The χ² functional of flux density F_all(p) in ( 5.4) depends on sixteen parameters, (c_k, γ_k, δ_k ), k=3,…,7 $k = 3, \dots, 7$ , and k_BT , which have to be varied simultaneously to find the minimum. Therefore, a good initial guess for the flux density is required, valid over the complete spectral range $1 GeV - 10^{11} GeV$ . As the spectral breaks are widely separated, c_k<<c_k+1 $c_{k} < < c_{k + 1}$ , see Table 1 and Fig. 3, Fig. 4, Fig. 5, Fig. 6 and Fig. 7, it is possible to iteratively find an initial guess for F_all(p) by varying six parameters in each step (seven in the last one).

In the analytic flux density F_all(p) in ( 5.4), we drop the exponential and the three factors defined by the last three breaks (c_k, γ_k, δ_k ), k=5,6,7 $k = 5, 6, 7$ , and denote this density by F_all(p; 3, 4). The χ²(3, 4) functional of F_all(p; 3, 4) is defined by using the Tibet-III and IceTop-73 flux points with energies below the spectral break at c₅ ≈ 1.2 × 10⁷ GeV, see Fig. 6. (The precise location of the spectral breaks c_k , k=3,…,7 $k = 3, \dots, 7$ , is not known at this stage, but can be estimated by plotting the properly rescaled flux points, see the remark after (2.4) and Fig. 8.) By minimizing χ²(3, 4), we obtain preliminary values for the parameters (c_k, γ_k, δ_k ), k=3,4. $k = 3, 4 .$ This is the first step of the iteration.

The second step is based on density F_all(p ; 3; 4, 5) defined in like manner by dropping the last two factors indexed k=6,7 $k = 6, 7$ and the exponential in (5.4). For the parameters (c₃, γ₃, δ₃), we substitute the values obtained from the previous χ²(3, 4) fit. The χ²(4, 5) functional of F_all(p; 3; 4, 5) thus depends on the parameters (c_k, γ_k, δ_k ), k=4,5 $k = 4, 5$ . The flux points included in χ²(4, 5) have energies up to the spectral break c₆ ≈ 1.1 × 10⁸ GeV, see Fig. 6. (That is, we extend the data set used in the χ²(3, 4) fit by adding flux points located between the spectral breaks c₅ and c₆.) As a consistency check, the values obtained for (c₄, γ₄, δ₄) in the χ²(3, 4) and χ²(4, 5) fits should be similar, which is likely to be the case if the breaks are sufficiently separated, c_k<<c_k+1 $c_{k} < < c_{k + 1}$ .

The third step of the iteration is performed with density F_all(p ; 3, 4; 5, 6) defined by dropping the factor indexed k=7 $k = 7$ and the exponential in (5.4). We substitute the parameters (c₃, γ₃, δ₃) obtained from the χ²(3, 4) fit and (c₄, γ₄, δ₄) from the χ²(4, 5) fit. The fitting parameters are now (c_k, γ_k, δ_k ), k=5,6 $k = 5, 6$ . In the χ²(5, 6) functional of density F_all(p; 3, 4; 5, 6), we include flux points with energies up to the spectral break c₇ ≈ 4.9 × 10⁹ GeV, see Fig. 7. For consistency, the parameters (c₅, γ₅, δ₅) obtained from the χ²(4, 5) and χ²(5, 6) fits should be similar.

In the fourth and last step, density F_all(p; 3, 4, 5; 6, 7, k_BT) is identical with F_all(p) in ( 5.4). We substitute the parameters (c₃, γ₃, δ₃) obtained from the χ²(3, 4) fit, (c₄, γ₄, δ₄) from the χ²(4, 5) fit and (c₅, γ₅, δ₅) from the χ²(5, 6) fit. The fitting parameters are now (c_k, γ_k, δ_k ), k=6,7 $k = 6, 7$ , and the temperature k_BT in the Boltzmann factor. The χ²(6, 7, k_BT) functional includes all high-energy flux points, see Fig. 5. As a consistency check, the parameters (c₆, γ₆, δ₆) obtained from the χ²(5, 6) and χ²(6, 7, k_BT) fits should be similar; the latter ones are used as initial guess for (c₆, γ₆, δ₆) in the final χ² fit.

In this way, we have found an initial guess for the fitting parameters (c_k, γ_k, δ_k ), k=3,…,7 $k = 3, \dots, 7$ and k_BT of flux density F_all(p) in ( 5.4), which can be employed in the χ² fit described after (5.4), where all parameters are simultaneously varied. This iterative method of wideband fitting over an extended spectral range covering several spectral breaks is applicable to flux densities based on multiply broken power laws such as (2.6).

References

- [1]
- M. Amenomori, et al.
- Astrophys. J., 678 (2008), p. 1165
- Full Text via CrossRef
  |
  View Record in Scopus
  |
  Citing articles (102)
- [2]
- IceCube Collaboration
- Phys. Rev. D, 88 (2013), Article 042004
- [3]
- IceCube Collaboration
- Proceedings of Science (2015) PoS(ICRC2015)334
- [4]
- KASCADE-Grande Collaboration
- Astropart. Phys., 36 (2012), p. 183
- [5]
- Pierre Auger Collaboration
- Proceedings of Science (2015) PoS(ICRC2015)271
- [6]
- Telescope Array Collaboration
- Proceedings of Science (2015) PoS(ICRC2015)035
- [7]
- AMS Collaboration
- Phys. Rev. Lett., 114 (2015), Article 171103
- [8]
- AMS Collaboration
- Phys. Rev. Lett., 115 (2015), Article 211101
- [9]
- Y.S. Yoon, et al.
- Astrophys. J., 728 (2011), p. 122
- Full Text via CrossRef
  |
  View Record in Scopus
  |
  Citing articles (104)
- [10]
- R. Tomaschitz
- Physica A, 385 (2007), p. 558
- Article
  |
  PDF (276 K)
  |
  View Record in Scopus
  |
  Citing articles (14)
- [11]
- R. Tomaschitz
- Physica A, 394 (2014), p. 110
- Article
  |
  PDF (858 K)
  |
  View Record in Scopus
  |
  Citing articles (4)
- [12]
- R. Tomaschitz
- Physica A, 451 (2016), p. 456
- Article
  |
  PDF (903 K)
  |
  View Record in Scopus
- [13]
- R. Tomaschitz
- Europhys. Lett., 104 (2013), p. 19001
- Full Text via CrossRef
  |
  View Record in Scopus
  |
  Citing articles (1)
- [14]
- R. Tomaschitz
- Physica A, 387 (2008), p. 3480
- Article
  |
  PDF (693 K)
  |
  View Record in Scopus
  |
  Citing articles (14)
- [15]
- R. Tomaschitz
- Physica B, 405 (2010), p. 1022
- Article
  |
  PDF (370 K)
  |
  View Record in Scopus
  |
  Citing articles (10)
- [16]
- A. Chandra, et al.
- Proceedings of Science (2015) PoS(ICRC2015)235
- [17]
- PAMELA Collaboration
- Astrophys. J., 818 (2016), p. 68
- [18]
- BESS-Polar Collaboration
- Astrophys. J., 822 (2016), p. 65
- [19]
- M. Amenomori, et al.
- Adv. Space Res., 47 (2011), p. 629
- Article
  |
  PDF (782 K)
  |
  View Record in Scopus
  |
  Citing articles (13)
- [20]
- M. Amenomori, et al.
- Astrophys. Space Sci. Trans., 7 (2011), p. 15
- Full Text via CrossRef
  |
  View Record in Scopus
  |
  Citing articles (7)

- [21]
- E.G. Berezhko, S.P. Knurenko, L.T. Ksenofontov
- Astropart. Phys., 36 (2012), p. 31
- Article
  |
  PDF (381 K)
  |
  View Record in Scopus
  |
  Citing articles (13)
- [22]
- KASCADE-Grande Collaboration
- Adv. Space Res., 53 (2014), p. 1456
- [23]
- KASCADE-Grande Collaboration
- Proceedings of Science (2015) PoS(ICRC2015)359
- [24]
- Pierre Auger Collaboration
- Phys. Rev. D, 90 (2014), Article 122006
- [25]
- A. Letessier-Selvon, T. Stanev
- Rev. Mod. Phys., 83 (2011), p. 907
- Full Text via CrossRef
  |
  View Record in Scopus
  |
  Citing articles (52)
- [26]
- L. O'C. Drury
- Astropart. Phys., 39 (2012), p. 52
- [27]
- F.A. Aharonian
- Astropart. Phys., 43 (2013), p. 71
- Article
  |
  PDF (1186 K)
  |
  View Record in Scopus
  |
  Citing articles (17)
- [28]
- P. Slane, et al.
- Space Sci. Rev., 188 (2015), p. 187
- Full Text via CrossRef
  |
  View Record in Scopus
  |
  Citing articles (7)
- [29]
- K.V. Ptitsyna, S.V. Troitsky
- Phys. Usp., 53 (2010), p. 691
- Full Text via CrossRef
  |
  View Record in Scopus
  |
  Citing articles (23)
- [30]
- R. Aloisio, V. Berezinsky, A. Gazizov
- Astropart. Phys., 39 (2012), p. 129
- Article
  |
  PDF (1055 K)
  |
  View Record in Scopus
  |
  Citing articles (16)
- [31]
- R. Tomaschitz
- Europhys. Lett., 106 (2014), p. 39001
- Full Text via CrossRef
  |
  View Record in Scopus
  |
  Citing articles (3)
- [32]
- R. Tomaschitz
- Phys. Lett., 378 (2014), p. 2337
- Article
  |
  PDF (966 K)
  |
  View Record in Scopus
  |
  Citing articles (3)
- [33]
- R. Tomaschitz
- Phys. Lett., 378 (2014), p. 2915
- Article
  |
  PDF (808 K)
  |
  View Record in Scopus
  |
  Citing articles (2)
- [34]
- R. Tomaschitz
- J. High Energy Astrophys., 8 (2015), p. 10
- Article
  |
  PDF (1480 K)
  |
  View Record in Scopus
  |
  Citing articles (1)

Partition function and thermodynamic parameters of the all-particle cosmic-ray flux

Highlights

Abstract

Keywords

1. Introduction

2. Flux densities with factorizing spectral kernel

3. Partition function and entropy

4. Quantization of stationary non-equilibrium densities

5. Proton and helium fluxes and the all-particle energy spectrum

5.1. Product representation of spectral flux densities over multiple spectral breaks

5.2. Decomposition of the all-particle spectrum into partial fluxes

6. All-particle spectrum as additive superposition of spectral peaks

6.1. Partial flux densities defined by spectral breaks

6.2. Thermodynamic parameters of cosmic-ray fluxes

7. Conclusion

Appendix A. Iterative wideband fitting

References