The high–and ultra-high–energy cosmic-ray flux (1014-1020 eV$10^{14}\text{-}10^{20}\ \text{eV}$) [1] has been measured by air shower arrays such as Tibet-III [2], KASCADE [3], KASCADE-Grande [4], IceCube [5], HiRes-I and II [6], the Telescope Array [7,8] and the Pierre Auger Observatory [9], as well as by the HEGRA [10], AKENO [11], GAMMA [12], Yakutsk [13] and AGASA [14] arrays. Here, the goal is to give a statistical description of the broadband spectrum by means of classical ultra-relativistic power-law distributions. We fit the all-particle spectrum with the flux density of a dilute three-component plasma (non-equilibrated, with components at different temperature) and locate the "knees" and "ankles" [1,15–17] of the spectrum. The knees are identified as the peaks of the power-law distributions of the plasma components, and the ankles turn out to be the minima in the cross-over regimes of the partial fluxes. These extrema clearly emerge in the E³-scaled flux representation employed in the wideband fit.

We sketch the basic thermodynamic formalism of classical ultra-relativistic power-law ensembles and nonthermal mixtures thereof, starting with the spectral number density parametrized with the Lorentz factor. We then define the specific number count and the internal energy of the gas components, as well as their entropy density. After this overview, we discuss the phase-space measure and probability density of nonthermal power-law ensembles, and derive the grand canonical partition function and the extensive entropy functional of multi-component power-law mixtures involving different particle species at different temperature.

We discuss the practical aspects of broadband spectral fitting with flux densities of classical unequilibrated (stationary non-equilibrium) power-law mixtures in the ultra-relativistic regime, and perform the spectral fit of the all-particle cosmic-ray spectrum above 1014 eV$10^{14}\ \text{eV}$. The spectral flux density of a power-law distribution admits a power-law ascent linear in the log-log flux representation, and a curved descending power-law slope. The curvature is caused by the exponential cutoff determined by the temperature parameter of the Boltzmann factor [18,19]. The broadband spectrum of the all-particle mixture consists of overlapping spectral peaks generated by the partial fluxes of the individual plasma components. By choosing a suitable flux representation (double-logarithmic E3F(E)$E^3F(E)$), we find secondary knees and ankles in the all-particle spectrum. We calculate the thermodynamic parameters of each gas component, and estimate the entropy production of high-energy cosmic rays.

An ultra-relativistic classical gas mixture composed of nonthermal power-law components is defined by the spectral number density

dρmix(γ)=∑i=1Mdρi(γ),dρi(γ)=m3isi2π2(ℏc)3e−βiγ−αiGi(γ)γ2−1−−−−−√γdγ,(1)

$$\begin{equation*} \begin{align*} \begin{equation} \begin{array}{c} \displaystyle\text{d}\rho_{\text{mix}}(\gamma)=\sum\limits_{i=1}^M \text{d}\rho_i (\gamma ),\\ \displaystyle\text{d}\rho_i (\gamma )=\frac{m_i^3 s_i }{2\pi^2( \hbar c)^3}\frac{\text{e}^{-\beta_i\gamma-\alpha_i}}{G_i(\gamma)}\sqrt{\gamma^2-1}\gamma\text{d}\gamma, \end{array}\end{equation} \end{align*} \tag{ 1 } \end{equation*}$$

parametrized by the Lorentz factor γ. The positive spectral functions Gi(γ)$G_i (\gamma )$ in the denominator of the power-law components dρi(γ)$\text{d}\rho_i (\gamma )$ are normalized as Gi(1)=1$G_i (1)=1$ and are to be determined from a spectral fit. An M-component mixture is composed of particle species (Zi,mi)$(Z_i ,m_i )$, i=1,…,M$i=1,\ldots ,M$, where m_i is the mass and Z_i the (positive or negative) integer charge number (in multiples of the electron charge e > 0) of the particles in the respective component (Zi,mi)$(Z_i ,m_i )$. The mass parameter m_i is a shortcut for m_i c², so that E=miγ$E=m_i \gamma$ is the particle energy in an ideal gas mixture. s_i is the spin multiplicity and βi=mi/(kBTi)$\beta_i =m_i /(k_{\text{B}} T_i )$ the dimensionless temperature parameter, with Boltzmann constant kB$k_{\text{B}}$. The mixture need not be equilibrated, admitting components at different temperature. The densities (1) are based on grand canonical partitions, allowing a fluctuating particle number via the dimensionless real fugacity parameter αi$\alpha_i$, which is related to the chemical potential μi$\mu_i$ of each gas component by αi=−βiμi/mi$\alpha_i =-\beta_i \mu_i /m_i$, and zi=e−αi$z_i=\text{e}^{-\alpha_i}$ is the fugacity. The ensemble average leading to (1) is discussed in the next section. The spectral function Gi(γ)$G_i (\gamma )$ is a linear combination of two power laws,

