Supernova remnants SNRs have long been suspected to be the primary sources of Galactic cosmic rays. Over the past decades great strides have been made in the modelling of particle acceleration magnetic field amplification and escape from SNRs. Yet while many SNRs have been observed in nonthermal emission in radio Xrays and gammarays there is no evidence for any individual object contributing to the locally observed flux. Here we propose a particular spectral signature from individual remnants that is due to the energydependent escape from SNRs. For young and nearby sources we predict fluxes enhanced by tens of percent in narrow rigidity intervals given the percentlevel flux uncertainties of contemporary cosmicray data such features should be readily detectable. We model the spatial and temporal distribution of sources and the resulting distribution of fluxes with a Monte Carlo approach. The decision tree that we have trained on simulated data is able to discriminate with very high significance between the null hypothesis of a smooth distribution of sources and the scenario with a stochastic distribution of individual sources. We suggest that this cosmicray energydependent injection time CREDIT scenario be considered in experimental searches to identify individual SNRs as cosmicray sources.
['Astrophysics - High Energy Astrophysical Phenomena', 'Astrophysics - Astrophysics of Galaxies']
Investigating the CREDIT history of supernova remnants as cosmicray sources
{'Investigating the CREDIT history of supernova remnants as cosmic-ray sources': 'Anton Stall, ∗ Chun Khai Loo, † and Philipp Mertsch ‡ Institute for Theoretical Particle Physics and Cosmology (TTK), RWTH Aachen University, 52056 Aachen, Germany (Dated: September 18, 2024)  \nSupernova remnants (SNRs) have long been suspected to be the primary sources of Galactic cosmic rays. Over the past decades, great strides have been made in the modelling of particle acceleration, magnetic field amplification, and escape from SNRs. Yet, while many SNRs have been observed in non-thermal emission in radio, X-rays, and gamma-rays, there is no evidence for any individual object contributing to the locally observed flux. Here, we propose a particular spectral signature from individual remnants that is due to the energy-dependent escape from SNRs. For young and nearby sources, we predict fluxes enhanced by tens of percent in narrow rigidity intervals; given the percent-level flux uncertainties of contemporary cosmic-ray data, such features should be readily detectable. We model the spatial and temporal distribution of sources and the resulting distribution of fluxes with a Monte Carlo approach. The decision tree that we have trained on simulated data is able to discriminate with very high significance between the null hypothesis of a smooth distribution of sources and the scenario with a stochastic distribution of individual sources. We suggest that this cosmic-ray energy-dependent injection time (CREDIT) scenario be considered in experimental searches to identify individual SNRs as cosmic-ray sources.', 'INTRODUCTION': "Cosmic rays with energies up to at least E knee ≃ 3 PeV are usually assumed to be of Galactic origin [1-3]. Supernova remnants (SNRs) [4] have long been considered the prime candidate for the origin of these Galactic cosmic rays (GCRs) for a number of reasons: First, observations of SNRs oftentimes show power-law spectra, e.g. in radio [5], X-rays [6] or gamma-rays [7]. Second, diffusive shock acceleration (e.g. [8]) can operate at the SNRs' blast waves, thus providing a likely acceleration mechanism. Finally, the locally observed energy density of cosmic rays can be explained if on average ∼ 10 % of the shock kinetic energy is converted to GCRs [9]. Yet, observational evidence for SNRs to accelerate cosmic rays up to E knee is scarce [2].  \nFrom a modelling perspective, it is not clear how SNRs can accelerate particles to a few PeV either. Adopting shock speeds of O (10 4 ) km s -1 and a (turbulent) magnetic field strength typical of the interstellar medium, i.e. ∼ 1 µ G, it can be shown [10] that the maximum particle energy is of the order of 10 TeV which falls short of E knee by some three orders of magnitude. Attaining energies of E knee instead would require significant amplification of the turbulent magnetic field. This amplification is believed to be achieved thanks to the non-resonant hybrid instability (also known as Bell instability) [11], which is driven by the macroscopic current of escaping cosmic rays.  \nIn the following, we lay out the concordance view of time-dependent acceleration at and escape from SNR shocks, e.g. [12-15]. The amplification of magnetic fields is time-dependent due to its dependence on environmental parameters, most importantly the shock speed. The non-resonant hybrid instability operates most efficiently at high shock speeds, which can be achieved in the ejecta-  \ndominated phase of the SNR evolution. Consequently, the maximum energy, or better rigidity R max 1 , is expected to be attained at t Sed , the start of the SedovTaylor phase, when the amount of matter accumulated by the shock equals the ejecta mass. For t > t Sed , the shock speed decreases and so does the maximum rigidity R max ( t ). Only particles of rigidities R < R max ( t ) can be confined by the magnetised turbulence across the shock, and so particles of rigidity R max ( t ) escape at time t . The trend continues down to a rigidity R b (e.g. 10 TV), below which particles can be confined without magnetic field amplification. All particles with R < R b escape at the same time t life when the supernova remnant shock dissipates. One could parametrise this time-dependence of cosmic-ray escape with a source term, that is the number of GCRs released from a source per unit time and rigidity, of the form Q ( R , t ) ≡ δ ( t -t esc ( R )) Q R ( R ) where t esc ( R ) is the inverse function of R max ( t ) and Q R ( R ) is an arbitrary function of rigidity, but likely close to powerlaw form. This concordance scenario has been widely adopted in the literature and is also supported by numerical simulations [16, 17]. We will refer to this as the Cosmic Ray Energy-Dependent Injection Time (CREDIT) scenario.  \nThe details of the cosmic-ray escape are usually not considered when trying to explain the locally observed GCR measurements. Instead, the rate of injection of GCRs into the interstellar medium is modelled as the product of a smooth source density n src ( ⃗r ) and a steady rate per unit rigidity ¯ Q ( R ). This rate can be considered the temporal average of the time-dependent source term  \n¯ Q ( R ) ≡ ν ∫ ∞ 0 d t δ ( t -t esc ( R )) Q R ( R ) = νQ R ( R ) . (1)  \nHere, ν is the Galactic rate of supernovae.  \nIt is hence natural to ask if there are observable consequences of the CREDIT scenario for locally measured cosmic rays. If one adopts a distribution of sources that is spatially smooth and steady, as is usually assumed, e.g. in numerical finite difference codes (see, however, [18]) then there are none, that is observationally the timedependent escape scenario cannot be distinguished from a time-independent scenario.  \nIn reality, SNRs are individual objects with spatial and temporal extents much smaller than the distances and times over which cosmic rays propagate. The smooth and steady assumption can be justified in situations where a large number of sources contribute to the flux of GCRs. In such a situation, the measured flux for a particular distribution of sources will differ little from the mean of an ensemble of such source distributions. This mean is identical to the prediction from the smooth source distribution due to the linearity of GCR transport. For instance, it is easy to estimate that for standard parameters of GCRs, the number of sources contributing to the flux at a 10 GV is about 100 times larger than at 100 TV, where only some tens of sources cause the bulk of the flux. This justifies the use of the smooth scenario at GV energies.  \nThere are other situations, however, when the assumption is not justified. Examples include the transport of GCR electrons and positrons above a few hundred GeV [19, 20] or of GCRs below a GV [21]; in both cases, the ranges are limited by energy losses. Consequently, the spectrum becomes very sensitive to the actual distribution of sources. For certain realisations of the source distribution, this can lead to deviations from the prediction of a smooth model and might be observable. Given our limited knowledge of the true source distances and ages, we have to resort to Monte Carlo simulations in which the total flux is simulated for different realisations of sources drawn from a smooth distribution [16, 22-24]. We refer to models of this sort as stochastic models.  \nIn the following, we argue that the CREDIT scenario put into the context of a stochastic model can give large and, more importantly, observable spectral features at TV rigidities. We illustrate this by computing the fluxes of GCR protons for a large ensemble of source distributions. We characterise the spectral features and employ machine learning to reliably discriminate the spectral fluctuations in a CREDIT scenario from the fluctuations due to statistical errors in a smooth model.", 'METHODOLOGY': "To study the influence of source stochasticity on the locally observed GCR spectrum, we consider a simplified model of the Galaxy. We assume that all sources lie in the Galactic disk ( z = 0) and cosmic rays diffuse through the Galactic halo that extends to z = ± H . We draw the source ages from a uniform distribution adopting the canonical supernova rate ν of 0 . 03 yr -1 [25, 26] and the source positions from an axisymmetric distribution just depending on galactocentric distance [27], assuming that the Sun is located at radius R ⊙ .  \nWe will only consider GCR protons with rigidities above some GV such that inelastic collision and advection can be ignored. We also assume an isotropic diffusion coefficient that follows a power-law in rigidity κ ( R ) = κ 0 β ( R ) R δ 2 . So, the transport equation is  \n∂ψ R ( t, x , R ) ∂t -κ ( R ) ∇ 2 ψ R ( t, x , R ) = Q ( R , x , t ) , (2)  \nwhere ψ R ( t, x , R ) = dn/d R denotes the isotropic CR density ( n is the number density). It is related to the differential flux Φ R = ( d 4 n ) / ( d R dAdtd Ω) = v/ (4 π ) ψ R and to the phase-space density f = f ( x , p, t ) through ψ R = (4 πp 2 Ze ) /c f .  \nWe assume that a single population of SNRs inject all GCR protons up to the knee rigidity R knee . In the CREDIT scenario, the particles' escape times are rigidity-dependent. This escape time t esc ( R ) = t 0 + ∆ t esc ( R ) depends on the time of the supernova explosion t 0 and the time it takes a proton of a certain rigidity to escape the source ∆ t esc ( R ). We assume that the latter follows a broken power-law (as indicated by [17]):  \n∆ t esc ( R ) = t Sed ( R R knee ) a , a = ln ( t life /t Sed ) ln ( R b / R knee ) , (3)  \nfor R > R b and t esc ( R ) = t life otherwise. Choosing R b = O (10 TV) is motivated by the maximum rigidity that can be confined in SNR environments according to [10]. Setting R b → ∞ lets all protons escape at the same time. We will refer to this scenario, which is usually considered as the standard injection in stochastic modelling [22-24], as burst-like injection.  \nThe source injection term on the RHS of equation (2) for a single source at ( t i , x i ) is given by  \nQ i ( R , x , t ) = Q R ( R ) δ (3) ( x -x i ) δ ( t -t esc ,i ( R )) , (4)  \nwhere Q R ( R ) is the time-integrated source spectrum as described in equation (1) and is the same for all sources. Motivated by diffusive shock acceleration [8], it has a  \nTABLE I. Fiducial values for the GCR proton spectrum  \npower-law dependence on rigidity given by Q R ( R ) ∝ R -2 . 2 . Note that the deviation from the canonical R -2 spectrum can be explained as due to the relative speed of the scattering centres [28, 29] or shock obliquity [30]. As the source spectrum normalisation only enters as an overall factor and we will only be interested in relative deviations of stochastic proton spectra from the ensemble mean, we do not need to fix it in this study.  \nWe solve the transport equation (2) under the assumption of two free escape boundary conditions, i.e. ψ ( z = ± H ) = 0, using the method of mirror charges. We do not consider a radial boundary condition, as we assume the observer to be sufficiently far away from the edge of the Galaxy such that the escape is dominantly in the z-direction. The Green's function of equation (2), that is the solution for a single source, is given by  \nG ( t, x ; t i , x i ) = Q R ( R ) (2 πσ 2 ) 3 / 2 e -( x i -x ) 2 2 σ 2 · ϑ ( z, σ 2 , z max ) (5)  \nwhere σ 2 ( R , t ; t i ) = 2 κ ( R ) ( t -t esc ,i ( R )). The Green's function contains the power-law source spectrum, a Gaussian, and a correction function ϑ , which accounts for the free escape boundary condition. This function is an infinite sum related to the Jacobi theta function [22]. It has values in [0 , 1] and gets much smaller than 1 if z ≈ ± H or if the diffusion length ( ≈ σ ) is comparable with H . We exclude source injections that lie outside the observer's past light-cone to remedy potential non-causal solutions of the diffusion equation [23].  \nPredictions for the local spectrum of GCR protons can be obtained by summing the Green's functions of tens of millions of sources, whose positions and ages are sampled from the respective distributions described above. This embarrassingly parallel problem can be efficiently solved on graphical processing units (GPUs). To this end, we have used the python package jax [31]. Pronounced jump-like features originating in the time-dependent escape can be seen in sizeable parts of the realisations. Locally measured GCR fluxes should thus reveal information about the escape of cosmic rays from sources.", 'RESULTS': "In the left panel of Fig. 1, we show a number of example spectra for a CREDIT scenario where we adopted the fiducial parameters presented in Tab. I. To better illustrate the fluctuations, we present the spectra for ten random realisations divided by the ensemble mean for R b = 1TV. One can see fluctuations around the average spectrum for all rigidities, but their character changes markedly at R b : For R < R b , the spectral features are very smooth and typically only extend a few percent above or below the average value. For R > R b instead, the spectral features are much more peaked, and the maxima of individual features can be as large as a few times the average spectrum.  \nOf course, given the limited event statistics, experimental bins are typically wider than these spectral features. In the middle panel of Fig. 1, we, therefore, present the same example spectra as in the left panel but sampled only at the central rigidities of bins employed by experimental collaborations. Specifically, for rigidities below and above ∼ 2 TV, we have adopted the published rigidity bins of the AMS-02 and DAMPE collaborations, respectively. While the features in the DAMPE range are not as high and as nicely resolved as before, the most dramatic spectral features are still easy to make out by eye. We checked that integrating a finely resolved spectrum over rigidity bins gives similar statistics of the fluxes.  \nStatistical errors also arise due to the finite number of registered events. For a large enough number of events per bin, the distribution of the measured fluxes will be Gaussian and the measured fluxes in different bins can be considered as independent. We use reported statistical errors of the calorimetric cosmic-ray experiments AMS-02 [32] and DAMPE [33] as standard deviations in the individual rigidity bins. We found that the relation σ stat [%] = 0 . 42 · ( R [TV]) 0 . 6 approximates the statistical errors of both experiments well. In the right panel of Fig. 1, we show the spectrum from the smooth model divided by the average, a flat line, plus ten random draws from a distribution of statistical errors only. Generally, the statistical fluctuations are much smaller than the expected features under the CREDIT scenario, rendering their detection possible.  \nHowever, it is conceivable that statistical fluctuations can produce a flux similar to a realisation with less dramatic features. We trained a decision tree (DT) to discriminate between the null hypothesis of a smooth source distribution (plus statistical errors) and the alternative hypotheses of a CREDIT ( R b ∈ [1 , 23] TV) or a burstlike ( R b → ∞ ) scenario 3 . The DT aims to split the  \nFIG. 1. GCR proton fluxes (normalised to the ensemble mean) predicted for the CREDIT scenario with R b = 1TV. The alternating white and grey vertical lines indicate the rigidity bins. Left panel: The coloured lines show 10 random realisations, adopting a resolution in rigidity much higher than currently realised by experiments. Middle panel: Same random realisation, but sampled at the centres of the rigidity bins employed by AMS-02 and DAMPE. Right panel: Fluxes in a smooth scenario plus statistical errors due to finite event numbers of the sizes published by AMS-02 and DAMPE.  \n<!-- image -->  \nFIG. 2. Confusion matrices (normalised to rows/ true labels for the percentages) with 20 distinct R b ∈ [1 , 23] TV for the CREDIT model  \n<!-- image -->  \ndata, the flux value in each rigidity bin, scattered in 81dimensional space, into hypercuboids that contain realisations of the same label [34]. This is done by continuously splitting one hypercuboid (and the dataset) into two smaller hypercuboids. This can be thought of as an if-statement, where if the value of a realisation at the bin is smaller than a fixed value, it is in one hypercube and otherwise in the other. Making a prediction for a new data point is just travelling through a number of ifstatements and obtaining the label at the end. We used  \nthe sklearn package [35] to train and validate our DT. No tuning of hyperparameters was required to improve the results. We estimate the significance of our hypothesis test and its power through the so-called confusion matrix. The confusion matrix shows how often an element of the validation set with a true label is classified with a predicted label. A diagonal matrix corresponds to a perfect classification. Using 10-fold cross-validation on our 4 · 10 6 realisations for each scenario, we find an almost diagonal confusion matrix, Fig. 2, which shows the strength of the DT to reliably distinguish the three models.  \nWe checked for the effect of the uncorrelated part of the systematic error by adding in quadrature an uncorrelated, relative systematic error of 1% to the statistical errors [36, 37]; the classifier's performance did not change significantly. We also considered a relaxation of the peaked injection of each rigidity at a single point in time and found that smoothing the injection by a Gaussian kernel with a conservative width of 1 kyr does not alter the results for R b ≳ 1 TV.", 'SUMMARY AND CONCLUSION': "We considered a stochastic source model for simulating the GCR proton spectrum. At rigidities above some TV, the spectrum is determined by relatively few sources. We discussed that rigidity-dependent escape times of protons from their sources (CREDIT scenario) enhance the observability of single source contributions in the spectrum. Such an escape model, motivated by the same physical mechanisms that are usually invoked to explain parti-  \ncle acceleration in Galactic SNRs beyond O (10 TV), predicts a high likelihood for source configurations that produce pronounced peaked features in the proton spectrum. The emergence of such peaks is characteristic of CREDIT models (see Fig. 1). Here, we used a decision tree to discriminate spectra produced by a CREDIT or burst-like source injection from a power-law spectrum and a smooth source distribution. The classifier works extremely well, reliably deciding whether a spectrum is compatible with a smooth source injection or a CREDIT/burst-like model.  \nWe suggest applying our classifier to experimental data in the future. If the data get classified as CREDIT-like, one can determine the time of the supernova event t i for an observed spectral feature, if other model parameters are assumed, i.e. diffusion coefficient κ ( R ), break rigidity R b , Sedov time t Sed , and lifetime of sources t life . Alternatively, for a candidate supernova time t i , one can constrain a combination of κ ( R ), R b , t Sed , and t life that fit a detected spectral feature. Instead, if the data get classified as burst-like, this corresponds to a lower bound on R b , which informs models of acceleration. Finally, if the classifier finds a smooth model, this should trigger a reckoning with the SNR paradigm for GCRs.  \nFor this, however, the cosmic-ray model parameters have to be chosen carefully such that no correlated spectral features attributed to the transport model emerge in the normalised spectra (like Fig. 1). This can be done using extensive fits like e.g. [38]. Deviations from the fiducial values given in Tab. I will change the spectral features, enhancing them whenever the effective number of contributing sources decreases (e.g. lower source rate, smaller halo height, or enhanced diffusion).  \nProducing extensive Monte Carlo simulation samples using state-of-the-art programming paradigms opens up the possibility to train classifiers to reliably distinguish features of different models in data. We expect this to be a powerful tool to constrain GCR models. The possibilities of this will be explored in future work, but are already demonstrated by the classification power to discriminate CREDIT against burst-like models.  \nAs the number of contributing sources falls with rigidity, extensions of high precision measurements beyond ≈ 100 TV can enhance the predictive power of the classifiers. One proposed space based direct detection experiment is AMS-100 [39] which proposedly will extend high precision proton measurements beyond the cosmicray knee. The leading high-precision direct detection experiment DAMPE is currently limited by statistics and suffers from increasingly limited power to discriminate proton and helium at rigidities beyond ≈ 100 TV [40]. The inclusion of helium spectra in our model and the use of the p+He spectra provided by DAMPE can also increase the predictive power of our model. In addition, an unbinned analysis would further improve the sensitivity to the spectral features. We postpone such a study to future work.  \n- ∗ stall@physik.rwth-aachen.de\n- † khai.loo@rwth-aachen.de\n- ‡ pmertsch@physik.rwth-aachen.de\n- [1] S. Navas et al. (Particle Data Group), Phys. Rev. D 110 , 030001 (2024).\n- [2] S. Gabici, C. Evoli, D. Gaggero, P. Lipari, P. Mertsch, E. Orlando, A. Strong, and A. Vittino, Int. J. Mod. Phys. D 28 , 1930022 (2019), arXiv:1903.11584 [astro-ph.HE].\n- [3] M. Kachelriess and D. V. Semikoz, Prog. Part. Nucl. Phys. 109 , 103710 (2019), arXiv:1904.08160 [astroph.HE].\n- [4] S. P. Reynolds, ARA&A 46 , 89 (2008).\n- [5] D. A. Green, Journal of Astrophysics and Astronomy 40 , 36 (2019), arXiv:1907.02638 [astro-ph.GA].\n- [6] F. Seward, P. Slane, R. Smith, B.-C. Koo, J.-J. Lee, J.-J. Lee, J. Hagler, L. Portolese, and R. Milaszewski, Chandra Supernova Remnant Catalog, https://hea-www. harvard.edu/ChandraSNR/ (2010), [Online; accessed 22 August 2024].\n- [7] F. Acero et al. (Fermi-LAT), ApJS 224 , 8 (2016), arXiv:1511.06778 [astro-ph.HE].\n- [8] M. A. Malkov and L. O. Drury, Reports on Progress in Physics 64 , 429 (2001).\n- [9] V. L. Ginzburg and S. I. Syrovatskii, The Origin of Cosmic Rays (Pergamon, 1964).\n- [10] P. O. Lagage and C. J. Cesarsky, A&A 125 , 249 (1983).\n- [11] A. R. Bell, MNRAS 353 , 550 (2004).\n- [12] D. Caprioli, P. Blasi, E. Amato, and M. Vietri, MNRAS 395 , 895 (2009), arXiv:0807.4261 [astro-ph].\n- [13] K. M. Schure and A. R. Bell, Mon. Not. Roy. Astron. Soc. 435 , 1174 (2013), arXiv:1307.6575 [astro-ph.HE].\n- [14] A. Bell, K. Schure, B. Reville, and G. Giacinti, Mon. Not. Roy. Astron. Soc. 431 , 415 (2013), arXiv:1301.7264 [astro-ph.HE].\n- [15] P. Cristofari, P. Blasi, and E. Amato, Astropart. Phys. 123 , 102492 (2020), arXiv:2007.04294 [astro-ph.HE].\n- [16] P. Blasi and E. Amato, J. Cosmology Astropart. Phys. 2012 , 010 (2012), arXiv:1105.4521 [astro-ph.HE].\n- [17] D. Caprioli, P. Blasi, and E. Amato, MNRAS 396 , 2065 (2009), arXiv:0807.4259 [astro-ph].\n- [18] T. A. Porter, G. Johannesson, and I. V. Moskalenko, Astrophys. J. 887 , 250 (2019), arXiv:1909.02223 [astroph.HE].\n- [19] M. Pohl and J. A. Esposito, Astrophys. J. 507 , 327 (1998), arXiv:astro-ph/9806160.\n- [20] P. Mertsch, J. Cosmology Astropart. Phys. 2018 , 045 (2018), arXiv:1809.05104 [astro-ph.HE].\n- [21] V. H. M. Phan, F. Schulze, P. Mertsch, S. Recchia, and S. Gabici, Phys. Rev. Lett. 127 , 141101 (2021), arXiv:2105.00311 [astro-ph.HE].\n- [22] P. Mertsch, J. Cosmology Astropart. Phys. 2011 , 031 (2011), arXiv:1012.0805 [astro-ph.HE].\n- [23] Y. Genolini, P. Salati, P. D. Serpico, and R. Taillet, A&A 600 , A68 (2017), arXiv:1610.02010 [astro-ph.HE].\n- [24] C. Evoli, E. Amato, P. Blasi, and R. Aloisio, Phys. Rev. D 104 , 123029 (2021), arXiv:2111.01171 [astro-ph.HE].\n- [25] S. van den Bergh and G. A. Tammann, ARA&A 29 , 363 (1991).\n- [26] G. A. Tammann, W. Loeffler, and A. Schroeder, ApJS 92 , 487 (1994).\n- [27] K. M. Ferri'ere, Reviews of Modern Physics 73 , 1031  \n- (2001), arXiv:astro-ph/0106359 [astro-ph].\n- [28] A. R. Bell, Mon. Not. Roy. Astron. Soc. 182 , 443 (1978).\n- [29] D. Caprioli, C. C. Haggerty, and P. Blasi, Astrophys. J. 905 , 2 (2020), arXiv:2009.00007 [astro-ph.HE].\n- [30] A. R. Bell, K. M. Schure, and B. Reville, Mon. Not. Roy. Astron. Soc. 418 , 1208 (2011), arXiv:1108.0582 [astroph.HE].\n- [31] J. Bradbury, R. Frostig, P. Hawkins, M. J. Johnson, C. Leary, D. Maclaurin, G. Necula, A. Paszke, J. VanderPlas, S. Wanderman-Milne, and Q. Zhang, JAX: composable transformations of Python+NumPy programs (2018).\n- [32] M. Aguilar et al. (AMS), Phys. Rep. 894 , 1 (2021).\n- [33] Q. An et al. (DAMPE), Science Advances 5 , eaax3793 (2019), arXiv:1909.12860 [astro-ph.HE].\n- [34] C. Bishop, Pattern Recognition and Machine Learning , Information Science and Statistics (Springer, 2006).\n- [35] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cour-  \n- napeau, M. Brucher, M. Perrot, and E. Duchesnay, Journal of Machine Learning Research 12 , 2825 (2011).\n- [36] L. A. Cavasonza, H. Gast, M. Kramer, M. Pellen, and S. Schael, Astrophys. J. 839 , 36 (2017), [Erratum: Astrophys.J. 869, 89 (2018)], arXiv:1612.06634 [hep-ph].\n- [37] M. Aguilar et al. (AMS), Phys. Rev. Lett. 114 , 171103 (2015).\n- [38] G. Schwefer, P. Mertsch, and C. Wiebusch, ApJ 949 , 16 (2023), arXiv:2211.15607 [astro-ph.HE].\n- [39] S. Schael, A. Atanasyan, J. Berdugo, T. Bretz, M. Czupalla, B. Dachwald, P. von Doetinchem, M. Duranti, H. Gast, W. Karpinski, T. Kirn, K. Lubelsmeyer, C. Ma˜na, P. S. Marrocchesi, P. Mertsch, I. V. Moskalenko, T. Schervan, M. Schluse, K. U. Schroder, A. Schultz von Dratzig, C. Senatore, L. Spies, S. P. Wakely, M. Wlochal, D. Uglietti, and J. Zimmermann, Nuclear Instruments and Methods in Physics Research A 944 , 162561 (2019), arXiv:1907.04168 [astro-ph.IM].\n- [40] F. Alemanno et al. (DAMPE), Phys. Rev. D 109 , L121101 (2024), arXiv:2304.00137 [astro-ph.HE]."}