What can cosmic-ray knees reveal about source populations?
Myrto Falalaki Vasiliki Pavlidou
 
	
	
	
	
	
	
	
		University of Crete, Department of Physics & Institute of
		Theoretical & Computational Physics, 70013 Herakleio, Greece 

		ph5690@edu.physics.uoc.gr; pavlidou@physics.uoc.gr
		Institute of Astrophysics,
		Foundation for Research and Technology-Hellas, 71110 Heraklion, Crete, Greece
		
	
	
	
	
	
	
	
	Breaks in the cosmic ray (CR) flux spectrum encode information on the properties of CR accelerator populations producing the observed flux. Spectral steepenings, known as knees, are  generally accompanied by a transition to heavier composition.
We seek generic features of CR source populations that imprint onto  knee observables in a manner that is robust enough and distinct enough to be discernible even in the presence of significant uncertainties and systematics in CR data. In particular, we explore how diversity among population members could imprint on the knee phenomenology, under the assumption that a knee is due to a fixed-rigidity cutoff in the source spectrum.

We use a simple theoretical model for a population of CR accelerators. Each member of the population accelerates CR stochastically to a power-law spectrum, up to a cutoff rigidity, resulting from source-confinement requirements. We allow for variance among members, in the cutoff rigidity and in the power-law slope.
We find that: (a) the slope step of the flux spectrum is ∼ 0.5, decreasing weakly with increasing spread in either property; (b) composition always breaks first; (c) the difference between the break energies in composition and flux increases with increasing diversity; (d) composition and flux break together only if population diversity is minimal. 
These trends are robust under our assumptions; deviations from them in observed data would indicate more complex physics than encoded in our simple model. 
Comparing these trends with observed CR knees, we conclude that: (i) the primary knee at ∼ 4×10^15 eV is consistent with a constant-rigidity cutoff according to KASCADE-Grande data processed with post-LHC hadronic models, but not according to other datasets; (ii) the second knee at ∼ 5 × 10^17 eV conclusively requires more complexity than the cutoff of a single CR source population; (iii) the spectral feature identified by Auger at ∼ 10^19 eV is consistent with a constant-rigidity source cutoff only if there is a substantial spread in both cutoff rigidity and slope in the parent source population. Interestingly, a significant spread in slope would also result in spectral curvature before the break, which would in turn be contributing to the ankle feature.
	
	

	Falalaki & Pavlidou
	
	
	What can cosmic-ray knees reveal about source populations? 
    Myrto Falalaki
		1,2
		Vasiliki Pavlidou1,2
		
    Received / accepted 
===========================================================


	
	

§ INTRODUCTION
 

    The energy spectrum of cosmic radiation is a critical observable for our understanding of the nature of cosmic ray (CR) accelerators.  This spectrum can  be described as a power law over many orders of magnitude (from ∼10^9 eV to ∼10^20 eV). The spectral power index, however, exhibits changes at characteristic energies, referred to as knees (when the spectrum steepens) and ankles (when the spectrum flattens). 
    
  
 Knees  have now been spectrally resolved in great detail at several different energies in the CR spectrum, ranging from few times 10^15 eV to 10^19 eV (e.g., ). These breaks are typically accompanied by transitions to heavier composition around the same energies as the spectral break (e.g., ). This behavior has led to the qualitative interpretation of the knee phenomenon as a constant-rigidity cutoff (also referred to as Peters cycle, ), either in Galactic confinement or in source acceleration / confinement: particles accelerated stochastically by a CR source cannot reach energies beyond the threshold where magnetic confinement is lost . This effect occurs roughly at the energy where a particle's gyroradius in the source's magnetic field becomes comparable to the size of the source (e.g., ). 
 
A rich literature exists on quantitative models of the various knee-like breaks of the CR spectrum, ranging from fits of specific datasets with complex phenomenological models, to comprehensive physical models, including specific hypotheses on source accelerator physics and source population properties, propagation effects, and simulations of atmospheric air showers to obtain direct observables (e.g., ). 