Gi(γ)=γδigi(γ)gi(1),gi(γ)=1+aiγσi,(2)

$$\begin{equation*} \begin{align*} \begin{equation} G_i(\gamma)=\gamma^{\delta_i}\frac{g_i(\gamma)}{g_i(1)},\quad g_i (\gamma )=1+a_i \gamma^{\sigma_i},\end{equation} \end{align*} \tag{ 2 } \end{equation*}$$

where δi$\delta_i$ and σi$\sigma_i$ are real power-law exponents and the amplitude a_i is positive. Spectral functions of this type are thermodynamically stable [18] and suitable to model wideband spectra [20–22]. Quantized power-law ensembles with ai=0$a_i =0$ have been studied in [23,24].

The total particle number Nmix=∑i=1MNi$N_{\text{mix}} =\sum\limits_{i=1}^M {N_i }$ and the internal energy Umix=∑i=1MUi$U_{\text{mix}} =\sum\limits_{i=1}^M {U_i}$ of a nonthermal gas mixture in a volume V are determined by the partial densities dρi(γ)$\text{d}\rho_i (\gamma)$ in (1),

Ni=V∫∞γcutidρi(γ),Ui=miV∫∞γcutiγdρi(γ).(3)

$$\begin{equation*} \begin{align*} \begin{equation} N_i =V\int_{\gamma_i^{\text{cut}}}^\infty {\text{d}\rho_i} (\gamma), \qquad U_i =m_i V\int_{\gamma_i^{\text{cut}}}^\infty {\gamma \text{d}\rho_i(\gamma)}.\end{equation} \end{align*} \tag{ 3 } \end{equation*}$$

The lower integration boundary is the cutoff Lorentz factor γcuti≥1$\gamma_i^{\text{cut}}\ge 1$, referring to particles of a given species (Zi,mi)$(Z_i ,m_i )$ with energies exceeding Ecuti=miγcuti$E_i^{\text{cut}} =m_i \gamma _i^{\text{cut}}$, and adjusted to the energy range of the available data sets. The ultra-relativistic limit of the power-law densities dρi(γ)$\text{d}\rho _i (\gamma)$ is obtained by replacing γ2−1−−−−−√→γ$\sqrt{\gamma^2-1} \to \gamma$ in (1) and is applicable if γcuti≫1$\gamma_i^{\text{cut}} \gg1$. The thermal relativistic Maxwell-Boltzmann distribution is recovered with δi=0$\delta_i =0$, gi(γ)=1$g_i (\gamma )=1$ and γcuti=1$\gamma_i^{\text{cut}} =1$ as lower integration boundary in (3).

The partition function of a gas component (Zi,mi)$(Z_i ,m_i)$ is denoted by Z(Zi,mi)$Z_{(Z_i ,m_i )}$. Since particle number and internal energy are related to the grand partition function by Ni=−(logZ(Zi,mi)),αi$N_i =-(\log Z_{(Z_i ,m_i)})_{,\alpha_i}$ and Ui=−mi(logZ(Zi,mi)),βi$U_i =-m_i (\log Z_{(Z_i ,m_i )} )_{,\beta_i }$, we can identify the partition function Zmix$Z_{\text{mix}}$ of a nonthermal mixture as logZmix=∑i=1Mξi$\log Z_{\text{mix}} =\sum\limits_{i=1}^M {\xi_i}$, where ξi=logZ(Zi,mi)$\xi_i =\log Z_{(Z_i,m_i)}$ and

