Calculation of conventional and prompt lepton fluxes at very high energy

An efficient method for calculating inclusive conventional and prompt atmospheric leptons fluxes is presented. The coupled cascade equations are solved numerically by formulating them as matrix equation. The presented approach is very flexible and allows the use of different hadronic interaction models, realistic parametrizations of the primary cosmic-ray flux and the Earth's atmosphere, and a detailed treatment of particle interactions and decays. The power of the developed method is illustrated by calculating lepton flux predictions for a number of different scenarios.


Introduction
Cosmic rays entering the Earth's atmosphere produce a multitude of secondary particles in interactions with air nuclei. Some of the secondary particles decay into muons and neutrinos, which are not absorbed in the atmosphere and can reach particle detectors at ground level. The spectra of these leptons contains not only information about the primary cosmic rays, but also about the particle physics of their production and the properties of the traversed atmosphere. Furthermore, searches for high-energy neutrinos from astrophysical sources have to cope with a large flux of atmospheric leptons as background. A better understanding of this flux, in particular its dependence on zenith angle and the changing properties of the atmosphere will help to develop improved methods to identify astrophysical neutrino fluxes and also contribute to a better understanding of hadronic interactions at high energy.
Many calculations of atmospheric lepton fluxes have been carried out since the early 1960's (e.g. [1], see [2] for a review). However, most of them incorporated approximations in solving the cascade equations that lead to increased uncertainties on the relation between the predicted fluxes and the physical input parameters, or the calculations could only be carried out in detail for just one or a few parameter/model combinations due to the large CPU time requirements. In this work we reduce the uncertainties related to the calculation method to a minimum by developing a numerical method with a level of detail comparable with Monte Carlo calculations [3][4][5]. In addition, the contributions of heavy flavor mesons and resonances to the flux of atmospheric leptons is accounted for in detail. The code is made publicly available . a e-mail: anatoli.fedynitch@cern.ch https://github.com/afedynitch/MCEq Our approach is based on the numerical solution of the coupled cascade equations, which have been rewritten into a matrix form to make use of modern implementations of linear algebra algorithms. While providing superior precision at very high energies where Monte Carlo methods are often statistically inefficient, the high performance of the algorithm allows us to perform calculations for many input parameter and model assumptions in a very short time. On an average portable computer it takes a few seconds to calculate lepton fluxes, while keeping most of relevant parameters accessible for users and easy to modify according to the current application.

Coupled cascade equation
The cascade equations for particle h can be written for one discrete energy bin E i It is part of a system of coupled ordinary differential equations, describing the evolution of the flux Φ of particles as a function of the atmospheric slant depth ing the mass density ρ, which is typically a function of the atmospheric height. The behavior of the particle cascade is driven by the competition of two source terms (1c), (1d) and two sink terms (1a), (1b). The interaction length λ h int,E i = m air /σ inel p−air (E i ) in units g/cm 2 [6] is independent of the slant depth and only varies slowly with energy due to the inelastic particle-air cross-section. The decay length λ h dec,E i (X) = cτ h E i ρ air (X)/m h is proportional to the lifetime τ h of particle h and can vary by orders of magnitude due to the relativistic time dilation. In our approximation new particles are created along the shower trajectory in hadronic interactions in (1c) or decays in (1d), with energy conservation restricting the range of possible source particles to energies E k equal or greater than E i . The interactions coefficients of particle l producing particle h, c l(E k )→h(E i ) , are obtained from hadronic interactions models by histogramming the particle yield as function of x (p) lab = E i /E k for l-air collisions. Suitable interaction models include, for example, QGSJET-II-04 [7] and EPOS LHC [8], and the upcoming versions of SIBYLL [9] and DPMJET-III [10]. Alternatively one can directly extrapolate results of fixed-target experiments on light nuclear targets using scaling arguments.
While in some cases it is possible to find analytical expressions for the decay coefficients d l(E k )→h(E i ) , see [6], numerical simulations have to be used for describing complex decays accurately. We have tabulated the decay coefficients as function of E k based on simulations with the Monte Carlo event generator Pythia 8 [11,12] that includes also rare decay channels and accounts for the effect of electroweak matrix elements where applicable.
The full system of equations of the hadronic cascade can be obtained by writing Eq. (1) for all possible types of hadrons and leptons. In the following we will concentrate on high lepton energies and neglect the interaction and/or decay terms of neutrinos and muons. To obtain, for example, the flux of atmospheric neutrinos at the surface, one needs to solve the full system taking into account the nonlinear X dependence of λ dec and the non-analytic forms of the particle spectra serving as input for the interaction coefficients c l→h .