Despite these intensive theoretical and experimental efforts, the exact characteristics, nature and origin of all observed CR knees is still under debate. One very important factor contributing to this continued uncertainty is that the data are not yet fully converged. High-level observables (i.e. spectra and moments of the log mass-number distribution) can differ between observatories using different observational techniques at the same energies (see e.g. review of data on the primary knee in ); between different observatories using similar techniques (see, e.g., differences in flux and composition at the highest energies reported by the Pierre Auger Observatory and Telescope Array, ); between events recorded by the same observatory when using different observables (see, e.g., differences in composition at ultrahigh energies derived when using shower depth or shower muon content, e.g., ); and even between identical recorded datasets when processed with different simulations of hadronic interactions (for example, with pre- versus post-LHC models, ; or with different hadronic interaction packages of the same generation, ). As a result, detailed fits to sophisticated source population models may fail or return parameters that appear astrophysically contrivied (e.g., ), not because the models are necessarily inconsistent with the actual source populations, but rather because our data from said source populations might be affected by systematics that are not adequately quantified or accounted for. 
 
The question then arises whether there exist any generic features of cosmic-ray source populations that imprint onto cosmic-ray observables in a manner that is robust enough and distinct enough to be discernible even in such a still-fluid experimental landscape. This is our aim in this paper. Using the simplest possible model for an underlying population of CR accelerators, we seek to build insight on how  the diversity between population members could imprint on the knee phenomenology, under the assumption that a knee is a fixed-rigidity-cutoff phenomenon. Such insight is necessary in order to identify any robust features and trends of simple population models and their most straightforward variations.  Absence of said features would then constitute evidence of more complex physics of either CR acceleration (complicated astrophysics) or CR interactions at the detection site (unexpected particle physics). 

 A knee is described phenomenologically by the characteristic energy where the break occurs in the CR flux spectrum, by the spectral slopes before and after this characteristic energy, and the (generally different) energy where the accompanying break in CR composition takes place.  We will investigate whether features of the population of CR sources responsible for the knee result in patterns in these observables that are simple and robust enough that they might be recognizable even in the presence of significant systematic uncertainties in the data. 
 
 To this end, we set up a simple, generic model for a population of CR accelerators with the following properties: (a) For energies well below its rigidity cutoff (equal to the proton cutoff energy,  ), each source contributes to the Galactic CR flux particles with a power-law spectrum of energies of slope γ (encoding both the source acceleration properties, and losses during propagation). (b) The rigidity cutoff can be described by an exponential suppression of the power-law spectrum. Different nuclei cutoff at different energies scaling as Z. (c) We allow diversity between members of the CR source population in both  and γ. 
A "knee" observed in the summed CR spectrum due to such a population of sources will consist of: a break in the spectrum, encoding, but not necessarily equal to, ; and a break in the composition at a similar, but not necessarily identical, energy. The location in energy of those breaks, and the slope difference between and after the knee will be modulated by the distribution of γ and  among the individual members of the population. 

We explore the features and trends of this simple model, and in particular we address the following specific questions: (1) Under what conditions do the flux spectrum and composition spectrum break together? (2) When they do not, which one breaks first? (3) How do the break energies relate to ? (4) How does the diversity of the population properties, manifesting as a spread in  and γ, affect the answers to these questions?

This paper is organized as follows. In <ref> we lay out the formulation of our model, emphasizing its free parameters and their impact on the observable flux and composition spectra. In <ref>
we explore the behavior of these observables as the  population diversity increases. We summarize our conclusions and discuss them in the context of observations of various CR knees in <ref>. 
 



§ THE MODEL
 


 §.§ Single Source

	We implement a fixed-rigidity cutoff, characteristic of CR sources relying on magnetic confinement
 
 <cit.>.





 We do so by assuming that  the differential flux F_z produced by a single cosmic ray source as a function of the energy E of primary particles of charge Z has the form:
    
    F_z(E) = F_0(Z)(EE_0)^-γexp[ -EZ ]

    where   is the rigidity cutoff, F_0(Z) is the flux of a specific species of  atomic number Z at some normalization energy E_0 ≪,  and γ is the low-energy power-law slope. In our model, we assume that γ is identical for all species in a single source, and that it encodes both acceleration and propagation physics (including losses and escape). 


The total flux from a single source at an energy E will be the sum of F_z(E) over all primary charges:
    
    F(E) =∑_ZF_z(E) .