ξi=V4πsim3i(2πℏc)3e−αi∫∞γcutie−βiγGi(γ)γ2−1−−−−−√γdγ.(4)

$$\begin{equation*} \begin{align*} \begin{equation} \xi_i =V\frac{4\pi s_i m_i^3 }{(2\pi \hbar c)^3}\text{e}^{-\alpha_i}\int_{\gamma_i^{\text{cut}}}^\infty \frac{\text{e}^{-\beta_i\gamma}}{G_i(\gamma)}\sqrt{\gamma^2-1}\gamma\text{d}\gamma.\end{equation} \end{align*} \tag{ 4 } \end{equation*}$$

Apparently ξi=Ni$\xi_i =N_i$ holds for this classical distribution, and Zmix$Z_{\text{mix}}$ is the product of the individual components Z(Zi,mi)$Z_{(Z_i ,m_i)}$. Zmix$Z_{\text{mix}}$ depends on the temperature variables βi$\beta_i$ and the fugacity parameters αi$\alpha_i$; the latter can be parametrized by particle number and temperature, αi(βi,Ni/V)$\alpha_i (\beta_i ,N_i /V)$, by solving (4). The entropy of a nonthermal gas mixture is an extensive quantity,

Smix=∑i=1MSi,Si=kB(logZ(Zi,mi)+βimiUi+αiNi),(5)

$$\begin{equation*} \begin{align*} \begin{equation} S_{\text{mix}} =\sum\limits_{i=1}^M {S_i}, \quad S_i =k_\text{B} \bigg(\log Z_{(Z_i ,m_i)} +\frac{\beta_i}{m_i}U_i +\alpha_iN_i\bigg),\end{equation} \end{align*} \tag{ 5 } \end{equation*}$$

obtained from an ensemble average of power-law partitions, see the next section.

We give a probabilistic derivation of the classical spectral number density dρmix(γ)$\text{d}\rho_{\text{mix}}(\gamma)$ of an ultra-relativistic gas mixture in stationary non-equilibrium, cf. (1). To this end, we consider the phase space of n particles of a gas component (Zi,mi)$(Z_i ,m_i )$ labeled by charge number and mass in an M-component mixture, i=1,…,M$i=1,\ldots ,M$. The phase-space measure of component (Zi,mi)$(Z_i ,m_i )$ is

d3n(Zi,mi)(p,q)=1(2πℏ)3nd3p1d3q1d3p2d3q2⋯d3pnd3qn,d3pk=4πm3ic3γ2k−1−−−−−√γkdγk,k=1,…,n.(6)

$$\begin{equation*} \begin{align*} \begin{eqnarray} &\text{d}_{(Z_i ,m_i )}^{3n} (p,q)=\frac{1}{(2\pi \hbar )^{3n}}\text{d}^3p_1 \text{d}^3q_1 \text{d}^3p_2 \text{d}^3q_2 \cdots \text{d}^3p_n \text{d}^3q_n,\nonumber\\ &\text{d}^3p_k =\frac{4\pi m_i^3}{c^3}\sqrt{\gamma_k^2 -1} \gamma_k\text{d}\gamma_k, \qquad k=1,\ldots ,n.\end{eqnarray} \end{align*} \tag{ 6 } \end{equation*}$$

The d3qk$\text{d}^3q_k$ integrations range over a volume V, and the dγk$\text{d}\gamma_k$ integrations refer to the interval γcuti≤γk≤∞$\gamma_i^{\text{cut}} \le \gamma_k \le \infty$, cf. after (3), where γcuti$\gamma_i^{\text{cut}}$ is the cutoff Lorentz factor of the respective gas component (Zi,mi)$(Z_i ,m_i )$. The angular integration over the solid angle dΩ$\text{d}\Omega$ has been carried out in (6) and gives the factor 4π$4\pi$. The scale factor (2πℏ)3n$(2\pi \hbar)^{3n}$ renders the measure dimensionless.

Probability density and grand partition function of mixtures with power-law components

The probability density on the n-particle states of a gas component (Zi,mi)$(Z_i,m_i)$ factorizes as

