Universality parameterization extension of electron-positron Angluar spectra in extensive air showers

ii May, 2020 Copyright ©2020 All Rights Reserved Abstract The angular spectrum of extensive electromagnetic atmospheric cascades initiated by cosmic rays may be parameterized by only particle number and momentum angle with respect to the shower axis as shown by Lafebre et al. [1] This work extends the parameterizations given by Lafebre et al. [1] in secondary energy, momentum angle theta, and normalized particle number n values. The parametrization limit in secondary electron/positron energy has been raised by a factor of 2.5 and the normalized n limits have been lowered two orders of magnitude for a range of secondary energies. The spectrum fits have been ex- tended in theta for secondary energies up to 840. Another point of investigation was the dependence of the angular spectrum on shower stage, where no clear e↵ects were found aside from phenomenological restrictions. While only very high energy primary cosmic ray energies were used in Lafebre et al. [1] (1018 eV), this work probed primary energies down to 1015 eV and found parame- terizations to model the data with this new degree of freedom. No significant deviations were found between the electron and positron distributions. 1 1 Introducton Cosmic rays are particles with energies of order 106 MeV and above that originate from extrasolar sources. Their velocities are generally a large fraction of the speed of light, primarily being composed of alpha particles or protons with heavier nuclei being more rare. The sources of these particles are still under investigation, but it is believed that they mostly eminate from supernovae and active galactic nuclei. The peak of the energy distribution for these particles lies at 0.3 GeV. When a cosmic ray collides with the atmosphere, it initiates an extensive air shower. An important characteristic of these showers is that some aspects may be modeled using data that is not unique to an individual cascade - a concept called shower Universality. Papers such as Nerling et al. [2] and Lafebre et al. [1] are based on shower universality, each giving parameterizations that aid in shower reconstruction which only depend on shower stage, primary and secondary energies, normalized particle counts, and particle momentum angle to shower axis. Physicists use these parameterizations to work backwards from their data to better understand the cascade event. There are a few large cosmic ray observatories: the Cosmic Ray Telescope Array based in Utah and the Pierre Auger observatory located in Argentina. Recent soft- ware developments have allowed better shower simulation which aids in reconstruction analysis. The data for this study was generated by CORSIKA (COsmic Ray SImu- lations for KAscade), developed at Karlsruhe Institute of Technology. This program simulates extensive air showers and gives measurements of the cascade particles as they cross specified observation levels. The data in this study was formatted in rows, each to represent a particle which passed through an observation plane including its species, momentum vector components, arrival time, and weight. It should be noted that weight here does not 1 mean the force of gravity on the particles but rather a number that loosely represents how many particles are assigned to that entry. 2 1.1 Angular spectrum The coordinate system used in CORSIKA is shown in figure 1 The azimuthal angle ¢ represents the rotation about a vertical axis. The elevation angle ✓ represents the rotation upward from the vertical axis. Observation levels were chosen at regular grammage intervals - that is to Figure 1: Coordinate system used in CORSIKA simulations [3] say that a particle will traverse an equal mass of atmo- sphere between each level but the physical separation will not be constant. This is be- cause the atmosphere does not have constant density as a function of height. The longitudinal cascade description used in this work is called relative evolution stage, or t. This is a parameter that gives a more universal description of each shower as it is defined in terms of the shower max, or the point at which the shower reaches its maximum particle number. At t = 0, the shower is at its maximum. At negative values of t, the shower has not yet reached its maximum and for positive t values the shower has progressed beyond its maximum. Relative evolution stage, t is defined as: t⌘ X - Xmax Xo (1) X represents the slant depth - the amount of matter the cascade has travelled through. Xmax represents the depth in the atmosphere of the shower max and Xo is the radi- ation length of electronscmin 2 air. These values all have units of g . 1.1 Angular spectrum This work focuses on the extending the electron-positron angular 1.1 Angular spectrum spectrum of ex- tensive air showers. The angular spectrum represents the distributional relationship 3 4 1.1 Angular spectrum between cascade particle (secondary particle) counts, energies, and momentum eleva- tion angle with respect to the shower axis. To make discussion easier, the following notation will be used: primary cosmic ray energy will be denoted by ⇧, secondary electron and positron energy will be denoted as ✏, momentum angle will be ✓, and normalized particle number will be n. We define the normalization condition as in Lafebre et al. [1] in terms of the generic variables µ and ⌫: Z ⌫max ⌫min n(t; µ, ⌫)d⌫ =1 (2) This condition ensures that there exists a common scale to compare the data and functional forms. In the case of the angular spectrum we consider µ to be ln ✏ and ⌫ to be ⌦, the solid angle that encompasses the desired particles. n(t; ln(✏), ⌦)= n(t; ln(✏), ✓) sin(✓) (3) It is more convenient to write the function in terms of ✓ in order for it to be integrable. Lafebre et al. [1] defines the functional form of the angular spectrum as follows: ⇥ ⇤n(t; ln(✏), ⌦)= Co (e 1 ✓↵1 )-1/ + (e 2 ✓↵2 )-1/ (4) This equation represents two joined power laws with a connection factor er, repre- senting the “smoothness” of the transition. The ↵1 and ↵2 factors control the slopes of the power laws and the /31 and /32 factors control the relative heights. Co is the normalization constant. According to Lafebre et al. [1] these parameters are only functions of ✏ and are given by: 5 Angular spectrum parameter functions ↵ 1 (✏)= -0.399 ↵ 2 (✏)= -8.35 + 0.440 ln(✏) 0.210 /3 1 (✏)= -3.73 + 0.92✏ 2 Method 2.1 Data /3 2 (✏)= 32.9 - 4.84 ln(✏) We define the five dimensional histogram M⇧t✓✏N with bins in ⇧, t, ✓, ✏, and N representing primary energy, stage, momentum angle, secondary energy, and run number respectively. The subscripts represent the index along which each variable is varied or more precisely, the index value represents the corresponding bin in that variable. If an operation is performed on the object, it is according to the last index. The entries in this object are sum of the weights that CORSIKA has assigned to each particle entry. Because computation time scales with primary energy, CORSIKA may take in a thinning fraction "th =⇧ ✏ such that if the following condition is met for i secondary particles: "th⇧ > X ✏i (5) i Then only kth secondary particle will be followed and attributed the weight: ⇥ ✏k ⇤-1 Mk = P i ✏i (6) For N individual simulations, we take the following sum: X N M⇧t✏✓N = M⇧t✏✓ (7) 2.2 6 Errors We then take the weight per ✓, where the widths of the bins in ✓ are defined as x✓, to be: = ⇥⇧t✏✓ M⇧t✏✓ x (8) ✓ Because we wish to satisfy the normalization condition over ✓ we define the normal- ization constant C as the following: ✓ ◆ 180 X ⇥⇧t✏ x = ⇡ C⇧t✏ ✓ ✓ (9) ✓ We now take the outer product of the matrix of normalization constants and the matrix of the sine of the bin centers in theta x✓ in order to compensate for the solid angle ⌦: C⇧t✏ ⌦ sin(x ✓ )= C⇧t✏✓ (10) We can finally define n to be the normalized four-dimensional histogram that will be fit as: = ⇧t✏✓ n⇧t✓✏ ⇥C (11) ⇧t✏✓ 2.2 Errors The process of particle detection follows Poisson statistics. In the case that the weight of each event is 1, the variance for N events is naturally N. In our case, the weights of the events are not 1. This necessitates the use of a weighted poisson distribution, where variance is given by: er2 = N hw2i. N here is the number of unweighted events and therefore the variance reduces to the sum of the weights squared. The errors are then: sX 2 er⇧t✓✏ = M ⇧t✓✏N (12) N 7 Which is normalized in the same manner as n⇧t✏✓: ⌃⇧t✓✏ (13) er⇧t✓ ✏ =C ⇧t✏✓ p Indeed, this method proved to be much better than strictly using M⇧t✓✏ in the sense that the ⌫x2 was closer to unity for each fit, but we introduce a new factor to compensate for inconsistent error size: R⇧t✓✏ = A⇧t✓✏ µ (14) ⇧t✓✏ Where A⇧t✓✏ and µ⇧t✓✏ represent the variance and the mean respectively, for entries across N. We then define the final errors to be: S⇧t✏✓ = R⇧t✏✓ ⇥ ⌃⇧t✓✏ 3 Analysis and Parameterization 3.1 Histogram ranges (15) One thousand simulations were performed. Ten sets of one hundred air showers were binned, each set of one hundred having the same primary energy, ⇧. Histogram data bins and ranges Data Index limit Bin ranges Spacing scheme ⇧ (eV) 10 1015 to 1018 Logarithmic N (event) 100 0 to 100 Linear t (stage) 10 -5 to 5 Linear ✏ (MeV) 44 100 to 104.4 Logarithmic 8 ✓ (Degree) 80 0 to 90 Logarithmic 3.2 3.2 9 Data selection Data selection Some ranges of data collected from each shower will be unique to that event. As such, this data needs to be cut from the fitting process such that the domain is restricted to capture only the data that can be modeled universally. Immediately it is understood that the hadronic component of the shower is not universal and can be seen in the angular spectrum data as a second peak for ✏ > 101.9MeV . The following function in ✏ gives the upper bound in ✓ for 101.9MeV < ✏ < 104.4MeV : ✓ (✏)= 5.0 · 103 ✏-1 (16) The non-universal portion of data that this function cut o↵ did not move with respect to stage. As such, it appeared to be a function of ✏ alone. Another cut for n < 10-6 was made as statistically, there aren’t enough particles in the shower to justify fitting below it. Figure 2 shows a typical normalized complete angular spectrum vs Figure 2: Raw normalized angular spectrum vs. cut normalized angular spectrum to be fit a spectrum with the applied cuto↵s. To encompass all of the limits on the domains of fitting, Figure 3 shows the ranges for ✓ and n vs. ✏ where they remain constant in ⇧. 10 3.3 Initial fitting Figure 3: Upper and lower bounds for ✓ and n with respect to ✏ 3.3 Initial fitting For each primary particle energy ⇧, there are ten three dimensional surfaces in n, over ✓ and ✏, for constant t. Using4eq. as contour lines across ✓ at each bin center x✏ we used a SciPy fitting routine to find the parameters for each distribution. Then we fit the parameter data as a function of ✏. According to Lafebre et al, these parameters depend only on ✏ but in this study, we seek to confirm t independence and investigate ⇧ dependence. The first results showed di↵erent functional forms for each parameter. Changes in the auxiliary parameters of the angular spectrum parameters are meant to capture the ⇧ and t dependence, if there exists any. In that sense one may think of the a, b, c and d parameters in each ↵ and /3 as possible functions of ⇧ and t giving the motivation for including them in the functional statements. ⌫ The x2 values remained between 10-1 and 101 for all ⇧ and t. Parameter forms ↵1(⇧, t, ✏ )= 0.61 + a ln (✏)+ b ln (✏+1.5) ↵2(⇧, t, ✏ )= a ln (0.00005✏)2 + b ln (✏)+ c 11 3.3 Initial fitting /31(⇧, t, ✏ )= a ln (✏ + b)2 + c ln (d✏) /32(⇧, t, ✏ )= a ln ✏b + c 3.4 12 Errors Figure 4: Data and fits for ⇧ = 3.16 ⇥ 1018 eV at t = 5.0 and t = 0.0 3.4 Errors There was behavior of note in the errors as a function of ✏. S⇧t✏✓ rises and falls as a function of ✏ for a particular ✓ as a result of the variance, A⇧t✏✓ over N. This causes the x2⌫ values to rise from ⇠1 at the high and low ✏ to larger values for median ✏ values but remain on the order of 10. 3.5 Stage dependence Stage dependence was investigated in the following manner. For each value of ⇧, one may imagine ten three dimensional matrices, or rather ten angular spectrum surfaces. Each of these ten surfaces represents the angular spectrum at a particular stage. We then consider for each value of ✏ ten contour lines of di↵ering stage value over n and ✓ . If the auxiliary parameters in the angular spectrum parameterizations are plotted against t we may see their behavior as a function of t. Fitting this data linearly in t we find that the slopes all remain very close to zero. This means we are unable to make any definitive statement on stage dependence except that there appears to be none to within the errors we have. This was an attempt to quantify the fact that the variations between each contour line in t for each ✏ appear to be negligible. With this in 3.4 Errors mind, we define a new M matrix and normalize it in the same manner as section 13 3.6 14 Primary energy dependence 2: XX N M⇧t✏✓N = M⇧✏✓ (17) t This aids in fitting as there are more entries in each bin - no longer discriminating based on t values. The tail behaviors on the distributions are now much better behaved. 3.6 Primary energy dependence Once the stage independence was confirmed, the work turned to ⇧ dependence. With the stages combined into one contour line per ✏ for each angular spectrum, it was easier to show the behavior of the parameters with respect to ⇧. We plot the angular spectrum parameters against ✏ again, but this time there are ten lines varied in ⇧ with combined t, instead of t for a constant ⇧. There is marked change in the functional form of the angular spectrum parameters with respect to ⇧. To capture this, we give the auxiliary parameters as functions of ⇧. All of the auxiliary variables seemed to follow either linear or quadratic forms in log space. The parameterizations of the auxiliary variables in ↵1, ↵2, /31, and /32 are outlined in the following sections. 3.6.1 ↵1 auxiliary parameterizations a(⇧)= 7.55 ln (⇧)2 · 10-6 - 4.0 ln (⇧) · 10-4 + 6.0 · 10-3 (18) b(⇧)= 3.24 ln (⇧)2 · 10-4 - 1.03 ln (⇧) · 10-2 - 5.07 · 10-1 (19) The x2⌫ values are 0.57 and 1.97 respectively. 3.6.2 ↵2 auxiliary parameterizations a(⇧)= -2.62 ln (⇧)2 · 10-5 - 9.28 ln (⇧) · 10-4 + 4.68 · 10-3 (20) 3.6 15 Primary energy dependence b(⇧)= 5.36 ln (⇧)2 · 10-4 - 5.31 ln (⇧) · 10-2 + 9.64 · 10-1 (21) c(⇧)= 1.18 ln (⇧)2 · 10-3 + 1.42 ln (⇧) · 10-1 - 9.98 (22) The x2⌫ values are 3.34, 4.09, and 4.43 respectively. 3.6.3 /31 auxiliary parameterizations a(⇧)= 2.27 ln (⇧)2 · 10-4 - 1.16 ln (⇧) · 10-2 + 1.85 · 10-1 (23) b(⇧)= -31.4 ln (⇧)+ 904.3 (24) c(⇧)= 7.85 ln (⇧) · 10-5 + 1.09 (25) d(⇧)= 3.50 ln (⇧) · 10-5 - 3.9 · 10-4 (26) The x2⌫ values are 0.41, 1.35, 0.39, and 1.16 respectively. 3.6.4 /32 auxiliary parameterizations a(⇧)= -3.46 ln (⇧)2 · 10-2 + 1.81 ln (⇧) - 31.0 (27) b(⇧)= -9.25 ln (⇧)2 · 10-4 + 5.46 ln (⇧) · 10-2 + 2.21 · 10-2 (28) c(⇧)= 3.87 ln (⇧)2 · 10-2 - 2.05 ln (⇧)+ 6.30 · 101 The x2⌫ values are 0.04, 0.04, and 0.07 respectively. (29) 3.6 Primary energy dependence 3.6.5 Primary dependence plots Figure 5: Angular spectrum parameters versus ✏. ⇧ values decrease from yellow to purple. The ↵2 and /31 lines are shifted upwards by factors of 1.0, ↵1 is shifted by factors of 0.1, and /32 is shifted upward by consecutive factors of 2.0 Figure 6: ↵1 auxiliary parameters and fit lines 16 3.6 Primary energy dependence Figure 7: ↵2 auxiliary parameters and fit lines Figure 8: /31 auxiliary parameters and fit lines 17 18 Figure 9: /32 auxiliary parameters and fit lines 4 Final fitting Using the auxiliary parameter functions within the angular spectrum equations, we reconstruct the fit portion of the angular spectrum in ⇧, ✏, and ✓. Figure 10: Simulation data versus complete reconstruction for highest and lowest ⇧ values The spectrum for lower ⇧ is more difficult to reconstruct. A possible solution to 19 this is to increase the number of simulations and find di↵erent functional forms for the auxiliary parameters that better capture the ⇧ dependence. values for ⌫ The x2 the final reconstruction are quite large, but decrease as ⇧ increases from order 105 to 102. 5 Conclusion In this work, the parameterizations for the angular spectrum in extensive air showers were extended upward in ✏ by a factor of 2.5. The n range has been lowered to 10-6 for values of ✏ between 101 and 102 MeV. New functional forms of the angular spectrum parameters have been determined. Stage independence has been verified by fitting the angular spectrum parameters versus t to a slope⌫ of zero and finding that within the errors, the fits are valid i.e. the x2 values are ⇠1. ⇧ dependence has been shown in the data and auxiliary parameterizations that capture this behavior have been given. References [1] Lafebre, S., et al. Universality of Electron-Positron Distributions in Extensive Air Showers. Astroparticle Physics, vol. 31, no. 3, 2009, DOI: 10.1016/j.astropartphys.2009.02.002 [2] Nerling, F., et al. Universality of Electron Distributions in High-Energy Air Showers- Description of Cherenkov Light Production. Astroparticle Physics, vol. 24, no. 6, 2006. DOI:10.1016/j.astropartphys.2005.09.002 [3] D. Heck, J. Knapp, J.N. Capdevielle, G. Schatz, T. Thouw CORSIKA: A Monte Carlo Code to Simulate Extensive Air Showers Institute for Nuclear Physics Forschungszentrum und Universitä t Karlsruhe, Karlsruhe College de 20 France, Paris