Matrix form
An efficient numerical computing scheme can be found by rewriting the cascade equations into matrix form. We group the different Φ h (E i ) into a column vector Φ by writing blocks for each particle type for the discrete energy spectra The energy grid E i = 50 GeV · 10 i/N is logarithmically spaced with roughly 8 bins per decade of energy across the energy range of the calculation between 50 GeV and 10 10 GeV. We explicitly include more than 50 types of mesons, baryons and leptons. Together with some additional technical groups the dimension of Φ is then ∼ 6000. The reciprocal coefficients 1/λ of interaction and decay lengths are arranged in diagonal matrices The decay length matrix Λ dec is constructed analogously usingλ h dec,E i = λ h dec,E i (X)/ρ air (X) to factorize out the dependence on the air density. Sub-matrices containing the interaction coefficients are defined as and sub-matrices for the decay D l→h are constructed in a similar way. The full interaction and decay matrices C and D are built from these sub-matrices according to the order of particle types in Φ Using the definitions above, the matrix form of the coupled cascade equations can be written as

Short-lived particles
One goal of this work is to accurately take into account contributions of heavy flavor mesons and resonances to the flux of atmospheric leptons. Their short decay lengths introduce quickly decaying modes in Eqs. (1b) and (1d), which appear as eigenvalues of the decay matrix D with a large modulus of the negative real part. Fig. 1 illustrates the large difference between the decay lengths of conventional and charmed mesons. The energy dependence originates from time dilation. If only pions and kaons would be considered as intermediate mesons, the cascade equations would become interaction dominated above 1 TeV. The decay would be a slow process and the equations could be easily integrated using a moderate step size of O(1 g/cm 2 ). If short-lived particles are included, the choice of the appropriate step size is driven by the smallest eigenvalue, avoiding oscillations of fluxes around zero that are unphysical (Φ ≥ 0). For this reasons the equation system becomes a stiff numerical problem involving step sizes of To reduce the stiffness of the equation system, we introduce the resonance approximation. In the resonance approximation very short-lived particles (resonances), e.g. η or ρ mesons, decay immediately after their creation at the vertex. For each particle h this approximation is valid in a regime, where The coefficients for the chained production of a long lived secondary particle l via production and immediate decay of the resonance η, neglecting its interactions, can be written as To understand this result, let η int n be the vector containing the fluxes of all resonances, which are created during the integration step n. Using the matrix notation and forward Euler integration we can write C res h→η denotes a (k × d Φ ) matrix similar to C defined in Eq. (7), which contains production coefficients for resonances in interactions of hadrons. According to the approximation, all created resonances have to decay into ordinary particles within the same integration step. By writing this condition for a single resonance type k dec,e f f η k,n · ∆X n , we obtain the effective decay length dec,e f f = ∆X n = ρ(X) λ η k dec,e f f .
Using λ η k dec,e f f instead of the true decay length for shortlived resonances we make sure that all particles decay after one integration step in X without having to change the numerical treatment of the cascade equations. It should be noted that λ η k dec,e f f does not depend on the properties of the resonance. The contribution to the ordinary particle flux due to resonance decay is then where D res η→h is a (d Φ × k) matrix, containing decay coefficients of resonances into hadrons and leptons. By inserting Eq. (11) in (13) we replicate the expression from Eq. (10) for the production of particles via intermediate resonances and define the square (d φ × d φ ) resonance matrix R.
Chained decays proceeding through two or more resonances are governed by additional left multiplications of decay matrices. Extending the matrix form of the cascade equations (8) with intermediate resonance production results in At high energies, where the interaction of resonances becomes important. We introduce the parameter t mix = λ dec (E)/λ int (E), which is a threshold value, separating the energy regime where the particle can be treated as resonance from a regime where it has to be a full member of the cascade and listed in Φ. A reasonable value is t mix = 0.05. This complicates somewhat the procedure to fill the C, D and R matrices, where for each particle the individual threshold has to be taken into account as cut in row and/or column.