f(Zi,mi)n(p,q)b(Zi,mi)nh(Zi,mi)n=snin!b(Zi,mi)n(p,q)h(Zi,mi)n(p,q),=exp(−ξi−βi∑k=1nγk−αin),=exp(−∑k=1nlogGi(γk)).(7)

$$\begin{equation*} \begin{align*} \begin{eqnarray} f_n^{(Z_i ,m_i )}(p,q)&=\frac{s_i^n}{n!}b_n^{(Z_i ,m_i )}(p,q)h_n^{(Z_i,m_i )}(p,q),\nonumber\\ b_n^{(Z_i ,m_i )} &=\exp \bigg(-\xi_i -\beta_i\sum\limits_{k=1}^n {\gamma_k} -\alpha_i n\bigg),\nonumber\\ h_n^{(Z_i ,m_i )} &=\exp \bigg(-\sum\limits_{k=1}^n\log G_i (\gamma_k)\bigg).\end{eqnarray} \end{align*} \tag{ 7 } \end{equation*}$$

The factor 1/n!$1/n!$ accounts for the indistinguishability within the particle species labeled (Zi,mi)$(Z_i ,m_i )$, and the positive integer s_i for the spin multiplicity of this component. βi=mi/(kBTi)$\beta_i=m_i/(k_\text{B} T_i )$ is the temperature parameter of the respective gas component, and ξi$\xi_i$ is a normalization constant.

The normalization condition for the phase-space probability of an M-component mixture is

∫∑n1,…,nM=0∞∏i=1Mf(Zi,mi)ni(p,q)d3ni(Zi,mi)(p,q)=∏i=1M∫∑n=0∞f(Zi,mi)n(p,q)d3n(Zi,mi)(p,q)=1,(8)

$$\begin{equation*} \begin{align*} \begin{eqnarray} &\int\sum\limits_{n_1 ,\ldots ,n_M =0}^\infty\prod\limits_{i=1}^M f_{n_i}^{(Z_i ,m_i )} (p,q)\text{d}_{(Z_i ,m_i )}^{3n_i}(p,q)=\nonumber\\ &\prod\limits_{i=1}^M\int\sum\limits_{n=0}^\infty f_n^{(Z_i ,m_i )}(p,q)\text{d}_{(Z_i ,m_i)}^{3n}(p,q)=1,\end{eqnarray} \end{align*} \tag{ 8 } \end{equation*}$$

where summation and integration sign can be interchanged, as well as the integration and product signs. The (p,q)$(p,q)$ variables of the density and measure will occasionally be omitted. The integrals in (8) factorize into one-particle states, and normalization is achieved by identifying ξi=logZ(Zi,mi)$\xi_i =\log Z_{(Z_i ,m_i )}$ in (7) with the partition function (4). Thus,

Z(Zi,mi)=∑n=0∞snin!∫exp(−βi∑k=1nγk−αin)h(Zi,mi)nd3n(Zi,mi),(9)

$$\begin{equation*} \begin{align*} \begin{eqnarray} &Z_{(Z_i ,m_i )}= &\sum\limits_{n=0}^\infty \frac{s_i^n}{n!}\int\exp \left(-\beta_i \sum\limits_{k=1}^n{\gamma_k } -\alpha_i n\right)h_n^{(Z_i,m_i)} \text{d}_{(Z_i ,m_i)}^{3n},\qquad\end{eqnarray} \end{align*} \tag{ 9 } \end{equation*}$$

for each gas component (Zi,mi)$(Z_i ,m_i )$. The partition function of an ideal mixture factorizes, Zmix=∏i=1MZ(Zi,mi)$Z_{\text{mix}} =\prod\limits_{i=1}^M {Z_{(Z_i ,m_i)}}$, or additively ξmix=logZmix$\xi_{\text{mix}} =\log Z_{\text{mix}}$.

↑ Close

Entropy of a nonthermal gas mixture

The expectation values of particle number and internal energy of a gas component (Zi,mi)$(Z_i ,m_i)$ read

∑n=0∞∫nf(Zi,mi)n(p,q)d3n(Zi,mi)(p,q)=−dξidαi=⟨n⟩(Zi,mi)=Ni,(10)

