University of Groningen
Reconstructing air shower parameters with LOFAR using event specific GDAS atmosphere
Mitra, P.; Bonardi, A.; Corstanje, A.; Buitink, S.; Krampah, G. K.; Falcke, H.; Hare, B. M.;
Horandel, J. R.; Huege, T.; Mulrey, K.
Published in:
Astroparticle Physics
DOI:
10.1016/j.astropartphys.2020.102470
IMPORTANT NOTE: You are advised to consult the publisher's version (publisher's PDF) if you wish to cite from
it. Please check the document version below.
Document Version
Publisher's PDF, also known as Version of record
Publication date:
2020
Link to publication in University of Groningen/UMCG research database
Citation for published version (APA):
Mitra, P., Bonardi, A., Corstanje, A., Buitink, S., Krampah, G. K., Falcke, H., Hare, B. M., Horandel, J. R.,
Huege, T., Mulrey, K., Nelles, A., Rachen, J. P., Rossetto, L., Scholten, O., ter Veen, S., Trinh, T. N. G.,
Winchen, T., & Pandya, H. (2020). Reconstructing air shower parameters with LOFAR using event specific
GDAS atmosphere. Astroparticle Physics, 123, [102470].
https://doi.org/10.1016/j.astropartphys.2020.102470
Copyright
Other than for strictly personal use, it is not permitted to download or to forward/distribute the text or part of it without the consent of the author(s) and/or copyright holder(s), unless the work is under an open content license (like Creative Commons).
Take-down policy
If you believe that this document breaches copyright please contact us providing details, and we will remove access to the work immediately and investigate your claim.
Downloaded from the University of Groningen/UMCG research database (Pure): http://www.rug.nl/research/portal. For technical reasons the number of authors shown on this cover page is limited to 10 maximum.
Astroparticle Physics 123 (2020) 102470
ContentslistsavailableatScienceDirect
Astroparticle
Physics
journalhomepage: www.elsevier.com/locate/astropartphys
Reconstructing
air
shower
parameters
with
LOFAR
using
event
specific
GDAS
atmosphere
P.
Mitra
a,∗,
A.
Bonardi
a,b,
A.
Corstanje
b,
S.
Buitink
a,b,
G.K.
Krampah
a,
H.
Falcke
b,c,d,k,
B.M.
Hare
e,
J.R.
Hörandel
a,b,c,
T.
Huege
h,a,
K.
Mulrey
a,
A.
Nelles
f,i,
J.P.
Rachen
a,
L.
Rossetto
b,
O.
Scholten
e,g,
S.
ter
Veen
d,
T.N.G.
Trinh
e,j,
T.
Winchen
a,k,
H.
Pandya
aa Astrophysical Institute, Vrije Universiteit Brussel, Pleinlaan 2, 1050 Brussels, Belgium
b Department of Astrophysics / IMAPP, Radboud University Nijmegen, P. O. Box 9010, 6500 GL, Nijmegen, The Netherlands c NIKHEF, Science Park Amsterdam, 1098 XG Amsterdam, The Netherlands
d Netherlands Institute of Radio Astronomy (ASTRON), Postbus 2, 7990 AA Dwingeloo, The Netherlands e KVI-CART, University Groningen, P. O. Box 72, 9700 AB Groningen, The Netherlands
f DESY, Platanenallee 6, 15738 Zeuthen, Germany
g Interuniversity Institute for High-Energy, Vrije Universiteit Brussel, Pleinlaan 2, 1050 Brussels, Belgium h Institut für Kernphysik, Karlsruhe Institute of Technology(KIT), P. O. Box 3640, 76021, Karlsruhe, Germany i Institut für Physik, Humboldt-Universität zu Berlin, 12489 Berlin, Germany
j Department of Physics, School of Education, Can Tho University Campus II, 3/2 Street, Ninh Kieu District, Can Tho City, Vietnam k Max-Planck Institute for Radio Astronomy, Bonn, Germany
a
r
t
i
c
l
e
i
n
f
o
Article history: Received 18 March 2019 Revised 28 March 2020 Accepted 28 May 2020 Available online 1 June 2020 Keywords:
LOFAR Cosmic Ray EAS
Radio detection technique Atmosphere GDAS Index of refraction Effects of humidity X max reconstruction
a
b
s
t
r
a
c
t
Thelimitedknowledgeofatmosphericparameterslikehumidity,pressure,temperature, andtheindex ofrefractionhasbeenoneoftheimportantsystematicuncertaintiesinreconstructingthedepthofthe showermaximum fromtheradioemissionofairshowers.Currentairshower MonteCarlosimulation codeslikeCORSIKAandtheradioplug-inCoREASusevariousaveragedparameterizedatmospheres. How-ever,time-dependentandlocation-specific atmosphericmodelsareneededforthe cosmicrayanalysis methodusedforLOFARdata.There,dedicatedsimulationsetsareusedforeachdetectedcosmicray,to takeintoaccounttheactualatmosphericconditionsatthetimeofthemeasurement.UsingtheGlobal DataAssimilationSystem(GDAS),aglobal atmosphericmodel,wehaveimplementedtime-dependent, realisticatmosphericprofilesinCORSIKAandCoREAS.Wehaveproducedrealisticevent-specific atmo-spheresforallairshowersmeasuredwithLOFAR,aneventsetspanningseveralyearsandmanydifferent weatherconditions.Acompletere-analysisofourdatasetshowsthatforthemajorityofdata,our pre-viouscorrectionfactorperformedratherwell;wefoundonlyasmallsystematicshiftof2g/cm2 inthe
reconstructedXmax.However,underextremeweatherconditions,forexample,verylowairpressure,the
shiftcanbeupto15g/cm2.Weprovideacorrectionformulatodetermine theshiftinX
max resulting
fromacomparisonofsimulationsdoneusingtheUS-StdatmosphereandtheGDAS-basedatmosphere. © 2020ElsevierB.V.Allrightsreserved.
1. Introduction
Inrecentyears, thefield ofradiodetectionofairshowershas advanced quite rapidly [1,2]. Estimating thedepth of the shower maximum, Xmax, with improved accuracy is of great interest for
the studyof theprimary particlecomposition [3,4].The develop-mentoftheairshowerinducedbyacosmicrayisgovernedbythe interactions anddecaysofthesecondary particles.The secondary electrons andpositronsin theairshower undergocharge
separa-∗ Corresponding author.
E-mail address: pmitra@vub.be (P. Mitra).
tion asthey travel through the magnetic field ofthe Earth. This leadstoa time-varyingtransversecurrent, producing radio emis-sion.Thereisanothersmallcontributiontotheradiationfromthe excessofnegativechargeaccumulatedattheshowerfront,known asthe “Askaryaneffect” [5]. The emission reachesthe groundas ashortpulseon theorderof10 to100 nswithaspecific lateral intensitydistribution,orfootprint, thatdependson Xmax;Xmax is
calculated interms of total atmospheric matter traversed by the airshowerfromthetopoftheatmospheretothepointwherethe particle numberreaches the maximum. It is therefore important toknow thealtitude-dependent airdensity. Anotheratmospheric parameterthatplaysacrucialroleintheradioemissionisthe
re-https://doi.org/10.1016/j.astropartphys.2020.102470 0927-6505/© 2020 Elsevier B.V. All rights reserved.
fractiveindexofair.Ifforagivenemissionregionalongtheshower axisan observerislocatedatthecorrespondingCherenkovangle, radiationemittedfromallalongthisregionarrivessimultaneously. Thisresultsinahighlycompressedsignalintime,forminga ring-likestructureontheground [6,7].Therefractiveindexdetermines thepropagation velocity of the radio signal at differentaltitudes andinfluences the time compression [8,9]. Forobservers located ontheCherenkovring,pulsesarecoherentuptoGHzfrequencies
[10].TheangleatwhichCherenkovemissionisemittedisinversely proportional to therefractive index. At higherfrequencies pulses are more sensitive to the refractive index. In general, at all fre-quencies,thevariationsintherefractive indexlead tochanges in theradio intensityfootprint [11].Boththedensityandthe refrac-tiveindexofairaredependent onairtemperature, humidityand pressure.Thus,havingagoodunderstandingoftheseatmospheric variablesiscrucial.
Theradiodetectiontechniquecanbeusedincombinationwith establishedtechniquessuch asfluorescencedetectionandsurface detectionwithscintillatorsandwaterCherenkovdetectors.Dense antennaarrayslikethecoreoftheLOFARradiotelescope [12] pro-videtheopportunity toinvestigatetheradiofootprint,i.e.the lat-eralintensitydistribution,inclosedetailandenablethe measure-mentofXmax up toa precisionof < 20g/cm2.The precision is
sensitive to the choice of an atmospheric model included inthe Monte Carlo air shower simulation codes. There are several pa-rameterizedatmosphericmodels incorporatedinthe CORSIKAair showersimulationcode,basedonaveraged profiles:U.S.standard atmosphereparameterizedaccordingtoJ.Linsley [13], parameter-izedatmospheresforthePierreAuger ObservatorynearMalargüe (Argentina) by M. Will and B. Keilhauer [14], South Pole atmo-spheresparameterized by P. Lipari andD. Chirkin etc. So far, the USstandardatmospherehasbeenusedinLOFARanalyses,through CORSIKAsimulations [13]andtheCoREASextension [13]whichis usedtocalculatetheradioemissionoftheairshowers.
A first order linear correction to the US standard atmosphere hasbeenappliedtoaccount forthefact thatthe US-standard at-mospheredoesnot reflectthe realistic atmosphericconditions at agiventime.Itispreferabletointegratearealisticatmosphere di-rectlyintothesimulations.Inparticular,thereconstructionofXmax
dependsontherefractiveindexofair,andsoarealisticrefractive indexprofileneedstobeincluded.
Theeffectsoftherefractiveindex,n,onthereconstructedXmax
havebeen previously reportedin Ref.[15] andRef.[11],using dif-ferentsimulationcodes.In Ref.[11],CoREAS wasused tosimulate two ensembles of showers,one with a globally higher refractiv-ityN=
(
n− 1)
106, another withstandard values. AMonte Carlobasedapproachwastakentostudythesystematicshiftin recon-structedXmaxbycomparingthesetofsimulationswithhigher
re-fractivitytothestandardones.TheshiftinthereconstructedXmax
fromthe defaultvalue wasfound tobe proportional tothe geo-metricdistance to Xmax.The effect wasstronger in thehigh
fre-quency band of 120–250 MHz than in the 30–80 MHz band. In Ref.[15],amorerealisticprofileoftherefractivitywasconstructed for one particular day using information from the Global Data Assimilation System, GDAS, a global weather database. The dif-ferencesbetweenthis atmosphereanddefaultatmospheres were studiedusingtheSELFASradioemissionsimulationcode [16].The resultsshowedthat correctingfortherealisticdensityisthemost important factor in the accurate reconstruction of Xmax, causing
about30g/cm2 biasinX
max.Andthesecondmostimportant
cor-rection was through the inclusion of the high frequency refrac-tivityformula, applicable at radio frequencies, contributingabout 5g/cm2 bias in X
max.The effects ofthe increasedrefractivity on
the time traces and the lateral distribution function (LDF) were alsoreported.Inthe 20–80MHzfrequencyband, relatively small differences in the amplitude of the electric field and LDF were
found, whereas considerabledifferenceswere found studying the high frequency band between 120–250 MHz. These results were inagreement withRef.[11].While both workspaved thewayfor the understanding of atmosphericeffects on radio simulations, a directapplicationto realdatausing simulations withrealistic at-mosphericconditionswasnotaddressed.
In thiswork, for thefirst time, GDAS-based atmospheric pro-files, automaticallyincludedinCoREAS simulationsareapplied to LOFAR data.The effects of atmospheric parameters like pressure and humidity on the reconstructed Xmax are studied and
com-paredto the results ofpreviously used linear corrections. Anew GDAS-based correction is introduced and compared to previous methods.Furthermore,atool isdevelopedtoextract GDAS atmo-spheric parameters whichare then interfacedwithCORSIKA. The utility ofthistool is not onlylimitedto LOFAR. Thiscode,called ”gdastool”, has been available for public use since the releaseof CORSIKAversion7.6300.It isflexibleandreadyto beadapted by the users to obtain parameterized atmospheric profiles for user-specifiedtimeandlocation. Sections2and3describethe process-ing ofGDAS data to extract the atmospheric state variables and examples of atmospheric profiles at the LOFAR site, respectively.
Section4coversthedetailsoftheimplementationofGDASin COR-SIKA. In sections 5 and 6, LOFAR cosmic ray data are evaluated withtheGDASatmosphericprofiles,theGDAS-correctionfactoris introducedandtheexpliciteffectsofhumidityonshower param-etersarediscussed.
2. ExtractingatmosphericvariablesfromGDASdata
The Global Data Assimilation System (GDAS) developed at NOAA’s1 NationalCenters forEnvironmental Prediction (NCEP) is
atoolusedtodescribetheglobalatmosphere.Itisrunfourtimes aday(0,6,12,and18UTC)andprovidesa3-,6-and9-hour fore-cast based on the interpolation of meteorological measurements fromallover theworld includingweatherstations onland,ships and airplanes as well as radiosondes and weather satellites [17]. Thethreehourlydataareavailableat23constantpressurelevels, from1000hPa(roughlysealevel)to20hPa(≈ 26km)onaglobal 1◦ spacedlatitude-longitude grid(180◦ by 360◦).Eachdata setis complementedbydataatthesurfacelevel.Thedataarestoredin weekly filesandmadeavailable online.In ordertomodel a real-isticatmosphereoneneedstoobtainthesuitableatmospheric pa-rametersfromGDAS. Parameterslike temperature(K),height (m) relativehumidity(H)andpressure(hPa)canbe directlyextracted fromthedatabase. IntheGDASdata,the altitudeisin geopoten-tialunitswithrespectto ageoid (meansealevel). Thisisan ad-justment to geometric height or elevation above mean sea level usingthevariationofgravity withlatitudeandelevation.To con-vertfromgeopotentialheighth(m)tostandardgeometricaltitude z(m)weusetheformula
z
(
h,)
=(
1+0.002644· cos(
2)
)
· h +(
1+0.0089· cos(
2))
h2 6245000 (1)where
isthegeometriclatitude [18].Tocalculatetheairdensity, therelativehumidityistobeconvertedintowatervaporpressure. The following approximation ofthe empirical Magnus formulais used tocalculate thewatervapor pressure (hPa)in termsof hu-midityandtemperature [18]:
e=100%H × 6.1064× exp
26521.88t .5 + tfor t≤ 0◦C
P. Mitra, A. Bonardi and A. Corstanje et al. / Astroparticle Physics 123 (2020) 102470 3 and e= H 100%× 6.1070× exp
17.15t 234.9 + t for t≥ 0◦C. (2)Thedensitycanbecalculatedfromtheidealgaslawas
ρ
=PMairRT (3)
whereP isthe atmosphericpressurein Pa,Tistemperaturein K andRistheuniversalgasconstant,havingavalueof8.31451JK−1 mol−1 andMairisthemolarmassofair.Moistaircanbe
decom-posedintothree componentstocalculateitsmolarmass: dryair, watervapor andcarbon dioxide.Themolar massofhumid airis the sumof the molarmassesof the components,weighted with thevolumepercentage
φ
iofthatcomponent [18],Mair=Mdry·
φ
dry+Mwater·φ
water+MCO2·φ
CO2. (4) The molar masses of dry air, water vapor and CO2 are 0.02897,0.04401 and 0.01802kg-mol−1 respectively. The volume percent-ageofCO2istakenas385ppmv,thepercentageofwater
φ
wateristhepartialpressureofwatervapor dividedbythepressureP;the dryairmakesuptherest.
Therefractivity,definedasN=
(
n− 1)
106,isafunctionofhu-midity,pressureandtemperaturecanbeexpressedas
N=77.6890KhPa−1pd T +71.2952KhPa −1pw T +375463K 2hPa−1pw T2 (5)
with pw, pd and T being the partial water vapor pressure
(
pw=e× 100Pa)
,partialdryairpressureandtemperaturerespec-tively [19]. The effect of humidity is important for our study as it tends to increase the refractivityin comparison to that of dry airattheradiofrequencies.Thereare differencesbetweenthe re-fractivities obtainedinradio andtheonesinthevisible,nearthe infraredandUVrangesasdescribedin [18].Toaccountforthe un-certainties in GDAS data one needs to perform insitu measure-mentswithweatherballoons.Sincethisisbeyondthescopeofthis work andwe referto [18],whichprovides acomparisonbetween GDASdataandweatherballoonmeasurementsinArgentina.Since globalatmosphericmodelsaretypicallymorepreciseinthe North-ern hemispherewheremoreweatherdataisavailablewe assume that the intrinsic uncertainty ofGDAS atthe LOFAR siteis simi-lartothatinArgentina.Variousrelevantuncertaintiesare: ± 0.5
◦C for temperature, 0.5 hPafor pressure, and0.05hPa for water
vapor pressure andlessthan 1g/cm2 inatmosphericdepth over
thealtituderangefrom3to6km.Theuncertaintyinwatervapor pressuretranslatesto2− 7%uncertaintyinhumidity.Theresulting relativeuncertaintyinNduetotheseparametersisaround0.5%at the same altitude range.The GDAS data have a resolution of 1◦ by 1◦ inlatitude longitude.Thiscan be roughly approximated as adistanceof100kmbetweentwoadjacentgridpoints.Forhighly inclinedshowersthedistancetotheregionofshowerdevelopment fromtheobservationsitecanbelargerthanthedistancebetween two gridpoints.Forairshowerscomingfrom70◦ zeniththis dis-tanceisaround70kmandforzenith > 75◦ itisabout100km. Inthesecases,thechoiceofanexactgrid pointbecomes compli-cated.Alsoatthispoint,forzenith angles > 70◦ thecorrection dueto curvedatmosphere becomesimportant.This doesnot oc-curforLOFARasthedetectedcosmicrays arelimitedtowithina <55◦zenithangleduetotheparticledetectorsusedfortriggering. InthisregimetheGDASmodelworkswell.
3. GDASatmosphericprofilesattheLOFARsite
In thissection severalGDASatmospheric profiles extractedat the LOFAR site are discussed. Fig. 1 (left) shows humidity as a
Table 1
Shift in X max for different zenith, energy and X max bins for dif-
ferent frequency bands.
Frequency band Zenith X max (g/cm 2 )
50–350 MHz low < 30 ◦ -6.24 ± 0.30
50–350 MHz high > 30 ◦ -6.19 ± 0.37
30–80 MHz low < 30 ◦ 0.10 ± 0.50
30–80 MHz high > 30 ◦ -0.05 ± 0.46
Frequency band True X max (g/cm 2 ) X max (g/cm 2 )
50–350 MHz low < 624 -6.78 ± 0.41 50–350 MHz high > 624 -6.30 ± 0.32 30–80 MHz low < 624 -0.61 ± 0.51 30–80 MHz high > 624 0.51 ± 0.46 Frequency band Energy(GeV) X max (g/cm 2 )
50–350 MHz low < 2.18 × 10 8 -6.86 ± 0.35
50–350 MHz high > 2.18 × 10 8 -6.92 ± 0.38
30–80 MHz low < 2.18 × 10 8 -0.48 ± 0.48
30–80 MHz high > 2.18 × 10 8 0. ± 0.49
functionofaltitudefor5arbitraryatmosphericprofiles for differ-entdaysintheyear2011,betweenJuneandNovember.A signifi-cantday-to-dayfluctuationisseen.Theredsolidandbluedashed linesindicatetwo verydifferentweatherconditions;theredsolid linehaving high saturating humiditybetween 5− 8km suggests highercloudcoverage andthebluedashed linewithlow humid-ityinthatrangeindicateslowcloudcoverage. Fig.1(right)shows thedifference inatmosphericdepthprofilebetweenthe US stan-dardatmosphereandthe GDASatmospheresatLOFAR for8 pro-filesovertheyears2011− 2016.TheGDASatmospheresvary sig-nificantlyfromtheUSatmosphere.Atmosphericprofileswith sim-ilaratmospheric depthatgroundcan evolve differentlyhigherin theatmosphere. This isimportantforcalculating thecorrect dis-tance to the shower maximum. Fig. 1 shows the mean profile forthe relative difference inrefractivity
Nrelative between GDAS
andtheUSstandardatmosphereasafunctionofaltitudeforover 3 years for 100 cosmic rays recorded at LOFAR. It is defined as
Nrelative=
(
NGDAS− NUS)
/NUS,whereNGDAS iscalculatedfromEq-5usingGDASatmospheresatLOFAR.NUSisobtainedfromthe
lin-ear relation NUS=ρρus
sealevelNsealevel, withNsealevel=292.This is the
defaultoptionforcalculatingrefractivityinCoREASaswell. The absolutevalue ofthe mean
Nrelative is around 10% near groundandaround3− 8%between3to10kmofaltitude,the re-gionimportantforshowerdevelopment.
Approximately 75% ofthe atmosphericmatter and99% of the total massof water vapor andaerosols are containedwithin the troposphere, the lowest layer of Earth’s atmosphere. Within the tropospherethetemperaturedrops withaltitude,reachinga con-stantvalueinthetropopause,theboundaryregionbetween tropo-sphereandstratosphere.In theU.S standardatmospherethe tro-posphereendsat11kmandtropopauseextendsto analtitudeof 20km.ForthelocalGDASatmospherestheseboundariesarenot sharplydefined. Theflatpartinthemean
Nrelative > 10 km in
Fig.2is the result ofconstant temperature in the tropopause. How-evercontribution fromthisregion tothe radio emission is mini-mal.Toconsidertheeffectsofrefractiveindexinthepropagation timeofradiosignalitisimportanttocalculatetheeffectiveN [1,8]. Thisisdefinedas
Neff=
N
(
h)
dh DwhereDisthedistancebetweenthelineofemissionandobserver. The valuesof relative effectiverefractivity
Nrelatie f f ve between the GDAS and US standard atmosphere are around 7− 10 % in the range of altitude mentioned above, for observers within < 100 moftheshoweraxis.
Fig. 1. Atmospheric profiles at LOFAR. Left : Example of 5 humidity profiles between June to November during the year 2011. Right : 8 profiles for the difference in atmo- spheric depth between US standard atmosphere and GDAS atmospheres as a function of altitude between the years 2011 − 2016 .
Fig. 2. Mean relative refractivity, defined as N relative = NGDASNUS−N ; profiles for 100 US
recorded cosmic rays at LOFAR spanning over the years 2011 to 2014. The black solid line denotes the mean profile and the blue dashed lines show the standard deviations.
4. ImplementationinCORSIKA/CoREAS
To incorporate the atmospheric parameters extracted from GDAS in CORSIKA and CoREAS we have developed a program named”gdastool” that downloads the requiredGDAS file given thetimeandlocation ofobservationoftheeventandreturns re-fractiveindicesbetweengroundandthehighestGDASlevel.Italso fits the density profile according to the standard 5 layer atmo-spheric model used in CORSIKA [13]. In this model the density
ρ
(h)hasanexponentialdependenceonthealtitudeleadingtothe functionalformofmassoverburdenT(h)whichisthedensity inte-gratedoverheight(km)asT
(
h)
=ai+bie−105h/c
i i=1,...,4. (6)
Thus,thedensityis
ρ
(
h)
=bi/cie−105h/c
i i=1,...,4. (7)
Inthe fifth layer the overburdenis assumedto decrease linearly withheight.The parameters ai,bi andci are obtainedina
man-nersuch that thefunction T(h) iscontinuousatthelayer bound-aries and can be differentiated continuously.The first three lay-ersconstituteof the24densitypointsobtainedfrom GDASdata. Thefirst layerconsistsof10 points,second layer of7pointsand
thethird layerof7points. SinceGDASprovidesdataon constant pressurelevels,not ofconstant heights,thelayerboundariesvary slightly betweendifferent atmospheric profiles. The mean values oftheboundariesfortheconditionsof100cosmic rayeventsare 3.56± 0.11km,9.09± 0.23km,26.27± 0.56kmfromboundary 1to3,respectively.
Next,wefitthedatatoEq- 7inthefollowingway:
For layer 1the density profile is fitted withtwo free param-eters. Then thedensity
ρ
1 atboundary 1is calculated usingEq-7withtheobtainedparametersb1,c1.Theconditionthatthe
den-sity hasto be continuousat the boundariesreducesthe number of free parameters to 1 which is the parameter c. Thus the pa-rameterb2 forsecond layercan beexpressed asafunctionof
ρ
1and c2 with c2 being the only free parameter. The same fitting
procedure isrepeatedforthe third layer.The fourthlayer ranges fromthehighestGDASaltitudeto100km.Atthesealtitudesthere arenophysicalGDASdata.Theparameterc4 isobtainedbyfitting
thelast GDASpointandthedensityat100kmfromUSstandard atmosphere. At thesealtitudes the mass overburdenis less than 0.1%ofthevalue atground.Theimportantfactoristosatisfy the boundaryconditionsthroughouttheatmosphere.Alongwith den-sitythecontinuityofmassoverburdenisalsopreserved.Forthat, onceasmoothprofileforthedensityisobtained,theparametera inEq- 6is solved foranalytically, using theboundary conditions forthe massoverburden. The parameterizationforthe fifth layer was adapted from the US standard atmosphere [13]. The ”gdas-tool” alsoreturnsa densityprofileplotwiththebestfit param-etersasafunction ofaltitudeandthe rmsoftherelativedensity differencebetweendataandfit.The relativedensityisdefinedas
ρfit−ρdata
ρfit . Fig.3(left)anditsrmsisusedasagoodnessoffit. Fig.2
(left) showsthe example of a density profile betweenthe fitted modelandGDAS. Themeanrelativeerrorindensityfor100 pro-filesasafunctionofaltitudeispresentedin Fig.3(right).Atlower altitudesthemodelfitsthedataverywell;deviations > 2%start at altitudes higher than 15 km which are not so important for the shower development. A bump inthe profile at10 km is ob-served, this can be explained by the change in the atmosphere atthetroposphereboundaryasdiscussedintheprevious section. Therewillbeanerrorontheatmosphericdepthintroducedbythe fitted model in Eq- 6. It is on the order of 2 g/cm2 on average
between the altitude range mentioned above with a variance of 4− 5g/cm2.
The”gdastool” canbeexecutedasastandalonescriptwithin CORSIKA.Giventhe coordinateandUTC time stamp asinput pa-rameters it downloads the required GDAS files and extracts at-mosphericdata.It thenreturnsan outputfile thatcontainsfitted
P. Mitra, A. Bonardi and A. Corstanje et al. / Astroparticle Physics 123 (2020) 102470 5
Fig. 3. Left : Example of one density profile, GDAS and the fitted 5-layered atmospheric model. The bottom panel shows the relative error defined as ρfit−ρdata
ρfit . Right : Mean
relative error in density for 100 different atmospheric profiles. The mean is calculated at each of the 24 GDAS points for all the profiles. The error bars indicate the standard deviation.
massoverburdenparametersandtabulatedrefractiveindices inter-polatedto1 mintervals. Thisoutputfile canbe invokedthrough theCORSIKAsteering file.Whencalled, itreplaces thedefault at-mosphericparametersinCORSIKAwiththenewonesandthe on-the-flyrefractiveindexcalculationinCoREASwiththelook-up ta-ble.
5. Effectsonthereconstructionofthedepthoftheshower maximum
ThehighestprecisionforthedeterminationofXmaxwiththe
ra-diotechniqueiscurrentlyachievedwiththeLOFARradiotelescope. SituatedinthenorthoftheNetherlands,thedensecoreofLOFAR consistsof288low-banddipoleantennaswithinacirclewitha di-ameterof320meters,knownastheSuperterp.Theradioemission from airshowersin the frequencyrange 30–80MHz isrecorded by theLOFARlow-bandantennas [12,20].Anarray ofparticle de-tectorsinstalledon theSuperterpprovides thetriggerforthe de-tectionoftheairshowers [21].
The Xmax reconstruction techniqueused atLOFAR isbased on
the productionofdedicated simulation setsfor eachdetected air shower. The number of simulations needed to reconstruct the showermaximumisoptimizedwithCONEX [22].Asetoffull COR-SIKA simulations withproton andiron primaries is produced for each detected cosmic ray. The radio emission is simulated in a star-shaped patternforantennapositionsintheshowerplane us-ingCoREAS.Anantennamodelisappliedtothesimulatedelectric fieldsandcomparedtothemeasuredsignalinthedipoleantennas
[23].Thetimeintegratedpulsepoweriscalculatedina55ns win-dowcenteredaroundthepulsemaximum,summedoverboth po-larizations.Finally,a two-dimensionalmap ofthetimeintegrated power iscreatedby interpolatingthestar-shaped pattern [24].In thepreviousanalysisahybridfittingtechniquewasusedinwhich boththeradioandparticledatawerefittedtothetwo-dimensional radiationmapandtheone-dimensionalparticlelateraldistribution function simultaneously. Inthiswork insteadofthe combinedfit we fitonly theradio datatothe radiosimulation. Theadvantage of switchingto theradio only fittingmethod isthat it resultsin reducedsystematicuncertainties.
Fig. 4 shows the fit quality for an air shower detected with LOFAR as a function ofXmax simulated withtwo different
atmo-spheres - one withthe corresponding GDAS atmosphereandthe otherwiththeUSstandardatmosphere.Thereconstructedvalueof Xmaxisfoundfromtheminimumofthefittedparabolaaroundthe
best fittedpoints.We chosea LOFAReventforwhichthe ground
pressure was much lower than the US standard atmosphere, by 20 hPa. The atmospheric profile for thisparticular eventis rep-resented by the blue line with circles in Fig. 2 (right). The re-constructedXmax withtheUSatmospherecorresponds toa much
higher mass overburden than the reconstructed Xmaxusing much
thinnerGDASatmosphere.Inthisexamplethistranslatestoa dif-ference of around 37.5 g/cm2 in thereconstructed X
max between
the two cases. This large deviation is attributed to the extreme weather conditionfor the shower chosen inthe example. In the previousLOFAR analysisacorrection factortotheUSatmosphere wasused to account forthe real atmosphere [3,24]. The simula-tions that are producedwith USstandard atmospherewould ap-proximately yield the correct geometrical altitude to the shower maximum.ThenthecorrectedXmaxiscalculatedbyintegratingthe
GDASdensityprofileobtainedatLOFAR,fromthetopofthe atmo-spheretothegeometricaltitudeofXmaxinthefollowingway: X
(
h)
= 1cos
θ
∞h
ρ
gdas(
h
)
dh. (8)The corrected Xmax for thisparticular example is658 g/cm2 and
the difference between the corrected and new Xmax is about 20
g/cm2.
Usingthesameapproachdescribedabovewehavestudied123 airshowersrecordedwithLOFARwiththreesimulationsets:
• Set A–theshowersweresimulatedwithCORSIKAv-7.6300and GDASatmosphere.
• Set B–theshowerswere simulatedwithCORSIKAv-7.4385and USstandardatmosphere.
• Set C–thisset isidenticaltoSet Bbutwiththeadditional at-mosphericcorrectionfactortoitasdescribedabove.
The effect of using different CORSIKA versions on the recon-structedXmax,irrespectiveoftheatmosphericmodel,wasprobed.
The difference inXmax found usingCORSIKAversions 7.6300and
7.4385 was found to be very small, around 1.4 g/cm2. This
con-firmsthatthedifferencesbetweenSet-A,Set-BandSet-Caredue todifferentatmosphericmodels,not anyartifactarisingfrom dif-ferentversionsofCORSIKA.
In Fig.5thedifferenceinmeanreconstructedXmaxbetweenthe
varioussimulationsetsmentionedaboveisplottedagainstground pressurebins obtained fromGDAS. Boththe blue circlesandred squaresconvergetozerowhereGDASpressureapproachestheUS standardpressureat1013hPa.The redsquareshavelarge
Xmax
Fig. 4. Quality of fit as a function of simulated X max for a LOFAR event of energy 1.4 × 10 8 GeV, with a zenith angle of 38 ◦. Left : simulated with default US standard atmo-
sphere, reconstructed X max = 675.8 g/cm 2 . Applying the linear first order atmospheric correction, the resulting X max = 658 g/cm 2 . Right : simulated with GDAS atmosphere,
reconstructed X max = 638.3 g/cm 2 , the reconstructed X max in both the cases is indicated by solid black lines.
Fig. 5. Difference in mean X max as a function of ground pressure. The total sample
contains 123 air showers recorded at LOFAR. The black line denotes the U.S standard atmospheric pressure.
involvedinSet-B. Theblue circleshowevershow a higher devia-tionboth atlowandhighpressure values.Thissuggeststhat the linearfirstordercorrection addedtothestandardUSatmosphere implementedinSet-Cis notsufficient.As therefractiveindex ef-fectscannot be included inthelinear firstorder correction,one needsfull GDAS-basedatmosphericprofiles formoreextreme at-mosphericconditions.
Here, we studythe possibility to introduce a newglobal cor-rectionfactor to the reconstructed Xmax withUS standard
atmo-spheretocorrectforrealistic atmospshereswithouthavingtorun full GDAS-based CoREAS simulations. Toachieve this we studied thecorrelationbetweenXmax,refractivity, andslanted mass
over-burdenwhichis definedastheintegrated densityfromtheedge oftheatmospheretoagivenheightattheslantofzenithangle,at differentaltitudes.It wasseen that boththe correlation between Xmax and refractivity and between Xmax and slanted mass
over-burdencorrelation are poor atgroundand atlower altitudes. At thehigheraltitudes,between4-6km,Xmaxandmassoverburden
showa higher correlationwhich isnot prominentin Xmax vs
re-fractivityprofilesat thesealtitudes.We havefound the strongest correlation at an altitude of 5 km. Fig. 6 (left) shows the
scat-ter plot of
Xmax defined as Xmaxgdas− Xmaxus and difference in the
slantedmassoverburden
X5km=X5kmgdas− X5kmus .Theprecise
corre-lationsuggeststheprofilecanbefitwithastraightlineandisused asa parameterizationofglobalcorrection factor,provided by the equation: Xcorr max− Xmaxus =0.9
Xus 5km− X gdas 5km +0.28. (9)Thehistogramin Fig.6(right)showstheresidualoftheXcorr maxfrom
Xmaxgdas.The profile is symmetricwithmean 0g/cm2 andstandard
deviation11.56g/cm2.Thefluctuationsare withinthetypical
sys-tematicuncertaintyofthe reconstructedXmax withLOFAR, which
is around 17g/cm2 [24]. This correction factorcan be used asa
rule ofthumb for theestimation of reconstructed Xmax withthe
followingcaveats.Itisspecific toLOFAR,assimulationswere per-formedinvolving weatherconditions,observationlevel,and mag-neticfieldparticulartoLOFAR.Correspondingcorrectionequations forotherexperimentscanbeconstructedinthesamemannerand canyielddifferentresultsdependingonatmosphericparameters.
However,whilethisglobalcorrectionfactorisveryusefulwhen a fastreconstruction is needed, we willuse thefull Monte Carlo approachinafuturecomposition analysis.Simulationswithevent specificGDASatmospheresarealwaysmoreaccuratethanthe cor-rectionfactor.Thecorrectionfactormightalsointroducebiases re-latedtothemassoftheprimaryparticles.Protonprimarieson av-eragegenerate showersthat reach maximumlower in the atmo-spherethaniron;thesekindofeffectsarenottakenintoaccount. 6. Effectsofhumidity
As described in section 2, inthe radio frequencyregime, hu-midity increases the refractive index. Forthis study, two sets of simulations were produced. In one set the showers were simu-latedwiththerespectiveGDASatmosphereandintheotherwith a GDAS atmosphere with vanishing humidity. This wasachieved byhard-codingthepartialwatervaporpressureinEq-2to negligi-blevalues.FortheGDASatmosphereanextremelyhumidweather conditionattheLOFARsitewaschosen.Thesameatmospheric pa-rametersareusedinbothcasestoensurethattheparticlesevolve in a similar way in the atmosphere and produce same shower maximum. Inthis waythe inclusion of humidity only influences the simulatedradio pulses. Thedifference in the refractiveindex manifestsintermsofpropagationeffectsonthepulsearrivaltime and power. The pulse propagating though an atmosphere with
P. Mitra, A. Bonardi and A. Corstanje et al. / Astroparticle Physics 123 (2020) 102470 7
Fig. 6. Left : scatter plot of X max = X maxgdas − X maxus vs difference in slanted mass overburden X 5km = X 5kmgdas − X 5kmus . The red line is a linear fit to the profile. Right : Histogram
shows the residual of fitted and actual X max ; residual = X maxcorr − X maxgdas .
Fig. 7. Unfiltered electric field components of a CoREAS pulse in time for two dif- ferent refractive index profiles for a 10 17 eV proton shower with a zenith angle of
45 ◦coming from east for an observer at 150 m from the axis. The solid and dashed
lines represent the profiles with lower and higher refractive indices respectively.
higherrefractiveindexwillhavealowervelocitycomparedtodry air.Thisresults inadelayedarrival timeof thesignal,asseenin
Fig.7.The differenceinpeakarrival timeislessthan1 nsforan observerat150m.Theeffectisfoundtobelessprominentfor ob-servers furtheraway fromthe axis.The lateraldistributionofthe energyfluence,thetime-integratedpowerperunitarea,for differ-entobserverpositionsisalsostudiedfordifferentfrequencybands forthesetwocases,asshownin Fig.8.Inthelowfrequencyband of30–80MHzrelevantforLOFARthedifferenceinthefluence be-tweenthetwosetsissmall,fromaround4%closertoshoweraxis to2% atadistanceof100mfromtheaxis.Inthehighfrequency bandof50–350MHzthevaluesarelarger,beingaround8%at100 mfromthecore.InthehigherfrequencybandtheCherenkov-like effects become stronger and the signal is compressed along the Cherenkovring [25].Aroughestimateoftheradiusoftheringcan be obtainedfromthe projectionofa conewithanopeningangle givenbytheCherenkovanglestartingfromtheshowermaximum. The opening angle is stronglydependent on theindex of refrac-tion. Thisexplains the higherdifference in power in Fig. 8. Sim-ilar effects in high andlow frequency bands were also reported in [15] by studying the LDF of the electric field profiles. Inside
theCherenkovradius pulsesarestretched dueto refractiveindex effects.For higherrefractive indicesthis willlead to lower pulse powerwhichexplainsthenegativesignintherelativefluencefor observerdistancesclosetothecore.
Theradiationenergyisthetotalenergycontainedintheradio signal.Itscalesquadraticallywiththecosmicrayenergy,thuscan beusedasacosmic rayenergyestimator [26,27].Thesurface in-tegralovertheradioLDFmentionedaboveyieldstheradiation en-ergy.TherelativedifferenceintheintegratedLDFbetweenthe hu-midandnon-humid profiles forboth thelowandhighfrequency regimesissmallerthan1%.Thisindicatesthathumidityhasalmost noeffectontheestimatedcosmicray energyasdeterminedfrom theradiationenergywhichwasalsoconcludedin [28].
Next, to investigate the effect of humidity on Xmax
measure-ments we have performed a Monte Carlo comparison study be-tweentwosetsofsimulationsthatdealswiththeatmospheresina similarwayasdescribedinthebeginningofthissection.Foreach ofthesescaseswehaveusedasetof40simulatedeventswith dif-ferentenergy,zenithandazimuthangles.Eachof thesesets con-sistofanensembleofprotonandironinitiatedshowersbasedon CONEXselectioncriteria.Oneshowerfromthesetwithhigher hu-midity is takenas referenceand all the simulated showersfrom thesetwithzerohumidityareusedtoperformthereconstruction. Thisyields areconstructedXrecothat canbe comparedto the
ac-tual Xreal ofthe reference shower.The samemethod is repeated foralltheshowersinthesetwithhigherhumidity.Showerswith extremevaluesofXmax were notincludedinthefit.Therangeof
thefitwastakenas ± 50g/cm2 oftheactual X
max forthe test
shower.
Thedifference Xreco− Xreal estimatestheeffectofhumidity on
thereconstructedXmax.Wedonotobserveanysignificantshiftin
Xmaxinthisstudy.Thisindicatesthattheseeffectsaremostlikely
smallerthantheoverallresolutioninreconstructedXmaxinthe
LO-FARfrequencyband.Wealsoperformedthesamestudyinahigher frequency band between 50 and350 MHz, corresponding to the SKA-low band. There, an overall shift of 6.8g/cm2 in the
recon-structedXmax wasobserved.Theseresults,shownin Fig.9,arein
linewiththeLDFstudiesdescribedearlierinthissection. In Ref.[11], larger shifts of about 10 to 22 g/cm2 in
recon-structedXmax inthehighfrequencyband of120–250MHzfor4%
higherrefractivityand3.5to11g/cm2 inthelow frequencyband
of30–80MHz were reported.A toy model wasused todescribe theeffects.Thetoymodelwasbasedontheassumptions thatthe sizeoftheradiofootprintonthegroundwouldbeproportionalto
Fig. 8. LDF profiles for a 10 17 eV proton shower coming from zenith 45 ◦with X max = 593 g / cm 2 . Observers are located to the west of the shower axis. Left : low frequency
band between 30–80 MHz, Right : high frequency band between 50–350 MHz. The upper panel shows the LDF of total fluence for the humid and non-humid sets, the lower panel shows the relative difference between these two.
Fig. 9. Histogram for the X max = X reco − X real between the reconstructed and true value of the X max obtained from the Monte Carlo study between the humid and non-humid
simulation sets. Left : for the low frequency band of 30–80 MHz. Right : for the high frequency band of 50–350 MHz. The shift in the X max is significant at 2 σ level.
thegeometric distancetoXmaxandtotheCherenkovangleatthe
altitudeof Xmax. The effect of constant higher refractivity would
correspondtoahigherCherenkovangleresultinginan underesti-mationofXmax.Thisthen leadstoa clearlinearrelationbetween
shiftin Xmax and distanceto Xmax. Without having prior
knowl-edgeofindividualatmosphericconditions,anoverallscalingofthe refractivityprofile had to suffice. However, the realistic scenario is quite different. There are strong interplays betweenhumidity, pressure,andtemperaturewhich are reflectedin refractivity.The relative refractivityprofile in Fig. 2 showsthat the shift isnot a constant, butis altitude dependent. From near ground to higher altitudesit switches frombeinga highervalue than USstandard atmosphereto alower value.Thismakes anone-to-one compari-sontoRef.[11]hard.However,wecanarguethatqualitativelysame traitinthe highandlow frequencyband hasbeenfoundinboth theworks.
The effects of different zenith angles, true Xmax and energy
wereprobedfortheshiftinXmaxforboththefrequencybins.The
simulationsetwasdividedintwogroups,eachgroupbelongingto highandlow valuesofthe parameters mentionedabove. No sig-nificanteffectwasseen.
7. Conclusionanddiscussion
Simulatingairshowerswithrealisticatmospheresisimportant for the precise reconstruction of Xmax with the radio technique.
TheGDASdatabaseisausefulplatformtoextractatmospheric pa-rametersforagiventime andlocation.Atmosphericeffectson ra-diosimulationswerepreviouslystudiedinRefs. [11]and [15].The studiesdemonstratedtheroleofcorrectdescriptionofatmospheric densityandrefractiveindexwhenincludedintheradiosimulation
P. Mitra, A. Bonardi and A. Corstanje et al. / Astroparticle Physics 123 (2020) 102470 9
codes.However,theapplicationofsimulationswithrealistic atmo-spherestorealdatawasnotaddressed.
We report, for the first time, the application of GDAS-based atmospheric profiles, automated in CoREAS simulation to cosmic ray data.Bysystematicallyperforming GDAS-basedCoREAS simu-lations fortheLOFAR dataset,wehavedone comparisonbetween GDAS-basedatmospheresandalineargeometricalfirst order cor-rection tothe USstandard atmosphereonXmax.While thelinear
correctionissufficientforthebulkoftheevents,itbecomes indis-pensabletousefullGDASbasedatmospheresforextremevaluesof theairpressure.When theairpressureatgroundlevel differsby lessthan 10hPafromthe USstandardatmospherevalue, the re-constructedXmax valueincludingthelinearcorrectionagreeswith
the full GDAS-based reconstruction value within 2 g/cm2.
How-ever,whenthegroundpressureismorethan10hPafromtheUS standard atmosphere, thisdifference grows significantly up to 15 g/cm2.
We have also introduced a GDAS-based correction factor for Xmax reconstructed withUSstandard atmosphere withouthaving
torunfullGDAS-basedCoREASsimulations.ItisspecifictoLOFAR, butsimilar relationscan be workedoutforother experiments as well. The uncertainty on thepredicted Xmax usingthe correction
factorisabout12g/cm2;thisiswithinthetypicalX
max
reconstruc-tionuncertaintywithLOFAR,around17g/cm2.
Wehaveprobedtheeffectsofhumidityonthelateral distribu-tionofradiopowerby comparingtwoprofileswithhighandlow humidity. Weperformedthisstudyfordifferentfrequencybands. IntheLOFARfrequencybandof30–80MHztherelativedifference inpowerissmall.Forahigherfrequencybandof50–350MHzthe sameeffectsarecomparativelylarger,upto10%.Wealsoestimated theradiationenergyfromtheLDFprofilestoseetheeffectsof hu-midity onthe reconstructedenergy.No significant difference was found for eitherfrequency regime which indicates that humidity doesnot influencetheestimatedenergy.AMonte Carlostudyon the reconstructed Xmax wasalso done forthesefrequency bands.
No significant effect of humidity is found on the reconstructed XmaxforthelowfrequencybandrelevantforLOFAR;forthehigher
frequencyband ameandifference ontheorderof7g/cm2 is
ob-served. Thiscouldbe importantforthehighprecision Xmax
mea-surements forthecosmic raydetection withtheSKAexperiment
[29].
IntheprocessofimplementingGDAS-basedparameterized den-sity andrefractiveindexprofilein CORSIKA/CoREAS,we have de-velopedatool,called”gdastool”,whichhasbeenavailablefor pub-licusesincethereleaseofCORSIKAversion7.6300,andisalready being used by other experiments in the community around the globe.
In the previous LOFAR analysis the effects of refractive index were included within the systematic uncertaintieson the recon-structed Xmax. The improved atmospheric correction will lead to
a reduced systematic uncertainty. An update on the mass com-position results is not within the scope of this study. It will be discussed in a futurepublication, which involves, along with at-mosphericcorrections,improvedcalibrationoftheradioantennas, energyscale,andnewXmax reconstructiontechniques.
Declarationofcomeptinginterest
Theauthorsdeclarenocompetingfinancialinterests. Acknowledgement
The LOFARcosmic raykey scienceprojectacknowledges fund-ing from an Advanced Grant of the European Research Council (FP/2007-2013)/ERC Grant Agreement no 227610. The projecthas also received funding fromthe European Research Council (ERC)
undertheEuropeanUnion’sHorizon2020researchandinnovation program(grantagreementNo 640130).We furthermore acknowl-edge financial support from FOM, (FOM-project 12PR304). AN is supported by the DFG (Emmy-Noether grant NE2031/2-1 ). LO-FAR, the Low Frequency Array designed and constructed by AS-TRON,hasfacilitiesinseveralcountries,thatareownedbyvarious parties (each withtheir own funding sources),and that are col-lectivelyoperatedbytheInternationalLOFARTelescopefoundation underajointscientificpolicy.WesincerelythanktheCORSIKA de-velopersfortheir assistanceregardingtheimplementation ofour workinCORSIKAmodules.
References
[1] T. Huege, Radio detection of cosmic ray air showers in the digital era, Physics Reports 620 (2016) 1–52, doi: 10.1016/j.physrep.2016.02.001 .
[2] F.G. Schrder, Radio detection of Cosmic-Ray Air Showers and High-Energy Neu- trinos, Prog. Part. Nucl. Phys. 93 (2017) 1–68, doi: 10.1016/j.ppnp.2016.12.002 . [3] S. Buitink, et al., A large light-mass component of cosmic rays at 10 17 - 10 17.5
eV from radio observations, Nature 531 (2016) 70, doi: 10.1038/nature16976 . [4] W.D. Apel, et al., Reconstruction of the energy and depth of maximum of
cosmic-ray air-showers from LOPES radio measurements, Phys. Rev. D90 (6) (2014) 062001, doi: 10.1103/PhysRevD.90.062001 .
[5] G.A. Askar’yan , Excess negative charge of an electron-photon shower and its coherent radio emission, Sov. Phys. JETP 14 (2) (1962) 4 41–4 43 . [Zh. Eksp. Teor. Fiz.41,616(1961)]
[6] A. Nelles, et al., Measuring a Cherenkov ring in the radio emission from air showers at 110–190 MHz with LOFAR, Astropart. Phys. 65 (2015) 11–21, doi: 10. 1016/j.astropartphys.2014.11.006 .
[7] H. Schoorlemmer, et al., Energy and Flux Measurements of Ultra-High Energy Cosmic Rays Observed During the First ANITA Flight, Astropart. Phys. 77 (2016) 32–43, doi: 10.1016/j.astropartphys.2016.01.001 .
[8] J. Alvarez-Muniz, W.R. Carvalho Jr., E. Zas, Monte Carlo simulations of radio pulses in atmospheric showers using ZHAireS, Astropart. Phys. 35 (2012) 325– 341, doi: 10.1016/j.astropartphys.2011.10.005 .
[9] K.D. de Vries, A.M. van den Berg, O. Scholten, K. Werner, Coherent Cherenkov Radiation from Cosmic-Ray-Induced Air Showers, Phys. Rev. Lett. 107 (2011) 061101, doi: 10.1103/PhysRevLett.107.061101 .
[10] R. da, et al., First Experimental Characterization of Microwave Emission from Cosmic Ray Air Showers, Phys. Rev. Lett. 113 (22) (2014) 221101, doi: 10.1103/ PhysRevLett.113.221101 .
[11] A. Corstanje, et al., The effect of the atmospheric refractive index on the radio signal of extensive air showers, Astropart. Phys. 89 (2017) 23–29, doi: 10.1016/ j.astropartphys.2017.01.009 .
[12] P.S. et al., Detecting cosmic rays with the LOFAR radio telescope, Astronomy and Astrophysics 560 (A98) (2013), doi: 10.1051/0 0 04-6361/201322683 . [13] D.H. et al. , CORSIKA: A Monte Carlo code to simulate extensive air showers,
Report FZKA 6019 (1998) .
[14] B. Keilhauer, J. Blumer, R. Engel, H.O. Klages, M. Risse, Impact of varying at- mospheric profiles on extensive air shower observation: - Atmospheric den- sity and primary mass reconstruction, Astropart. Phys. 22 (2004) 249–261, doi: 10.1016/j.astropartphys.20 04.08.0 04 .
[15] F. Gaté, B. Revenu, D. García-Fernández, V. Marin, R. Dallier, A. Escudié, L. Mar- tin, Computing the electric field from extensive air showers using a realistic description of the atmosphere, Astropart. Phys. 98 (2018) 38–51, doi: 10.1016/j. astropartphys.2018.01.007 .
[16] V. Marin, B. Revenu, Simulation of radio emission from cosmic ray air shower with SELFAS2, Astropart. Phys. 35 (2012) 733–741, doi: 10.1016/j.astropartphys. 2012.03.007 .
[17] National oceanic and atmospheric administration, global data assimilation system, ( https://www.ncdc.noaa.gov/data- access/model- data/model- datasets/ global- data- assimilation- system- gdas ).
[18] P. Abreu, et al., Description of Atmospheric Conditions at the Pierre Auger Ob- servatory using the Global Data Assimilation System (GDAS), Astropart. Phys. 35 (2012) 591–607, doi: 10.1016/j.astropartphys.2011.12.002 .
[19] J. Rueger , Refractive index formulae for radio waves, Proceedings of FIG XXII International Congress (2002) .
[20] M.P. van Haarlem et al., LOFAR: The LOw-Frequency ARray, Astronomy and As- trophysics 556 (2013) 56, doi: 10.1051/0 0 04-6361/201220873 .
[21] S.T. et al., LORA: A scintillator array for LOFAR to measure extensive air show- ers, Nucl.Instrum.Meth A767 (2014) 339–346, doi: 10.1016/j.nima.2014.08.021 . [22] S. Buitink, et al., Cosmic ray mass composition with LOFAR, PoS ICRC2017
(2018) 499, doi: 10.22323/1.301.0499 .
[23] K. Mulrey, et al., Calibration of the LOFAR low-band antennas using the Galaxy and a model of the signal chain, Astropart. Phys. 111 (2019) 1–11, doi: 10.1016/ j.astropartphys.2019.03.004 .
[24] S.B. et al., Method for high precision reconstruction of air shower X max using
two-dimensional radio intensity profiles, Phys. Rev. D 90 (8) (2014), doi: 10. 1103/PhysRevD.90.082003 .
[25] A.N. et al., A parameterization for the radio emission of air showers as predicted by CoREAS simulations and applied to LOFAR measurements, As- tropart.Phys. 60 (2015) 13–24, doi: 10.1016/j.astropartphys.2014.05.001 .
[26] A. Aab, et al., Energy Estimation of Cosmic Rays with the Engineering Radio Array of the Pierre Auger Observatory, Phys. Rev. D93 (12) (2016a) 122005, doi: 10.1103/PhysRevD.93.122005 .
[27] A. Aab, et al., Measurement of the Radiation Energy in the Radio Signal of Extensive Air Showers as a Universal Estimator of Cosmic-Ray Energy, Phys. Rev. Lett. 116 (24) (2016b) 241101, doi: 10.1103/PhysRevLett.116.241101 .
[28] C. Glaser, M. Erdmann, J.R. Hrandel, T. Huege, J. Schulz, Simulation of Radia- tion Energy Release in Air Showers, JCAP 1609 (09) (2016) 024, doi: 10.1088/ 1475-7516/2016/09/024 .
[29] T. Huege, et al., High-precision measurements of extensive air showers with the SKA, PoS ICRC2015 (2016) 309, doi: 10.22323/1.236.0309 .