Initial state
The initial state of the cascade equation is the flux of cosmic rays at the top of the atmosphere. Since the calculation relies on the properties of the average air shower, the superposition theorem is sufficient to model the flux and composition. A cosmic ray nucleus is modeled as Z protons and A − Z neutrons, with each nucleon carrying the fraction E/A of the total kinetic energy. The relevant input for the initial condition is, therefore, the all-nucleon spectrum separated in the proton and neutron components. Alternatively one can use a single particle per energy bin and calculate particle yields at the surface for comparisons with full Monte Carlo methods such as [3][4][5]. Fig. 2 shows the spectra of three models we typically use for calculations. We emphasize, that it is crucial to include the knee and ankle in the flux calculations with respect to the  [15] has been often used for calculation of the prompt flux in the past. The poly-gonato model [16] focuses on the flux below and at the knee and it is not applicable at very high energies.  . Atmospheric density dependence on X, calculated using parameterizations for the US Standard Atmosphere [19] and the South Pole as implemented in CORSIKA [20], and the NRLMSISE-00 model. Solid lines represent a trajectory for θ = 0 • and dashed for θ = 70 • .
at tens of PeV. Improving the knowledge of the spectrum and composition as measured by air shower experiments would help to disentangle these ambiguities.

Geometry and Atmosphere
Treating the atmosphere in planar approximation, the relation between the height, slant depth, and local density can be taken directly from measurements. An often used approach, based on the idea by Linsley [17], is a parametrization of the relation between height and mass overburden X v (h) (slant depth for vertical trajectory) using 5 piecewise defined exponential functions, representing layers of the atmosphere. A higher flexibility is achieved if tabulated atmospheric data (e.g. from satellites) or detailed numerical models, such as NRLMSISE-00 [18], are used. At large zenith angles also the curvature of the surface of the Earth has to be accounted for. Therefore we compute and tabulate the relation ρ(h atm (X)) for each provided zenith angle θ and parametrization of the atmosphere. As shown in Fig. 4, this results in a linear smooth curve which is well suited for interpolation with splines. In Fig. 5 the ratio of the flux Φ, calculated with different models of the atmosphere to the flux calculated using the US Standard Atmosphere [19] is shown. Seasonal variations are of comparable magnitude in the CORSIKA parameterizations as in NRLSMSISE-00. Both models predict a seasonal variation of muon and neutrino rates at the order of ±10% in agreement with what IceCube has observed [21,22]. Another important feature is the nontrivial zenith angle dependence on atmospheric variations. A more detailed modeling of the atmosphere can be used to improve the experimental investigation of the prompt flux [23]. Furthermore, the flux at energies > PeV, which is dominated by prompt leptons, exhibits also a ∼ 15% dependence on the atmosphere.