In this work, we have taken the relative abundances F_0(Z) from <cit.> (see their table 1, where they present their fitted parameters for a power-law like model for every nuclei). While <cit.> fit a different slope γ_z for different species, we have adopted here an effective value of γ≈ 2.66 for all species. However, our results are independent of the exact choice for the value of γ, since we always plot deviations (differences) of fitted slopes from the underlying source γ.  We simulate fluxes for the following individual species: H, He, C, O, Ne, Mg, Si and Fe. These abundances are observationally motivated for energies below the primary CR knee; however our qualitative conclusions do not depend sensitively on this choice, and so they hold for any source population accelerating particles of mixed composition roughly comparable to that of Galactic CR, even if the pre-break abundances differ in their details from the ones we have adopted here. 

We quantify the resulting composition spectrum by the average of the logarithm of the mass number as a function of energy, ⟨ln A⟩ (E), since this metric is frequently used to summarize CR composition observations. For a single source, this will be given by
  
    ⟨lnA⟩ (E) = ∑_Z F_Z(E) ln (A_Z)/∑_Z F_Z(E) .





 §.§ Source population
 

The CR observables on Earth at a given energy are produced by particles accelerated by a population of sources with a distribution of sizes, B-fields and other properties. Variations in these  physical conditions among population members will result in corresponding variations of source model parameters <cit.>. In our simple source model, these parameters are F_0(Z), γ, and . In this work, we are interested in any signatures of 
the diversity in  and γ imprinted on the observables of a  knee produced by a constant-rigidity cutoff in the spectra of sources. To isolate the effects of each of these source properties, we explore the effect of increasing spread in one of , γ, while keeping the distribution of the other fixed to a delta function. We also make the simplifying assumption that  relative values of f̃_0(Z)=F_0(Z)/F_0(Z=1) (the accelerated particles relative abundances at energies well below the proton cutoff) are identical among different population members, so that the only quantities that may vary between individual CR sources are F_0(Z=1),  , and γ.  

Mathematically, we model this picture as follows. Choosing a normalization energy E_0 well below the lowest  encountered in the specific population, the exponential suppression factor in Eq. (<ref>) at E_0 is equal to 1 for all sources. The resulting total flux from all population sources at that energy 
can be calculated as:

    F_Z,tot(E_0≪min) =∑_i F_0,i(Z=1) 
       f̃_0(Z) ≡ F_Z,tot,0 ,

where the summation is over different population member sources. This equation defines the flux normalization of species Z for the population, F_Z,tot,0. 

To calculate the spectrum at higher energies, we introduce
the probability distributions  p_E() and p_γ(γ) of the cutoff rigidity, and the low-energy power-law slope, respectively, in the population. Formally, p_E()d is the fraction of particles at energy E_0 that were accelerated by sources with rigidity cutoffs between  and +d; and p_γ(γ)dγ is the fraction of particles accelerated by sources that, had they been responsible for the entire CR spectrum at low energies, would have produced  (after acceleration and propagation) a  spectrum which would have been a superposition of power laws with slopes between γ and γ + dγ. 
Then, the total flux of species Z due to the entire population can be calculated through 
    
    F_Z,  tot(E)
            = F_Z, tot,0∫_0^∞ d p_E()
             e^-EZ∫_0^∞ dγ p_γ(γ)
     (E/E_0)^-γ.

The total CR flux due to all species  will be 

    F_ tot(E)
        = ∑_Z     F_z,  tot(E) .

Finally, the composition spectrum due to the population will be given by 
  
    ⟨lnA⟩_ tot (E) = ∑_Z F_Z,  tot(E) ln (A_Z)/∑_Z F_Z,  tot(E) .

    
 

 §.§ Varying 

 We will estimate the effect of a spread in  among population members by assuming p_E() is lognormal: 
    
    p_E() = 1σ√(2π)exp[ -(ln-μ)^22σ^2] .

The parameters of the distribution μ and σ are related to the mean and the standard deviation of the distribution through 
   
    μ = ln( ^2/√(^2+^2))

    and 
    
    σ^2 = ln[1+()^2] .

The advantage of the lognormal distribution is that  is positive definite and as a result arbitrarily large values of  can be accomodated, allowing us to examine asymptotic behaviors for very diverse populations. The distribution p_γ remains fixed to a delta function. 




 §.§ Varying γ