$$\begin{equation*} \begin{align*} \begin{eqnarray} &\sum\limits_{n=0}^\infty\int nf_n^{(Z_i ,m_i )} (p,q)\text{d}_{(Z_i ,m_i)}^{3n}(p,q)=\nonumber\\ &-\frac{\text{d}\xi_i}{\text{d}\alpha_i} =\langle n\rangle_{(Z_i ,m_i)} =N_i,\end{eqnarray} \end{align*} \tag{ 10 } \end{equation*}$$

∑n=0∞∫(mi∑k=1nγk)f(Zi,mi)n(p,q)d3n(Zi,mi)(p,q)=−midξidβi=⟨u⟩(Zi,mi)=Ui,(11)

$$\begin{equation*} \begin{align*} \begin{eqnarray} &\sum\limits_{n=0}^\infty \int \bigg(m_i \sum\limits_{k=1}^n {\gamma_k }\bigg)f_n^{(Z_i ,m_i)}(p,q)\text{d}_{(Z_i ,m_i)}^{3n}(p,q)=\nonumber\\ &-m_i \frac{\text{d}\xi_i}{\text{d}\beta_i} =\langle u\rangle_{(Z_i ,m_i )} =U_i,\end{eqnarray} \end{align*} \tag{ 11 } \end{equation*}$$

where u=mi∑k=1nγk$u=m_i \sum\limits_{k=1}^n{\gamma_k}$, and Ni=ξi$N_i =\xi_i$, cf. (4). The spectral number density (1) follows from the phase-space average (10). The total particle count Nmix$N_{\text{mix}}$ and internal energy Umix$U_{\text{mix}}$ are obtained by summing over the particle species, cf. (3). The total entropy Smix=∑i=1MSi$S_{\text{mix}} =\sum\limits_{i=1}^M {S_i}$ of the mixture is found by adding the partial entropies defined by the average

Si=−kB∫∑n=0∞snin!h(Zi,mi)nb(Zi,mi)nlogb(Zi,mi)nd3n(Zi,mi),(12)

$$\begin{equation*} \begin{align*} \begin{equation} S_i =-k_\text{B} \int\sum\limits_{n=0}^\infty\frac{s_i^n}{n!}h_n^{(Z_i,m_i )} b_n^{(Z_i ,m_i )} \log b_n^{(Z_i ,m_i )} \text{d}_{(Z_i ,m_i )}^{3n},\end{equation} \end{align*} \tag{ 12 } \end{equation*}$$

which is the expectation value ⟨s⟩(Zi,mi)$\langle s\rangle_{(Z_i ,m_i )}$ of the phase-space random variable s=kB(ξi+βi∑k=1nγk+αin)$s=k_\text{B} (\xi_i +\beta_i \sum\limits_{k=1}^n {\gamma_k } +\alpha_i n)$,

Si=kB∫∑n=0∞(ξi+βi∑k=1nγk+αin)f(Zi,mi)nd3n(Zi,mi).(13)

$$\begin{equation*} \begin{align*} \begin{equation} S_i =k_\text{B} \int\sum\limits_{n=0}^\infty\bigg(\xi_i +\beta_i \sum\limits_{k=1}^n{\gamma_k} +\alpha_i n\bigg)f_n^{(Z_i ,m_i )} \text{d}_{(Z_i ,m_i )}^{3n}.\end{equation} \end{align*} \tag{ 13 } \end{equation*}$$

The partial entropies S_i stated in (5) are recovered from (13) by substitution of the averages (10) and (11).

↑ Close

We write the components of the total spectral number density dρmix$\text{d}\rho_{\text{mix}}$ in (1) as

dρi(E)dE=4πm2isie−αi(2πℏc)3gi(1)gi(γ)γ1−δie−βiγγ2−1−−−−−√,(14)

$$\begin{equation*} \begin{align*} \begin{equation} \frac{\text{d}\rho_i (E)}{\text{d}E}=4\pi \frac{m_i^2 s_i \text{e}^{-\alpha_i}}{(2\pi\hbar c)^3}\frac{g_i(1)}{g_i(\gamma)}\gamma^{1-\delta_i}\text{e}^{-\beta_i\gamma}\sqrt{\gamma^2-1},\end{equation} \end{align*} \tag{ 14 } \end{equation*}$$