Calculation of the prompt flux
In a related contribution [24], we discuss a model of charmed hadron production as it is implemented in the Monte Carlo model SIBYLL-2.3 RC1. For the calculation of the prompt flux using the method of this work, it is sufficient to take the interaction cross sections and Feynmanx F or the x Lab = E secondary /E projectile distributions from Monte Carlo. Two alternative models, where it is possible to extract x Lab distributions, are the Martin-Ryskin-Stasto (MRS) [25] model and the charm model of DPMJET-II.55 [26]. MRS considers perturbative production of charm quarks, based on a saturation model, and DPMJET includes contributions from non-perturbative, perturbative and fragmentation mechanisms, similar to SIBYLL-2.3. The comparison in Fig. 6 shows large differences in shape and cross-section between the models. MRS predicts a softer spectrum with a moderate forward crosssection. Due to a hard spectrum and the largest crosssection, we can expect that calculations using DPMJET will result in the highest fluxes. SIBYLL should produce more inclusive leptons compared to MRS, since the harder spectrum will yield charmed mesons at higher x (see discussion of spectrum weighted moments in [24]).
To examine the validity of charm production in DPM-JET, we compare the total cc cross-sections with ALICE data in Fig. 6 and the differential D ± -meson cross-section with forward data from LHCb in Fig. 7. The comparisons show that the charm model in DPMJET overestimates the cross-sections in the forward phase-space. DPMJET predictions are disfavored by LHC data. The recent measurements reduce the uncertainty of charm production for pp collisions at PeV (Lab) energies. However, a larger fraction of uncertainty comes also from nuclear effects. To assess this uncertainty in a quantitative way, we make the els of charm production 14 roduction duction based Inclusive cc cross-section in pp collsions. The ALICE measurement is corrected for the invisible part of the cross-section and extrapolated to full phase-space [27].
EPJ Web of Conferences Figure 7. Comparison of the differential D ± cross-section with data, measured in pp collsions by LHCb [28], with SIBYLL and DPMJET calculations.
assumption for the cc nuclear modification factor R p−air = dN cc p−air /dp T N coll dN cc pp /dp T ≡ 1.
In other words, there are no screening effects and the production of charmed quarks is a point-like process. Quantitatively this means for the inclusive cc cross-section in proton-air collisions σ cc,p−air = A air σ cc,p−p = 14.5 σ cc,p−p .
In calculations labeled SIBYLL-2.3 PL (PL for pointlike) we use the charm cross sections and distributions of SIBYLL-2.3 RC1 as predicted for pp interactions and scale the yields according to Eq. 17. In Fig. 8 the results of the different charm production models are compared with the Enberg-Reno-Sarcevic (ERS) [29] and the Thunman-Ingelman-Gondolo (TIG) [15] calculations of the prompt muon neutrino flux. SIBYLL-2.3 RC1 performs similarly to ERS and it produces notably higher fluxes than the MRS saturation model. All models predict different spectral shapes. Using DPMJET-II results in an order of magnitude higher fluxes. We consider the SIBYLL-2.3 PL predictions as upper boundary for nuclear uncertainties of charm production.

Partial contributions of intermediate particles
Since fluxes of all possible intermediate mesons and baryons are stored in the state vector Φ, we can easily trace back the mother particles of leptons at the surface. Fig. 9 is a break down of the different contributions of intermediate particles to the conventional and prompt flux. We define a lepton as prompt, if its mother particle of the last decay has a cτ < cτ(K 0 S ) = 2.68 cm. To improve the clarity of the graph, the simple broken power-law primary spectrum of TIG is employed as initial state. The dominant contributions to conventional muons are decays of charged pions and kaons, while prompt muons are originating from decays of charged and neutral D mesons and unflavored mesons, such as η, ω and φ. The latter resonances break the correlation between muon and neutrino fluxes at very high energies [30] as shown in the flavor ratios of Fig. 10. Prompt muon neutrino and electron neutrino fluxes are roughly equal and originate from decays of D-mesons and Λ + C baryons. The fractional contribution of D and Λ + C becomes equal at several hundreds of PeV. Decays of K ± are responsible for the largest fraction of conventional muon neutrinos at energies above a few TeV. For electron neutrinos, other channels are important, such as decays of K 0 L . At several hundreds of TeV there is an additional contribution from K 0 S as recently discussed in [31]. In any case, the flux from charmed particles is expected to be significantly higher than that due to these other channels.

Summary and Outlook
An efficient numerical treatment of cascade equations has been developed for the calculation of atmospheric lep-ISVHECRI 2014  ton fluxes at very high energy, with particular attention put on the transition from conventional to prompt production processes. As a first application we have shown calculations to illustrate the importance of the primary flux parametrization, the model of the atmosphere, the role of short-lived particles, and the model of charmed hadron production. More details about the numerical solution and the code will be published elsewhere.