To estimate the effect of a spread in γ, we keep p_E() fixed to a delta function, 
and we implement 
a Gaussian  p_γ(γ), with mean γ̅ and spread σ_γ. In this case, Eq.(<ref>) has an analytic solution,
     
    F_Z,  tot(E)
            =  
    F_z,tot,0exp[ -E/Z]
     (E/E_0)^-γ̅exp[ -σ^2_γ/2ln^2( EE_0) ].




 §.§ Observables


Even though the flux suppression of individual CR species is exponential, the fact that the suppression sets in at increasing energy for increasing Z  results to a flattening of the knee. This will be true for individual sources (Eq. <ref>), and more so for a population (Eq. <ref>) where the spread in  and/or γ will result in a smoother break. Observed knees can be fitted well by broken power laws, and we attempt the same for our model results. In particular, we test whether 
a broken power law of the form:
    
    F_tot(E) = C ×(EE_b)^-γ_1 ,     E ≤ E_b 
    (EE_b)^-γ_2 ,     E > E_b

    can describe adequately the total flux around the knee. From this fit, we exctract "observables" γ_1 (slope before the break), γ_2 (slope after the break), and E_b (break energy). 

    We additionally evaluate a fourth "observable", the composition break energy, E_A. We do so by:
    (a) fitting the composition spectrum of Eq. (<ref>) for E≪ by a constant (horizontal line); and
    (b) fitting the composition spectrum  by a logarithmic increase (linear increase with log E) for  E≫. Then, E_A is the energy at which the two lines intersect (see lower-right panel of Fig. <ref>). 
    
	

§ RESULTS
 

We start by examining the trends induced on the "observables" by a gradually increasing spread in  among population members. 


  
 
  In the upper left panel of Fig. <ref> we show the all-particle spectrum (flattened by E^γ) as we evaluate it from Eqs. (<ref>) and (<ref>) and for a lognormal p(), for ∈[0,3]. For broader p(), the all-particle spectrum deviates sooner (at lower energies) from its low-energy asymptotic behavior, as a result of the property of the lognormal distribution to peak around values of <. At the same time, the spectrum falls less steeply at high energies, owing to the tail of the longnormal towards high  values.  


In the lower-left panel of Fig. <ref>, we show an example of fitting the all-particle flux spectrum around the knee with a broken power law.  The case depicted here is produced by setting = in the lognormal p().  It is through such fits that we obtain the break point E_b and the spectral power indices γ_1,γ_2 (the "observables" discussed in the context of Eq. <ref>), which we present and discuss in Figs. 
<ref> and <ref> as functions of .

The upper-right panel of Fig. <ref> shows the composition spectrum for the same models as in the upper-left panel. The break in flux is accompanied by a break in composition. The impact of increasing  on the composition spectrum is more pronounced than that on the flux spectrum, both towards lower and higher energies, as well as in terms of the steepness of the break. To quantify this behavior, we again fit the low-energy and high-energy trends with logarithmic functions (⟨ln A ⟩ linear in log E), as shown in the lower-right panel of Fig. <ref>. The intersection of the two fits defines the composition break energy, E_A.


    

The effect of  on the "observables" is shown in Figs. <ref> and <ref>. 
Figure <ref> shows the exponents γ_1 and γ_2 (slopes before and after the break, red and blue points respectively), evaluated as difference from the single-source slope γ.  The error budget is completely dominated by systematic uncertainties, primarily driven by the choice of points to include in the fit. Here, error bars  correspond to the difference between including or dropping an extra point at the high-energy tail of the fit (see lower-left panel of Fig. <ref>).
 
    
    
The difference between the two slopes starts at ∼ 0.6 for a population comprised of members with identical  (=0), and decreases slowly as  increases, reaching ∼ 0.5 for = 3. This mild trend is a result of two compounding effects. First, as  increases, the number of sources with low  also increases, so the low-energy branch of the broken power law is also affected (becomes steeper, γ_1 increses) as some sources have already started becoming suppressed at low energies. Second, a high  also results in a larger number of sources with high . The flux of these sources is not suppressed until higher energies, resulting in a smoother decline of the all-particle spectrum. The high-energy branch of the power law thus becomes shallower (γ_2 decreases). The overall conclusion is that a fixed-rigidity cutoff in a CR source population with identical power-law slopes γ and a pre-break composition roughly comparable to that of Galactic CR produces a knee with a slope change in the range of ∼ 0.5-0.6, regardless of spread in . 