with spectral function gi(γ)$g_i (\gamma)$ in (2), and replace the Lorentz factor by the particle energy γ=E/mi$\gamma =E/m_i$. The mass m_i stands for m_i c², s_i is the spin degeneracy of component (Zi,mi)$(Z_i ,m_i )$, cf. after (1), and zi=e−αi$z_i =\text{e}^{-\alpha_i}$ the fugacity. βi=mi/(kBTi)$\beta_i =m_i/(k_\text{B}T_i)$, and δi$\delta_i$ is a real power-law exponent. The amplitude of the second power law in gi(γ)=1+aiγσi$g_i (\gamma )=1+a_i \gamma^{\sigma_i}$ is positive, ai≥0$a_i \ge 0$, and we will use positive exponents σi$\sigma_i$. Hence,

dρi(E)dE=4πsie−αˆi(2πℏc)3E2−δie−βˆiE1+aˆiEσi1−m2iE2−−−−−−−√,(15)

$$\begin{equation*} \begin{align*} \begin{equation} \frac{\text{d}\rho_i (E)}{\text{d}E}=\frac{4\pi s_i \text{e}^{-\hat{\alpha}_i}}{(2\pi\hbar c)^3} \frac{E^{2-\delta_i}\text{e}^{-\hat{\beta}_iE}}{1+\hat{a}_i E^{\sigma_i}}\sqrt{1-\frac{m_i^2}{E^2}},\end{equation} \end{align*} \tag{ 15 } \end{equation*}$$

where we introduced the shortcuts βˆi=1/(kBTi)$\hat{\beta}_i =1/(k_\text{B}T_i)$, aˆi=ai/mσii$\hat{a}_i =a_i /m_i^{\sigma_i}$ and

αˆi=αi−log(1+ai)−δilogmi.(16)

$$\begin{equation*} \begin{align*} \begin{equation} \hat{\alpha}_i =\alpha_i -\log (1+a_i)-\delta_i \log m_i.\end{equation} \end{align*} \tag{ 16 } \end{equation*}$$

We also note e−αˆi=zimδii(1+ai)$\text{e}^{-\hat{\alpha}_i}=z_im_i^{\delta_i}(1+a_i)$. We will use dimensionless quantities: E and m_i in eV units, and the units chosen for length, time and temperature are, respectively, meter, second and kelvin.

The spectral flux density of a gas component is found by multiplying the number density (15) with the particle velocity υi=c1−m2i/E2−−−−−−−−−√$\upsilon_i =c\sqrt{1-m_i^2 /E^2}$,

Fi(E)=c4π1−m2iE2−−−−−−−√dρidE=csie−αˆi(2πℏc)3E2−δie−βˆiE1+aˆiEσi(1−m2iE2).(17)

$$\begin{equation*} \begin{align*} \begin{eqnarray} F_i (E)&=\frac{c}{4\pi}\sqrt{1-\frac{m_i^2}{E^2}} \frac{\text{d}\rho_i}{\text{d}E}=\nonumber\\ &\quad\frac{cs_i \text{e}^{-\hat{\alpha}_i}}{(2\pi\hbar c)^3}\frac{E^{2-\delta_i}\text{e}^{-\hat{\beta}_iE}}{1+\hat{a}_iE^{\sigma_i}}\bigg(1-\frac{m_i^2}{E^2}\bigg).\end{eqnarray} \end{align*} \tag{ 17 } \end{equation*}$$

Density Fi(E)$F_i (E)$ is the number flux density per steradian; the energy flux is obtained by rescaling Fi(E)$F_i (E)$ with E. Wideband spectra are assembled by adding the partial flux densities of the gas components, Fmix(E)=∑i=1MFi(E)$F_{\text{mix}} (E)=\sum\limits_{i=1}^M {F_i(E)}$. To locate the "knees" and "ankles" of the spectrum in log-log plots, it is convenient to rescale Fmix(E)$F_{\text{mix}}(E)$ with a suitable power E^k. The rescaled density EkFmix(E)$E^kF_{\text{mix}} (E)$ is thus obtained by adding the partial fluxes