The atmospheric lepton fluxes play a crucial role in many particle and astroparticle physics experiments e.g. in establishing the neutrino signal and the muon background for neutrino oscillation measurements or the atmospheric background for astrophysical neutrino searches. The Matrix Cascade Equations MCEQ code is a numerical tool used to model the atmospheric lepton fluxes by solving a system of coupled differential equations for particle production interaction and decay at extremely low computational costs. Previously the MCEQ framework accommodated only longitudinal development of air showers an approximation that works well for neutrino and muon fluxes at high energies O 10 GeV  and above. However for accurate calculations of atmospheric lepton angular distributions at lower energies the lateral component of hadronic cascades becomes significant. We introduce 2D MCEQ an efficient numerical approach for combined longitudinal and angular evolution of air showers that retains the low computational complexity. The accuracy of the 2D MCEQ is affirmed by its benchmark comparison with the standard Monte Carlo code CORSIKA. Our method can be used for twodimensional evolution of hadronic cascades in arbitrary media and paves the way for efficient threedimensional calculations of atmospheric neutrino fluxes.
['Astrophysics - High Energy Astrophysical Phenomena', 'High Energy Physics - Phenomenology']
Atmospheric lepton fluxes via twodimensional matrix cascade equations
{'Atmospheric Lepton Fluxes via Two-Dimensional Matrix Cascade Equations': "Tetiana Kozynets, 1, ∗ Anatoli Fedynitch, 2, 3 and D. Jason Koskinen 1  \n1 Niels Bohr Institute, University of Copenhagen, Blegdamsvej 17, 2100 Copenhagen, Denmark 2 Institute of Physics, Academia Sinica, Taipei City, 11529, Taiwan 3 Institute for Cosmic Ray Research, University of Tokyo,  \n5-1-5 Kashiwa-no-ha, Kashiwa, Chiba 277-8582, Japan  \n(Dated: April 4, 2024)  \nThe atmospheric lepton fluxes play a crucial role in many particle and astroparticle physics experiments, e.g. in establishing the neutrino signal and the muon background for neutrino oscillation measurements, or the atmospheric background for astrophysical neutrino searches. The Matrix Cascade Equations ( MCEq ) code is a numerical tool used to model the atmospheric lepton fluxes by solving a system of coupled differential equations for particle production, interaction, and decay at extremely low computational costs. Previously, the MCEq framework accommodated only longitudinal development of air showers, an approximation that works well for neutrino and muon fluxes at high energies ( O (10 GeV) and above). However, for accurate calculations of atmospheric lepton angular distributions at lower energies, the lateral component of hadronic cascades becomes significant. We introduce '2D MCEq ', an efficient numerical approach for combined longitudinal and angular evolution of air showers that retains the low computational complexity. The accuracy of the 2D MCEq is affirmed by its benchmark comparison with the standard Monte Carlo code corsika . Our method can be used for two-dimensional evolution of hadronic cascades in arbitrary media and paves the way for efficient three-dimensional calculations of atmospheric neutrino fluxes.", 'I. INTRODUCTION': "Interactions of cosmic rays with the atomic nuclei in the Earth's atmosphere produce cascades of secondary particles, referred to as the extensive air showers. These cascades have two components - electromagnetic (encompassing production and subsequent reinteraction of energetic electrons and photons) and hadronic (including production and subsequent reinteraction/decay of unstable mesons). One of the byproducts of hadronic cascades are neutrinos, which span the entire energy range from MeV to PeV and thereby form a broad landscape for probing fundamental physics. In particular, the GeV-scale neutrinos produced in the air showers constitute the main signal for atmospheric neutrino oscillation studies, including e.g. muon neutrino disappearance [1-4], tau neutrino appearance [5-9], and searches for physics beyond the Standard Model [10-13]. In addition, atmospheric neutrinos are an 'irreducible background' for astrophysical neutrino searches (e.g. [14-17]), which further strengthens the motivation for accurate modelling of neutrino production in the Earth's atmosphere.  \nThe unoscillated neutrino fluxes depend on several main inputs: the primary cosmic ray flux (including composition and spectrum [18-21]); the hadronic interaction model (prescribing the probabilities of secondary particle yields in hadron-nucleus collisions; e.g. [22-27]); decay probabilities and branching ratios of unstable particles [18, 28]; and model for the atmospheric density as a function of altitude for specific geographical locations [29-32]. At O (GeV) and sub-GeV energies, the angular  \ndistributions of the air shower secondary particles ('secondaries') are further affected by the Earth's magnetic field, which curves the trajectories of the charged cosmic ray primary particles ('primaries') and the secondary muons. In addition, the angular spread of the low-energy secondaries with respect to the primary particle axis becomes non-negligible: as the transverse momentum is Lorentz-invariant, the deflection angle grows with decreasing energy and can vary from a few degrees to tens of degrees at GeV-scale energies. Both of these effects are necessary to properly describe the angular evolution of individual air showers and the resulting full-sky angular distribution of the O (GeV) atmospheric neutrinos [33].  \nMonte Carlo simulations are the most natural approach to incorporate the many stages of the air shower modelling into a single computational framework [34, 35]. The Monte Carlo treatment implies that the generation and propagation of the cosmic ray primaries and the interactions and decays of the secondaries are executed on an event-by-event basis. These processes are stochastic and follow the probabilistic particle yield prescriptions of a given event generator . The most widely used realizations of this method to date include the general-purpose codes such as Geant4 [36], fluka [37], mcnp [38] and phits [39] for particle propagation in matter, as well as corsika [32] and aires [40] codes specialized in air shower evolution. The closed-source HKKMS [34] and Bartol [35] atmospheric neutrino flux models are also based on Monte Carlo simulations, employing event generators Jam + DPMJet [26, 41-43] and Target [44], respectively. The HKKMS model has been tuned to reproduce the muon flux data [41, 45] and is set out to incorporate the fixed-target experiment measurements of hadronic interaction yields in the future [46]. It is therefore commonly used as the baseline atmospheric neutrino  \nflux model in experimental analyses by e.g. the SuperKamiokande and the IceCube collaborations, as well as in projections for the upcoming Hyper-Kamiokande, JUNO, and DUNE experiments [47, 48]. Both the HKKMS and the Bartol flux models include the geomagnetic effects and the deflection of the secondaries from the primary axis, and therefore are the standard reference for O (GeV) atmospheric lepton flux calculations. While the Monte Carlo approach for the inclusive atmospheric lepton flux calculations provides high level of detail as an inherent advantage, it is computationally expensive, fairly complex, and lacks sufficient flexibility for extraction of systematic uncertainties (e.g. those related to the cosmic ray flux model and the hadronic interaction model parameters [20, 49, 50]).  \nAnother natural path towards the inclusive atmospheric neutrino flux modelling is via a solution to the cascade equations describing particle production, interaction, and decay in the atmosphere (see e.g. [30] for a broad introduction into this topic). Numerous studies have tackled these equations semi-analytically, with [51-53] being the latest developments. The semi-analytical method was further overtaken by the high-precision numerical solutions provided by the MCEq software [22, 54]. The MCEq cascade equation solver relies on the probabilities of the secondary particle yields in the interaction and decay processes extracted from event generators and stored as matrices. Avoiding the need to repeatedly run event generators within the user interface, MCEq enables computation of inclusive secondary particle fluxes on millisecond timescales, compared to several CPU-hours typically required by the Monte Carlo calculations.  \nDespite the significant speedup over the Monte Carlo approaches and the flexibility to study the impact of the systematic parameters, the MCEq code could not be readily used to predict the angular distributions of the O (GeV) cascade secondaries. The reason for this constraint is that MCEq was originally written in the 1D approximation of the air shower development, i.e., under the assumption of strictly collinear (with respect to the primary cosmic ray axis) secondary particle production and propagation.  \nIn this study, we are seeking to extend the MCEq framework with the angular evolution of the individual air showers. We develop the numerical technique and the practical implementation of a two-dimensional cascade equation solver, where the secondaries are allowed to deviate from the primary particle trajectory. Our code, '2D MCEq ' 1 , enables numerical computation of the resulting angular distributions of secondary particles. This advancement has broad applications in the analyses involving O (GeV) atmospheric leptons (and more generally, any hadronic cascade secondaries in arbitrary media) and contributes to future development of fully numerical  \nor hybrid three-dimensional calculations of atmospheric neutrino fluxes and air showers.  \nThis paper is structured as follows. In Sec. II A, we review the analytical cascade equations in the onedimensional approximation and further show how to incorporate the second (angular) dimension in Secs. II B and II C. The numerical (matrix) form of the 1D equations forming the basis of the MCEq code is reviewed in Sec. III A. We then derive the matrix form for the 2D equations by reformulating them in the frequency domain (Sec. III B). The pipeline of the 2D MCEq software is described in Sec. IV, which includes the steps to prepare the interaction/decay probability matrices (Secs. IV A and IVB) and the principles of the 2D cascade equation integrator (Secs. IV C and IV D). Finally, we compare the 2D MCEq angular distributions to those obtained with the corsika Monte Carlo in Sec. V, focusing on 2D cascades induced by a single cosmic ray primary.", 'A. One-dimensional cascade equations': "In this work, we employ the cascade theory to characterize the spatial development of the secondary particle showers induced by a single cosmic ray projectile. The mathematical basis of this theory is the system of coupled partial integro-differential cascade equations , which are a form of the Boltzmann transport equations for multiple particle species. For the secondary particle species h , we define the single-differential particle density n h with respect to the kinetic energy E :  \nn h ( E ) = d N h d E , (1)  \nwhich represents the number of particles N h per energy interval. In the one-dimensional cascade theory, this single-differential density evolves as a function of the atmospheric slant depth X :  \nX ( h o ) = w h o 0 d lρ air ( l ) , (2)  \nwhere h o is the observation altitude above the surface of the Earth, ρ air is the depth-dependent air density, and the integral is evaluated along the trajectory l of the shower core. With ρ air given in g cm -3 , and l taken in cm, the unit of X is g cm -2 . Then, the one-dimensional coupled  \ncascade equations [22, 30, 54] read  \nd n h ( E,X ) d X = -n h ( E,X ) λ int ,h ( E ) -n h ( E,X ) λ dec ,h ( E,X ) (3a)  \n-∂ ∂E ( µ E n h ( E,X )) (3b)  \n+ ∑ ℓ w ∞ E d E ℓ d N dec l ( E ℓ ) → h ( E ) d E n ℓ ( E ℓ , X ) λ int ,ℓ ( E ℓ ) (3c)  \n+ ∑ ℓ w ∞ E d E ℓ d N int ℓ ( E ℓ ) → h ( E ) d E n l ( E ℓ , X ) λ dec ,ℓ ( E ℓ , X ) . (3d)  \nThe 'sink' terms in Eq. 3a represent the decrease in the density of particle type h as the result of its interactions in the atmosphere after travelling the interaction length λ int ,h , or its decay after travelling the decay length λ dec ,h . Another sink term in Eq. 3b stands for the energy losses of the charged particles due to ionization, where µ E = -⟨ d E d X ⟩ is the average stopping power per unit length. The 'source' terms in Eq. 3c and Eq. 3d describe the increase of n h due to the interactions and decays of other particle species ℓ with energy E ℓ . The respective yields of the particle h are reflected in the differential production cross sections d N ℓ ( E ℓ ) → h ( E ) d E . The energy conservation constraint is given in the integral bounds ( r ∞ E ) of Eqs. 3c and 3d: it requires that the total energy E ℓ of the primary particle must be greater than, or equal to, the total energy E of the secondary particle.", 'B. Incorporating the second (angular) dimension': "In high-energy inelastic collisions or decays, the angular deflection θ ℓ → h of the secondary particles h from the primaries ℓ is minor ( ≪ 1 · at energies ≫ 10 GeV), justifying the use of the 1D approximation in the evolution of highenergy hadronic cascades [33]. In this regime, the velocity unit vector ˆ u ℓ → h of h translates to (0 , 0 , 1) ⊤ in a Cartesian coordinate system where the z axis aligns with ℓ . Lower energies necessitate consideration of the x and y components of ˆ u ℓ → h and explicit inclusion of the azimuthal angle φ ℓ → h . Then, the velocity vector becomes ˆ u ℓ → h = (sin θ ℓ → h cos φ ℓ → h , sin θ ℓ → h sin φ ℓ → h , cos θ ℓ → h ) ⊤ . To second order in θ , this can be approximated as  \nˆ u ℓ → h =    θ ℓ → h cos φ ℓ → h θ ℓ → h sin φ ℓ → h 1 -( θ ℓ → h ) 2 2    . (4)  \nSimilarly, the initial particle ℓ can be assigned a unit velocity vector ˆ u ℓ = ( θ ℓ cos φ ℓ , θ ℓ sin φ ℓ , 1 -( θ ℓ ) 2 2 ) ⊤ in a fixed-frame Cartesian coordinate system, where θ ℓ is the angle between the direction of ℓ and the z axis of this system, and φ ℓ is the respective azimuthal angle.  \nIn a Monte Carlo simulation, where the interactions or decays would be treated on an event-by-event basis, the direction ˆ u h of the secondary particle h in the fixed (lab) frame could be found via a simple addition of ˆ u ℓ and ˆ u ℓ → h . However, to incorporate angular evolution into the semi-analytical cascade theory, the distributions of the particle travel directions have to be formulated in terms of angular densities . Invoking the azimuthal symmetry, i.e., the invariance with respect to φ , we define the double-differential particle density η with respect to the energy and the polar angle as  \nη h ( E,θ ) = 1 θ d 2 N h ( θ ) d θ d E , (5)  \nwhich is normalized to the single-differential density:  \nn h ( E ) = w θ max 0 η h ( E,θ ) θ d θ. (6)  \nThroughout this study, we assume θ max = π/ 2, i.e. consider only forward-going particles, as well as the delta function-like angular distributions of the primaries. As the cascade develops, more secondaries will be produced off-axis and the angular distribution η h of the secondaries will evolve as a function of slant depth. On the other hand, the distribution of the relative angles between the primaries and the secondaries, θ l → h , is defined by the allowed phase space in a given interaction or decay process and constitutes a fixed convolution kernel . The angular distribution of the secondary particle h can then be obtained as a two-dimensional convolution of the primary angular density with the convolution kernel in the plane orthogonal to the primary particle direction. For the case of secondaries obtained in interactions, we denote this kernel as ς l → h and write  \nη h ( E,θ ) = η ℓ ( E,θ ℓ ) ∗ ∗ ς ℓ → h ( E ℓ , E, θ ℓ → h ) , (7a)  \n= w θ max 0 η ℓ ( θ ℓ ) ς ℓ ( E ℓ ,θ ℓ ) → h ( E,θ ) θ ℓ d θ ℓ , (7b)  \nwhere the ' ∗∗ ' operator represents two-dimensional convolution. The convolution kernels for decays will be denoted as δ l → h throughout this study. Following the formalism of [55] for the convolution of two azimuthally symmetric functions, we have absorbed the integration over the azimuthal variables into the definition of ς ℓ → h and δ ℓ → h . To illustrate the 2D convolution principle, we consider as an example a proton-induced hadronic cascade, as shown in Fig. 1. In this simplified setup, a beam of protons with the energy density n enters the atmosphere at the slant depth X 0 aligned with the downward-pointing z axis, hence θ primary = 0 . The direction of this proton beam is represented by the unit vector ˆ u primary . In the 1D geometry, the velocity unit vector ˆ u secondary of ν µ is aligned with ˆ u primary , while in the 2D geometry, this does not hold beyond X 0 . As the proton interacts with the atmospheric nuclei at X 1 , the secondary products of the interaction (including the π + ) gain transverse momentum,  \nand their angular distribution widens. This is represented by the convolution with the kernel ς p → π + . The angular distribution of muon neutrinos at X 2 further widens due to the convolution of the pion angular density with the decay kernel δ π + → ν µ .  \nMathematically, the production of the particle h with the energy E by the interactions of the primary ℓ with the energy E ℓ leads to the following change in the angular density of h :  \nd η int h ( θ ) d X = 1 λ int ,ℓ w π/ 2 0 η ℓ ( θ ℓ ) ς ℓ ( E ℓ ,θ ℓ ) → h ( E,θ ) θ ℓ d θ ℓ . (8)  \nAn equivalent expression can be formulated for decays by replacing ς ℓ → h with δ ℓ → h . The appearance of the θ ℓ factor in the integrals of Eq. 8 is an important feature of the 2D convolution in the xy plane, where θ ℓ and θ are  \ninterpreted as the radii of the ˆ u l and ˆ u h velocity vectors projected onto xy .", 'C. Two-dimensional cascade equations in the angular domain': 'Equation 8 and its equivalent for decays are the source terms in the two-dimensional cascade equations; they directly modify the angular densities of the secondaries, which evolve longitudinally as a function of the slant depth X in the atmosphere. At the same time, the sink terms in Eq. 3a and Eq. 3b do not change the angular distribution of the primaries and only contribute to the change in the overall normalization. With these two observations combined, we can write down the 2D version of Eq. 3 as follows:  \nd η h ( E,X,θ ) d X = -η h ( E,X,θ ) λ int ,h ( E ) -η h ( E,X,θ ) λ dec ,h ( E,X ) -∂ ∂E ( µ E η h ( E,X,θ )) + ∑ ℓ w π/ 2 0 θ ℓ d θ ℓ w ∞ E d E ℓ ς ℓ ( E ℓ ,θ ℓ ) → h ( E,θ ) λ int ,ℓ ( E ℓ ) η ℓ ( E ℓ , X, θ ℓ ) (9) + ∑ ℓ w π/ 2 0 θ ℓ d θ ℓ w ∞ E d E ℓ δ ℓ ( E ℓ ,θ ℓ ) → h ( E,θ ) λ dec ,ℓ ( E ℓ , X ) η ℓ ( E ℓ , X, θ ℓ ) .  \nThe longitudinal development of the secondary particle cascades is computed through the forward difference integration of Eq. 9. The new component in Eq. 9 compared to Eq. 3 is the angular development of the secondaries, which is taken care of via the 2D convolutions of the angular densities of the primaries with the interaction/decay convolution kernels.', 'III. MATRIX CASCADE EQUATIONS AND THE MCEQ CODE': 'The basic principle of the MCEq code is to evolve the hadronic and electromagnetic cascades in the atmosphere, given a cosmic ray primary flux and the probabilities of interactions and decays of all primary and secondary particles. With these inputs, MCEq solves the cascade equations for the densities of the secondaries of interest. We begin this section with a review of the numerical form of the one-dimensional equations (Eq. 3) based on the formalism derived in [22, 50, 54]. We further extend this numerical framework to two dimensions, building on [56].', 'A. Review of the matrix cascade equations in 1D': 'For a shower particle h which can interact with nuclei in the atmosphere (e.g. 14 N or 16 O), the interaction cross section is energy-dependent, as is the yield of the interaction products in an inelastic collision. Additionally, if the particle is unstable, the energy spectra of its decay products depend on the boost of the parent particle. It is therefore natural to discretize the transport equation in energy, i.e. to represent the particle densities in discrete energy bins E i , i ∈ [0 , N E -1]. The discrete one-dimensional cascade equation reads:  \nd n h E i ( X ) d X = -n h E i ( X ) λ h -n h E i ( X ) λ h ( X )  \nint ,E i dec ,E i (10a) -∇ i [ µ h E i n h E i ( X )] (10b) + ∑ ℓ ∑ E ∗ k ≥ E ∗ i c ℓ ( E k ) → h ( E i ) λ ℓ int ,E k n ℓ E k ( X ) (10c) + ∑ ℓ ∑ E ∗ k ≥ E ∗ i d ℓ ( E k ) → h ( E i ) λ ℓ dec ,E k ( X ) n ℓ E k ( X ) , (10d)  \nwhere we arranged the terms in the same order as in Eq. 3 to clarify the correspondences between the continuous and  \n/uni0302  \n/uni0302  \n/uni0302  \n/uni03B8  \nFIG. 1. Schematic development of a hadronic cascade ( p → π + → ν µ ) in the 1D (longitudinal-only) and the 2D (longitudinal + angular) geometries. In this diagram, the longitudinal propagation is performed over three discrete steps along the slant depth X for illustrative purposes. At each step in X , the angular distribution of the primaries from the previous step is shown as the dotted line, and the current angular distribution of the specified particle as the solid line. The distributions of secondaries get wider further down the chain due to the convolution with the kernels ς p → π + and δ π + → ν µ (see text for details).  \n<!-- image -->  \nsin  \n/uni0302  \n≈  \n/uni03B8  \nthe discrete equation versions. In Eq. 10c, we defined the coefficient c for the yield of particle h in interactions as  \nc ℓ ( E k ) → h ( E i ) = d N ℓ ( E k ) → h ( E i ) d E ∣ ∣ ∣ E = E i ∆ E k , (11)  \nwhich translates as the energy density of particles h with energy E i generated per primary ℓ within the energy bin E k . The decay coefficients d in Eq. 10d are defined in the same way. The equations for the different particle species are coupled through the yield coefficients c ℓ → h and d ℓ → h . The solution is obtained by solving Eq. 10 in X [22, 54] iteratively in the matrix form. The yield coefficients are derived from the event generators (e.g. UrQMD [25], DPMJet [26, 27], Sibyll [22, 23], or epos-lhc [24] for hadron-nucleus collisions, and Pythia [28] for decays) by histogramming the secondary particle yields as a function of the secondary particle kinetic energy, E i . In the 1D approximation, all secondary particle angles with respect to the primary particle direction of motion are contributing to the yield coefficient, thereby resulting in an angle-integrated interaction/decay probability.', 'B. 2D matrix cascade equations in the Hankel frequency domain': "Depending on the energy scales of hadronic interactions and unstable particle decays in the atmosphere, the widths of the angular distributions of the secondary particles can vary by orders of magnitude. For example, the average emission angle of 2 GeV pions produced in a collision of a 100 GeV is 10 · with respect to the primary proton direction, while 20 GeV pions deflect only by ∼ 1 · from the proton axis and produce ν µ at angles as small as 0 . 01 · relative to the pion direction. Due to the evolution of hadronic cascades over multiple generations (shower age), the total deflection is amplified by the number of generations even if the angular deflection in a single interaction/decay is small. As a result, Eq. 9 requires a 'universal' θ grid which could accommodate both large and small angular deflections. Making such a grid linear would imply an extremely fine discretization, and the numerical evaluation of the 2D convolution integrals would become prohibitively expensive. If the θ grid was logarithmic, the computation of the convolution integral would become more complicated due to the mis-alignment of the input and the output grids. While the techniques for convolving functions defined on logarithmic grids exist, they usually come with hyperparameters to be tuned  \nangular densities  \nby the user in order to keep the numerical errors to the minimum [57-59]. This extra freedom in the choice of hyperparameters could lead to unpredictable numerical behaviour in the integration of Eq. 9 over thousands of steps in X .  \nTo avoid the complications of the 2D convolutions in the θ space (which we will also refer to as the 'real' space), we choose to operate in the spectral ('frequency') domain instead. This is motivated by the existence of the convolution theorem , which transforms the convolutions in the real space into multiplications in the frequency space. For the specific case of the 2D convolution of the azimuthally symmetric functions ς ℓ → h , δ ℓ → h , and η ( X,θ ℓ ), the correct transform enabling the use of the convolution theorem is the zeroth-order Hankel transform H [55]:  \nH [ f ( θ )]( κ ) = w ∞ 0 f ( θ ) J 0 ( κθ ) θ d θ, (12)  \nwhere f ( θ ) is a function of the continuous variable θ , κ is the spectral frequency mode ( κ ≥ 0), and J 0 is the zeroth-order Bessel function of the first kind. In the formal definition of H , the upper limit of the θ integral in Eq. 12 is ∞ , however we only consider the forward-going particles with θ ≤ π/ 2.  \nThe convolution theorem states that, for the azimuthally symmetric functions f ( θ ) and g ( θ ),  \nH [ f ( θ ) ∗ ∗ g ( θ )] = H [ f ( θ )]( κ ) · H [ g ( θ )]( κ ) , (13)  \ni.e., the Hankel transform of the convolution result is a product of the Hankel transforms of the input functions in the frequency space [55]. Equation 13 is fully applicable to the two-dimensional cascade equations in Eq. 9. We therefore bring the convolution kernels and the angular densities of the cascade particles to the Hankel frequency space by defining their zeroth-order Hankel transforms as follows:  \n˜ η h E i ( X,κ ) ≡ H [ η h E i ( X,θ )]( κ ); (14a)  \n˜ ς ℓ ( E k ) → h ( E i ) ( κ ) ≡ H [ ς ℓ ( E k ) → h ( E i ) ( θ )]( κ ); (14b)  \n˜ δ ℓ ( E k ) → h ( E i ) ( κ ) ≡ H [ δ ℓ ( E k ) → h ( E i ) ( θ )]( κ ) . (14c)  \nThen, we can reformulate Eq. 9 as  \nd˜ η h E i ( X,κ ) d X = -˜ η h E i ( X,κ ) λ h int ,E i -˜ η h E i ( X,κ ) λ h dec ,E i ( X ) -∇ i [ µ h E i ˜ η h E i ( X,κ )] + ∑ E ∗ k ≥ E ∗ i ∑ ℓ [˜ ς ℓ ( E k ) → h ( E i ) · ˜ η l E k ]( κ ) λ ℓ int ,E k (15) + ∑ E ∗ k ≥ E ∗ i ∑ ℓ [ ˜ δ ℓ ( E k ) → h ( E i ) · ˜ η ℓ E k ]( κ ) λ ℓ dec ,E k ( X ) ,  \nwhich is the main equation to be solved in '2D MCEq '. For practical applications, the κ grid is made discrete and integer-valued. This implies that in Eq. 15, the  \nmultiplication of the Hankel-transformed convolution kernels and the Hankel-transformed angular densities of the primaries is performed elementwise with respect to the discrete frequency modes κ . The 1D MCEq equation (Eq. 10) is a special case of Eq. 15 for κ = 0, as J 0 (0) = 1 and Eq. 12 becomes equivalent to our earlier definition of the angular density normalization from Eq. 6. Therefore, Eq. 15 retains the computational complexity of Eq. 10, up to a linear scaling by the number of the frequency modes ( N κ ). One can then either choose to solve the N κ equations (one for each κ ) sequentially or in parallel, or to assemble the Hankel-transformed yield coefficients and angular densities into a more complex matrix structure. Our current implementation relies on the sequential solution of the N κ equations but can easily be adapted to the user's preference. In Secs. IV A and IV B, we explain how to arrive at the Hankel-transformed yield coefficients for the relevant interaction and decay channels, as well as provide further details on our choice of the κ grid where these coefficients are stored.", 'IV. 2D MCEQ PIPELINE AND SOLUTION SCHEME': 'The main computational advantage of the MCEq code compared to the Monte Carlo simulations comes from the pre-calculation of the particle yields, which is done outside of the user interface. The pre-tabulated interaction and decay coefficients are used to build the matrices for Eq. 10 (1D MCEq ) or Eq. 15 (2D MCEq ) during the code initialization. These coefficients are derived from the histogrammed kinematic properties of the secondary particles in hadronic interactions or decays, as evaluated via the chromo code [60]. The histograms for 1D MCEq include primary and secondary particle energies ( E prim and E sec ). For 2D MCEq , an additional Hankel mode κ is incorporated as a third dimension, which stores the angular densities of the secondaries in a compact form. Below, we describe how the output of any given event generator is used to populate the ( E prim , E sec , κ ) grid.', 'A. Compact representation of event information': "For the low-energy atmospheric neutrino flux calculations ( O (GeV) and below), we deploy a logarithmicallyspaced kinetic energy grid ranging from 10 MeV to 10 TeV 2 . This grid utilizes a bin width of ∆ log 10 E kin = 0 . 1, resulting in N E = 60 energy bins. The 2D MCEq code currently excludes electromagnetic cascades, deemed unnecessary for the evolution of hadronic cascades that  \nproduce O (GeV) atmospheric leptons 3 . Thus, we apply Eq. 15 only to hadrons and leptons (excluding the τ lepton), totaling H = 21 particle species: 6 baryons ( p/ ¯ p, n/ ¯ n , and Λ 0 / ¯ Λ 0 ), 5 mesons ( π ± , K ± , K 0 L , ), and 10 leptons ( µ ± R / L , µ ± , ν e / ¯ ν e , and ν µ / ¯ ν µ ). Muons contribute 6 species, where each µ ± includes two polarizations: lefthanded 'L', right-handed 'R', and an unpolarized component (denoted as µ ± without a subscript).  \nFor hadronic interactions, the chromo code [60] runs the following models: UrQMD [25], epos-lhc [24], Sibyll-2.3d [22, 23], and DPMJet-III 19.1 [26, 27]. Pythia 8.306 [28] is used for unstable particle decays. However, it cannot simulate the production of polarized muons in π ± and K ± two-body decays or the three-body decays of polarized muons, instead generating events in the spin-averaged phase spaces. Hence, muon polarization modeling is done separately, as outlined in Sec. IV E and detailed in Appendix A 3.  \nOur event generation and histogramming scheme is consistent across all interaction/decay channels. For every primary in the kinetic energy bin k , we assign the logarithmic center ( E k = √ E k ' · E k '' for [ E k ' , E k '' ]) as the primary's initial energy. This primary particle enters the chosen event generator with four-momentum p µ prim = ( E k , 0 , 0 , √ E 2 k +2 E k m prim ) ⊤ , moving along the positive z axis. A stationary nitrogen nucleus ( 14 N) is the target for simulating hadronic interactions. Using other atmospheric nuclei like 16 O has little effect on secondary particle yields. For decays, we set the unstable particle at rest and boost its decay products to the lab frame, recording the daughter energies as in 1D MCEq .  \nTo solve the 2D cascade equation (Eq. 15), we additionally need to compute the Hankel-transformed angular densities of secondary particles, i.e., ˜ ς l ( E k ) → h ( E i ) ( κ ). We illustrate this process using the muon neutrino production chain from Fig. 1, i.e., p + 14 N → π + + X ∗ → ν µ + µ + , where X ∗ denotes all other secondary particles and nuclear remnants. As angular distributions of all particles fill a discrete energy grid, we select the primary proton, secondary pion, and tertiary neutrino energy bins around 100GeV, 10GeV, and 4GeV, respectively. This illustrative choice reflects the characteristic energies for the low-energy neutrino production in air showers.  \nThe simulation chain begins with protons incident on 14 N. We choose the epos-lhc hadronic model and compute the yield of the secondary π + . The top left panel of Fig. 2 shows the distribution of angles θ p → π + that the secondary pions make with the primary proton axis. The number of entries in the histogram,  \nn π + ≡ n p · c p → π + , is equal to the total secondary pion yield in this interaction. Each pion contributes a delta function δ ( θ -θ j p → π + ) , j ∈ [0 , n π + -1] , to the angular density. The Hankel transform of the delta function has an analytical representation, H [ 1 a δ ( θ -a ) ] ( κ ) = J 0 ( κa ) 4 , which can populate the κ -grid immediately after the event generation. Summation of these Hankel-transformed delta functions approaches the Hankel transform of the secondary pions' underlying angular density. This density can be expressed via inverse Hankel transform as:  \nς p → π + ( θ ) = 1 n p H -1 [ ˜ ς p → π + ( κ ) ] ,  \n= 1 n p H -1   n π + ∑ j =1 J 0 ( κθ j p → π + )   ( θ ) .  \n(16a) (16b)  \nEquations 15 and 16 could theoretically extend to an infinite number of modes κ . In practice, a truncated κ -grid with 24 logarithmically spaced integer modes between 0 and 2000 suffices to accurately represent the angular distributions of GeV-scale atmospheric leptons. In the example shown in Fig. 2, two key observations validate this approach. Firstly, the inverse Hankel transform from Eq. 16 effectively represents the original pion angular distribution, thus demonstrating the utility of the Hankel transform for compacting angular densities of secondary particles. Secondly, the amplitudes of the higher-frequency modes with κ ≥ 100 are negligible, indicating a sufficiently broad angular distribution of the GeV-scale pions produced in proton14 N interactions. For sharper-edged distributions, like in the pion decay to muon neutrinos (middle panel of Fig. 2), the truncated κ grid may not sufficiently reconstruct the angular density, resulting in a characteristic 'ringing'. However, in a realistic air shower, the effect of this artefact is minimal due to the wider pion angular distribution.  \nThe final ν µ angular distribution from the Monte Carlo simulation chain is successfully reconstructed through the inverse transform of the convolution result in the Hankel space:  \nς p → ν µ ( θ ) = H -1 [ ˜ ς p → π + ( κ ) · ˜ δ π + → ν µ ( κ ) ] ( θ ) . (17)  \nThis demonstrates the application of the convolution theorem from Eq. 13. Notably, after successive convolutions, the final θ refers to the angle of the secondary particle with respect to the original shower-inducing primary.  \nglyph[triangleleft]  \nglyph[triangleleft]  \n6  \nglyph[triangleleft]  \nglyph[triangleleft]  \n6  \nglyph[triangleleft]  \nk  \nglyph[triangleleft]  \n6  \nglyph[triangleleft]  \nh  \nh  \nh  \ndms ne sgd rdbnmcGLYPH<246>qx oGLYPH<246>qshbkd  \ne  \nglyph[triangleleft]  \nθ↪  \n1  \n,  \n,  \nh  \n)  \n0  \n)  \nΓ  \nadenqd+ sgqntfg sgd enqvGLYPH<246>qc  \nhr  \nΛ  \n0GLYPH<246>  \ncn mns bgGLYPH<246>mfd sgd GLYPH<246>mftkGLYPH<246>q  \nh  \nGLYPH<213>  \nN  \n(  \nΓ  \nc  \nΓ  \nc  \nglyph[triangleleft]  \n×  \nπ  \nc  \nE  \nglyph[triangleleft]  \n-  \nglyph[triangleleft]  \n-  \nΓ  \nT  \nπ  \nΓ  \n(  \nc  \nglyph[triangleleft]  \ni  \ni  \n-  \nθ↪  \n1  \nGLYPH<0>  \n·  \nh  \nglyph[triangleleft]  \n-  \nglyph[triangleleft]  \n-  \nk  \ndms ne sgd rdbnmcGLYPH<246>qx oGLYPH<246>qshbkd  \np- 5- Sgd mdv bnlonmdms hm  \nqhdr GLYPH<246>mc  \nrhlokx bnmsqhatsd sn  \nhr sgd GLYPH<246>mftkGLYPH<246>q  \ncdudknoldms ne  \n(  \n∣  \n·  \nh  \n∣  \nadenqd+ sgqntfg sgd enqvGLYPH<246>qc  \nGLYPH<246>jdm  \n(  \nk  \n(  \n)  \nne uhGLYPH<246>  \nsgd 1C bnmun,  \nGLYPH<0>  \n,  \nbGLYPH<246>qd p- 5- Sgd mdv bnlonmdms hm cn mns bgGLYPH<246>mfd sgd GLYPH<246>mftkGLYPH<246>q  \n0GLYPH<246>  \n∣  \nglyph[triangleleft]  \n-  \ni  \n·  \nglyph[triangleleft]  \n-  \nh  \n∣  \nhshdr ne sgd oqhlGLYPH<246>qhdr vhsg sgd  \nhr sgd GLYPH<246>mftkGLYPH<246>q  \nqhdr GLYPH<246>mc  \njdqmdkr-  \nc sgd  \nf  \nhmcdw hm  \nGLYPH<246>jdm  \nbGLYPH<246>qd  \nne uhGLYPH<246>  \n+ sghr dptGLYPH<246>shnm fhudr sgd dun,  \nhshdr ne sgd oqhlGLYPH<246>qhdr vhsg sgd  \ndmdqfx cdmrhsx ne GLYPH<246>kk oGLYPH<246>qshbkd  \ne  \nc sgd  \nf  \nhmcdw hm  \nc  \nπ  \njdqmdkr-  \nglyph[triangleleft]  \nh  \nGLYPH<213>  \nNglyph[triangleright]  \nδ  \n(  \nd rkGLYPH<246>ms  \ncdosg  \n+ sghr dptGLYPH<246>shnm fhudr sgd dun,  \nq ne oGLYPH<246>qshbkdr ne sxod  \nd dptXshnmr hm sgd GXmidk  \nc  \nN  \ndmdqfx cdmrhsx ne GLYPH<246>kk oGLYPH<246>qshbkd  \nd rkGLYPH<246>ms  \ncdosg  \nbx cnlXhm  \ne  \nglyph[triangleleft]  \nc  \nπ  \nh  \nGLYPH<213>  \nNglyph[triangleright]  \nδ  \n(  \nπ  \nGLYPH<213>  \nNglyph[triangleright]  \nδ  \n(  \nq ne oGLYPH<246>qshbkdr ne sxod  \ne  \nglyph[triangleleft]  \nh  \nd dptXshnmr hm sgd GXmidk  \ne  \nhr  \nΛ  \nglyph[triangleleft]  \nc  \nN  \n×  \nbx cnlXhm  \nh  \nrbGLYPH<246>kdr  \nne gGLYPH<246>cqnmhb hmsdqGLYPH<246>bshnmr  \ne  \nglyph[triangleleft]  \nh  \nGLYPH<213>  \nN  \n(  \nπ  \nGLYPH<213>  \nNglyph[triangleright]  \nδ  \n(  \nr hm sgd GLYPH<246>slnrogdqd+ sgd vhcsg  \nhr  \nΛ  \ndms ne sgd rdbnmcGLYPH<246>qx oGLYPH<246>qshbkd  \ne  \nglyph[triangleleft]  \n×  \nm ne sgd rdbnmcGLYPH<246>qx  \nrbGLYPH<246>kdr  \nne gGLYPH<246>cqnmhb hmsdqGLYPH<246>bshnmr  \nadenqd+ sgqntfg sgd enqvGLYPH<246>qc  \nh  \nGLYPH<213>  \nN  \n(  \nd knv,dmdqfx rdbnmcGLYPH<246>qhdr vhkk  \nr hm sgd GLYPH<246>slnrogdqd+ sgd vhcsg  \np- 5- Sgd mdv bnlonmdms hm  \ndms ne sgd rdbnmcGLYPH<246>qx oGLYPH<246>qshbkd  \nGLYPH<246>c  \nGLYPH<246>ants  \nm ne sgd rdbnmcGLYPH<246>qx  \nhr sgd GLYPH<246>mftkGLYPH<246>q  \nadenqd+ sgqntfg sgd enqvGLYPH<246>qc  \nqfx nmdr vhkk mns cdGLYPH<222>dbs eqnl  \nd knv,dmdqfx rdbnmcGLYPH<246>qhdr vhkk  \nGLYPH<246>jdm  \nbGLYPH<246>qd  \nne uhGLYPH<246>  \np- 5- Sgd mdv bnlonmdms hm  \nkx- Etqsgdqlnqd+ sgd dunktshnm  \nGLYPH<246>c  \nGLYPH<246>ants  \nhshdr ne sgd oqhlGLYPH<246>qhdr vhsg sgd  \nhr sgd GLYPH<246>mftkGLYPH<246>q  \nhydGLYPH<246>akd  \nrkGLYPH<246>ms  \nqfx nmdr vhkk mns cdGLYPH<222>dbs eqnl  \nc sgd  \nf  \nhmcdw hm  \n+ sghr dptGLYPH<246>shnm fhudr sgd dun,  \njdqmdkr-  \nGLYPH<246>jdm  \nbGLYPH<246>qd  \nne uhGLYPH<246>  \nkGLYPH<246>q  \nchrsqhatshnmr ne rdbnmcGLYPH<246>qhdr  \nc sgd  \nf  \nhmcdw hm  \nkx- Etqsgdqlnqd+ sgd dunktshnm  \nhshdr ne sgd oqhlGLYPH<246>qhdr vhsg sgd  \nm hm GLYPH<246>  \nrhmfkd hmsdqGLYPH<246>bshnm.cdbGLYPH<246>x  \ndmdqfx cdmrhsx ne GLYPH<246>kk oGLYPH<246>qshbkd  \nd rkGLYPH<246>ms  \ncdosg  \nhydGLYPH<246>akd  \nrkGLYPH<246>ms  \njdqmdkr-  \nDp- 5 cdlGLYPH<246>mcr  \n+ sghr dptGLYPH<246>shnm fhudr sgd dun,  \nq ne oGLYPH<246>qshbkdr ne sxod  \nkGLYPH<246>q  \nchrsqhatshnmr ne rdbnmcGLYPH<246>qhdr  \nlncGLYPH<246>sd  \nansg kGLYPH<246>qfd  \ndmdqfx cdmrhsx ne GLYPH<246>kk oGLYPH<246>qshbkd  \nd dptXshnmr hm sgd GXmidk  \nbx cnlXhm  \nm hm GLYPH<246>  \nrhmfkd hmsdqGLYPH<246>bshnm.cdbGLYPH<246>x  \nπ  \nhmf rtbg GLYPH<246>  \ne  \nglyph[triangleleft]  \nDp- 5 cdlGLYPH<246>mcr  \nh  \nGLYPH<213>  \nNglyph[triangleright]  \nδ  \n(  \nd rkGLYPH<246>ms  \ncdosg  \nd dptXshnmr hm sgd GXmidk  \n×  \nbx cnlXhm  \nq ne oGLYPH<246>qshbkdr ne sxod  \nrbqdshyGLYPH<246>shnm+ GLYPH<246>mc  \nlncGLYPH<246>sd  \nansg kGLYPH<246>qfd  \nπ  \nktshnm hmsdfqGLYPH<246>kr vntkc adbnld  \ne  \nglyph[triangleleft]  \nhr  \nΛ  \ne  \nglyph[triangleleft]  \nhmf rtbg GLYPH<246>  \nh  \nGLYPH<213>  \nN  \n(  \nrbGLYPH<246>kdr  \nne gGLYPH<246>cqnmhb hmsdqGLYPH<246>bshnmr  \nh  \nGLYPH<213>  \nNglyph[triangleright]  \nδ  \n(  \nm sgd nsgdq gGLYPH<246>mc+ he sgd  \nr hm sgd GLYPH<246>slnrogdqd+ sgd vhcsg  \n×  \nm ne sgd rdbnmcGLYPH<246>qx  \nrbqdshyGLYPH<246>shnm+ GLYPH<246>mc  \ndms ne sgd rdbnmcGLYPH<246>qx oGLYPH<246>qshbkd  \noqnctbs ne bnmunktshnm vntkc  \nhr  \nΛ  \ne  \nglyph[triangleleft]  \nGLYPH<213>  \nN  \n(  \nrbGLYPH<246>kdr  \nne gGLYPH<246>cqnmhb hmsdqGLYPH<246>bshnmr  \nd knv,dmdqfx rdbnmcGLYPH<246>qhdr vhkk  \nktshnm hmsdfqGLYPH<246>kr vntkc adbnld  \nadenqd+ sgqntfg sgd enqvGLYPH<246>qc  \nGLYPH<246>r  \nsgd ch  \nΦ  \nr hm sgd GLYPH<246>slnrogdqd+ sgd vhcsg  \nGLYPH<246>c  \nGLYPH<246>ants  \nh  \np- 5- Sgd mdv bnlonmdms hm  \nm sgd nsgdq gGLYPH<246>mc+ he sgd  \ndmrhsx+ dudm he sgd fqhcr ne sgd  \ndms ne sgd rdbnmcGLYPH<246>qx oGLYPH<246>qshbkd  \nhr sgd GLYPH<246>mftkGLYPH<246>q  \noqnctbs ne bnmunktshnm vntkc  \nm ne sgd rdbnmcGLYPH<246>qx  \nqfx nmdr vhkk mns cdGLYPH<222>dbs eqnl  \nadenqd+ sgqntfg sgd enqvGLYPH<246>qc  \nGLYPH<246>jdm  \nbGLYPH<246>qd  \nne uhGLYPH<246>  \nsgd sdbgmhptdr enq bnmunkuhmf  \nd knv,dmdqfx rdbnmcGLYPH<246>qhdr vhkk  \nkx- Etqsgdqlnqd+ sgd dunktshnm  \nGLYPH<246>r  \nsgd ch  \nΦ  \np- 5- Sgd mdv bnlonmdms hm  \nhshdr ne sgd oqhlGLYPH<246>qhdr vhsg sgd  \ndmrhsx+ dudm he sgd fqhcr ne sgd  \nhr sgd GLYPH<246>mftkGLYPH<246>q  \nGLYPH<246>c  \nGLYPH<246>ants  \nhydGLYPH<246>akd  \nrkGLYPH<246>ms  \nqfx nmdr vhkk mns cdGLYPH<222>dbs eqnl  \nkGLYPH<246>q  \nchrsqhatshnmr ne rdbnmcGLYPH<246>qhdr  \nsgd sdbgmhptdr enq bnmunkuhmf  \njdqmdkr-  \nGLYPH<246>jdm  \nbGLYPH<246>qd  \nne uhGLYPH<246>  \nhshdr ne sgd oqhlGLYPH<246>qhdr vhsg sgd  \nkx- Etqsgdqlnqd+ sgd dunktshnm  \nm hm GLYPH<246>  \nrhmfkd hmsdqGLYPH<246>bshnm.cdbGLYPH<246>x  \nhydGLYPH<246>akd  \nrkGLYPH<246>ms  \nDp- 5 cdlGLYPH<246>mcr  \njdqmdkr-  \nd dptXshnmr hm sgd GXmidk  \nkGLYPH<246>q  \nchrsqhatshnmr ne rdbnmcGLYPH<246>qhdr  \nlncGLYPH<246>sd  \nansg kGLYPH<246>qfd  \nm hm GLYPH<246>  \nrhmfkd hmsdqGLYPH<246>bshnm.cdbGLYPH<246>x  \nhmf rtbg GLYPH<246>  \nbx cnlXhm  \nDp- 5 cdlGLYPH<246>mcr  \nd dptXshnmr hm sgd GXmidk  \nrbqdshyGLYPH<246>shnm+ GLYPH<246>mc  \nbx cnlXhm  \nlncGLYPH<246>sd  \nansg kGLYPH<246>qfd  \nktshnm hmsdfqGLYPH<246>kr vntkc adbnld  \nrbGLYPH<246>kdr  \nne gGLYPH<246>cqnmhb hmsdqGLYPH<246>bshnmr  \nhmf rtbg GLYPH<246>  \nm sgd nsgdq gGLYPH<246>mc+ he sgd  \nr hm sgd GLYPH<246>slnrogdqd+ sgd vhcsg  \nrbqdshyGLYPH<246>shnm+ GLYPH<246>mc  \noqnctbs ne bnmunktshnm vntkc  \nm ne sgd rdbnmcGLYPH<246>qx  \nrbGLYPH<246>kdr  \nne gGLYPH<246>cqnmhb hmsdqGLYPH<246>bshnmr  \nktshnm hmsdfqGLYPH<246>kr vntkc adbnld  \nGLYPH<246>r  \nsgd ch  \nΦ  \nd knv,dmdqfx rdbnmcGLYPH<246>qhdr vhkk  \nr hm sgd GLYPH<246>slnrogdqd+ sgd vhcsg  \nm sgd nsgdq gGLYPH<246>mc+ he sgd  \ndmrhsx+ dudm he sgd fqhcr ne sgd  \nGLYPH<246>c  \nGLYPH<246>ants  \nm ne sgd rdbnmcGLYPH<246>qx  \noqnctbs ne bnmunktshnm vntkc  \nGLYPH<246>r  \nsgd ch  \nvd bgnnrd sn nodqGLYPH<246>sd hm sgd rodbsqGLYPH<246>k GLYPH<213>[eqdptdmbxGLYPH<248>( cnlGLYPH<246>hm  \nSn GLYPH<246>unhc sgd bnlokhbGLYPH<246>shnmr ne sgd 1C bnmunktshnmr hm  \nπ  \nΓ  \nk  \nΓ  \nglyph[triangleleft]  \nk  \nΓ  \nglyph[triangleleft]  \ni  \nk  \nh  \ni  \nk  \nhmrsdGLYPH<246>c-  \nHm oGLYPH<246>qshbtkGLYPH<246>q+  \nvd vhrg sn sqGLYPH<246>mrenql sgd ch  \nΦ  \ntrhnm  \nκ  \n(  \nk  \nems  \nπ  \nglyph[triangleleft]  \ncda  \nπ  \nglyph[triangleleft]  \nsgd  \nδ  \nroGLYPH<246>bd GLYPH<213>vghbg vd vhkk GLYPH<246>krn qdedq sn GLYPH<246>r sgd [qdGLYPH<246>kGLYPH<248> roGLYPH<246>bd(+  \nπ  \nΓ  \n≥  \ni  \ne  \nΓ  \nglyph[triangleleft]  \ni  \nk  \nglyph[triangleleft]  \nk  \ni  \nglyph[triangleleft]  \nπ  \ni  \nGLYPH<213>  \nNglyph[triangleright]  \nδ  \nk  \n(  \nδ  \nk  \n(  \nδ  \nk  \nc  \nδ  \nk  \nc  \nδ  \nk  \nk  \nglyph[triangleright]  \njdqmdkr GLYPH<246>mc sgd oqhlGLYPH<246>qx oGLYPH<246>qshbkd GLYPH<246>mftkGLYPH<246>q  \ncdmrhshdr hm rtbg  \nvd bgnnrd sn nodqGLYPH<246>sd hm sgd rodbsqGLYPH<246>k GLYPH<213>[eqdptdmbxGLYPH<248>( cnlGLYPH<246>hm  \nκ  \nGLYPH<213>  \nN  \n(  \nsgd cdbqdGLYPH<246>rd hm sgd nudqGLYPH<246>kk mnqlGLYPH<246>khyGLYPH<246>shnm-  \nBnlahmhmf sghr  \nhmrsdGLYPH<246>c-  \nHm oGLYPH<246>qshbtkGLYPH<246>q+  \nvd vhrg sn sqGLYPH<246>mrenql sgd ch  \nΦ  \ntrhnm  \nGLYPH<246>  \nvGLYPH<246>x  \nsgGLYPH<246>s  \nsgdhq bnmunktshnm hm sgd qdGLYPH<246>k roGLYPH<246>bd hr qdokGLYPH<246>bdc  \ni  \nk  \ncda  \nπ  \nglyph[triangleleft]  \nκ  \n4  \nk  \n(  \ne  \nΓ  \nglyph[triangleleft]  \nh  \nπ  \nΓ  \n(  \nk  \nglyph[triangleleft]  \n≥  \nnardquGLYPH<246>shnm vhsg Dp- 3+ vd bGLYPH<246>m vqhsd cnvm sgd 1C udqrhnm  \ni  \nGLYPH<213>  \nN  \n(  \nπ  \ni  \nGLYPH<213>  \nNglyph[triangleright]  \nδ  \nk  \n(  \nδ  \nk  \nc  \nδ  \nk  \nglyph[triangleright]  \nax ltkshokhbGLYPH<246>shnm hm sgd eqdptdmbx roGLYPH<246>bd- Hm nsgdq vnqcr+  \n<!-- image -->  \nrbqdshyGLYPH<246>shnm+ GLYPH<246>mc  \nsgd mtldqhbGLYPH<246>k  \nqfx nmdr vhkk mns cdGLYPH<222>dbs eqnl kx- Etqsgdqlnqd+ sgd dunktshnm hydGLYPH<246>akd rkGLYPH<246>ms cdosgr bGLYPH<246>m gGLYPH<246>ud GLYPH<246> kGLYPH<246>q chrsqhatshnmr ne rdbnmcGLYPH<246>qhdr m hm GLYPH<246> rhmfkd hmsdqGLYPH<246>bshnm.cdbGLYPH<246>x Dp- 5 cdlGLYPH<246>mcr GLYPH<246> [tmhudqrGLYPH<246>kGLYPH<248> lncGLYPH<246>sd ansg kGLYPH<246>qfd GLYPH<246>mc rlGLYPH<246>kk hmf rtbg GLYPH<246> fqhc khmdGLYPH<246>q vntkc hm oqGLYPH<246>bshbd GLYPH<246>mc sgdqdenqd qdrsqhbs ntq GLYPH<246>ssdmshnm sn δ ≡ ν ↪ 1=ookxhmf sgd GGLYPH<246>mjdk sqGLYPH<246>mrenql sn ansg φ k ≥ e ↪Γ k ≥ e GLYPH<246>mc π GLYPH<213> Nglyph[triangleright] δ k (+ vd bGLYPH<246>m cd ff md GLYPH<208> π e glyph[triangleleft] h GLYPH<213> Nglyph[triangleright] θ ( × ∈ W π e glyph[triangleleft] h GLYPH<213> Nglyph[triangleright] δ (GLYPH<210>GLYPH<213> θ ( GLYPH<246>mc GLYPH<208> φ k Γ glyph[triangleleft] i ( ≥ e Γ glyph[triangleleft] h ( GLYPH<213> θ ( × ∈ W φ k Γ glyph[triangleleft] i ( ≥ e Γ glyph[triangleleft] h ( GLYPH<213> δ (GLYPH<210>GLYPH<213> θ (- Sgdm+ trhmf sgd bnmunktshnm sgdnqdl+ vd qdenqltkGLYPH<246>sd Dp- 5 GLYPH<246>r c GLYPH<208> π e glyph[triangleleft] h GLYPH<213> Nglyph[triangleright] θ ( c N ; -GLYPH<208> π e glyph[triangleleft] h GLYPH<213> Nglyph[triangleright] θ ( κ e -GLYPH<208> π e glyph[triangleleft] h GLYPH<213> Nglyph[triangleright] θ ( κ e GLYPH<213> N ( d knv,dmdqfx rdbnmcGLYPH<246>qhdr vhkk GLYPH<246>c GLYPH<246>ants sgd oqhlGLYPH<246>qx oGLYPH<246>qshbkd qfx nmdr vhkk mns cdGLYPH<222>dbs eqnl kx- Etqsgdqlnqd+ sgd dunktshnm hydGLYPH<246>akd rkGLYPH<246>ms cdosgr bGLYPH<246>m gGLYPH<246>ud GLYPH<246> kGLYPH<246>q chrsqhatshnmr ne rdbnmcGLYPH<246>qhdr m hm GLYPH<246> rhmfkd hmsdqGLYPH<246>bshnm.cdbGLYPH<246>x Dp- 5 cdlGLYPH<246>mcr GLYPH<246> [tmhudqrGLYPH<246>kGLYPH<248> cd ff mhshnm ne ∈ + sgd toodq khlhs ne sgd δ hmsdfqGLYPH<246>k hm Dp- 6 hr /similarequal + gnvdudq vd nmkx bnmrhcdq sgd enqvGLYPH<246>qc,fnhmf oGLYPH<246>qshbkdr hm oqGLYPH<246>bshbd GLYPH<246>mc sgdqdenqd qdrsqhbs ntq GLYPH<246>ssdmshnm sn δ ≡ ν ↪ 1=ookxhmf sgd GGLYPH<246>mjdk sqGLYPH<246>mrenql sn ansg φ k ≥ e ↪Γ k ≥ e GLYPH<246>mc π GLYPH<213> Nglyph[triangleright] δ k (+ vd bGLYPH<246>m cd ff md GLYPH<208> π e glyph[triangleleft] h GLYPH<213> Nglyph[triangleright] θ ( × ∈ W π e glyph[triangleleft] h GLYPH<213> Nglyph[triangleright] δ (GLYPH<210>GLYPH<213> θ ( GLYPH<246>mc GLYPH<208> φ k Γ glyph[triangleleft] i ( ≥ e Γ glyph[triangleleft] h ( GLYPH<213> θ ( × ∈ W φ k Γ glyph[triangleleft] i ( ≥ e Γ glyph[triangleleft] h ( GLYPH<213> δ (GLYPH<210>GLYPH<213> θ (- Sgdm+ trhmf sgd bnmunktshnm sgdnqdl+ vd qdenqltkGLYPH<246>sd Dp- 5 GLYPH<246>r sgd sdbgmhptdr enq bnmunkuhmf glyph[triangleleft] -i · glyph[triangleleft] -h k cda π glyph[triangleleft] i Φ trhnm jdqmdk GLYPH<246>mc sgd dmrhsx+ dudm he sgd fqhcr ne sgd sgd sdbgmhptdr enq bnmunkuhmf ) ∣ glyph[triangleleft] -i · glyph[triangleleft] -h ∣ k W GLYPH<208> Γ k Γ glyph[triangleleft] i ( ≥ e Γ glyph[triangleleft] h ( -GLYPH<208> π k glyph[triangleleft] i GLYPH<210>GLYPH<213> θ ( κ k cda π glyph[triangleleft] i GLYPH<213> N ( glyph[triangleright] FIG. 2. Top left : Angular distribution of secondary pions in p + 14 N → π + + X , derived from the epos-lhc event generator (solid blue line), and the inverse Hankel transform result (dashed red line) of the top right panel. Top right : Hankel transform of the secondary pions' angular distribution, calculated using Eq. 16. Middle : Similar to the top panels but evaluated for daughter ν µ in π + → µ + + ν µ decay. Bottom : Demonstrating the convolution theorem on the p → π + → ν µ chain: the inverse transform of the product of the Hankel transforms (top right and middle right) reproduces the angular density of the tertiary neutrinos (bottom left). The dotted black line in the bottom left panel indicates the properly normalized angular distribution of secondary pions for comparison. Energy settings for all considered particles are found in Table A1.  \nlncGLYPH<246>sd  \nansg kGLYPH<246>qfd  \nGLYPH<246>mc  \nrlGLYPH<246>kk  \nktshnm hmsdfqGLYPH<246>kr vntkc adbnld  \nhmf rtbg GLYPH<246>  \nfqhc khmdGLYPH<246>q  \nvntkc  \nm sgd nsgdq gGLYPH<246>mc+ he sgd  \nδ  \nfqhc  \nrbqdshyGLYPH<246>shnm+ GLYPH<246>mc  \nsgd mtldqhbGLYPH<246>k  \noqnctbs ne bnmunktshnm vntkc  \nktshnm hmsdfqGLYPH<246>kr vntkc adbnld  \nGLYPH<246>r  \nsgd ch  \nΦ  \ntrhnm jdqmdk GLYPH<246>mc sgd  \nm sgd nsgdq gGLYPH<246>mc+ he sgd  \nδ  \nfqhc  \nc  \nπ  \njdqmdkr GLYPH<246>mc sgd oqhlGLYPH<246>qx oGLYPH<246>qshbkd GLYPH<246>mftkGLYPH<246>q  \ncdmrhshdr hm rtbg  \nne Dp- 0 GLYPH<246>r enkknvr9  \nsgd cdbqdGLYPH<246>rd hm sgd nudqGLYPH<246>kk mnqlGLYPH<246>khyGLYPH<246>shnm-  \nBnlahmhmf sghr  \ne  \nglyph[triangleleft]  \nh  \nGLYPH<213>  \nNglyph[triangleright]  \nθ  \n(  \nc  \nN  \nGLYPH<208>  \ne  \nglyph[triangleleft]  \nems  \nπ  \nglyph[triangleleft]  \nh  \nκ  \ni  \nGLYPH<213>  \nNglyph[triangleright]  \nθ  \n(  \ne  \nems  \nπ  \nglyph[triangleleft]  \n·  \nglyph[triangleleft]  \n-  \nh  \nh  \nh  \nk  \n∣  \n∣  \nGLYPH<208>  \nπ  \nglyph[triangleleft]  \n-  \ncda  \nπ  \nglyph[triangleleft]  \ne  \nglyph[triangleleft]  \nGLYPH<208>  \nπ  \nκ  \ni  \nk  \nΓ  \nglyph[triangleleft]  \nWGLYPH<208>  \nφ  \n-  \nWGLYPH<208>  \nW  \nφ  \nk  \nΓ  \nglyph[triangleleft]  \nΓ  \nk  \nΓ  \nglyph[triangleleft]  \nGLYPH<208>  \n(  \ne  \ncda  \nπ  \nglyph[triangleleft]  \ni  \n(  \ni  \nh  \nGLYPH<213>  \nNglyph[triangleright]  \nθ  \n(  \n≥  \nh  \ne  \nΓ  \nglyph[triangleleft]  \nh  \nk  \nGLYPH<213>  \nN  \n(  \nems  \nπ  \nglyph[triangleleft]  \nh  \n(  \nh  \n(  \n-  \ni  \n-  \nκ  \n≥  \n(  \ne  \nΓ  \nglyph[triangleleft]  \ne  \nΓ  \nglyph[triangleleft]  \nκ  \nk  \nh  \n(  \n;  \n-  \n)  \n)  \nπ  \nGLYPH<208>  \nk  \nglyph[triangleleft]  \nk  \nglyph[triangleleft]  \nGLYPH<208>  \nGLYPH<208>  \nπ  \nπ  \nk  \ni  \ni  \nglyph[triangleleft]  \ni  \nGLYPH<210>GLYPH<213>  \nθ  \n(  \nGLYPH<210>GLYPH<213>  \nθ  \n(  \nGLYPH<210>GLYPH<213>  \nθ  \n(  \nk  \nh  \nΓ  \nΓ  \nπ  \nΓ  \n(  \nπ  \nGLYPH<213>  \nNglyph[triangleright]  \nδ  \nGLYPH<213>5(  \n4  \nGLYPH<213>7(  \nGLYPH<213>7(", 'B. 2D MCEq matrix production': 'We generate 10 million events per primary species and energy bin on the MCEq grid, following the methodology in Sec. IV A. Since different hadronic models have different valid energy ranges, we create interpolated matrices which allow to smoothly transition between the hadronic model choices. This emulates the approach used in codes such as corsika and is explained in Appendix A 2.  \nThe generation of secondary particle yields involves H = 21 primary particles in 2D MCEq . The yield of the secondaries is stored as a function of the secondary and primary kinetic energies (on a grid of length N E ) and the Hankel frequency mode κ . The resulting dimension of the 2D MCEq matrices is N κ × ( N E · H ) × ( N E · H ), equating to 24 × (1260 × 1260). These are treated as 24 separate sparse matrices when solving the 2D cascade equation in the Hankel frequency domain.', 'C. Solution in the Hankel space': "̸  \nThe convolution theorem transforms the 2D convolution into multiplication in the Hankel frequency domain. The amplitude of the primary angular distribution corresponding to the mode κ is multiplied by the amplitude of the interaction/decay kernel corresponding to exactly the same mode, i.e., the modes κ 1 and κ 2 are not coupled if κ 1 = κ 2 . This allows us to treat each of the N κ = 24 equations of the 2D MCEq completely independently and solve them using the strategy analogous to that of the 1D MCEq [22, 54]. This involves applying the forward Euler integrator to Eq. 15, i.e., performing the longitudinal evolution of the Hankel-transformed angular densities of the cascade secondaries in discrete slant depth steps ∆ X . This approach is best summarized in the matrix form:  \n˜ η ( X t +1 , κ ) = ˜ η ( X t , κ ) -∇ E [diag( µ ) · ˜ η ( X t , κ )] +∆ X t [ ( -I + C k ) Λ int + 1 ρ ( X t ) ( -I + D k ) Λ dec ] ˜ η ( X t , κ ) , (18)  \nwhere X t and X t +1 = X t +∆ X t are two consecutive slant depth values, and C κ and D κ are the slices of the yield coefficient 'cubes' ς l → h ( κ ) and δ l → h ( κ ) at the frequency mode κ . Following [22, 54], we also construct the diagonal matrices Λ int and Λ dec from interaction and decay lengths. Each diagonal entry corresponds to a specific particle h and energy bin E i , arranged similarly to the particle density vector ˜ η ( X,κ ) that we seek to evolve. After ˜ η ( X final , κ ) is computed, the final step is to reconstruct the angular densities of η ( X final , θ ) via the inverse Hankel transform.", 'D. Reconstruction of the real-space solutions': "After the final integration step, the 2D MCEq solver returns the state vector ˜ η ( X final , κ ). This includes the Hankel frequency space amplitudes for all cascade particles across the MCEq energy grid. The inverse Hankel transform enables the reconstruction of the secondary particle densities as a function of the angle θ relative to the primary particle axis:  \nη ( X , θ ) = H -1 [ ˜ η ( X , κ )] ( θ ) ,  \nfinal final (19a) ≡ w ∞ 0 ˜ η ( X final , κ ) J 0 ( κθ ) κ d κ, (19b) ≃ w κ max 0 ˜ η ( X final , κ ) J 0 ( κθ ) κ d κ. (19c)  \nAlthough the κ grid in 2D MCEq is discrete with logarithmically spaced modes, accurate computation of the integral in Eq. 19c requires a quasicontinuous κ range, which is achieved via spline interpolation of ˜ η ( X final , κ ). This can be done at any X t (0 ≤ X t ≤ X final ) along the integration path, should the solution at a particular slant depth/altitude be of interest to the user.  \nThe reconstruction of angular densities of high-energy secondaries (with kinetic energies of 10 GeV and above), as well as those created very early in the cascade evolution, must be treated with care if the starting angular distribution of the primaries is narrow (e.g. delta function-like). Direct application of Eq. 19 may lead to 'ringing' in the reconstructed angular densities for such secondary particles. Thus, it is recommended to apply Eq. 19 for secondary particles with energy < ∼ 10GeV at several kilometers into the atmosphere, and to use the 1D approximation at high energies/altitudes.", 'E. Modeling of muon transport': "Muons play a crucial role in air shower development and atmospheric neutrino flux calculations. They either decay in-flight, generating muon and electron (anti)neutrinos, or survive until the surface level, creating a source of background in neutrino detection.  \nMuon polarization occurs as atmospheric muons are produced in the two-body decays of π ± and K ± . In the parent meson rest frame, the muons are fully polarized, with their momenta perfectly aligned or antialigned with their spin direction. This affects the angular distributions and the energy spectra of neutrinos originating from muon decays [53]. A simplified representation of the muon population across the entire continuous range of helicities is achieved by including only six muon species into the 2D MCEq cascade equations (see Sec. IV A).  \nAnother phenomenon affecting muon transport is multiple scattering, which modifies the trajectories of the muons as they scatter with atmospheric nuclei. We implement a Gaussian approximation to the Moli'ere formalism for describing this effect [32, 61]. This results in an overall  \nwidening of the muon angular distributions compared to the case without multiple scattering.  \nAppendix A3 details the implementation of muon polarization and multiple scattering, as well as their impact on the results of 2D MCEq .", 'A. Simulation setup': "To validate our solutions to the two-dimensional matrix cascade equations via 2D MCEq , we use the corsika v7.75 Monte Carlo code [32] as a benchmark. We aim to compare the angular distributions of the GeV-scale atmospheric neutrinos and muons generated in the cosmicray induced air showers. All of our simulations are run for a single angle of incidence of the cosmic ray primary, and the secondary particle angular distributions are computed with respect to the primary particle axis. In the terminology of Fig. 1, we are comparing the distributions of arccos(ˆ u primary · ˆ u secondary ) between 2D MCEq and corsika .  \nTo ensure a fair comparison, we equalize the configuration settings between MCEq and corsika to the extent possible. Most importantly, we match the choice of hadronic interaction models by using UrQMD [25] as the low-energy model 5 and epos-lhc [24] as the highenergy model in both corsika and 2D MCEq . The transition energy between models is set to 150 GeV 6 . In Appendix B, we provide a comprehensive list of other relevant physics settings in MCEq and corsika .  \nOur setup for the benchmarking of lepton densities consists of a proton primary incident onto the Earth's atmosphere at an inclination angle θ 0 with respect to the negative z axis (downward direction). We test both vertical ( θ 0 = 0 · ) and inclined showers (30 · ≤ θ 0 ≤ θ max = 80 · ). The energy of the proton either is fixed or follows a spectrum with a power-law dependence (e.g. ∝ E -2 . 7 ). For each considered initial condition, we simulate ∼ 1 million events in corsika with different random seeds. This lets us gather enough statistics for the low-energy muons and neutrinos at several observation altitudes. The corresponding binned angular distributions are compared directly to the angular densities obtained with 2D MCEq by solving Eq. 15.", 'B. Results': "To provide a representative example of how the 2D MCEq solutions compare to the corsika Monte Carlo outputs, we choose the case of a 100 GeV proton primary incident at θ 0 = 30 · . In Fig. 3, we show the angular densities of neutrinos ( ν µ + ¯ ν µ ; ν e + ¯ ν e ) and muons ( µ -+ µ + ) at low energies (up to 5 GeV).  \nFrom Fig. 3, we find that the angular distributions of O (few GeV) leptons with respect to the proton primary axis are in a very good agreement between 2D MCEq and corsika . For neutrinos, the differences between the two codes are predominantly statistical and reach at most 10% in the tails of the distributions. This level of agreement holds across all altitudes and energy bins considered. For muons, a characteristic tilt of the corsika -toMCEq angular distribution ratio is observed at all altitudes, reaching ∼ 20% in the distribution tails. This is indicative of a bias of the corsika angular distribution towards smaller angles or the MCEq angular distributions towards larger angles. One possible explanation for this pattern is that all particles in 2D MCEq travel exactly the same distance to reach a specific depth, including those that deflect by as much as 20-30 · from the primary axis. However, muons with such large deflection angles naturally travel longer distances than those at 0 · , introducing a factor of 1 / cos θ increase in the integration step length. This means that muons with large deflections from the primary direction must lose more energy than currently modelled in 2D MCEq and migrate from the energy bins shown in Fig. 3 to the bins of lower energy. Qualitatively, this explains why corsika could have fewer muons at large angles. However, since the amount of energy lost is directly proportional to the distance travelled (∆ E ≃ ⟨ d E d X ⟩ ∆ X ), any discrepancies related to the energy loss are expected to accumulate with distance. This shows to a small degree in Fig. 3, where the tilt in the corsika -toMCEq ratio of the muon angular densities develops mildly as a function of altitude. While we cannot definitively connect the angle-dependent discrepancy between 2D MCEq and corsika with the simplified treatment of the angle-dependent propagation step length in 2D MCEq , the latter remains a relevant feature to be implemented in future iterations of the MCEq code. At present, we point out that the 2D MCEq angular densities still provide a very good overall representation of the corsika angular distributions, which can be seen from the agreement of the distribution moments ( ⟨ θ ⟩ and ⟨ θ 2 ⟩ ) in Fig. 4. The sub-degree level of difference observed in the angular distribution moments will not be possible to resolve under any realistic experimental resolution at GeV-scale energies, implying a negligible impact on experimental analyses.  \nWhile the main objective of the 2D MCEq code is to evolve angular distributions of the secondaries in addition to the energy spectra already provided by 1D MCEq , we compare the energy spectra from 2D MCEq to those from corsika to provide further validation to our approach.  \nN  \n2  \nd  \nFIG. 3. Angular distributions of atmospheric leptons in a proton-induced air shower ( E 0 = 100GeV, θ 0 = 30 · ), as computed numerically in 2D MCEq (solid lines) and simulated in the corsika Monte Carlo (filled histograms with errorbars). The angle θ on the x axis is the angle a given secondary makes with the direction of the primary proton. The different colors correspond to the different energy bands, and the bottom sub-panel in each plot shows the ratio of corsika ('C') to MCEq ('M').  \n<!-- image -->  \nTo obtain the energy spectra from 2D MCEq , we extract the κ = 0 mode from the Hankel-space solutions, which is equivalent to the angle-integrated particle densities as per Sec. III B. As seen in Fig. 5, the spectra from MCEq and corsika agree within a few % in the 1-10 GeV region, which is the main energy range of interest in this study. Above 10 GeV, the difference between the two codes grows as a function of energy, reaching ∼ 20% at the maximum neutrino energies available from 100 GeV primary showers. The energy dependence of the corsika -toMCEq ratio  \ncould point to the difference in the treatment of hadronic interactions, e.g. the hadron yields between the different UrQMD model versions. The same level of agreement is observed when the primaries have power law-like spectra, as demonstrated in Figs. C.1 and C.2 in Appendix C.  \nFor shower inclinations less than 60 · , our tests demonstrate the same level of agreement as in Figs. 3 and 5. In Figs. D.1 and D.2 in Appendix D, a comparison between 2D MCEq and corsika is further presented for highly inclined showers (80 · ), revealing up to 25% differences in  \n+  \nFIG. 4. Comparison of the first ( ⟨ θ ⟩ , or the distribution mean) and the second ( √ ⟨ θ 2 ⟩ , or the distribution width) moments of the angular distributions of atmospheric leptons at the Earth's surface as computed in corsika ('C') and MCEq ('M') for the same initial conditions as in Fig. 3. The bottom panel shows the difference between the corsika and the MCEq estimates.  \n<!-- image -->  \nthe angular distributions at large angles of deflection from the primary particle axis. While it is possible that the different implementations of muon energy losses or muon propagation geometries are contributing to the observed mismatch, the precise impact of these factors has not been quantified in this study. However, we emphasize that the level of agreement of the angular distributions and spectra in the 1-10 GeV energy region is still satisfactory even for highly inclined showers, and the mild angle-dependent discrepancy observed for single-primary showers will be smeared out by the integration of the secondary particle fluxes across the full sky.  \nFinally, we note that both the 2D MCEq solutions and the corsika Monte Carlo outputs are subject to systematic uncertainties due to the choice of the hadronic interaction model used to describe the particle yields in the hadron-nucleus inelastic collisions. For the energy spectra and angular distributions of the low-energy leptons, the choice of the low-energy hadronic interaction model ( E primary ≤ 150 GeV) has the largest impact. We test two hadronic interaction models internally within 2D MCEq in Figs. E.1 and E.2 in Appendix E, finding up to 20% differences between the models at lepton energies below 10 GeV.", 'VI. SUMMARY AND OUTLOOK': "In this work, we have detailed the development and application of 2D MCEq , an extended version of the 1D MCEq software. The 2D MCEq code provides an efficient numerical approach to angular evolution of hadronic cascades with broad particle physics applications - in particular, in atmospheric lepton flux modelling. This  \ntool considers all crucial aspects of hadronic and leptonic physics, such as inelastic interactions of hadrons with atmospheric nuclei, decays of unstable particles, energy losses, muon polarization, and muon multiple scattering.  \nValidation of 2D MCEq was performed against the standard Monte Carlo code, corsika . The results display agreement within 1-10 % for neutrino angular distributions in air showers up to medium inclinations. Larger differences between the two codes were observed in the distribution tails (corresponding to large angles and high energies).  \nGiven the very high level of agreement with corsika and a significant computational superiority over the Monte Carlo approach, 2D MCEq provides a very appealing option for atmospheric lepton flux calculations. The computational cost of the 2D MCEq calculations is between several CPU-seconds for vertical showers and 1 CPUminute for the near-horizontal showers, compared to multiple CPU-hours to gather sufficient statistics for inclusive flux calculations via the Monte Carlo simulations. Our tool therefore opens the pathway to fast exploration of the systematic uncertainties on the angular distributions of atmospheric leptons, including those associated with the hadronic interaction models and the cosmic ray primary flux. The 2D MCEq code can further be utilized within hybrid air-shower calculation frameworks, such as the integration of corsika and conex [64], with the added feature of explicit angular dependence.  \nFuture enhancements will involve the integration of three-dimensional calculations, accounting for factors such as the Earth's spherical geometry, the initial angular distribution of cosmic ray primaries, the geomagnetic cutoff for these primaries, and the deflection of cascade secondaries within the geomagnetic field.", 'p ( E =100GeV) at 30  inclination': "FIG. 5. Energy spectra of atmospheric leptons in a proton-induced air shower ( E 0 = 100GeV, θ 0 = 30 · ), as computed numerically in 1 MCEq (solid lines) and simulated in the corsika Monte Carlo (filled markers). Here, ' MCEq ' corresponds to the κ = 0 slice of the 2D MCEq solution. The bottom sub-panel in each plot shows the ratio of corsika ('C') to MCEq ('M'). The shaded gray band represents the region where the MCEq solution is numerically unstable (see the caption of Fig. A.3 for details).  \n<!-- image -->", 'ACKNOWLEDGEMENTS': "The authors, T.K. and D.J.K., acknowledge the support from the Carlsberg Foundation (project no. 117238). The author, A.F., is grateful for the supportive environment provided by Prof. Hiroyuki Sagawa's group at the ICRR as a recipient of the JSPS International Research Fellowship (JSPS KAKENHI Grant Number 19F19750).  \nThe authors acknowledge the invaluable computational resources provided by the Academia Sinica Grid-Computing Center (ASGC), which is supported by Academia Sinica.", '1. Event generation chain example': 'TABLE A1. Energy grid settings for the event generation chain example in Fig. 2. The three columns (from left to right) correspond to the left bin edge, the bin center, and the right bin edge of the respective particles on the MCEq kinetic energy grid.', '2. Hadronic model interpolation': "The Sibyll-2.3d and epos-lhc generators are the common choice for the high-energy hadronic interaction model in the air shower codes such as corsika . They are nominally valid down to the primary energies of E thresh = 80GeV [32]. At E primary < E thresh , the UrQMD or DPMJet-III 19.1 models are recommended. In our tests, the transition threshold is lifted to E thresh = 150GeV, where the low-energy UrQMD model is still valid [62, 63]. To transition between the two energy regimes, we run both 'high' and 'low' energy models in chromo and interpolate across the two energy bins adjacent to E thresh (one on each side of E thresh ). In this case, 'linear spline' refers to the order-1 spline interpolation on the logarithmic MCEq energy grid. An example result of the interpolated angular distributions of the secondary pions in the p + 14 N → π + + X ∗ process is shown in Fig. A.1.", 'p (141GeV) + (4GeV)': "1  \nFIG. A.1. Angular density of the secondary pions obtained in the p + 14 N → π + + X ∗ process, as simulated in the UrQMD (orange) and the epos-lhc event generators. Assuming the threshold of 150 GeV in the low-/high-energy model transition, the protons in the MCEq energy bin centered at ∼ 140 GeV fall into the 'intermediate' energy regime, where we use a model linearly interpolated between UrQMD and epos-lhc .  \n<!-- image -->", '3. Muon production and propagation': 'Along with the generic hadronic cascade development, the cascade equations need to account for muon-specific phenomena - namely, muon polarization and muon multiple scattering - for accurate GeV and sub-GeV neutrino flux predictions. The implementation of these phenomena in 2D MCEq is outlined below.', 'a. Muon polarization': "Most atmospheric muons are produced in the two-body decays of π ± and K ± . In the rest frame of the decaying mesons, the muons are completely polarized, i.e. have their momenta perfectly aligned or antialigned with their spin direction. This is a direct consequence of the angular momentum conservation. For example, µ -always has a negative helicity (projection of the muon spin onto its momentum) in the π -/K -rest frame, since ¯ ν µ generated in the same decay must necessarily be right-handed (i.e., have a positive helicity). A similar argument holds for µ + and ν µ , with the helicity assignments flipped. In the lab frame, the µ -helicity in the two-body decay of a meson M -reads [53]:  \nP ( β M , θ ∗ ) = 1 β µ · (1 -r M ) + (1 + r M ) cos θ ∗ β M (1 + r M ) + (1 -r M ) cos θ ∗ β M , (A1)  \nwhere r M = ( m µ /m M ) 2 , β is the Lorentz velocity factor ( β X ≡ v X /c ), and θ ∗ is the angle of the muon emission in the decaying meson rest frame (defined e.g. with respect to the z axis). For each individual muon, the helicity defined in Eq. A1 is a continuous quantity spanning the range between -1..1. To simplify the solution of the cascade equations, where the helicity expectation values of the population of muons rather than the spin states of individual muons are of interest, we switch to the helicity basis , i.e. the basis of purely right-handed muons µ ± R and purely left-handed muons µ ± L . Then, on average, the probability of finding µ -in the right-/left-handed state is  \nP R , L ( β M , θ ∗ ) = 1 2 [1 ±P ( β M , θ ∗ )] , (A2)  \nwhere '+' corresponds to the right-handed state and '-' - to the left-handed state. For µ + , the correspondence between the right-/left-handedness and the +/- sign in Eq. A2 is flipped. In practice, to find the probabilities of the polarized muon production, we let π ± and K ± decay at rest in Pythia and select the two-body decay events. This immediately gives us the muon emission angle θ ∗ and, after boosting to the lab frame, the velocity factors β µ and β M . Then, we compute the probabilities of the right-/left-handedness according to Eq. A2 and assign the respective helicity to each generated muon. We also keep track of the muon lab frame energies and their emission angles with respect to the primary meson boost axis. Finally, we fill in the 'cubes' of the polarized muon yield coefficients in the Hankel frequency space as outlined before in Sec. IV A.  \nTo compute and pre-histogram the neutrino yields in the decays of polarized muons, we use the whizard Monte Carlo code [65, 66], which can take the spin direction of the parent muon as an input and sample from the final three-particle phase space simultaneously. This is in contrast with the analytical expressions for the daughter lepton momenta in the polarized muon decay provided in [53] and [67], which are given separately for each daughter and marginalized over momenta of the other two decay products. Thus, we prefer the whizard Monte Carlo approach for simplicity and efficiency of implementation. For each of µ -and µ + decaying at rest, and each of the two choices of spin configurations (ˆ s µ = ⟨ 0 , 0 , ± 1 ⟩ ), we generate 10,000,000 three-body decay events. We then boost the daughter neutrinos to the muon energies corresponding to the kinetic energy grid of MCEq . This gives us the angular distributions and the energy spectra of neutrinos originating from the decays of the left-helical and the right-helical muon states µ ± L and µ ± R . The superposition of the latter, as prescribed by Eq. A2, gives an accurate representation of the muon population across the entire continuous range of helicities accessible to the polarized muons via Eq. A1. This justifies the inclusion of only six muon species ( µ ± L , µ ± R , as well as unpolarized muons µ ± originating from the three-body decays of K ± ) into the MCEq cascade equations.  \nIn Fig. A.2, we show example angular distributions the electron antineutrinos resulting from the µ -→ ν µ + ¯ ν e + e -decay as computed in whizard . For example, at E µ = 5GeV and E ¯ ν e ≤ 2GeV, we find that both the shape of the neutrino angular distributions and their normalization differs depending on muon polarization. While [22, 54] already included the muon polarization effects in the one-dimensional approximation via the analytical prescription of [53], this simple example further illustrates the importance of the muon polarization treatment for the two-dimensional solver.  \nIn Figs. A.3 and A.4, we show the impact of the muon polarization on the energy spectra and angular distributions of atmospheric neutrinos produced in a full proton-induced air shower. In this example, the primary proton has a fixed energy of 100 GeV and a 30 · inclination. The neutrino fluxes with muon polarization enabled/disabled are obtained via 2D MCEq as the solutions to Eq. 15 at the sea level. This corresponds to X ≈ 1196 g cm -2 in the US Standard  \n1  \nFIG. A.2. Comparison of the ¯ ν e angular distributions with E ¯ ν e ≤ 2GeV in the decay of a 5 GeV muon for three muon polarization cases: left-handed (yellow), right-handed (purple), and unpolarized (red). The displayed events were generated with the whizard Monte Carlo code [65, 66].  \n<!-- image -->  \nFIG. A.3. Impact of muon polarization on the energy spectra of the sea-level atmospheric neutrino fluxes in a 100 GeV proton air shower with 30 · inclination. The dashed lines refer to the case of all atmospheric muons being treated as unpolarized. The solid lines correspond to the muon polarization treatment as described in the main text of the present section. The shaded gray band represents the region where the MCEq solution is unstable due to the numerical implementation of the delta function-like initial condition on the discrete MCEq energy grid (see [68] for the discretization details).  \n<!-- image -->  \natmosphere [32]. We find that muon polarization has the largest impact on the ν e fluxes and energy spectra. Assuming that all muons are unpolarized can lead to 10-30% error in the ν e spectrum normalization and the nearly the same bias in the angular density. The spectrum and the angular distribution of ν µ are affected at the level of a few percent; ¯ ν µ experience an up 10% effect growing towards higher energies due to the energy-dependent π + /π -and µ + /µ -ratios. The muon normalizations and angular distributions remain unchanged as expected.", 'b. Muon multiple scattering': "An additional effect modifying the muon angular distributions is their Coulomb scattering with atmospheric nuclei, e.g. 14 N or 16 O. The effect of multiple scattering is described by the Moli'ere theory [61] if the number of scatters is  \nFIG. A.4. Impact of muon polarization on the angular distributions of the sea-level atmospheric neutrinos with E ν ≃ 2GeV. The initial conditions of the air shower are the same as in Fig. A.3.  \n<!-- image -->  \nlarge ( ≫ 20 across a layer of matter with a given thickness), and by the Poissonian probability of scattering with the Rutherford cross section if the number of scatters is small [32, 69]. We implement a Gaussian approximation to the Moli'ere formalism, which is one of the multiple scattering handling options provided in the corsika Monte Carlo [32]. In the Gaussian approximation, the probability P of a muon deflecting by the space (unprojected) angle θ after traversing ∆ X of the atmospheric slant depth is defined as follows:  \nP ( θ, ∆ X ) = 1 πθ 2 s ∆ X · exp [ -θ 2 ∆ Xθ 2 s ] , (A3)  \nwhere θ 2 s = 1 λ s ( E s E µ, lab β µ ) 2 , E s = 0 . 021GeV, λ s = 37 . 7gcm -2 , E µ, lab is the total muon energy in the lab frame, and β µ - its lab-frame Lorentz velocity factor [32, 69, 70]. In Fig. A.5, we show several representative angular densities computed according to Eq. A3. For illustration, we use ∆ X = 1gcm -2 . In general, however, ∆ X varies with X in response to the longitudinal atmospheric density variations, and the width of the muon multiple scattering kernel is variable. While in just 1 g cm -2 the expected muon deflection is small ( O (0 . 1 · ) at GeV energies), this effect accumulates with the slant depth and results in a noticeable shift of the sea-level muon angular distribution, especially in horizontal showers.  \nTo incorporate this additional convolution kernel into the solution of Eq. 15, we need to Hankel-transform Eq. A3. In the Hankel frequency space, the muon multiple scattering kernel reads:  \n˜ P ( κ, ∆ X ) = exp [ -κ 2 ∆ Xθ 2 s 4 ] , (A4)  \nwhich is scaled so that the overall normalization of the muon angular distribution (represented by the κ = 0 mode) remains unchanged. We can then directly multiply Eq. A4 by the Hankel amplitudes of the muon angular distributions after each integration step ∆ X . This way, the simplified muon multiple scattering model becomes a natural part of the matrix cascade equations. The treatment of multiple scattering is identical for all of the muon species, i.e. µ ± L , µ ± R , and µ ± .  \nIn Fig. A.6, we show the impact of the muon multiple scattering on the sea-level muon angular distributions in a proton-induced air shower, given the same initial conditions as in Figs. A.3 and A.4. For a representative example, we focus on the muons with E µ ≃ 2GeV (i.e., the parents of O (GeV) electron and muon neutrinos). We find that the cumulative effect of muon multiple scattering is a ∼ 1 · shift of the angular density peak, compared to the air shower evolved without muon multiple scattering. We confirmed that the 'tilt' seen in the lower panel of Fig. A.6 grows with the distance traversed by the muons. The effect on the neutrino angular distributions was, however, found to be negligible, introducing at most O (1%) bias at the sea level if muon multiple scattering was not included in the cascade equations.  \n1  \n2  \nFIG. A.5. Gaussian approximation of the muon deflection angles due to their multiple scattering on atmospheric nuclei (see [32, 61, 69] and Eq. A3). The assumed slant depth traversed by the muons is ∆ X = 1gcm -2 . As expected, lower-energy muons (e.g. 1 GeV; yellow line) get deflected more than the higher-energy muons (e.g. red and purple lines for 2 GeV and 5 GeV, respectively).  \n<!-- image -->  \n1e 3  \nFIG. A.6. Impact of muon multiple scattering on the angular distributions of the sea-level atmospheric muons with E µ ≃ 2 GeV. The initial conditions of the air shower are the same as in Fig. A.3.  \n<!-- image -->", 'Appendix B: Comparison of the configuration settings between CORSIKA and MCEq': "The geomagnetic field and the respective curving of the charged particle trajectories are not implemented in 2D MCEq at the time of writing. We therefore effectively disable the geomagnetic field in corsika by setting B x = B z = 10 -5 µ T ( corsika requires | B | > 0). In addition, since the 2D MCEq code currently excludes electromagnetic cascades, we also disable the electromagnetic interactions in corsika by setting all electromagnetic flags (' elmflg ') to false. While the choice of the hadronic models is matched between corsika in MCEq , the switch between the low- and the high-energy regimes means a sharp transition between the two hadronic interaction models in corsika and a smooth interpolation between the models in 2D MCEq (see Appendix A 2 for details). Thus, when both energy regimes are covered in a simulated case, small discrepancies between the two codes are possible due to the different implementations of the low-energy/high-energy model transition. When comparing the simulation outputs at the energies below the transition threshold, one further has to be mindful of the different low-energy model versions. We expect the differences due to the low-/high-energy transition implementation and the internal model versions to be smaller than due to a full change of the low-energy interaction model to a different one, which is investigated in Appendix E.  \nIn corsika , the azimuthal angle of the primary particle incidence is fixed at φ 0 = 0. The height of the first possible interaction of the proton with atmospheric nuclei is set to 112 . 8 km in both MCEq and corsika . The atmospheric density as a function of the slant depth X is modelled according to the Linsley parametrization of the US Standard atmosphere [32]. In MCEq , the average stopping power of the charged particles due to ionization, bremsstrahlung and pair production is taken from tables provided by the Particle Data Group [18, 22], whereas the energy derivative ∂ ∂E is approximated as a five-point stencil. In corsika , the average stopping power is calculated analytically via the Bethe-Bloch prescription [69, 71], and is directly used to reduce the energy of the charged particles between two propagation steps. The Gauss approximation is employed in both codes for muon angular deflections due to multiple scattering.", 'Appendix C: MCEq-CORSIKA benchmarking for a power law cosmic ray spectrum': 'FIG. C.1. Angular distributions of atmospheric leptons in proton-induced air showers with a power law starting spectrum (see title).  \n<!-- image -->', 'p ( np E 2.7 ) at 30  inclination': 'FIG. C.2. Energy spectra of atmospheric leptons in proton-induced air showers with a power law starting spectrum (see title).  \n<!-- image -->  \np ( E =100GeV) at 80  inclination  \nFIG. D.1. Angular distributions of atmospheric leptons in a proton-induced air shower ( E 0 = 100GeV, θ 0 = 80 · ), as computed numerically in 2D MCEq (solid lines) and simulated in the corsika Monte Carlo (filled histograms with errorbars).  \n<!-- image -->', 'p ( E =100GeV) at 80  inclination': "N  \nd  \nFIG. D.2. Energy spectra of atmospheric leptons in a proton-induced air shower ( E 0 = 100GeV, θ 0 = 80 · ), as computed numerically in 1D MCEq (solid lines) and simulated in the corsika Monte Carlo (filled markers).  \n<!-- image -->  \n2  \nFIG. E.1. Angular distributions of low-energy atmospheric leptons ( E ≃ 2 GeV) in a proton-induced air shower, as computed numerically in 2D MCEq using two different low-energy hadronic interaction models ( DPMJet -III 19.1 ('D') and UrQMD 3.4 ('U')). The bottom sub-panel in each plot shows the DPMJet / UrQMD solution ratio.  \n<!-- image -->  \nD / U  \nFIG. E.2. Energy spectra of atmospheric leptons in a proton-induced air shower ( E 0 = 100GeV, θ 0 = 30 · ), as computed numerically in 2D MCEq using two different low-energy hadronic interaction models. The legend follows that of Fig. E.1. The shaded gray band represents the region where the MCEq solution is numerically unstable (see the caption of Fig. A.3 for details).  \n<!-- image -->  \n- [1] M. G. Aartsen et al. (IceCube Collaboration), Determining neutrino oscillation parameters from atmospheric muon neutrino disappearance with three years of IceCube DeepCore data, Phys. Rev. D 91 , 072004 (2015).\n- [2] M. G. Aartsen et al. (IceCube Collaboration), Measurement of atmospheric neutrino oscillations at 6-56 GeV with IceCube DeepCore, Phys. Rev. Lett. 120 , 071801 (2018).\n- [3] K. Abe et al. (Super-Kamiokande Collaboration), Atmospheric neutrino oscillation analysis with external constraints in Super-Kamiokande I-IV, Phys. Rev. D 97 , 072001 (2018).\n- [4] P. Adamson et al. (MINOS Collaboration), Combined analysis of ν µ disappearance and ν µ → ν e appearance in MINOS using accelerator and atmospheric neutrinos, Phys. Rev. Lett. 112 , 191801 (2014).\n- [5] M. G. Aartsen et al. (IceCube Collaboration), Measurement of atmospheric tau neutrino appearance with IceCube DeepCore, Phys. Rev. D 99 , 032007 (2019).\n- [6] Z. Li et al. (Super-Kamiokande Collaboration), Measurement of the tau neutrino cross section in atmospheric neutrino oscillations with Super-Kamiokande, Phys. Rev. D 98 , 052006 (2018).\n- [7] K. Abe et al. (Super-Kamiokande Collaboration), A Measurement of atmospheric neutrino flux consistent with tau neutrino appearance, Phys. Rev. Lett. 97 , 171801 (2006).\n- [8] T. Eberl, S. Hallmann, and J. Hofestadt (KM3NeT Collaboration), Tau neutrino appearance with KM3NeT/ORCA, Proc. Sci., ICRC2017 ( 2018 ) 1025.\n- [9] S. Hallmann, Sensitivity to atmospheric tauneutrino appearance and all-flavour search for neutrinos from the Fermi Bubbles with the deepsea telescopes KM3NeT/ORCA and ANTARES, Ph.D. thesis, Erlangen-Nuremberg U., 2021, https://ecap.nat.fau.de/wp-content/uploads/ 2021/02/Dissertation\\_HallmannSteffen\\_Opus.pdf .\n- [10] T. Stuttard and M. Jensen, Neutrino decoherence from quantum gravitational stochastic perturbations, Phys. Rev. D 102 , 115003 (2020).\n- [11] P. Coloma, J. Lopez-Pavon, I. Martinez-Soler, and H. Nunokawa, Decoherence in neutrino propagation through matter, and bounds from IceCube/DeepCore, Eur. Phys. J. C 78 , 614 (2018).\n- [12] A. Schneider, B. Skrzypek, C. A. Arguelles, and J. M. Conrad, Closing the neutrino BSM gap: Physics potential of atmospheric through-going muons at DUNE, Phys. Rev. D 104 , 092015 (2021).\n- [13] L. Fischer (IceCube Collaboration), Search for heavy neutral lepton production and decay with the IceCube Neutrino Observatory, Proc. Sci., ICHEP2022 ( 2022 ) 190.\n- [14] M. G. Aartsen et al. (IceCube Collaboration), Evidence for high-energy extraterrestrial neutrinos at the IceCube detector, Science 342 , 1242856 (2013).\n- [15] M. G. Aartsen et al. (IceCube Collaboration), Observation and characterization of a cosmic muon neutrino flux from the Northern hemisphere using six years of IceCube data, Astrophys. J. 833 , 3 (2016).\n- [16] R. Abbasi et al. (IceCube Collaboration), First all-flavor search for transient neutrino emission using 3-years of IceCube DeepCore data, JCAP 01 , 027 (2022).  \n- [17] R. Abbasi et al. (IceCube Collaboration), Search for sub-TeV neutrino emission from novae with IceCubeDeepCore, Astrophys. J. 953 , 160 (2023).\n- [18] R. L. Workman (Particle Data Group Collaboration), Review of particle physics, Prog. Theor. Exp. Phys. 2022 , 083C01 (2022).\n- [19] T. K. Gaisser, Spectrum of cosmic-ray nucleons, kaon production, and the atmospheric muon charge ratio, Astropart. Phys. 35 , 801 (2012).\n- [20] H. P. Dembinski, R. Engel, A. Fedynitch, T. Gaisser, F. Riehn, and T. Stanev, Data-driven model of the cosmicray flux and mass composition from 10 GeV to 10 11 GeV, Proc. Sci., ICRC2017 ( 2018 ) 533, [arXiv:1711.11432].\n- [21] T. K. Gaisser, T. Stanev, and S. Tilav, Cosmic ray energy spectrum from measurements of air showers, Front. Phys. (Beijing) 8 , 748 (2013).\n- [22] A. Fedynitch, F. Riehn, R. Engel, T. K. Gaisser, and T. Stanev, Hadronic interaction model Sibyll 2.3c and inclusive lepton fluxes, Phys. Rev. D 100 , 103018 (2019).\n- [23] F. Riehn, R. Engel, A. Fedynitch, T. K. Gaisser, and T. Stanev, Hadronic interaction model Sibyll 2.3d and extensive air showers, Phys. Rev. D 102 , 063002 (2020).\n- [24] T. Pierog, I. Karpenko, J. M. Katzy, E. Yatsenko, and K. Werner, EPOS LHC: Test of collective hadronization with data measured at the CERN Large Hadron Collider, Phys. Rev. C 92 , 034906 (2015).\n- [25] S. A. Bass et al. , Microscopic models for ultrarelativistic heavy ion collisions, Prog. Part. Nucl. Phys. 41 , 255 (1998).\n- [26] S. Roesler, R. Engel, and J. Ranft, The Monte Carlo event generator DPMJET-III, in International Conference on Advanced Monte Carlo for Radiation Physics, Particle Transport Simulation and Applications (MC 2000) (2000), pp. 1033-1038, arXiv:hep-ph/0012252.\n- [27] A. Fedynitch, Cascade equations and hadronic interactions at very high energies, Ph.D. thesis, KIT, Karlsruhe, Dept. Phys., 2015.\n- [28] C. Bierlich et al. , A comprehensive guide to the physics and usage of PYTHIA 8.3 (2022), arXiv:2203.11601.\n- [29] J. M. Picone, A. E. Hedin, D. P. Drob, and A. C. Aikin, NRLMSISE-00 empirical model of the atmosphere: Statistical comparisons and scientific issues, J. Geophys. Res. Space Phys. 107 , 1468 (2002).\n- [30] T. K. Gaisser, R. Engel, and E. Resconi, Cosmic Rays and Particle Physics (Cambridge University Press, 2016).\n- [31] M. Honda, M. Sajjad Athar, T. Kajita, K. Kasahara, and S. Midorikawa, Atmospheric neutrino flux calculation using the NRLMSISE-00 atmospheric model, Phys. Rev. D 92 , 023004 (2015).\n- [32] D. Heck, J. Knapp, J. N. Capdevielle, G. Schatz, and T. Thouw, CORSIKA: A Monte Carlo code to simulate extensive air showers, Technical Report No. FZKA-6019, Forschungszentrum Karlsruhe (1998), https://publikationen.bibliothek.kit.edu/ 270043064/3813660 .\n- [33] P. Lipari, The geometry of atmospheric neutrino production, Astropart. Phys. 14 , 153 (2000).\n- [34] M. Honda, T. Kajita, K. Kasahara, and S. Midorikawa, A new calculation of the atmospheric neutrino flux in a 3-dimensional scheme, Phys. Rev. D 70 , 043008 (2004).\n- [35] G. D. Barr, T. K. Gaisser, P. Lipari, S. Robbins, and\n- T. Stanev, A three-dimensional calculation of atmospheric neutrinos, Phys. Rev. D 70 , 023006 (2004).\n- [36] S. Agostinelli et al. (GEANT4 Collaboration), GEANT4a simulation toolkit, Nucl. Insrum. Methods Phys. Res., Sect. A 506 , 250 (2003).\n- [37] G. Battistoni, S. Muraro, P. R. Sala, F. Cerutti, A. Ferrari, S. Roesler, A. Fasso, and J. Ranft, The FLUKA code: Description and benchmarking, AIP Conf. Proc. 896 , 31 (2007).\n- [38] C. J. Werner, J. S. Bull, C. J. Solomon, F. B. Brown, G. W. McKinney, M. E. Rising, D. A. Dixon, R. L. Martz, H. G. Hughes, L. J. Cox, et al. , MCNP version 6.2 release notes (2018), https://permalink.lanl.gov/object/tr? what=info:lanl-repo/lareport/LA-UR-18-20808 .\n- [39] T. Sato, Y. Iwamoto, S. Hashimoto, T. Ogawa, T. Furuta, S. ichiro Abe, T. Kai, P.-E. Tsai, N. Matsuda, H. Iwase, et al. , Features of particle and heavy ion transport code system (PHITS) version 3.02, J. Nucl. Sci. Technol. 55 , 684 (2018).\n- [40] S. J. Sciutto, AIRES: A system for air shower simulations, arXiv:astro-ph/9911331.\n- [41] M. Honda, T. Kajita, K. Kasahara, S. Midorikawa, and T. Sanuki, Calculation of atmospheric neutrino flux using the interaction model calibrated with atmospheric muon data, Phys. Rev. D 75 , 043006 (2007).\n- [42] M. Honda, T. Kajita, K. Kasahara, and S. Midorikawa, Improvement of low energy atmospheric neutrino flux calculation using the JAM nuclear interaction model, Phys. Rev. D 83 , 123001 (2011).\n- [43] Y. Nara, JAM: an event generator for high energy nuclear collisions, Eur. Phys. J. Web Conf., 208, 11004 (2019).\n- [44] R. Engel, G. Barr, T. K. Gaisser, S. Robbins, and T. Stanev, TARGET 2.2: A hadronic interaction model for studying inclusive muon and neutrino fluxes, in 28th International Cosmic Ray Conference (2003), pp. 16031606.\n- [45] M. Honda, M. Sajjad Athar, T. Kajita, K. Kasahara, and S. Midorikawa, Reduction of the uncertainty in the atmospheric neutrino flux prediction below 1 GeV using accurately measured atmospheric muon flux, Phys. Rev. D 100 , 123022 (2019).\n- [46] K. Sato, H. Menjo, Y. Itow, and M. Honda, Upgrade of Honda atmospheric neutrino flux calculation with implementing recent hadron interaction measurements, Proc. Sci., ICRC2021 ( 2021 ) 1210 (2021).\n- [47] A. Abusleme et al. (JUNO Collaboration), JUNO sensitivity to low energy atmospheric neutrino spectra, Eur. Phys. J. C 81 , 10 (2021).\n- [48] S. Fukasawa, M. Ghosh, and O. Yasuda, Complementarity between Hyperkamiokande and DUNE in determining neutrino oscillation parameters, Nucl. Phys. B 918 , 337 (2017).\n- [49] G. D. Barr, T. K. Gaisser, S. Robbins, and T. Stanev, Uncertainties in atmospheric neutrino fluxes, Phys. Rev. D 74 , 094009 (2006).\n- [50] A. Fedynitch and M. Huber, Data-driven hadronic interaction model for atmospheric lepton flux calculations, Phys. Rev. D 106 , 083018 (2022).\n- [51] V. A. Naumov and T. S. Sinegovskaya, Simple method for solving transport equations describing the propagation of cosmic ray nucleons in the atmosphere, Phys. At. Nucl. 63 , 1927 (2000).  \n- [52] A. A. Kochanov, T. S. Sinegovskaya, and S. I. Sinegovsky, bibinfo title High-energy cosmic ray fluxes in the Earth atmosphere: calculations vs experiments, Astropart. Phys. 30 , 219 (2008).\n- [53] P. Lipari, Lepton spectra in the Earth's atmosphere, Astropart. Phys. 1 , 195 (1993).\n- [54] A. Fedynitch, R. Engel, T. K. Gaisser, F. Riehn, and T. Stanev, Calculation of conventional and prompt lepton fluxes at very high energy, EPJ Web Conf. 99 , 08001 (2015).\n- [55] N. Baddour, Operational and convolution properties of two-dimensional fourier transforms in polar coordinates, J. Opt. Soc. Am. A 26 , 1767 (2009).\n- [56] T. Kozynets, A. Fedynitch, and D. J. Koskinen, A numerical approach to angular distributions in hadronic cascades, Proc. Sci. ICRC2021 ( 2021 ) 1209.\n- [57] G. V. Haines and A. G. Jones, Logarithmic Fourier transformation, Geophys. J. Int. 92 , 171 (1988).\n- [58] A. J. S. Hamilton, Uncorrelated modes of the nonlinear power spectrum, Mon. Not. R. Astron. Soc. 312 , 257 (2000).\n- [59] J. Lang and B. Frank, Fast logarithmic fourier-laplace transform of nonintegrable functions, Phys. Rev. E 100 , 053302 (2019).\n- [60] chromo : Cosmic ray and HadROnic interactiOn MOnte-carlo frontend (2022), https://github.com/ impy-project/chromo .\n- [61] H. A. Bethe, Moliere's theory of multiple scattering, Phys. Rev. 89 , 1256 (1953).\n- [62] J. Steinheimer, V. Dexheimer, H. Petersen, M. Bleicher, S. Schramm, and H. Stoecker, Hydrodynamics with a chiral hadronic equation of state including quark degrees of freedom, Phys. Rev. C 81 , 044913 (2010).\n- [63] UrQMD User Guide, Frankfurt Institute for Advanced Studies (2014), https://vfs.fias.science/d/ cd8fa2d938/files/?p=/urqmd-3.4.pdf .\n- [64] T. Bergmann, R. Engel, D. Heck, N. N. Kalmykov, S. Ostapchenko, T. Pierog, T. Thouw, and K. Werner, One-dimensional hybrid approach to extensive air shower simulation, Astropart. Phys. 26 , 420 (2007).\n- [65] W. Kilian, T. Ohl, and J. Reuter, WHIZARD: Simulating multi-particle processes at LHC and ILC, Eur. Phys. J. C 71 , 1742 (2011).\n- [66] M. Moretti, T. Ohl, and J. Reuter, O'Mega: An optimizing matrix element generator, arXiv:hep-ph/0102195.\n- [67] S. Mine, Systematic measurement of the spinpolarization of the cosmic-ray muons, Ph.D. thesis, Institute for Nuclear Study, University of Tokyo, 1996, https://repository.dl.itc.u-tokyo.ac. jp/record/2005339/files/A12180.pdf .\n- [68] J. Heinze, A. Fedynitch, D. Boncioli, and W. Winter, A new view on Auger data and cosmogenic neutrinos in light of different nuclear disintegration and air-shower models, Astrophys. J. 873 , 88 (2019).\n- [69] D. E. Groom and S. R. Klein, Passage of particles through matter, Eur. Phys. J. C 15 , 163 (2000).\n- [70] B. B. Rossi, High-Energy Particles , Prentice-Hall Physics Series (Prentice-Hall, Englewood Cliffs, N.J, 1952).\n- [71] E. Segr'e, H. Staub, H. A. Bethe, and J. Ashkin, Experimental Nuclear Physics. Volume I (John Wiley & Sons, Chapman & Hall, New York, London, 1953)."}