In Fig. <ref> we overplot the break energies of the spectrum (E_b, blue points) and of the composition (E_A, red points) as a function of the spread in . Error bars are again dominated by systematics. For E_b, they are obtained in the same way as the error bars in γ_1 and γ_2. In the case of E_A, they correspond to the difference resulting from retaining or dropping one point towards the break in the rising part of the composition spectrum. We observe that in the case of a population with members of identical, or very similar, , composition and flux spectrum break together, at an energy about half of . As  increases both the composition and the flux spectrum break earlier, however E_A decreases faster than E_b, and as a result composition always breaks first. 

 


        


        


We now turn to trends in "observables" that result from a spread in γ among population members. For these calculations, we keep p_E() fixed to a delta function. In Fig.  <ref> we plot the all-particle spectrum (flattened by E^γ̅) as we evaluate it from Eqs. (<ref>) and (<ref>) and for a Gaussian  p_γ(γ) .  Different line colors correspond to different values of σ_γ. Here, the range of σ_γ we consider is much narrower than the range of . The reason is that a large spread in power law indices results in significant spectral curvature (see e.g. black dashed line in Fig. <ref>), which is not generally seen in the CR spectrum. Interestingly, small spreads in γ appear to generate  diversity in the flux spectrum comparable to that produced by very substantial spreads in . 

We investigate the qualitative direction of the  trends seen in the "observables" with increasing σ_γ, and quantify them, in Figs. <ref> and <ref>. The trend of the sharpness of the break with σ_γ is shown in  Fig. <ref>. Colors and error bars are as in Fig. <ref>. Here again the slope step decreases with increasing population spread, from about 0.6 for a population with no spread, to about 0.45 for σ_γ=0.2, driven primarily by the post-break slope becoming shallower. Still, the effect is very mild: in a standard fixed-rigidity-cutoff knee, the slope change does not become very different from a 0.5 step, even if the underlying accelerator population exhibits significant spread in its properties. 

The trends in energy breaks however are now different, as is shown in Fig. <ref>. 
The location of the composition break, E_A, is unaffected by any spread in γ. 
The flux spectrum break, E_b, on the other hand increases with increasing population spread: the more diverse in γ the population, the later the flux spectrum breaks. Importantly, the overall result goes in the same direction as in the case of a spread in : composition always breaks first. 
 
    
    

§ CONCLUSIONS AND DISCUSSION 
 


 We have used the simplest possible model of a knee induced by a fixed-rigidity CR-source cutoff, in order to explore how diversity among the members of the underlying cosmic-ray accelerator population affects the knee phenomenology. In particular, we explored how diversity in (a) the rigidity cutoff  and (b) in the pre-knee cosmic-ray slope γ (encoding both  acceleration and loss/propagation physics) is imprinted in (i) the difference in pre- and post-break flux spectrum slopes, and (ii) in the break energies of the flux and composition spectra. We have identified the following robust trends.

 1. A knee induced by a fixed-rigidity cutoff in the CR source population exhibits a steep break in the all-particle spectrum, with a slope step around 0.5. Diversity in the population (in either  or γ) tends to somewhat reduce the steepness of the break, but the effect is weak.   

 2. Diversity in  moves both the flux break energy, E_b, and the composition break energy, E_A, to lower values, with the composition being more strongly affected. In contrast, diversity in γ leaves E_A practically unaffected, but moves E_b to higher energies.

 3. Any diversity in either  or γ works to separate E_b from E_A, in the same direction: composition breaks first. The difference between the two can be as large as a factor of several, especially if both  and γ vary substantially among population members.  

 4. An observation of the composition and flux breaking around the same energy would constitute an indication of minimal spread in both  and γ in the underlying CR source population.

 Any qualitative deviation from these trends would require a model of significantly higher astrophysical complexity than the one we discussed here. Additions to the physics generating the observed flux and composition spectra that could break these trends include, for example, the presence of multiple, distinct populations of cosmic ray accelerators contributing at different flux levels, or intergalactic propagation effects if relevant at the energies under consideration. 

 We can use these insights to discuss qualitatively the likelihood of each of the steepening features in the broadband cosmic-ray spectrum being a simple fixed-rigidity CR-source-cutoff knee, even without any detailed fits to a specific model or the underlying accelerator populations. 

