Lookup tables to compute high energy cosmic ray induced atmospheric ionization and changes in atmospheric chemistry
LLookup tables to compute high energy cosmic rayinduced atmospheric ionization and changes inatmospheric chemistry.
Dimitra Atri , Adrian L. Melott and Brian C. Thomas . University of Kansas, Department of Physics and Astronomy, 1251 Wescoe Dr. Washburn University, Department of Physics and Astronomy, Topeka, KS 66621,United States of AmericaE-mail: [email protected]
Abstract.
A variety of events such as gamma-ray bursts and supernovae may exposethe Earth to an increased flux of high-energy cosmic rays, with potentially importanteffects on the biosphere. Existing atmospheric chemistry software does not have thecapability of incorporating the effects of substantial cosmic ray flux above 10
GeV . Anatmospheric code, the NASA-Goddard Space Flight Center two-dimensional (latitude,altitude) time-dependent atmospheric model (NGSFC), is used to study atmosphericchemistry changes. Using CORSIKA, we have created tables that can be used tocompute high energy cosmic ray (10
GeV - 1
P eV ) induced atmospheric ionizationand also, with the use of the NGSFC code, can be used to simulate the resultingatmospheric chemistry changes. We discuss the tables, their uses, weaknesses, andstrengths.
Keywords : high energy cosmic rays, cosmic ray theory a r X i v : . [ a s t r o - ph ] M a y igh energy cosmic ray induced atmospheric ionization
1. Introduction
Nearby supernovae [1], gamma ray bursts [2, 3] and possibly galactic shocks [4] maybathe the Earth in cosmic rays (CRs) of much higher than usual incident energies. It isof considerable interest to investigate the effect of such events on the Earth’s atmosphere,and consequent possible connections to mass extinctions and other events in the fossilrecord. Until now, cosmic rays of this high energy have been included only by meansof a simple phenomenological approximation, enhancing the existing background CRionization, or not at all. Computations of the effects of gamma-ray bursts [2] have beenbased on only the effects of photons. However, a sufficiently strong background of veryhigh-energy cosmic rays can punch through the galactic and terrestrial magnetic fields[3] and irradiate the Earth, with effects potentially competitive with those of photons.Other potential scenarios which may account for a 62 My periodicity in biodiversity [4]are based solely on CR, and to date the terrestrial effects have only been approximated.Thus there is considerable value in developing software to model the effects of a spectrumof CR with energies above those normally included in atmospheric computations.We have developed a method to model changes in atmospheric chemistry whenhigh energy cosmic rays (HECRs) ionize the atmosphere. When HECRs hit theatmosphere, like other CRs they will interact with atmospheric constituents, primarilythe molecules of N and O . This interaction will either result in a nuclear reaction orelectromagnetic interaction ionizing the atmosphere, the latter being the interaction ofprimary importance for chemistry [5]. To study this interaction we created a normalizedHECR ionization spectrum data table that contains ionization energy due to HECRs atvarious energies and altitudes. We will discuss the methods used, how the data tablewas generated, and its use in a widely used time-dependent atmospheric ionizationand chemistry code. This result naturally complements and extends basic work on thelower-energy CRs that normally dominate atmospheric ionization [6].In order to obtain ionization energy from HECRs, we used CORSIKA (COsmic RaySImulations for KAscade), a high energy cosmic ray extensive air shower simulator [7,8]. The code follows the interactions of a primary cosmic ray and its secondary particlesthrough the atmosphere to the ground. CORSIKA uses Monte Carlo calculations toaccount for high energy strong and electromagnetic interactions using a number ofextensive air shower simulators (of which we used UrQMD-low energy and EPOS-highenergy) [9].Using CORSIKA we created tables so that atmospheric ionization can be calculatedfor an arbitrary CR spectrum between 10 GeV and 1
P eV . We intend to continuouslyupdate this table to higher energies, eventually reaching the highest energies observed.
2. Methods
The tables we have generated give energy deposition binned in two different ways, andcan be used to either simply compute atmospheric ionization or can be used with the igh energy cosmic ray induced atmospheric ionization
Figure 1.
Fractional energy deposition for a 10 TeV primary at zenith angles 5 o (solid), 45 o (dotted), and 85 o (dashed) in bins of the NGSFC code in log pressure,proportional to total column density traversed by the shower, which is approximatelylinear in altitude. CRs of 10 − eV are of much larger primary energy than those that dominatenormal galactic CRs, so one should not simply turn up the usual background, as ina previous supernova study [11]. We used CORSIKA, which is designed to performdetailed simulations of extensive air showers initiated by high energy cosmic ray igh energy cosmic ray induced atmospheric ionization intervals of primary energy between 10 GeV and 1
P eV , i.e. at 50 different primaryenergies.We created lookup tables using data from CORSIKA runs that containsatmospheric ionization energy deposition per primary due to CRs in the range of 10
GeV - 1
P eV . The sum of deposition over altitude is less than the total of the primaryenergy, as not all energy is deposited in the the atmosphere by electromagnetic processes.Nuclear interactions also occur between HECRs and atmospheric particles, but nuclearenergy is dumped into nuclei, mostly not into atmospheric chemistry. Energy that goesinto nuclear interactions or reaches the ground is not included in the deposition table.An arbitrary spectrum can be convolved with this data and the results used in theNGSFC code (or other similar codes) to simulate the effects this energy flux depositionwill have on atmospheric chemistry. We investigated the effect of zenith angle by running100 shower ensembles (1500 primary particle simulation runs) at zenith angles of 5 o , 45 o ,and 85 o at 10 eV . In Figure 1 we show the fractional energy deposition for each of thesezenith angles (excluding nuclear interactions) per interval (intervals of the NGSFC code)in log pressure, proportional to total column density traversed by the shower, which isapproximately linear in altitude. Note that the lateral displacement in the lines, and thelocation of their maxima, are reasonably approximated by a ( cos θ ) − factor, confirmingthe simplest thing one would expect from a column density factor. For a general-purposecode at the energies we consider, it is appropriate to assume the flux is isotropic. Foreach shower, we recorded the fractional energy deposition in each of 1000 bins of 1 gcm − column density. We used their mean to construct a lookup table describing theenergy deposition for a particle of given primary energy as a function of pressure. Belowwe describe in detail the construction and use of the table.The greatest deposition of energy per bin corresponds to the first interactionbetween the incoming CR and the atmosphere, which occurs very high in theatmosphere. Since log column density is nearly linear with altitude, we analyze themaximum energy deposition per bin of log column density. This has weak trends withprimary energy, as shown in Figure 2. This gives about 400 g cm − as the site of thegreatest deposition of energy per unit distance, corresponding to an altitude of about5 km . This can be compared with (a) about 13 km as the mean altitude of maximalenergy deposition density for the normal CR spectrum as implemented in the NGSFCcode, strongly biased toward latitudes greater than above about 60 degrees, and (b) 22to 35 km as the peak deposition for keV − M eV photons, depending upon energy [2,15].The table was originally created to compute ozone ( O ) depletion and otheratmospheric chemistry changes [5] resulting from energy deposition by HECRs. O lives at altitudes of 10-35 km [16] with considerable latitude dependence [15]. Ourdata is inaccurate from 46-90 km because CORSIKA runs on a linear column densityscale starting with 1 g cm − , but main effects of atmospheric chemistry changes onthe biosphere occur at lower altitudes. A 46-90 km altitude is less significant becausehigh-energy CRs rarely have their first interaction that high in the atmosphere. igh energy cosmic ray induced atmospheric ionization Figure 2.
Electromagnetic energy deposition from a 100
GeV primary (solid), a 10
T eV primary (dotted) and an
P eV primary (dashed) per g cm − , in bins of NGSFCcode . Energy deposition for the higher energy primaries is deeper in the atmosphere.Energy going into nuclear interactions or hitting the ground is not included in thistable
3. Construction of the Lookup tables
The CORSIKA model 6.900 is used in all simulations. When installing CORSIKA,UrQMD was used for the low energy hadronic interaction model and EPOS as thehigh energy hadronic interaction model [22]. CORSIKA was installed with CURVED,SLANT and UPWARD options to simulate primaries at large zenith angles and trackupward going particles in a curved atmosphere. The following are the importantvariables used in the input file for CORSIKA. The Longitudinal Shower Developmentvariable (LONGI) is set as LONGI T 1 T F giving longitudinal ionization for every 1 g cm − bin to facilitate interpolation to the bins of GSFC photochemical code. Thevariable THETAP gives the angle of incidence of particles and was set as THETAP z z,where z is the zenith angle. The code ran for zenith angles 5 o , 15 o ...upto 85 o , a total of9 angles. The energy range variable (ERANGE) is set as ERANGE x x, where x is thevariable energy of the primary particle. For our preliminary runs (10 GeV - 100
T eV ),CORSIKA was run 100 times at each of the energies stated above using different SEED igh energy cosmic ray induced atmospheric ionization g cm − is the only data we are interested in order to evaluate the atmosphericionization. In the future we plan to extend this table above the P eV range.CORSIKA ran at 51 different energy levels 100 times each at intervals of 0.1 in logenergy for 9 angles. Each time it was run CORSIKA simulated 15 protons entering theatmosphere. Therefore, we compiled 13500 particles for each of the 51 energy levels.The energy was averaged for each of the 900 runs at each longitudinal pressure bin.Minor problems arose when transferring the CORSIKA data to the NGSFC code.First, the smallest atmospheric depth bin at which energy deposition is given is 1 g cm − in CORSIKA, while 22 of the 46 pressure levels are below 1 g cm − in the NGSFC code.Secondly, CORSIKA outputs deposition in a linear pressure scale. We interpolated toa logarithmic scale for the NGSFC code. We converted CORSIKA atmospheric depthunits of g cm − to millibars (needed for NGSFC) by a conversion factor of 0.98.To resolve the problem of CORSIKA’s highest interaction altitude being below 22of the NGSFC codes altitudes we linearly interpolated the data using the ionizationfrom the 1 g cm − bin. Because the first interaction point of primaries is rarely withinthe first g cm − , most of the CORSIKA output files have no energy deposition for thisaltitude. Also, the data we used from CORSIKA is binned in 0.1 GeV intervals, meaningfor the lower energy runs many of the 1036 altitude levels are zero. For this reason ouroutput file may have no energy in the 22 highest altitude bins for certain energy levels.For our purpose of looking at the ozone these higher altitudes, above ∼ km , are notimportant, since the amount of ionization in them from HECR would be quite small.We use the same linear interpolation method used for each bin with less than 1 g cm − , depositing the energy from ground level to the highest logarithmic altitudelevel in the last bin. The logarithmic altitude bins are centered on the bins used in theNGSFC code (Table 2).The geomagnetic field should be mentioned with regard to the lookup table. Itsimpact on CRs is minor for protons with energies greater than 17 GeV [17]. Becausethis only impacts a small region of the CR range of importance (10
GeV - 1
EeV ) effectsdue to the magnetic field are small. If desired, alterations to the latitude distributionof energy deposition can be made in the 10 to 17
GeV range to simulate the effects ofdeflection by the geomagnetic field.
4. Comparison with existing simulations
We have compared our data with existing data generated by CORSIKA [6] as well asanalytical methods proposed by Velinov et al. (2008) [21].We found very good agreement with Usoskin et al. data at 1
T eV (figure 3), 100
GeV (figure 4) and 10
GeV (figure 5). However, there is about 5% higher ionizationin our data at lower altitudes. This can be attributed to the UrQMD model we haveused, which is more efficient in tracking muons compared to FLUKA used in the earlier igh energy cosmic ray induced atmospheric ionization
P eV (figure 6) with Tanguy Pierog’s(CORSIKA developer) simulation generated by the same hadronic interaction modelsand found very good agreement. This simulation provides the energy deposition profileof a 1 PeV primary using the EPOS-UrQMD hadronic interaction models producedwith CORSIKA. Fluctuations in the energy deposition profile were removed using theCONEX code.The analytical method [21] can be used for a thin target, which is above 50 kmaltitude. Also, the equation used to calculate ionization limits primary energy range of600
M eV T eV . We applied the analytical solution to the energy range 17
GeV T eV in order to ignore the effects of solar modulation. Our values of ionization are about 5%higher than predicted by the approximation (Table 3). As discussed earlier, CORSIKAdoes not allow altitude bins smaller than 1 g cm − , which contains the entire atmosphereabove 50 km, and therefore, it is not possible to resolve different levels within 1 g cm − bin. Atmospheric ionization produced by the primaries in our energy range is a very tinyfraction of the total ionization from the GCR flux, and therefore we can not compareour results with certain models [20], but it will be possible to compare for a given CRspectrum. igh energy cosmic ray induced atmospheric ionization Figure 3.
Comparison of our data (dash) at 1 TeV with Usoskin et al. (solid) [6]. igh energy cosmic ray induced atmospheric ionization Figure 4.
Comparison of our data (solid) at 100 GeV with Usoskin et al. (dash) [6] igh energy cosmic ray induced atmospheric ionization Figure 5.
Comparison of our data (dash) at 10 GeV with Usoskin et al. (solid) [6] igh energy cosmic ray induced atmospheric ionization Figure 6.
Comparison of our data (dots) at 1 PeV with T. Pierog’s data (solid)produced with CORSIKA/CONEX in bins of 10 g cm −
5. How to use the lookup tables
The lookup tables are formatted into 51 columns corresponding to 51 energies but thenumber of rows for both tables are different. For the atmospheric chemistry model,the table has 46 rows corresponding to the altitude bins of the NGSFC code shown intable 2. As is, it displays ionization energy deposition for a spectrum of 1 particle persteradian per square meter per second for every 0.1 logarithmic energy bin, in unitsof
GeV /m .sr.s . Trivially, one must multiply by the number of particles per unit areaper steradian per second in their bin to get the total energy flux at each altitude binfor a given spectral form. To input this data into the NGSFC code the energy at eachaltitude must be added over all energy levels, creating a 1-dimensional data set of totalenergy deposition flux at each altitude.The general purpose table has 51 columns as mentioned above and 21 rows of binsize 50 g cm − each. The table displays the ionization yield function in units of ionpairs per meter square per steradian per second. Total number of ion pairs for a givenspectrum can be obtained by the method mentioned above. igh energy cosmic ray induced atmospheric ionization o − o , o − o , o − o etc). For an isotropic flux, the sameflux is entered for each latitude bin. The final data file to be input into the NGSFCcode should now be a 2-dimensional set of ionization energy flux deposition at the 18latitude bins for 46 altitudes, a total of 828 data points. For a point source, the inputinto latitude bins may be adjusted by the appropriate factor, including a correction forthe cos θ factor as in [2]. This may result in hemispheric differences in the results.The energy flux data as a function of altitude and latitude, generated as describedabove, may be used in the NGSFC code by way of a simple read-in subroutine.Depending on the units of the spectrum used, conversion to cm − s − may be necessarysince the NGSFC code uses cgs units. In order to use the input as a source of N O y , theenergy deposition rate must first be converted to ionization rate. This is accomplishedusing 35 eV per ion pair [18], which finally gives values in units of ions cm − s − .Constituents in the code are stored with units of number cm − s − . Therefore, thearea ionization rate is converted to a volume rate using the height of the altitude bins,which depends on the current density of each bin at read in (the density depends ontemperature, which depends on the presence of sunlight, etc.). This ionization rate isthen used as a source for N O y , assuming 1.25 molecules are created for each ion pair[16]. The model then runs as usual, incorporating this source of N O y in the relevantchemistry computations. The general procedure here is the same as that used in previouswork with both photon and solar proton ionization sources [2, 10].It must be noted that alpha particles were not considered to produce these tables.At such high energies, an alpha particle can be approximated by four protons with thesame energy per nucleon and the resulting ionization can be computed using the datagenerated with protons.
6. Discussion
The inclusion of CR effects will make possible a more accurate investigation of theeffects of gamma-ray bursts on the Earth, which have the potential to explain certainmass extinction events [19]. A quantitative treatment of possible time-varying flux ofhigh-energy CRs becomes possible [5]. It may have other applications which we cannotanticipate, and should be made generally available.This lookup table will be made freely available via ftp, and upgraded in the future asappropriate. Look for a link at http : //kusmos.phsx.ku.edu/ melott/Astrobiology.htm upon publication of this paper. igh energy cosmic ray induced atmospheric ionization Acknowledgments
We thank Tanguy Pierog and Dieter Heck for their help and advice in using CORSIKA.This research was supported in part by the National Science Foundation throughTeraGrid resources provided by the National Center for Supercomputing Applicationsand by NASA grant NNX09AM85G. B.T. acknowledges a Small Research Grant fromWashburn University.
References [1] B.D. Fields, T. Athanassiadou, and S.R. Johnson,
Supernova Collisions with the Heliosphere , Astrophys. J. , 549 (2008).[2] B.C. Thomas, A.L. Melott, C.H. Jackman, C.M. Laird, M.V. Medvedev, R.S. Stolarski, N. Gehrels,J.K. Cannizzo, D.P. Hogan, and L.M. Ejzak,
Gamma-Ray Bursts and the Earth: Exploration ofAtmospheric, Biological, Climatic, and Biogeochemical Effects , Astrophys. J. , 509 (2005).[3] C.D. Dermer and J.M. Holmes,
Cosmic Rays from Gamma-Ray Bursts in the Galaxy , Astrophys.J. , L21 (2005).[4] M.M. Medvedev and A.L. Melott,
Do Extragalactic Cosmic Rays Induce Cycles in FossilDiversity? , Astrophys. J . , 879 (2007).[5] A.L. Melott, A.J. Krejci, B.C. Thomas, M.V. Medvedev, G.W. Wilson, and M.J. Murray, Atmospheric consequences of cosmic-ray variability in the extragalactic shock model , J. Geophys.Res. , E10007, doi:10.1029/2008JE003206 (2008).[6] I.G. Usoskin and G.A. Kovaltsov,
Cosmic ray induced ionization in the atmosphere: Full modelingand practical applications , J. Geophys. Res. , D21206, doi:10.1029/2006JD007150 (2006).[7] T. Djemil, R. Attallah, and J. N Capdevielle,
Simulation of the Atmospheric Muon Flux withCorsika , Int. J. Mod. Phys.
A 20 .[9] T. Pierog, R. Engel, D. Heck, S. Ostapchenko,and K. Werner, Prepared for 30th InternationalCosmic Ray Conference,
Latest Results from the Air Shower Simulation Programs CORSIKAand CONEX , (ICRC 2007) , Merida, Yucatan, Mexico, arXiv:0802.1262v1 (2007).[10] C.H. Jackman and R.D. McPeters, The Effect of Solar Proton Events on Ozone and OtherConstituents , Geophys. Monogr. Ser. edited by J.M. Pap and P. Fox, , 305 (2004).[11] N. Gehrels, C. M. Laird, C. H. Jackman, J. K. Cannizzo, B. J. Mattson, C. Wan,
Ozone depletionfrom nearby supernovae , Astrophys J. , 1169 (2003).[12] C.H. Jackman, M.T. DeLand, G.J. Labow, E.L. Fleming, D.K. Weisenstein, M.K.W. Ko, M.Sinnhuber, and J.M. Russell,
Neutral atmospheric influences of the solar proton events inOctoberNovember 2003 , J. Geophys. Res. , A09S27, doi: 10.0129/2004JA010888 (2005).[13] B.C. Thomas, C.H. Jackman, and A.L. Melott,
Modeling atmospheric effects of the September1859 solar flare , Geophys. Res. Lett. , L06810, doi:10.1029/2006GL029174, (2007)[14] C.H. Jackman, E.L. Fleming, S. Chandra, D.B. Considine, J.E. Rosenfield, Past, present, andfuture modeled ozone trends with comparisons to observed trends , J. Geophys. Res. , 28,753(1996).[15] L.M. Ejzak, A.L. Melott, M.V. Medvedev, and B.C. Thomas,
Terrestrial Consequences of Spectraland Temporal Variability in Ionizing Photon Events , Astrophys. J. , 373 (2007).[16] M.B. Harfoot, D. J. Beerling, B. H. Lomax, and J. A. Pyle,
A two-dimensional atmosphericchemistry modeling investigation of Earth’s Phanerozoic O and near-surface ultravioletradiation history , J. Geophys. Res.
D07308 (2007), doi:10.1029/2006JD007372.[17] http : − cutof f − rigidity.php [18] H.S. Porter, C.H. Jackman, A.E.S. Green, Efficiencies for production of atomic nitrogen and oxygenby relativistic proton impact in air , J. Chem. Phys. , 154 (1976). igh energy cosmic ray induced atmospheric ionization [19] A.L. Melott and B.C. Thomas, Late Ordovician geographic patterns of extinction comparedwith simulations of astrophysical ionizing radiation damage , Paleobiology , in press (2009),arXiv:0809.0899[20] F.M. Vitt and C.H. Jackman,
A comparison of sources of odd nitrogen production from 1974through 193 in the Earth’s middle atmosphere as calculated using a two-dimensional model , J.Geophys. Res.
D3, 6729 (1996).[21] P.I.Y. Velinov and L. Mateev,
Analytical approach to cosmic ray ionization by nuclei with chargeZ in the middle atmosphere: Distribution of galactic CR effects , ASR igh energy cosmic ray induced atmospheric ionization igh energy cosmic ray induced atmospheric ionization Table 1.
Input File for CORSIKA - Complete Description of Variables at http : − ik.f zk.de/corsikausersguide/corsika t ech.html igh energy cosmic ray induced atmospheric ionization Table 2.