Detecting cosmic rays with the LOFAR radio telescope
P. Schellart, A. Nelles, S. Buitink, A. Corstanje, J.E. Enriquez, H. Falcke, W. Frieswijk, J.R. Hörandel, A. Horneffer, C.W. James, M. Krause, M. Mevius, O. Scholten, S. ter Veen, S. Thoudam, M. van den Akker, A. Alexov, J. Anderson, I.M. Avruch, L. Bähren, R. Beck, M.E. Bell, P. Bennema, M.J. Bentum, G. Bernardi, P. Best, J. Bregman, F. Breitling, M. Brentjens, J. Broderick, M. Brüggen, B. Ciardi, A. Coolen, F. de Gasperin, E. de Geus, A. de Jong, M. de Vos, S. Duscha, J. Eislöffel, R.A. Fallows, C. Ferrari, M.A. Garrett, J. Grießmeier, T. Grit, J.P. Hamaker, T.E. Hassall, G. Heald, J.W.T. Hessels, M. Hoeft, H.A. Holties, M. Iacobelli, E. Juette, A. Karastergiou, W. Klijn, J. Kohler, V.I. Kondratiev, M. Kramer, M. Kuniyoshi, G. Kuper, P. Maat, G. Macario, G. Mann, S. Markoff, D. McKay-Bukowski, J.P. McKean, J.C.A. Miller-Jones, J.D. Mol, D.D. Mulcahy, H. Munk, R. Nijboer, M.J. Norden, E. Orru, R. Overeem, H. Paas, M. Pandey-Pommier, R. Pizzo, A.G. Polatidis, A. Renting, J.W. Romein, H. Röttgering, A. Schoenmakers, D. Schwarz, J. Sluman, O. Smirnov, C. Sobey, B.W. Stappers, M. Steinmetz, J. Swinbank, Y. Tang, C. Tasse, C. Toribio, J. van Leeuwen, R. van Nieuwpoort, R.J. van Weeren, N. Vermaas, R. Vermeulen, C. Vocks, C. Vogt, R.A.M.J. Wijers, S.J. Wijnholds, et al. (5 additional authors not shown)
AAstronomy & Astrophysics manuscript no. Detecting_cosmic_rays_with_the_LOFAR_radio_telescope c (cid:13)
ESO 2013November 7, 2013
Detecting cosmic rays with the LOFAR radio telescope
P. Schellart , A. Nelles , , S. Buitink , , A. Corstanje , J. E. Enriquez , H. Falcke , , , , W. Frieswijk ,J. R. Hörandel , , A. Horneffer , C. W. James , M. Krause , M. Mevius , , O. Scholten , S. ter Veen ,S. Thoudam , M. van den Akker , A. Alexov , , J. Anderson , I. M. Avruch , L. Bähren , R. Beck ,M. E. Bell , P. Bennema , M. J. Bentum , G. Bernardi , , P. Best , J. Bregman , F. Breitling ,M. Brentjens , J. Broderick , M. Brüggen , B. Ciardi , A. Coolen , F. de Gasperin , E. de Geus , A. deJong , M. de Vos , S. Duscha , J. Eislöffel , R. A. Fallows , C. Ferrari , M. A. Garrett , , J. Grießmeier ,T. Grit , J. P. Hamaker , T. E. Hassall , , G. Heald , J. W. T. Hessels , , M. Hoeft , H. A. Holties ,M. Iacobelli , E. Juette , A. Karastergiou , W. Klijn , J. Kohler , V. I. Kondratiev , , M. Kramer , ,M. Kuniyoshi , G. Kuper , P. Maat , G. Macario , G. Mann , S. Markoff , D. McKay-Bukowski , ,J. P. McKean , J. C. A. Miller-Jones , , J. D. Mol , D. D. Mulcahy , H. Munk , R. Nijboer , M. J. Norden ,E. Orru , R. Overeem , H. Paas , M. Pandey-Pommier , R. Pizzo , A. G. Polatidis , A. Renting ,J. W. Romein , H. Röttgering , A. Schoenmakers , D. Schwarz , J. Sluman , O. Smirnov , , C. Sobey ,B. W. Stappers , M. Steinmetz , J. Swinbank , Y. Tang , C. Tasse , C. Toribio , J. van Leeuwen , , R. vanNieuwpoort , , R. J. van Weeren , N. Vermaas , R. Vermeulen , C. Vocks , C. Vogt , R. A. M. J. Wijers ,S. J. Wijnholds , M. W. Wise , , O. Wucknitz , , S. Yatawatta , P. Zarka , and A. Zensus (Affiliations can be found after the references) Received 16 September 2013 / Accepted 4 November 2013
ABSTRACT
The low frequency array (LOFAR), is the first radio telescope designed with the capability to measure radio emissionfrom cosmic-ray induced air showers in parallel with interferometric observations. In the first ∼ of observing,405 cosmic-ray events in the energy range of − eV have been detected in the band from −
80 MHz . Each ofthese air showers is registered with up to ∼ independent antennas resulting in measurements of the radio emissionwith unprecedented detail. This article describes the dataset, as well as the analysis pipeline, and serves as a referencefor future papers based on these data. All steps necessary to achieve a full reconstruction of the electric field at everyantenna position are explained, including removal of radio frequency interference, correcting for the antenna responseand identification of the pulsed signal. Key words. astroparticle physics – methods: data analysis – instrumentation: interferometers
1. Introduction
With the development of ever faster electronics and the in-crease in computational power, the construction of radiotelescopes as large interferometric arrays of rather simpleantennas opens a new window for observations. The lowfrequency array (LOFAR; van Haarlem et al. 2013), is thefirst large-scale implementation of this technique. In addi-tion to producing the first high quality images at these lowfrequencies of −
240 MHz , LOFAR was designed to studyshort, pulsed signals in the time-domain. With a vast arrayof antennas observing the whole sky simultaneously, obser-vations are not limited to a predefined direction, thereforeproviding optimal conditions for cosmic-ray detection.Cosmic rays, accelerated charged particles from astro-physical sources, can be observed over several decades ofenergy. When cosmic rays of high energies reach the Earth,they do not reach the surface as primary particles, but in-stead interact with atmospheric nuclei. Thereby, a cascadeof particles is created, consisting mostly of photons and asignificant fraction of charged particles. While propagat- ing through the atmosphere, the charged particles of thisextensive air shower emit electromagnetic radiation, whichadds up coherently for wavelengths comparable to the di-mensions of the shower front (Huege 2013).Already in the 1960s it was proven that cosmic ray-induced air showers emit nanosecond duration pulses withsignificant power in the MHz radio frequency range (Al-lan & Jones 1966; Jelley et al. 1965), but due to lack ofsufficiently sophisticated and fast electronics the techniquewas not pursued further. Only in the past decade interestin the detection technique was rekindled and successfullyapplied (Falcke 2008). The proof of principle and largeprogress in the understanding of the emission was made atthe LOFAR Prototype Experimental Station (LOPES; Fal-cke et al. 2005; Huege et al. 2012) and further refined bymeasurements at the CODALEMA experiment (Ardouinet al. 2005).Similar to optical measurements of the fluorescenceemission from atoms excited by interaction with the airshower, radio emission directly traces the longitudinal
Article number, page 1 of 16 a r X i v : . [ a s t r o - ph . I M ] N ov hower development, which is closely related to the type ofthe primary particle. Unlike optical fluorescence measure-ments, radio emission measurements are less dependent onobserving conditions and can operate day and night match-ing the duty cycle of particle detector measurements.Due to the very steep energy spectrum, measuring thehighest-energy cosmic rays requires vast detector areas.Cost constraints therefore limit the density of detectorswithin this area giving a wide spacing between the individ-ual antennas. Theoretical models describing the differentemission mechanisms at play point to a very detailed andnon-symmetrical emission pattern at ground level (Werneret al. 2012; Alvarez-Muñiz et al. 2012; Marin & Revenu2012; Huege et al. 2013). Testing these models therefore re-quires dense sampling of the electric field over a sufficientlylarge area.LOFAR offers a high number of antennas clustered onan irregular grid, with increasingly large spacing betweenantenna clusters further away from the center. In the coreof the array about 2300 antennas are installed within about , which allows air showers to be measured with un-precedented spatial resolution. These measurements willcontribute significantly to conclusively confirm theoreticalmodels for the radio emission on a shower by shower basis,a goal previously unattainable due to lack of sufficientlyhigh quality data.Measurements and converging theoretical predictions ofthe expected radio signal from a cosmic-ray induced airshower give a short, nanosecond time-scale bi-polar pulse,which is mostly linearly polarized. This article describes thedetection set-up and automated processing pipeline used atLOFAR to measure and identify these signals.Starting with a description of the instrumental set-upat LOFAR in Sect. 2, an overview of the data reductionpipeline is given in Sect. 3. Finally, Sect. 4 describes thecharacteristics of the dataset obtained between June 2011and April 2013.The LOFAR dataset will be used in forthcoming publi-cations to verify existing models for radio emission from airshowers and to develop new techniques that use radio emis-sion to measure important characteristics of the incomingparticle, such as energy and mass.
2. LOFAR
LOFAR is a distributed radio telescope. Its antennas aredistributed over northern Europe with the densest concen-tration in the north of the Netherlands, in the Provinceof Drenthe. The observation support center and process-ing facilities are also located near this central core. Theantennas of LOFAR are grouped into stations , each sta-tion taking the role of a single dish in a traditional radiointerferometer array. A station consist of a number of low-band antennas (LBAs, −
90 MHz ) and high-band anten-nas (HBAs, −
240 MHz ). The 24 stations within the ∼ wide core are distributed in an irregular patternthat maximizes uv -coverage, or spatial frequencies for stan-dard interferometric observations. The 16 additional Dutchremote stations are distributed with increasing distance tothe core. International stations are currently located inGermany, France, the United Kingdom, and Sweden, givingLOFAR a maximum baseline of for interferomet-ric observations. Core stations and remote stations consistof 96 LBAs plus 48 HBAs. International stations have 96 x [m]-200 -100 0 100 200 300 y [ m ] -200-150-100-50050100150 Fig. 1:
Layout of the center of LOFAR. The six stationsto the left form the Superterp. The crosses indicate theLBA inner and outer antenna sets, respectively. The opensquares show the positions of the HBA tiles, which are splitinto two groups per station. The filled squares indicate thepositions of the LORA particle detectors.
Fig. 2:
Low-band antennas at the central core of LOFAR,the Superterp. In the background the black box of a LORAparticle detector can be seen.LBAs and 96 HBAs. At the center of the LOFAR coresix stations are located in a roughly
320 m diameter area,called the
Superterp , providing both the shortest baselinesfor interferometric observations and the densest populationof antennas ideal for cosmic-ray observations. While everyLOFAR station is equipped with the necessary electronicsto observe cosmic rays, the current data set is taken withthe central 24 stations, where additional information fromparticle detectors is available (see Sect. 2.3). The positionsof the antennas of the seven most central LOFAR stationsare shown in Fig. 1.
The LBAs are the main tool for cosmic-ray detection. AnLBA consists of two orthogonal inverted V-shaped dipoles,each with a length of .
38 m . These are supported by acentral polyvinyl chloride pole, which holds the low-noise
Article number, page 2 of 16. Schellart et al.: Detecting cosmic rays with the LOFAR radio telescope
W ESN
X Y ˆ e x ˆ e y ◦ Fig. 3:
Geometry of the LBA. The X and Y dipoles areoriented NE-SW and NW-SE respectively. This is rotatedby 225 degrees with respect to the standard local Cartesiancoordinate system used in Sect. 3.4.amplifier and guides the signal cables, as shown in Fig. 2.The dipoles X and Y , that make up each antenna, areoriented southwest to northeast (SW-NE) and southeast tonorthwest (SE-NW), as can be seen in Fig. 3.The low-noise amplifier has an intentional impedancemismatch with the antenna. This mismatch, combined withthe characteristic length of the dipoles, makes the systemsensitive in a broad band from −
90 MHz . In principle,this allows observations from the ionospheric cutoff up tothe start of the commercial FM radio band. For most ob-servations the frequency range is limited by a combinationof selectable hardware and software filters to −
80 MHz to suppress strong Radio Frequency Interference (RFI) inthe outer bands. The LBAs are designed to be sky noiselimited after RFI has been removed (van Cappellen et al.2007). After amplification the signals from the individualdipoles are transmitted through coaxial cables to the elec-tronics cabinet located at every station.The HBAs have been optimized for a frequency band of −
240 MHz . The design clusters 16 antenna elementsinto a “tile”, the signals from these elements are amplifiedand combined in an analog beam-former. This means thatwhile the LBAs are sensitive to the whole sky the HBAs aremost sensitive within the ∼ ◦ of the tile-beam, of whichthe direction is chosen at the start of every observation.This results in a smaller effective area for cosmic-ray ob-servations, as the measurement will only be optimal if thedirection of the cosmic ray happens to coincide with thebeam direction of the observation. Therefore, the analysisof HBA data and their interesting higher frequency rangerequires a different approach for cosmic-ray studies. Re-sults of these measurements will be described in a laterpublication. After being forwarded to the electronics cabinet the sig-nals of the LBAs are again amplified, filtered, and digitizedby a
12 bit
A/D converter with a sampling frequency of
200 MHz . Due to signal path limitations in the Dutch sta- A 160 MHz clock is also available.
Number of detectors triggered6 8 10 12 14 16 18 20 E n e r g y t h r es ho l d ( P e V ) Energy thresholdEvent rate E ve n t r a t e p e r d ay . / h r , P e V . / h r , P e V Fig. 4:
Energy threshold in PeV (left) and the event rateper day (right) are shown as a function of the number oftriggered particle detectors. Two possible trigger conditionsare indicated with the dotted lines.tions only 48 dual-polarized or 96 single-polarized antennascan be processed at a given time. For the dual-polarizedoption the antennas are grouped into an inner and an outerset, which has to be chosen before an observation.For astronomical observations the data are then beam-formed and sent to the central processing facility. In ad-dition, there is the possibility to store a snapshot of theoriginal data. Every station is equipped with ring-buffers,the so called Transient Buffer Boards (TBBs). These con-tinuously store the last . of raw data (an extension to is currently being deployed). When triggered, the contentsof the TBBs are frozen, read out via the Wide Area Net-work and stored on disk for further analysis. The triggercan be generated based on various parameters in an FPGA at the local receiver unit. Alternatively, the trigger can begenerated by an array of particle detectors (see Sect. 2.3) orreceived from outside of LOFAR. Currently, the main trig-ger for cosmic-ray observation is provided by the particledetectors. Later, a radio self-trigger will be implemented,using the current dataset as a training set to deduce triggercriteria, so that the FPGA trigger can be run independentlyat every LOFAR station. These criteria have to reduce falsetriggers to limit the data rate. Using every LOFAR stationindividually will dramatically increase the effective area.Essential for measuring cosmic rays with LOFAR as aradio telescope is that the whole process of triggering andstoring radio-pulse data can take place without interferingwith the ongoing observations. LORA, the LOFAR Radboud Air Shower Array, is an arrayof particle detectors co-located with the center of LOFAR.The array provides a reconstruction of basic parameters ofrecorded air showers, such as the direction and the positionof impact, as well as the energy of the incoming cosmicray (Thoudam et al. in prep.). It also provides the time ofarrival, which is used to trigger the read-out of the radioantennas. Field Programmable Gate Array.Article number, page 3 of 16
ORA consists of 20 detector units distributed on theSuperterp, as shown in Fig. 1. Each detector contains twoscintillators ( .
45 m , type: NE 114), which are individuallyread out through a photomultiplier tube. The detectors areinside weatherproof shelters and have been tested to notcreate any interference at radio frequencies.Conditions at which triggers are sent to LOFAR can beadjusted to match the desired energy threshold. There aretwo constraints on the desired rate: the rate of events inter-esting for radio observations has to be maximized, while thenetwork load on the LOFAR system has to be kept low inorder to avoid interfering with the primary observation. Atrigger in a single detector is generated when a particle sig-nal of more than σ above the noise is registered. In orderto only detect air showers a coincidence of several detectorsis needed. Events of less than eV have a very low prob-ability to be observable in radio above the sky-noise level.The energy threshold and the corresponding event rate areshown in Fig. 4 as the function of the number of triggereddetectors. Requiring triggers in 13 detectors yields a thresh-old energy of . · eV , with an average trigger-rate of . / hour . This trigger rate has been selected as theoptimal setting for the observations. After the commissioning phase LOFAR is to be used ona proposal-based schedule. Proposals are open to thecommunity for imaging or beam-formed observations, aswell as TBB observations. Some fraction of the observ-ing time is reserved for participating consortia and key sci-ence projects. The LOFAR cosmic ray key science project(CRKSP) is one of six LOFAR key science projects.To maximize the duty cycle TBB observations can berun in the background of all other observations that do notneed the full network bandwidth. This does however meanthat the array configuration is determined by the primaryobservation, therefore the amount of data in a specific arrayconfiguration (such as the selection of LBA or HBA antennatype) available for analysis is not determined by the cosmic-ray project itself, except when LOFAR is otherwise idle andthe observing configuration can be chosen freely.During the observation, triggers from LORA are re-ceived by the LOFAR control system. The system checkswhether a dump from the TBBs is allowed. If so, the ring-buffers are frozen and a specified block of data around thetrigger time is dumped to disk. For each cosmic-ray event . of radio data are stored, which corresponds to
77 MB per station. This provides sufficient frequency resolution forhigh quality RFI cleaning while minimizing data transferand storage requirements.Every evening, the data-files are archived at LOFARand compressed for transport. They are stored in the LongTerm Archive (van Haarlem et al. 2013), from where theycan be retrieved for data analysis.
3. Reconstruction of cosmic-ray data
All newly recorded data are processed every evening, afterhaving been copied via the network to the processing clusterof the Astrophysics department at the Radboud UniversityNijmegen. In addition to the hdf5 files, containing the A tree-like file format (Alexov et al. 2012). data of one LOFAR station each, the recorded data fromthe particle detectors and a trigger log file are transferred.With this information an automated pipeline is run. Thepipeline is based on the task oriented PyCRTools frame-work consisting of fast low-level C++ routines embedded inPython for maximum flexibility. All results are stored in aPostgreSQL database for subsequent data mining analysis.The goal of the processing pipeline is to autonomously iden-tify a full set of physics quantities for each air shower de-tected with LOFAR. The pipeline is optimized to identifythose nanosecond pulses that are not generated by terres-trial sources.All data are first processed per station, i.e. per file. Theset of files received for a single trigger form an event . Whenthe data from one station pass the criteria for containinga cosmic-ray signal (see Sect. 3.3), the corresponding eventis called a cosmic-ray event . It is not necessary to observea pulse in all stations, only the stations with a significantsignal are used in a combined analysis.
The reconstruction pipeline comprises a number of stepsthat will be individually explained in the following sections.An overview of the steps and the overall structure is de-picted in Fig. 5.
Before proceeding to extract the cosmic-ray signal fromthe data, some preparatory steps have to be performed.Knowledge about the system is applied in the form of cal-ibration procedures, the data are cleaned of narrowband-transmitters, and antennas that show malfunctions areflagged.
There are known signal path differences between the LO-FAR antennas. Measured differences of cable lengthsbetween the antennas are corrected for up to the sam-ple level already at the stations before the data are writ-ten to disk. Additionally, relative time offsets between theantennas are corrected for at sub-sample accuracy usingstandard LOFAR calibration tables. These tables are gen-erated by phase-calibrating on the strongest astronomicalradio sources and are regularly tested and updated if neces-sary (van Haarlem et al. 2013). Sub-sample corrections areapplied as phase offsets to the Fourier transformed signalin the cosmic-ray pipeline, before processing it in the dataanalysis.
Narrow-band RFI in the time series signal can be revealedby making an average power spectrum. An example isshown in the top panel of Fig. 6, where most of the strongRFI is visible outside the −
80 MHz range. The averagepower spectrum is created by averaging the square of theabsolute value of the Fourier transform over several blocksof data. The block size can be freely chosen within thefull data length to obtain a desired frequency resolution;here samples are used, giving a resolution of ∼ , Article number, page 4 of 16. Schellart et al.: Detecting cosmic rays with the LOFAR radio telescope
Data fileBlockselection LORA timestampPhasecalibration LOFAR tablesRFIcleaningRFI iden-tification Gaincalibration(relative) Galactic noiseBeamformingUnfoldingantennaresponse Antenna modelPulsedetectionDirectionfittingCoordinatetransfor-mationUpdatedirection Extractingpulseparameters
Fig. 5:
General structure of the analysis pipeline. Rect-angles represent input and rounded squares are processingsteps.enough to resolve most RFI lines. A reasonable data lengthis needed for this procedure to produce a stable average,which sets the limit for the chosen block length to be storedfrom the TBBs, as mentioned in Sect. 2.4. In order to min-imize artificial side lobes a half-Hann window is applied tothe first and last 10% of each trace prior to the Fouriertransformation.The standard approach to RFI cleaning (or RFI flag-ging) is to identify peaks sticking out significantly above L o g - Sp e c t r a l P o w e r [ A D U ] L o g - Sp e c t r a l P o w e r [ A D U ] L o g - Sp e c t r a l P o w e r [ A D U ] Fig. 6:
The average spectrum of a typical LBA event.The raw data ( top ), with flagged contaminated channels( middle ), cleaned and clipped to −
80 MHz ( bottom ).the overall spectral shape, also called the baseline , and setthe corresponding Fourier component amplitudes to zero.However, this requires ‘a priori’ knowledge of the baseline.While the baseline can be obtained through a smoothing orfitting procedure, this is often not stable in the presence ofstrong RFI, requiring an iterative approach.An alternative approach to RFI cleaning uses the phaseinformation in the complex-valued spectrum instead. If anRFI transmitter is measured in all antennas, the phase dif-ference, or relative phase, between each pair of antennas Article number, page 5 of 16 ill be a constant value as function of time with a smallnon constant random noise contribution. Note that theexact value of the constant, which only depends on the ge-ometric delay between antennas, is not relevant, only itsnon time-varying nature. When no transmitter is present,the relative phase is expected to be both random and timevarying, as the signal then consists of the added signalsfrom many incoherent sources on the sky with additionalrandom noise. Therefore, RFI can be identified by lookingat the stability of phase differences between antennas overtime. For each antenna-dipole j = 0 , , . . . , in a stationand data block k , the phase spectrum is calculated as φ j,k ( ω ) = arg( x j,k ( ω )) , (1)where x j,k ( ω ) is the complex frequency component ω of thespectrum.Subtracting the phase of one of the antennas as referenceantenna gives the relative phases and results in a set ofphases for every frequency channel, one for each block ofdata. Only one reference antenna is used and this is takento be the one with median power to avoid selecting a brokenantenna.The average phase is defined as ¯ φ j ( ω ) = arg (cid:32) N − (cid:88) k =0 exp( iφ j,k ( ω )) (cid:33) , (2)and the phase variance as s j ( ω ) = 1 − N (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) N − (cid:88) k =0 exp( iφ j,k ( ω )) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (3)where N is the number of data blocks.For completely random phases one expects s j ( ω ) ≈ as opposed to s j ( ω ) = 0 when all phases are equal. Thephase variance per frequency channel will now either be ata value close to 1, including some random ‘noise’, or at asignificantly lower level. The latter reveals the presence of aradio transmitter, as shown in Fig. 7, where a contaminatedpart of the spectrum is shown with the corresponding phasevariance.Since RFI lines will result in peaks toward smaller val-ues of the phase variance, and noise has no preferred peakdirection, calculating the standard deviation σ in this plotonly for values above the median will ensure a stable re-sult. All frequencies that have a phase variance of at least σ below the median are flagged as containing RFI. Ad-ditionally a −
80 MHz bandpass filter is applied, flag-ging the most heavily RFI polluted low and high frequencyparts of the bandwidth by default. To prevent pulse-ringing the −
80 MHz block filter is first convolved with a, σ tapering = 2 . , Gaussian . After removing the flaggedchannels, the resulting cleaned spectrum is shown in thebottom panel of Fig. 6.In general, there is very little RFI at the LOFAR Su-perterp. A lot of effort has been made to remove local Assuming a Gaussian distribution, σ can be estimated by sort-ing the data points, and comparing the value at 95 percentile tothe median. This difference amounts to ∼ . σ . This effect also occurs when flagging large blocks of RFI butthis does not happen in practice and so no tapering window isapplied for this case.
85 90 95Frequency [MHz]1110987654 L o g - Sp e c t r a l P o w e r [ A D U ] P h a s e V a r i a n c e Fig. 7:
Average LBA spectrum (bottom, left axis) withthe corresponding phase variance (top, right axis). RFIlines can clearly be identified in the phase variance withpeaks toward lower values, representing more stable phasedifferences between antennas over time.sources that could disturb the LOFAR measurements anda protected zone has even been established (Offringa et al.2013). This relative quietness is illustrated in Fig. 8. Itshows the result of the RFI cleaning for all events for allfrequencies. While every event has some RFI, no single RFIline is present in every event. Within the −
80 MHz band,there are only two lines that are present in more than 40%of the events. In total there are rarely events with morethan 2% flagged channels out of the more than 32000 fre-quency channels in a block of data. This is is shown inFig. 9, where the total fraction of events is plotted againstthe number of flagged channels.
Occasionally, one or more antennas give invalid signals, e.g.due to hardware malfunction. To identify these bad anten-nas the integrated spectral power is calculated P = (cid:90)
80 MHz30 MHz | x ( ω ) | dω, (4)where x ( ω ) is the ω frequency component of the cleanedspectrum. The power in every antenna is required to be inthe range of one half to two times the median power fromall antennas. Antennas outside this range are marked asbad and excluded from further analysis. There are ongoing efforts for an absolute calibration of thevoltage traces of LOFAR and therefore the reconstructedelectric field. Those efforts will be described in a forth-coming publication and include calibration on astronomi-cal sources, terrestrial transmitters, and already conducteddedicated measurement campaigns, similar to those per-formed at other experiments, e.g. Nehls et al. (2008). Onceimplemented, the reconstruction pipeline will deliver cali-brated electric field strengths and their polarization compo-
Article number, page 6 of 16. Schellart et al.: Detecting cosmic rays with the LOFAR radio telescope
Frequency [MHz]0 20 40 60 80 100 F r ac t i on o f eve n t s f l a gg e d Fig. 8:
Fraction of events that is affected by narrow-bandRFI in each of the ∼ frequency channels as functionof frequency. Number of flagged channels0 500 1000 1500 2000 2500 3000 F r ac t i on o f eve n t s -3 -2 -1 Fig. 9:
Relative fraction of events with a certain number offlagged channels. Over 60% of the events have less than 100channels ( ≈
300 kHz) flagged out of the full used bandwidthof more than 16000 channels.nents for all events. However, significant progress in under-standing the mechanisms of radio emission in air showerscan already be made with a relative calibration.
The LBA measurement is dominated by sky noise, which inturn is dominated by the Galaxy moving through the an-tenna beam pattern. Therefore, the noise as seen by eachantenna is a function of the Local Sidereal Time (LST) andcan be used to correct for differences in gain between an-tennas. Instead of correcting all antennas at all times toa fixed value, which would be over- or underestimating thenoise at certain times, the received power can be normal-ized to a LST-dependent reference value. In Fig. 10 theintegrated spectral power (equation 4), after RFI cleaning,is given as a function of LST for the instrumental polar-ization X and Y . The data have been retrieved from allcosmic-ray events measured within the first year of data-taking. One can define a reference value for the integrated : : : : : : : Local Sidereal Time0.0050.0060.0070.0080.0090.0100.0110.012 g a l a c t i c p o w e r [ a . u . ] : : : : : : : Local Sidereal Time0.0050.0060.0070.0080.0090.0100.0110.012 g a l a c t i c p o w e r [ a . u . ] Fig. 10:
Integrated spectral power normalized to the band-width, after RFI cleaning, as a function of local siderealtime for the X (NE-SW) ( top ) and Y (NW-SE) ( bottom ) in-strumental antenna polarizations. Also shown is the fittedsecond order Fourier transform (solid line). The uncertain-ties on the data still include systematic effects due to theset-up itself, as well as possible artifacts of the RFI clean-ing, when having certain frequencies that are contaminatedin a significant fraction of the data.spectral power as a function of LST by fitting a functionto these data points. Since the movement of the Galaxythrough the antenna beam pattern is periodic by nature itis fitted with the 2nd order Fourier series P ref ( t ) = a (cid:88) n =1 a n sin( nt ) + b n cos( nt ) , (5)thereby avoiding artificial jumps in the fit at 0:00 LST. Thetime t is given in units of radian here. This results in a gaincorrection for each antenna as x (cid:48) ( ω ) = (cid:115) P ref ( t ) P ( t ) x ( ω ) , (6)where the square root is needed, because the correction isapplied to the amplitude spectrum. Article number, page 7 of 16 [ns]
LOFAR - t
LORA (t-1500 -1000 -500 0 500
Fig. 11:
Difference in time between the time of a pulseidentified in the radio signal and the trigger time set bythe signal in the particle detectors. This plot shows thedistribution summed over all Superterp stations.
After cleaning and calibration of the data, the central ele-ment of the pipeline is the identification and characteriza-tion of the radio pulse as the signal of the air shower.
In order to restrict the search for the radio pulse to a smallerregion in the trace, the information from the trigger time ofthe particle detectors is used. Figure 11 shows the differencein time between the trigger from the particle detectors andthe pulse location in the radio data obtained from a searchwith a large window. The distribution shows a clear peakat the region of the coincidences at an offset of ±
168 ns .In absolute timing the offset between LORA and LOFARis , of which are already accounted for inthe triggering system.Average offset is obtained by fitting a Gaussian to thedistribution of pulse positions with respect to the triggertime. This is only an approximation, as the real offset perevent depends on the position of the core and the incom-ing direction of the air shower. Also, effects due to thepropagation of particles and radiation in the atmospherecan play a role. The overall difference is due to the factthat both detectors operate independently on different tim-ing systems. Both are based on GPS timing, but correctfor drifts ( <
20 ns ) in different ways and have a differingabsolute time. The spread on the differences is howeversufficiently small for Superterp stations to not require ad-ditional synchronization of the two systems. Stations fur-ther away can have larger offsets due to the signal traveltime, which can be corrected for after a reconstruction ofthe shower.These measurements allow for the pulse search to berestricted to a small fraction of the full time trace, limitingthe chance to pick up random noise fluctuations.
The trigger threshold of the scintillator array is chosen tobe lower than the threshold to detect a radio signal. Thisensures a full sample, but also makes it necessary to identifyin a first quality check whether there is a detectable signalpresent. Therefore, per antenna polarization, the signalsare first beamformed in the direction reconstructed fromthe data of the particle detectors. This direction is givenin the local Cartesian coordinate frame of the station by n and the position of each antenna j is given by r j . A planarwavefront arriving at the phase center (0 , , at time t = 0 will arrive at antenna j with a delay given by ∆ t j = − c n · r j | n | = − c ˆe n · r j , (7)where c is the speed of light. The beamformed signal, infrequency space, in this direction is then given by x bf ( ω ) = N a (cid:88) j =0 x j ( ω ) e πiω ∆ t j , (8)where x j ( ω ) is frequency component ω of the Fourier trans-form of the signal from antenna j and N a is the number ofantennas. The inverse Fourier transform gives the beam-formed time series signal. Due to beamforming any signalcoming from the direction of the air shower is amplified by afactor N a in amplitude while uncorrelated noise is only am-plified by a factor √ N a . Therefore, if no significant signal isdetected in the beamformed trace, the event very unlikelycontains a cosmic-ray signal strong enough to be detectedat single dipole level by the rest of the pipeline. Thus, theanalysis of the data of that station is aborted.To test this assumption, Fig. 12 shows the distributionof the peak amplitude in the beamformed signal per sta-tion, distinguishing between events in which ultimately acosmic-ray was identified and those in which there was not.The peak amplitude is normalized by the root mean squareof the trace, as a proxy for the noise contribution. Fromthis it can be seen that the fraction of events where a strongsignal is observed in the beamformed trace is significantlyhigher for stations where eventually a cosmic-ray signal isdetected. All events in the tail of the non-detected distri-bution were visually inspected and identified as broad-bandRFI, with pulses differing significantly in shape from thoseof cosmic rays and directions ultimately deviating signif-icantly from the direction as measured with the particledetectors. This distribution shows that an initial filteringbased on a moderate signal-to-noise of beamformed pulsesis a quick and effective way to filter out those events thatare potentially interesting, as well as further narrowing thesearch window per antenna reducing false positives for pulsedetection. The sensitivity of the LOFAR LBA is a complex function ofboth frequency and direction. Correcting for this antennapattern , i.e. unfolding, requires an initial guess for the pulsedirection and in turn may influence the position of the pulsein time and thus the direction by changing the phase atwhich each frequency arrives. Therefore the correction hasto be done in an iterative loop as indicated in Fig. 5. Each
Article number, page 8 of 16. Schellart et al.: Detecting cosmic rays with the LOFAR radio telescope
S/N in beamformed trace5 10 15 20 25 30 35 40 45 50 F r ac t i on o f eve n t s -4 -3 -2 -1
101 Pulse detectedPulse not detected
Fig. 12:
Distribution of the signal-to-noise ratio (S/N) ininitial beamforming. The S/N is defined the ratio of thepeak amplitude of the beam-formed trace and the RMS ofthis trace. Two cases are separated: a cosmic-ray eventwas ultimately detected by the pipeline (solid line) or not(dashed line). The initial cut, which is applied in thepipeline, is indicated by the dotted line.iteration starts with an increasingly accurate signal direc-tion and proceeds by unfolding the antenna pattern, pulsedetection, and direction fitting. The loop is concluded whenthe direction no longer significantly changes, which usuallyhappens in less than ∼ iterations.For the antenna pattern of the LBA a simulation isused, which is made using the software WIPL-D (Kolundz-ija 2011) and a customized software model of the electronicschain.From the impedance and radiation pattern in a trans-mitting situation the open circuit voltage is calculated as afunction of frequency and direction for an incoming planewave with an electric field strength of / m . The equiv-alent circuit of the antenna in a receiving situation is avoltage source with an internal resistance equal to the an-tenna impedance. This is combined with measured dataof the amplifier directly behind the antenna. The result ofthe model is the output voltage of the amplifier over a
75 Ω resistor .Any wave coming from a direction ˆe n can be seen asa linear superposition of monochromatic plane waves, po-larized in the ˆe φ and ˆe θ direction. Here φ and θ are thestandard spherical coordinate angles with the x and z axisrespectively, e.g. E ( t ) = (cid:88) ω ( E θ,ω ˆe θ + E φ,ω ˆe φ ) e − i ( k n · x + ωt ) . (9)This geometry can be seen in Fig. 13.These terms are related to the output voltage of theamplifier for each dipole, and for each frequency, via theJones matrix J (Jones 1941; Hamaker et al. 1996) of theantenna model (cid:18) V X V Y (cid:19) = (cid:18) J Xθ J Xφ J Y θ J Y φ (cid:19) (cid:18) E θ E φ (cid:19) , (10) Matched to the impedance of the coaxial cables connectingthe antenna to the station electronics cabinet.
East, x X North, yY Zenith, z φ θ ˆe θ ˆe φ ˆe n Fig. 13:
On-sky polarization coordinate frame ( ˆe θ , ˆe φ , ˆe n ).Also depicted is the (north, east, zenith) coordinate frameof the simulations, where the unit vectors ( ˆe x , ˆe y , ˆe z ) cor-respond to the x , y and z -axis, respectively. Furthermorethe dipole antennas X and Y are shown.where J Xθ is the complex response of the antenna and am-plifier of the X -dipole to a wave purely polarized in the ˆe θ direction.Therefore, in order to both correct for the antenna re-sponse and convert from output voltage to electric fieldstrength in the on-sky frame (see Sect. 3.4), each pair ofFourier components from the signal in the two instrumen-tal polarizations ( X , Y ) is multiplied by the inverse Jonesmatrix, followed by an inverse Fourier transform back tothe time domain.The components of the Jones matrix of the antennamodel are simulated on a grid with steps of in fre-quency, ◦ in θ and ◦ in φ . In order to obtain the com-ponents at the frequency and direction of observation, tri-linear interpolation is performed on the real and imaginaryparts of the complex table when needed. Examples of theresponse are depicted as a function of frequency in Fig. 14and as a function of direction in Fig. 15. Estimating the direction of the incoming air shower, seeSect. 3.3.5, can either be done using beamforming orthrough pulse timing. Beamforming was found to be verysensitive to the optimization algorithm used, essentially re-quiring a grid search to avoid getting stuck in a local mini-mum. This is computationally very expensive, moreover itonly provides relative time differences between any two an-tennas rather than an absolute time needed for extractionof relevant physical parameters (see Sect. 3.4).In order to use pulse timing, individual pulses have to beidentified. This can be done by using the cross-correlationmethod, where one looks for the maximum in the cross cor-relation of the signals between all antennas. This howeverhas the same drawback as beamforming, as only relativetiming is calculated. A method to retrieve the absolutepulse timing is through the use of the
Hilbert envelope ,which is used in this pipeline. A detailed comparison ofthe methods is given in Sect. 3.3.5.
Article number, page 9 of 16 -1 | ∆ V |
10 20 30 40 50 60 70 80 90Frequency [MHz] − π − π π + π P h a s e Fig. 14:
Jones matrix components of the antenna modelamplitudes ( top ) and phases ( bottom ) for a dipole receivinga wave polarized in the ˆe θ direction (circles) and a wavepolarized in the ˆe φ direction (stars) for an arrival directionof φ = 345 ◦ , θ = 50 ◦ . Also plotted, as the dotted line, arethe interpolated values.A sensible definition of the pulse arrival time is themeasured arrival time of the maximum of the electric fieldstrength. In practice, however, using directly max( | x ( t ) | ) is highly dependent on the filter characteristics of the re-ceiving system and the sampling used. Therefore, the ar-rival time is defined as the position of the maximum in theamplitude envelope of the analytic signal, also called theHilbert envelope A ( t ) = (cid:112) x ( t ) + ˆ x ( t ) . (11)Where ˆ x ( t ) is the Hilbert transform , or imaginary propaga-tion, of the signal x ( t ) defined by F (ˆ x ( t ))( ω ) = − i · sgn( ω ) · F ( x ( t ))( ω ) (12)where F denotes the Fourier transform.In order to find the pulse maximum with subsampleprecision, the signal is first up-sampled by a factor , such ∆ V Fig. 15:
An example Jones matrix component describ-ing the dipole response, at
60 MHz , | J X,θ | in the form ofthe output Voltage ( ∆ V ) as a function of direction for anincoming wave that is purely linearly polarized in the ˆe θ direction. µs )6040200204060 A m p li t u d e ( A D U ) pulse maximum SignalEnvelopeRMS
Fig. 16:
The solid light line shows the up-sampled signal.Overlaid is the Hilbert envelope and the RMS noise in blackdashes. A pulse is accepted whenever the signal to noiseratio exceeds three.that the maximum search will not be the limiting factorin the timing resolution. Subsequently, a simple maximumsearch is performed on the envelope. In addition, the signal-to-noise ratio is calculated where the signal is defined asthe maximum and the noise as the root mean square of theenvelope. An example can be seen in Fig. 16.This maximum search is performed on each of the on-sky polarizations E θ ( t ) and E φ ( t ) separately and any pulsewith a signal to noise greater than three is marked as apossible cosmic-ray signal to be used for direction fitting.Because the pulse is expected to be intrinsically strongerin one of the two polarizations, depending on the anglebetween the shower axis and the geomagnetic field, the po-larization with the highest average signal to noise (over all Article number, page 10 of 16. Schellart et al.: Detecting cosmic rays with the LOFAR radio telescope
Average Residuals [ns] N u m b e r o f eve n t s Monte Carlo Residuals of Random Pulse [ns] N u m b e r o f eve n t s Fig. 17:
Average residual delays derived from a plane-fit to data ( top ) and from random samples in the searchwindow with respect to a horizontal shower front ( bottom ).The vertical line indicates the cut value derived from thesimulated distribution, which is applied to the data.antennas) is first identified and only its maximum positionsare used for the subsequent direction fit.
As described above, every station is processed separately,meaning that the data do not provide a large lever arm fordirection fitting. However, it also means that the actualshape of the shower front is an insignificant factor in thedirection fitting. For a measurement with a single station,which has a maximum baseline of
80 m , the shower frontcan be approximated by a plane wave. Thus, to determinethe arrival direction of the cosmic ray a planar wavefront isfit to the arrival times of the pulses.This method assumes that essentially all antennas areon a single plane, which certainly holds for all LOFARstations as the ground was flattened during construction.Given a vector of arrival times t , and the vectors x and y for the coordinates of the antennas, the best fitting solutionfor a plane wave: ct = Ax + By + C, (13) can be found using a standard least squares approach. From A and B the Cartesian directions φ, θ can be extracted as: A = sin( θ ) sin( φ ) , (14) B = sin( θ ) cos( φ ) . (15)The plane wave fit itself is done in several iterations. Aftera fit is performed the residual delays are investigated andthose antennas that have residual delays larger than 3 timesthe standard deviation on the residual delays, are removedfrom the set and the data are refitted. The fit is termi-nated when there are less than four antennas left in the setor if no further antennas need to be removed. For this bestdirection all residual delays, including those of removed an-tennas, are calculated again and used for quality cuts later.There are several quality criteria in the pipeline relatedto the plane wave fit. If the fit fails, a station is not consid-ered further. In addition, a cut is made on the remainingaverage residual delays with respect to the expectation ofthe best fit. This cut can be derived from the distributionof all occurring plane wave residual delays, as shown in thetop panel of Fig. 17. The first peak with events of an aver-age residual delay of less than
10 ns corresponds to excellentevents, in which a clear cosmic-ray pulse can be identifiedin all antennas. The largest peak corresponds to all thoseevents in which random noise fluctuations are identified asa pulse. This can be illustrated by a small Monte Carlosimulation. A random sample is picked from the range ofthe search window and its residual to the middle of thesearch window (corresponding to a vertical shower) is cal-culated. This results in the distribution in the bottom panelof Fig. 17. The peak in the distribution obtained from dataand the Monte Carlo distribution are centered around thesame value and can therefore be identified with each other.Second order effects, being the directions of the air showersand non-infinite sampling, can influence the shape of thepeak. The longer tail of the first peak (up to about
50 ns )corresponds to events that have some antennas with cor-rectly identified pulses and varying numbers of outliers, i.e.antennas where a random pulse is identified.Therefore one can safely choose the value
90 ns as a firstcut for good cosmic-ray events. Further cuts for higherquality events or stations can be applied in later analyses.The plane wave fit results now also allow for a justifi-cation of the choice of the Hilbert envelope as the methodfor pulse timing, as opposed to cross correlation. Figure18 shows the ratio of the number of antennas in which apulse has been identified by either method with respectto the remaining residuals on a test-set of randomly cho-sen events that contain a cosmic-ray signal. The distri-bution clearly shows that the Hilbert envelope finds sig-nificantly more signals in the first bin, i.e. in the correctbin with small residuals. In general, cross correlating isexpected to be better for pulses with lower signal-to-noiseratio. For pulses with a high signal-to-noise, however, theHilbert transform performs more accurately. When usingthe Hilbert envelope, the position of the maximum is onlydetermined by the recorded individual pulse, whereas thepeak of the cross correlation is determined by the degree towhich two signals correlate. This degree of correlation canbe influenced by correlations in the noise (for instance resid-ual RFI) or lacking similarity of the pulse shape betweenantennas, thereby making the cross-correlation less accu-rate for timing of pulses with a high signal-to-noise ratio.
Article number, page 11 of 16 esiduals from planewavefit [ns] C r o ss C o rr e l a t i on / N H il b e r t E n ve l op e N Fig. 18:
Difference in reconstruction between Hilbert En-velope and Cross Correlation. The different quality of thereconstruction is illustrated by plotting the fraction of thenumbers of antennas N , identified by each method, withrespect to the residual that was found in the plane wavereconstruction. For values above one the Hilbert Envelopeidentified more antennas, which is the case for the desiredcorrectly identified signals, which can be found below
20 ns . After the antenna pattern unfolding cycle completes with asuccessful direction fit for a given station, the electric fieldcomponents in on-sky polarizations E θ ( t ) , E φ ( t ) , and theshower arrival direction n are known. However, to com-pare measured data to air-shower simulations the three-dimensional electric field at ground level E ( t ) = E x ( t ) ˆe x + E y ( t ) ˆe y + E z ( t ) ˆe z (16)is needed, where ˆe x , ˆe y and ˆe z form the right handed co-ordinate system pointing east, north and up, respectively.This geometry can also be seen in Fig. 13.Assuming the signal has no electric field component inthe propagation direction − ˆe n , this follows from a simplerotation ( E x , E y , E z ) T = R · ( E θ , E φ , T , with the rotationmatrix R = (cid:32) cos θ cos φ − sin φ sin θ cos φ cos θ sin φ cos φ sin θ sin φ − sin θ θ (cid:33) . (17)Note that this assumption is only an approximation,since the signal is measured in the near field of the showerand the source is moving. However, these are second ordereffects and the E θ ( t ) and E φ ( t ) components are expected todominate over the E n ( t ) component Huege (2013). More-over, since LOFAR uses a dual polarization set-up it is notpossible to extract the E n ( t ) component of a linearly po-larized signal.The pipeline concludes by storing pulse parameters foreach antenna in the projected directions. In addition to the shower arrival direction, obtained frompulse timing, two more parameters are extracted: for each
Fig. 19:
Footprint of an air shower measured with LO-FAR. The signal strength (peak amplitude of the radio sig-nal) is encoded logarithmically in the size of the markerand the color shows the time of arrival. The pentagonsrepresent the positions of the particle detectors, their sizeis proportional to the number of recorded particles. Thereconstructed shower axis is indicated by the blue cross forthe core position and the line for the projected arrival di-rection.antenna the peak amplitude and integrated power of thepulse are calculated.Without multiplicative unit conversion factors, ignoredfor current lack of absolute calibration, the integrated pulsepower is defined through the integration of the instanta-neous Poynting vector and the electric field strength as: P = (cid:88) k P k ∝ (cid:88) k (cid:90) ∆ t | E k ( t ) | dt, (18)where k = ( x, y, z ) are the polarization components of theelectric field and ∆ t is taken as a symmetric window aroundthe pulse maximum.This is calculated in discrete sampling x i as P k = 1 f (cid:88) signal | x i | − N signal N noise (cid:88) noise | x i | , (19)where f = 200 MHz is the sampling frequency and N signal and N noise are the number of samples in the signal and noisewindows respectively. The noise window consists of the full block excluding the pulse window. The result of the reconstruction pipeline is a full three-dimensional electric field vector per antenna position as afunction of time. There are various ways in which this resultcan be visualized. The shower footprint, Fig. 19, shows thesignal strength (peak amplitude of the radio signal) at themeasured antenna locations as well as the time of arrival.Here, one can see that both the radio signal strength andthe arrival times are consistent with the air-shower directionand core position as determined by the scintillator array.
Article number, page 12 of 16. Schellart et al.: Detecting cosmic rays with the LOFAR radio telescope
Distance to shower axis [m] S i gn a l [ ADU ] -1 xyz Fig. 20:
Distribution of radio signals (peak amplitudein arbitrary units) with respect to the distance from theshower axis as reconstructed from the scintillator data.Shown are the three components of the reconstructed elec-tric field.
Numer of LORA detectors per event E ve n t s d e t ec t e d / E ve n t s t r i gg e r e d Fig. 21:
Fraction of air showers with a detectable radiosignal over the number of air showers triggered with a scin-tillator signal is plotted against the number of particle de-tectors above threshold in an event. The red straight lineis a fit to the data.Both effects are distinctive properties of radio emission fromair showers and are not produced by RFI.Another common way to visualize the result is in theform of the lateral distribution, shown in Fig. 20. Herethe radio signal strength, in all three polarization compo-nents, is shown as a function of projected distance to theshower axis. This projection retains the spatial distributionof the antennas (i.e. stations can be seen as groups), butazimuthal symmetry in the shower plane is assumed. Thisrather complicated looking distribution can be explainedusing detailed models of the radio emission, which also in-clude non-rotational symmetrical effects. Further details ofevent by event characteristics will be reported in forthcom-ing publications. ] LORA s Angular Difference [0 2 4 6 8 10 12 14 16 18 20 N u m b e r o f eve n t s Fig. 22:
Angular difference between the shower axis recon-structed from the particle data and the direction estimatefrom the radio signal. To make the events comparable, thedifference is scaled with the uncertainty of the individualreconstruction σ LORA .
4. Properties of reconstructed air showers
In order to verify the data quality and the method of re-construction a short overview of the first data taken withLOFAR is given. The data set used here (June 2011 untilApril 2013) contains 3341 recorded triggers, of which 1597pass the strict quality cut for a good data reconstructionof the particle measurement. Of all triggers, 405 eventscontain signals of cosmic rays as identified by the pipeline,with a threshold energy of · eV . On the reconstruction of air showers from the particle dataquality cuts are applied. The reconstruction is consideredreliable, when the reconstructed shower core is containedwithin the array, the shower is not too horizontal ( θ < ◦ )and the reconstructed Molière radius falls in the range of −
100 m . After cuts, the lowest energy of a shower thattriggered a read-out of the LOFAR buffers is . · eV and the highest is . · eV . The LORA scintillator arraybecomes fully efficient above · eV .All triggers sent by the scintillator array follow anearly uniform distribution in azimuth and a sin( θ ) cos( θ ) -distribution in zenith angle as it is expected from the ge-ometry for a horizontal array with flat detectors.The number of events with a detectable radio signalincreases with the number of triggered particle detectors,as can be seen in Fig. 21, where the fraction of triggeredevents, with and without a detected radio signal, is plottedagainst the number of particle detectors per event. Thefraction is clearly increasing with the number of triggereddetectors, as shown by a fitted straight line. According tothis fit, at a threshold of 13 detectors about 10% of theevents contain a cosmic-ray signal. Characteristic transverse size of an air shower.Article number, page 13 of 16 °45°90°135°180° 225° 270° 315°103060
Fig. 23:
Arrival directions of the cosmic-ray events de-tected with LOFAR from June 2011 until April 2013. Eastis ◦ and north corresponds to ◦ . Also indicated (cross)is the direction of the magnetic field at LOFAR. For a first estimate all reconstructed triggers are consid-ered valid events which show radio pulses coming from adirection that agrees to ◦ angular distance with the di-rection that was reconstructed from the arrival times mea-sured with the particle detectors. This choice is based onthe results shown in Fig. 22. This figure shows the an-gular difference between the two reconstructed axes for allevents. A steep fall-off in number of events with an increas-ing angular difference can be seen. Any event that deviatesmore than σ LORA certainly lies outside the correct dis-tribution. The shower axis is on average reconstructed withan uncertainty σ LORA ∼ ◦ from the data of the particledetectors. Thus, a quality cut of ◦ is chosen.Figure 23 shows all 405 cosmic-ray events successfullydetected with the LBAs as distributed on the local sky.Visible is a clear north-south asymmetry, where 276 eventsarrive from the northern hemisphere. This corresponds toa probability p = 0 . ± . for a detected event to arrivefrom the north. As the magnetic field at LOFAR is parallelto the north-south axis this is expected, if the main contri-bution to the signal is of geomagnetic origin (Falcke et al.2005; Ardouin et al. 2009).The effect is also illustrated in Fig. 24, which shows thefraction of detected air showers as a function of azimuth an-gle for the events with radio signal, as well as for all LORAtriggers sent. While the events registered with the LORAdetectors are uniformly distributed in azimuth, the radioevents show a clear deficit from the south. Due to the ori-entation of the LOFAR antennas and thereby the reducedsensitivity for purely east-west polarized signals, events ar-riving directly form the north are not necessarily preferred,as their signal is expected to be mainly polarized in theeast-west direction (Huege 2013). The detection efficiencyas a function of direction follows from a deconvolution ofthe expected emission strength with the antenna patternand will not be discussed in detail here.The energies of the air showers with a detectable radiosignal are shown in Fig. 25. The depicted energy is the one ] (cid:176) Azimuth [0 50 100 150 200 250 300 350 R e l a t i ve f r ac t i on o f eve n t s [ % ] Radio (LOFAR)Particle (LORA)
Fig. 24:
Binned distribution of the azimuth angles of allevents measured with the particle detectors (black squares)and those in coincidence of particle detectors and radio an-tennas (red triangles). The best fit of a straight line to theparticle data is also shown. The fit has a χ / nDoF = 0 . . log(Energy [eV]) 15.5 16 16.5 17 17.5 18 18.5 N u m b e r o f eve n t s Fig. 25:
Distribution of the energies of the cosmic rayswhich had a measurable radio signal in the LOFAR data.The depicted energy is the one reconstructed from the cor-responding particle data. The quality cuts, as described inSect. 4.1, are applied.reconstructed from the corresponding particle data. Thisreconstruction has an overall systematic uncertainty of %and varying event by event uncertainties (Thoudam et al. inprep.). One clearly sees that below ∼ eV the detectionof air showers through their radio signal is not fully efficient,as the strength of the radio signal scales with the energyof the shower. Higher energies in this distribution are con-strained by the steeply falling cosmic-ray energy spectrumand limited size of the detector array, which leads to limitedevent statistics at the highest energies. There are significanthints that showers of higher energies have been measuredwith LOFAR (especially when including the stations out-side the Superterp), but these events are not well enoughconstrained by the data from the particle detectors in orderto have a reference energy of the necessary accuracy. Aftera calibration of the energy of the radio measurements, thoseevents will be used in a radio-stand-alone reconstruction. Article number, page 14 of 16. Schellart et al.: Detecting cosmic rays with the LOFAR radio telescope
5. Conclusions
At LOFAR cosmic-ray induced air showers are regularlymeasured with an array of particle detectors, LORA, anda large array of radio antennas. The cosmic-ray pipeline isroutinely finding their distinctive radio signatures in themeasurements and a full three-dimensional electric fieldvector is reconstructed for every antenna position.A large dataset has been gathered with hundreds ofidentified cosmic-ray events in data from the LBAs. Withup to a thousand antennas per events, these are the firsthighly detailed measurement of the radio signal of air show-ers. These measurements will be used for a detailed charac-terization of the shower shape and will be the benchmarkdata for comparison with models of radio emission in airshowers.
Acknowledgements
The LOFAR cosmic ray key science project very much ac-knowledges the scientific and technical support from AS-TRON, especially in constructing the LORA particle detec-tors. We thank the KASCADE Collaboration for providingthe scintillator detectors. Furthermore, we acknowledge fi-nancial support from the Netherlands Research School forAstronomy (NOVA), the Samenwerkingsverband Noord-Nederland (SNN) and the Foundation for FundamentalResearch on Matter (FOM) as well as support from theNetherlands Organization for Scientific Research (NWO),VENI grant 639-041-130. We acknowledge funding froman Advanced Grant of the European Research Council un-der the European Union’s Seventh Framework Program(FP/2007-2013) / ERC Grant Agreement n. 227610. LO-FAR, the Low Frequency Array designed and constructedby ASTRON, has facilities in several countries, that areowned by various parties (each with their own fundingsources), and that are collectively operated by the Inter-national LOFAR Telescope (ILT) foundation under a jointscientific policy. The authors would like to thank boththe internal and external referees for carefully reading themanuscript. Chiara Ferrari and Giulia Macario acknowl-edge financial support by the "Agence Nationale de laRecherche" through grant ANR-09-JCJC-0001-01.
References
Alexov, A., Schellart, P., ter Veen, S., et al. 2012, in Astronomical So-ciety of the Pacific Conference Series, Vol. 461, Astronomical DataAnalysis Software and Systems XXI, ed. P. Ballester, D. Egret, &N. P. F. Lorente, 283Allan, H. & Jones, J. 1966, Nature, 212, 129Alvarez-Muñiz, J., Carvalho, W. R., & Zas, E. 2012, AstroparticlePhysics, 35, 325Ardouin, D., Belletoile, A., Berat, C., et al. 2009, AstroparticlePhysics, 31, 192Ardouin, D., Belletoile, A., Charrier, D., et al. 2005, InternationalJournal of Modern Physics A, 20, 6869Falcke, H., Apel, W. D., Badea, A. F., et al. 2005, Nature, 435, 313Falcke, H. e. a. 2008, in ICRC Merida, Vol. Rapporteur, 30th Inter-national Cosmic Ray ConferenceHamaker, J. P., Bregman, J. D., & Sault, R. J. 1996, A&AS, 117, 137Huege, T. 2013, ARENA 2012, AIP Conf. Proc. 1535, 121Huege, T., Apel, W. D., Arteaga, J. C., et al. 2012, Nuclear Instru-ments and Methods in Physics Research A, 662, 72Huege, T., Ludwig, M., & James, C. W. 2013, ARENA 2012, AIPConf. Proc. 1535, 128Jelley, J., Fruin, J., Porter, N., et al. 1965, Nature, 205, 327 Jones, R. C. 1941, Journal of the Optical Society of America (1917-1983), 31, 488Kolundzija, B. 2011, in Proceedings of the 5th European Conferenceon Antennas and Propagation (EUCAP)Marin, V. & Revenu, B. 2012, Astroparticle Physics, 35, 733Nehls, S., Hakenjos, A., Arts, M. J., et al. 2008, Nucl. Inst. Meth. A,589, 350Offringa, A. R., de Bruyn, A. G., Zaroubi, S., et al. 2013, A&A, 549,A11Thoudam, S., Buitink, S., Corstanje, A., et al. in prep., Nuclear In-strument and Methods Avan Cappellen, W., Ruiter, M., & Kant, G. 2007, LOFAR-ASTRON-ADD-009van Haarlem, M. P., Wise, M. W., Gunst, A. W., et al. 2013, A&A,556, A2Werner, K., de Vries, K. D., & Scholten, O. 2012, AstroparticlePhysics, 37, 5 Department of Astrophysics/IMAPP, Radboud UniversityNijmegen, P.O. Box 9010, 6500 GL Nijmegen, The Nether-lands Nikhef, Science Park Amsterdam, 1098 XG Amsterdam, TheNetherlands Netherlands Institute for Radio Astronomy (ASTRON),Postbus 2, 7990 AA Dwingeloo, The Netherlands Max-Planck-Institut für Radioastronomie, Auf dem Hügel69, 53121 Bonn, Germany KVI, University Groningen, 9747 AA Groningen, TheNetherlands ECAP, University of Erlangen-Nuremberg, 91058 Erlangen,Germany Astronomical Institute ’Anton Pannekoek’, University ofAmsterdam, Postbus 94249, 1090 GE Amsterdam, TheNetherlands Kapteyn Astronomical Institute, PO Box 800, 9700 AVGroningen, The Netherlands Leiden Observatory, Leiden University, PO Box 9513, 2300RA Leiden, The Netherlands Jodrell Bank Center for Astrophysics, School of Physics andAstronomy, The University of Manchester, Manchester M139PL,UK Astrophysics, University of Oxford, Denys Wilkinson Build-ing, Keble Road, Oxford OX1 3RH School of Physics and Astronomy, University of Southamp-ton, Southampton, SO17 1BJ, UK Max Planck Institute for Astrophysics, Karl SchwarzschildStr. 1, 85741 Garching, Germany International Centre for Radio Astronomy Research - CurtinUniversity, GPO Box U1987, Perth, WA 6845, Australia STFC Rutherford Appleton Laboratory, Harwell Science andInnovation Campus, Didcot OX11 0QX, UK Institute for Astronomy, University of Edinburgh, Royal Ob-servatory of Edinburgh, Blackford Hill, Edinburgh EH9 3HJ,UK LESIA, UMR CNRS 8109, Observatoire de Paris, 92195Meudon, France Argelander-Institut für Astronomie, University of Bonn, Aufdem Hügel 71, 53121, Bonn, Germany Leibniz-Institut für Astrophysik Potsdam (AIP), An derSternwarte 16, 14482 Potsdam, Germany Thüringer Landessternwarte, Sternwarte 5, D-07778 Tauten-burg, Germany Astronomisches Institut der Ruhr-Universität Bochum, Uni-versitaetsstrasse 150, 44780 Bochum, Germany Laboratoire de Physique et Chimie de l’ Environnement et del’ Espace, LPC2E UMR 7328 CNRS, 45071 Orléans Cedex02, France SRON Netherlands Insitute for Space Research, PO Box 800,9700 AV Groningen, The Netherlands Center for Information Technology (CIT), University ofGroningen, The Netherlands Article number, page 15 of 16 Centre de Recherche Astrophysique de Lyon, Observatoire deLyon, 9 av Charles André, 69561 Saint Genis Laval Cedex,France ARC Centre of Excellence for All-sky astrophysics (CAAS-TRO), Sydney Institute of Astronomy, University of SydneyAustralia University of Hamburg, Gojenbergsweg 112, 21029 Hamburg,Germany Astro Space Center of the Lebedev Physical Institute, Prof-soyuznaya str. 84/32, Moscow 117997, Russia Centre for Radio Astronomy Techniques & Technologies(RATT), Department of Physics and Electronics, RhodesUniversity, PO Box 94, Grahamstown 6140, South Africa SKA South Africa, 3rd Floor, The Park, Park Road,Pinelands, 7405, South Africa Harvard-Smithsonian Center for Astrophysics, 60 GardenStreet, Cambridge, MA 02138, USA Laboratoire Lagrange, UMR7293, Universitè de Nice Sophia-Antipolis, CNRS, Observatoire de la Cóte d’Azur, 06300Nice, France Space Telescope Science Institute, 3700 San Martin Drive,Baltimore, MD 21218, USA Sodankylä Geophysical Observatory, University of Oulu,Tähteläntie 62, 99600 Sodankylä, Finland Netherlands eScience Center, Science Park 140, 1098 XG,Amsterdam, The Netherlands36