For the primary CR "knee" at ∼ 4× 10^15 eV <cit.>, the situation is still unclear, as observations from different experiments and using different techniques have unfortunately not yet fully converged, even at the very coarse level needed for the type of comparison we are seeking to make in this work. For example, early results presented in the <cit.> review, but also very recent results from the Large High Altitude Air Shower Observatory (LHAASO) experiment <cit.>, indicate a small slope step (between 0.2 and 0.4), and the composition breaking after the spectrum. In the context of a 
simple fixed-rigidity source-cutoff knee described here, a slope step on the low side could be conceivably achieved with a combination of variations in both  and γ in the population of contributing sources. However, such a physical picture cannot accommodate the composition breaking after the spectrum.  Even if we were to accept that, within uncertainties, composition and spectrum could be breaking together, this would point towards a population with very little spread in both γ and , which would then result to a sharper slope step than the one observed. This combination of phenomenological observables then hint towards for a more complex picture - plausibly, for example, towards more than one source populations contributing cosmic rays around the primary knee, or additional physical effects (e.g., ).  On the other hand, in the latest, post-LHC, reanalysis of data from the KASCADE-Grande experiment <cit.>, the composition appears to be already getting heavier before the break in spectrum, while the slope step is reported at ∼ 0.5, as might expected from the simplest version of fixed-rigidity knee explored here.  

For the "second knee" around 5×10^17 eV (see e.g.  for a review), the observational situation is similarly unconverged.  Different datasets disagree over the location of the composition break, and on whether, before the break, the composition was getting heavier <cit.> or lighter <cit.>.  However there is consensus that the break is soft (slope step between 0.2 and 0.3), while the composition across the second knee is becoming lighter. Most likely then in this case there is a second, light (i.e., still efficiently accelerating) population contributing (e.g., ), so the simple physics we explored is not adequate to model this transition. 


The knee-like "new feature" identified by the Pierre Auger Observatory in the cosmic ray spectrum around 10^19 eV <cit.> features a slope step of ∼ 0.5, right in the expected range for a fixed-rigidity source-cutoff knee. This lends support to the interpretation that the break seen in composition-sensitive observables before the break in the spectrum <cit.> indeed indicates a transition to heavier primaries. This is further reinforced by the reconstructed shape of the ⟨ln A ⟩ spectrum. Although the overall normalization of the composition spectrum is very dependent on the choice for hadronic interaction modelling, the relative change of ⟨ln A ⟩ with energy does feature the overall shape seen in the upper-right panel of Fig. <ref>: ⟨ln A ⟩ flattens off between 10^18 and 10^18.5 eV, before starting to increase logarithmically around 10^18.6 eV. Our interpretation is then in overall agreement with the one proposed by <cit.>. Interestingly, however, the difference between the locations of the composition break (around 10^18.6 eV if we follow the same procedure we have used here to determine E_A) and the spectrum break (reported by  at 10^19.1 eV) is higher that the differences we have seen produced by variations in  alone. The implication is that γ may also be significantly varying in the dominant CR source population at these energies. This in turn could produce a non-negligible curvature in the spectrum before the break (see black dashed line in Fig. <ref>), that may be contributing to the curvature of the ankle, which is also located at the same energy as the composition break. 

 


	
	
 
MF acknowledges support 
by the European Research Council under the European Union's Horizon 2020 research and innovation programme, grant agreement No. 771282 (PASIPHAE). VP acknowledges support by the Hellenic Foundation for Research and Innovation under the “First Call for H.F.R.I. Research Projects to support Faculty members and Researchers and the procurement of high-cost research equipment grant”,  Project 1552 CIRCE, and by the  Foundation of Research and Technology - Hellas Synergy Grants Program (project MagMASim).  

	
	
	
	
	
	






















        




