© ESO, 2013

1. Introduction

With the development of ever faster electronics and the increase in computational power, the construction of radio telescopes as large interferometric arrays of rather simple antennas opens a new window for observations. The low frequency array (LOFAR; van Haarlem et al. 2013), is the first large-scale implementation of this technique. In addition to producing the first high quality images at these low frequencies of 10−240 MHz, LOFAR was designed to study short, pulsed signals in the time-domain. With a vast array of antennas observing the whole sky simultaneously, observations are not limited to a predefined direction, therefore providing optimal conditions for cosmic-ray detection.

Cosmic rays, accelerated charged particles from astrophysical sources, can be observed over several decades of energy. When cosmic rays of high energies reach the Earth, they do not reach the surface as primary particles, but instead interact with atmospheric nuclei. Thereby, a cascade of particles is created, consisting mostly of photons and a significant fraction of charged particles. While propagating through the atmosphere, the charged particles of this extensive air shower emit electromagnetic radiation, which adds up coherently for wavelengths comparable to the dimensions of the shower front (Huege 2013).

Already in the 1960s it was proven that cosmic ray-induced air showers emit nanosecond duration pulses with significant power in the MHz radio frequency range (Jelley et al. 1965; Allan & Jones 1966), but due to lack of sufficiently sophisticated and fast electronics the technique was not pursued further. Only in the past decade interest in the detection technique was rekindled and successfully applied (Falcke et al. 2008). The proof of principle and large progress in the understanding of the emission was made at the LOFAR Prototype Experimental Station (LOPES; Falcke et al. 2005; Huege et al. 2012) and further refined by measurements at the CODALEMA experiment (Ardouin et al. 2005).

Similar to optical measurements of the fluorescence emission from atoms excited by interaction with the air shower, radio emission directly traces the longitudinal shower development, which is closely related to the type of the primary particle. Unlike optical fluorescence measurements, radio emission measurements are less dependent on observing conditions and can operate day and night matching the duty cycle of particle detector measurements.

Due to the very steep energy spectrum, measuring the highest-energy cosmic rays requires vast detector areas. Cost constraints therefore limit the density of detectors within this area giving a wide spacing between the individual antennas. Theoretical models describing the different emission mechanisms at play point to a very detailed and non-symmetrical emission pattern at ground level (Werner et al. 2012; Alvarez-Muñiz et al. 2012; Marin & Revenu 2012; Huege et al. 2013). Testing these models therefore requires dense sampling of the electric field over a sufficiently large area.

LOFAR offers a high number of antennas clustered on an irregular grid, with increasingly large spacing between antenna clusters further away from the center. In the core of the array about 2300 antennas are installed within about 4km², which allows air showers to be measured with unprecedented spatial resolution. These measurements will contribute significantly to conclusively confirm theoretical models for the radio emission on a shower by shower basis, a goal previously unattainable due to lack of sufficiently high quality data.

Measurements and converging theoretical predictions of the expected radio signal from a cosmic-ray induced air shower give a short, nanosecond time-scale bi-polar pulse, which is mostly linearly polarized. This article describes the detection set-up and automated processing pipeline used at LOFAR to measure and identify these signals.

Starting with a description of the instrumental set-up at LOFAR in Sect. 2, an overview of the data reduction pipeline is given in Sect. 3. Finally, Sect. 4 describes the characteristics of the dataset obtained between June 2011 and April 2013.

The LOFAR dataset will be used in forthcoming publications to verify existing models for radio emission from air showers and to develop new techniques that use radio emission to measure important characteristics of the incoming particle, such as energy and mass.

2. LOFAR

LOFAR is a distributed radio telescope. Its antennas are distributed over northern Europe with the densest concentration in the north of the Netherlands, in the Province of Drenthe. The observation support center and processing facilities are also located near this central core. The antennas of LOFAR are grouped into stations, each station taking the role of a single dish in a traditional radio interferometer array. A station consist of a number of low-band antennas (LBAs, 10−90 MHz) and high-band antennas (HBAs, 110−240 MHz). The 24 stations within the ~2 km wide core are distributed in an irregular pattern that maximizes uv-coverage, or spatial frequencies for standard interferometric observations. The 16 additional Dutch remote stations are distributed with increasing distance to the core. International stations are currently located in Germany, France, the United Kingdom, and Sweden, giving LOFAR a maximum baseline of 1292km for interferometric observations. Core stations and remote stations consist of 96 LBAs plus 48 HBAs. International stations have 96 LBAs and 96 HBAs. At the center of the LOFAR core six stations are located in a roughly 320m diameter area, called the Superterp, providing both the shortest baselines for interferometric observations and the densest population of antennas ideal for cosmic-ray observations. While every LOFAR station is equipped with the necessary electronics to observe cosmic rays, the current data set is taken with the central 24 stations, where additional information from particle detectors is available (see Sect. 2.3). The positions of the antennas of the seven most central LOFAR stations are shown in Fig. 1.

2.1. The antennas

The LBAs are the main tool for cosmic-ray detection. An LBA consists of two orthogonal inverted V-shaped dipoles, each with a length of 1.38 m. These are supported by a central PVC¹ pole, which holds the low-noise amplifier and guides the signal cables, as shown in Fig. 2. The dipoles X and Y, that make up each antenna, are oriented northeast to southwest (NE-SW) and northwest to southeast (NW-SE), as can be seen in Fig. 3.

The low-noise amplifier has an intentional impedance mismatch with the antenna. This mismatch, combined with the characteristic length of the dipoles, makes the system sensitive in a broad band from 10−90 MHz. In principle, this allows observations from the ionospheric cutoff up to the start of the commercial FM radio band. For most observations the frequency range is limited by a combination of selectable hardware and software filters to 30−80 MHz to suppress strong Radio Frequency Interference (RFI) in the outer bands. The LBAs are designed to be sky noise limited after RFI has been removed (van Cappellen et al. 2007). After amplification the signals from the individual dipoles are transmitted through coaxial cables to the electronics cabinet located at every station.

The HBAs have been optimized for a frequency band of 110−240 MHz. The design clusters 16 antenna elements into a tile, the signals from these elements are amplified and combined in an analog beam-former. This means that while the LBAs are sensitive to the whole sky the HBAs are most sensitive within the ~20° of the tile-beam, of which the direction is chosen at the start of every observation. This results in a smaller effective area for cosmic-ray observations, as the measurement will only be optimal if the direction of the cosmic ray happens to coincide with the beam direction of the observation. Therefore, the analysis of HBA data and their interesting higher frequency range requires a different approach for cosmic-ray studies. Results of these measurements will be described in a later publication.

2.2. The transient buffer boards

After being forwarded to the electronics cabinet the signals of the LBAs are again amplified, filtered, and digitized by a 12bit A/D converter with a sampling frequency of 200MHz². Due to signal path limitations in the Dutch stations only 48 dual-polarized or 96 single-polarized antennas can be processed at a given time. For the dual-polarized option the antennas are grouped into an inner and an outer set, 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 addition, there is the possibility to store a snapshot of the original data. Every station is equipped with ring-buffers, the so called Transient Buffer Boards (TBBs). These continuously store the last 1.3 s of raw data (an extension to 5s is currently being deployed). When triggered, the contents of the TBBs are frozen, read out via the Wide Area Network and stored on disk for further analysis. The trigger can be generated based on various parameters in an FPGA³ at the local receiver unit. Alternatively, the trigger can be generated by an array of particle detectors (see Sect. 2.3) or received from outside of LOFAR. Currently, the main trigger for cosmic-ray observation is provided by the particle detectors. Later, a radio self-trigger will be implemented, using the current dataset as a training set to deduce trigger criteria, so that the FPGA trigger can be run independently at every LOFAR station. These criteria have to reduce false triggers to limit the data rate. Using every LOFAR station individually will dramatically increase the effective area.

Essential for measuring cosmic rays with LOFAR as a radio telescope is that the whole process of triggering and storing radio-pulse data can take place without interfering with the ongoing observations.

2.3. The LOFAR Radboud Air Shower Array

LORA, the LOFAR Radboud Air Shower Array, is an array of particle detectors co-located with the center of LOFAR. The array provides a reconstruction of basic parameters of recorded air showers, such as the direction and the position of impact, as well as the energy of the incoming cosmic ray (Thoudam et al., in prep.). It also provides the time of arrival, which is used to trigger the read-out of the radio antennas.

LORA consists of 20 detector units distributed on the Superterp, as shown in Fig. 1. Each detector contains two scintillators (0.45 m², type: NE 114), which are individually read out through a photomultiplier tube. The detectors are inside weatherproof shelters and have been tested to not create any interference at radio frequencies.

Conditions at which triggers are sent to LOFAR can be adjusted to match the desired energy threshold. There are two constraints on the desired rate: the rate of events interesting for radio observations has to be maximized, while the network load on the LOFAR system has to be kept low in order to avoid interfering with the primary observation. A trigger in a single detector is generated when a particle signal of more than 4σ above the noise is registered. In order to only detect air showers a coincidence of several detectors is needed. Events of less than 10¹⁶ eV have a very low probability to be observable in radio above the sky-noise level. The energy threshold and the corresponding event rate are shown in Fig. 4 as the function of the number of triggered detectors. Requiring triggers in 13 detectors yields a threshold energy of 2.4 × 10¹⁶ eV, with an average trigger-rate of 0.8 events/h. This trigger rate has been selected as the optimal setting for the observations.

2.4. Observations

After the commissioning phase LOFAR is to be used on a proposal-based schedule. Proposals are open to the community for imaging or beamformed observations, as well as TBB observations. Some fraction of the observing time is reserved for participating consortia and key science 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 be run in the background of all other observations that do not need the full network bandwidth. This does however mean that the array configuration is determined by the primary observation, therefore the amount of data in a specific array configuration (such as the selection of LBA or HBA antenna type) available for analysis is not determined by the cosmic-ray project itself, except when LOFAR is otherwise idle and the observing configuration can be chosen freely.

During the observation, triggers from LORA are received by the LOFAR control system. The system checks whether a dump from the TBBs is allowed. If so, the ring-buffers are frozen and a specified block of data around the trigger time is dumped to disk. For each cosmic-ray event 2.1 ms of radio data are stored, which corresponds to 77MB per station. This provides sufficient frequency resolution for high quality RFI cleaning while minimizing data transfer and storage requirements.

Every evening, the data-files are archived at LOFAR and compressed for transport. They are stored in the Long Term Archive (van Haarlem et al. 2013), from where they can be retrieved for data analysis.

3. Reconstruction of cosmic-ray data

All newly recorded data are processed every evening, after having been copied via the network to the processing cluster of the Astrophysics department at the Radboud University Nijmegen. In addition to the hdf5⁴ files, containing the data of one LOFAR station each, the recorded data from the particle detectors and a trigger log file are transferred. With this information an automated pipeline is run. The pipeline is based on the task oriented PyCRTools framework consisting of fast low-level C++ routines embedded in Python for maximum flexibility. All results are stored in a PostgreSQL database for subsequent data mining analysis. The goal of the processing pipeline is to autonomously identify a full set of physics quantities for each air shower detected with LOFAR. The pipeline is optimized to identify those nanosecond pulses that are not generated by terrestrial sources.

All data are first processed per station, i.e. per file. The set of files received for a single trigger form an event. When the data from one station pass the criteria for containing a cosmic-ray signal (see Sect. 3.3), the corresponding event is called a cosmic-ray event. It is not necessary to observe a pulse in all stations, only the stations with a significant signal are used in a combined analysis.

3.1. Pipeline structure

The reconstruction pipeline comprises a number of steps that will be individually explained in the following sections. An overview of the steps and the overall structure is depicted in Fig. 5.

3.2. Preparing the data

Before proceeding to extract the cosmic-ray signal from the data, some preparatory steps have to be performed. Knowledge about the system is applied in the form of calibration procedures, the data are cleaned of narrowband-transmitters, and antennas that show malfunctions are flagged.

3.2.1. Timing offsets and phase calibration

There are known signal path differences between the LOFAR antennas. Measured differences of cable lengths between the antennas are corrected for up to the 5ns sample level already at the stations before the data are written to disk. Additionally, relative time offsets between the antennas are corrected for at sub-sample accuracy using standard LOFAR calibration tables. These tables are generated by phase-calibrating on the strongest astronomical radio sources and are regularly tested and updated if necessary (van Haarlem et al. 2013). Sub-sample corrections are applied as phase offsets to the Fourier transformed signal in the cosmic-ray pipeline, before processing it in the data analysis.

3.2.2. RFI cleaning

Narrow-band RFI in the time series signal can be revealed by making an average power spectrum. An example is shown in the top panel of Fig. 6, where most of the strong RFI is visible outside the 30−80 MHz range. The average power spectrum is created by averaging the square of the absolute value of the Fourier transform over several blocks of data. The block size can be freely chosen within the full data length to obtain a desired frequency resolution; here 2¹⁶ samples are used, giving a resolution of ~3kHz, enough to resolve most RFI lines. A reasonable data length is needed for this procedure to produce a stable average, which sets the limit for the chosen block length to be stored from the TBBs, as mentioned in Sect. 2.4. In order to minimize artificial side lobes a half-Hann window is applied to the first and last 10% of each trace prior to the Fourier transformation.

The standard approach to RFI cleaning (or RFI flagging) is to identify peaks sticking out significantly above the overall spectral shape, also called the baseline, and set the corresponding Fourier component amplitudes to zero. However, this requires a priori knowledge of the baseline. While the baseline can be obtained through a smoothing or fitting procedure, this is often not stable in the presence of strong RFI, requiring an iterative approach.

An alternative approach to RFI cleaning uses the phase information in the complex-valued spectrum instead. If an RFI transmitter is measured in all antennas, the phase difference, or relative phase, between each pair of antennas will be a constant value as function of time with a small non constant random noise contribution. Note that the exact value of the constant, which only depends on the geometric delay between antennas, is not relevant, only its non time-varying nature. When no transmitter is present, the relative phase is expected to be both random and time varying, as the signal then consists of the added signals from many incoherent sources on the sky with additional random noise. Therefore, RFI can be identified by looking at the stability of phase differences between antennas over time. For each antenna-dipole j = 0, 1, ..., 95 in a station and data block k, the phase spectrum is calculated as (1)where x_j,k(ω) is the complex frequency component ω of the spectrum.

Subtracting the phase of one of the antennas as reference antenna gives the relative phases and results in a set of phases for every frequency channel, one for each block of data. Only one reference antenna is used and this is taken to be the one with median power to avoid selecting a broken antenna.

The average phase is defined as (2)and the phase variance as (3)where N is the number of data blocks.

For completely random phases one expects s_j(ω) ≈ 1 as opposed to s_j(ω) = 0 when all phases are equal. The phase variance per frequency channel will now either be at a value close to 1, including some random noise, or at a significantly lower level. The latter reveals the presence of a radio transmitter, as shown in Fig. 7, where a contaminated part of the spectrum is shown with the corresponding phase variance.

Since RFI lines will result in peaks toward smaller values of the phase variance, and noise has no preferred peak direction, calculating the standard deviation σ in this plot only for values above the median will ensure a stable result. All frequencies that have a phase variance of at least 5σ below the median⁵ are flagged as containing RFI. Additionally a 30−80 MHz bandpass filter is applied, flagging the most heavily RFI polluted low and high frequency parts of the bandwidth by default. To prevent pulse-ringing the 30−80 MHz block filter is first convolved with a, σ_tapering = 2.5 MHz, Gaussian⁶. After removing the flagged channels, the resulting cleaned spectrum is shown in the bottom panel of Fig. 6.

In general, there is very little RFI at the LOFAR Superterp. A lot of effort has been made to remove local sources that could disturb the LOFAR measurements and a protected zone has even been established (Offringa et al. 2013). This relative quietness is illustrated in Fig. 8. It shows the result of the RFI cleaning for all events for all frequencies. While every event has some RFI, no single RFI line is present in every event. Within the 30−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 more than 2% flagged channels out of the more than 32 000 frequency channels in a block of data. This is is shown in Fig. 9, where the total fraction of events is plotted against the number of flagged channels.

3.2.3. Flagging bad antennas

Occasionally, one or more antennas give invalid signals, e.g. due to hardware malfunction. To identify these bad antennas the integrated spectral power is calculated (4)where x(ω) is the ω frequency component of the cleaned spectrum. The power in every antenna is required to be in the range of one half to two times the median power from all antennas. Antennas outside this range are marked as bad and excluded from further analysis.

3.2.4. Absolute gain calibration

There are ongoing efforts for an absolute calibration of the voltage traces of LOFAR and therefore the reconstructed electric field. Those efforts will be described in a forthcoming publication and include calibration on astronomical sources, terrestrial transmitters, and already conducted dedicated measurement campaigns, similar to those performed at other experiments, e.g. Nehls et al. (2008). Once implemented, the reconstruction pipeline will deliver calibrated electric field strengths and their polarization components for all events. However, significant progress in understanding the mechanisms of radio emission in air showers can already be made with a relative calibration.

3.2.5. Relative gain calibration

The LBA measurement is dominated by sky noise, which in turn is dominated by the Galaxy moving through the antenna beam pattern. Therefore, the noise as seen by each antenna is a function of the Local Sidereal Time (LST) and can be used to correct for differences in gain between antennas. Instead of correcting all antennas at all times to a fixed value, which would be over- or underestimating the noise at certain times, the received power can be normalized to a LST-dependent reference value. In Fig. 10 the integrated spectral power (Eq. (4)), after RFI cleaning, is given as a function of LST for the instrumental polarization X and Y. The data have been retrieved from all cosmic-ray events measured within the first year of data-taking. One can define a reference value for the integrated spectral power as a function of LST by fitting a function to these data points. Since the movement of the Galaxy through the antenna beam pattern is periodic by nature it is fitted with the 2nd order Fourier series (5)thereby avoiding artificial jumps in the fit at 0:00 LST. The time t is given in units of radian here. This results in a gain correction for each antenna as (6)where the square root is needed, because the correction is applied to the amplitude spectrum.

3.3. Identifying cosmic-ray signals

After cleaning and calibration of the data, the central element of the pipeline is the identification and characterization of the radio pulse as the signal of the air shower.

3.3.1. Using information from the particle detectors

In order to restrict the search for the radio pulse to a smaller region in the trace, the information from the trigger time of the particle detectors is used. Figure 11 shows the difference in time between the trigger from the particle detectors and the pulse location in the radio data obtained from a search with a large window. The distribution shows a clear peak at the region of the coincidences at an offset of 253 ± 168 ns. In absolute timing the offset between LORA and LOFAR is 10 253 ns, of which 10 000 ns are already accounted for in the triggering system.

Average offset is obtained by fitting a Gaussian to the distribution of pulse positions with respect to the trigger time. This is only an approximation, as the real offset per event depends on the position of the core and the incoming direction of the air shower. Also, effects due to the propagation of particles and radiation in the atmosphere can play a role. The overall difference is due to the fact that both detectors operate independently on different timing systems. Both are based on GPS timing, but correct for drifts (<20ns) in different ways and have a differing absolute time. The spread on the differences is however sufficiently small for Superterp stations to not require additional synchronization of the two systems. Stations further away can have larger offsets due to the signal travel time, which can be corrected for after a reconstruction of the shower.

These measurements allow for the pulse search to be restricted to a small fraction of the full time trace, limiting the chance to pick up random noise fluctuations.

3.3.2. Finding candidate events

The trigger threshold of the scintillator array is chosen to be lower than the threshold to detect a radio signal. This ensures a full sample, but also makes it necessary to identify in a first quality check whether there is a detectable signal present. Therefore, per antenna polarization, the signals are first beamformed in the direction reconstructed from the data of the particle detectors. This direction is given in the local Cartesian coordinate frame of the station by n and the position of each antenna j is given by r_j. A planar wavefront arriving at the phase center (0,0,0) at time t = 0 will arrive at antenna j with a delay given by (7)where c is the speed of light. The beamformed signal, in frequency space, in this direction is then given by (8)where x_j(ω) is frequency component ω of the Fourier transform of the signal from antenna j and N_a is the number of antennas. The inverse Fourier transform gives the beamformed time series signal. Due to beamforming any signal coming from the direction of the air shower is amplified by a factor N_a in amplitude while uncorrelated noise is only amplified by a factor . Therefore, if no significant signal is detected in the beamformed trace, the event very unlikely contains a cosmic-ray signal strong enough to be detected at single dipole level by the rest of the pipeline. Thus, the analysis of the data of that station is aborted.

To test this assumption, Fig. 12 shows the distribution of the peak amplitude in the beamformed signal per station, distinguishing between events in which ultimately a cosmic-ray was identified and those in which there was not. The peak amplitude is normalized by the root mean square of the trace, as a proxy for the noise contribution. From this it can be seen that the fraction of events where a strong signal is observed in the beamformed trace is significantly higher for stations where eventually a cosmic-ray signal is detected. All events in the tail of the non-detected distribution were visually inspected and identified as broad-band RFI, with pulses differing significantly in shape from those of cosmic rays and directions ultimately deviating significantly from the direction as measured with the particle detectors. This distribution shows that an initial filtering based on a moderate signal-to-noise of beamformed pulses is a quick and effective way to filter out those events that are potentially interesting, as well as further narrowing the search window per antenna reducing false positives for pulse detection.

3.3.3. Correction for the antenna response

The sensitivity of the LOFAR LBA is a complex function of both frequency and direction. Correcting for this antenna pattern, i.e. unfolding, requires an initial guess for the pulse direction and in turn may influence the position of the pulse in time and thus the direction by changing the phase at which each frequency arrives. Therefore the correction has to be done in an iterative loop as indicated in Fig. 5. Each iteration starts with an increasingly accurate signal direction and proceeds by unfolding the antenna pattern, pulse detection, and direction fitting. The loop is concluded when the direction no longer significantly changes, which usually happens in less than ~5 iterations.

For the antenna pattern of the LBA a simulation is used, which is made using the software WIPL-D (Kolundzija 2011) and a customized software model of the electronics chain.

From the impedance and radiation pattern in a transmitting situation the open circuit voltage is calculated as a function of frequency and direction for an incoming plane wave with an electric field strength of 1V/m. The equivalent circuit of the antenna in a receiving situation is a voltage source with an internal resistance equal to the antenna impedance. This is combined with measured data of the amplifier directly behind the antenna. The result of the model is the output voltage of the amplifier over a 75Ω resistor⁷.

Any wave coming from a direction ê_n can be seen as a linear superposition of monochromatic plane waves, polarized in the ê_φ and ê_θ direction. Here φ and θ are the standard spherical coordinate angles with the x and z axis respectively, e.g. (9)This geometry can be seen in Fig. 13.

These terms are related to the output voltage of the amplifier for each dipole, and for each frequency, via the Jones matrix J (Jones 1941; Hamaker et al. 1996) of the antenna model (10)where J_Xθ is the complex response of the antenna and amplifier of the X-dipole to a wave purely polarized in the ê_θ direction.

Therefore, in order to both correct for the antenna response and convert from output voltage to electric field strength in the on-sky frame (see Sect. 3.4), each pair of Fourier components from the signal in the two instrumental polarizations (X,Y) is multiplied by the inverse Jones matrix, followed by an inverse Fourier transform back to the time domain.

The components of the Jones matrix of the antenna model are simulated on a grid with steps of 1MHz in frequency, 5° in θ and 10° in φ. In order to obtain the components at the frequency and direction of observation, trilinear interpolation is performed on the real and imaginary parts of the complex table when needed. Examples of the response are depicted as a function of frequency in Fig. 14 and as a function of direction in Fig. 15.

3.3.4. Pulse detection

Estimating the direction of the incoming air shower, see Sect. 3.3.5, can either be done using beamforming or through pulse timing. Beamforming was found to be very sensitive to the optimization algorithm used, essentially requiring a grid search to avoid getting stuck in a local minimum. This is computationally very expensive, moreover it only provides relative time differences between any two antennas rather than an absolute time needed for extraction of relevant physical parameters (see Sect. 3.4).

In order to use pulse timing, individual pulses have to be identified. This can be done by using the cross-correlation method, where one looks for the maximum in the cross correlation of the signals between all antennas. This however has the same drawback as beamforming, as only relative timing is calculated. A method to retrieve the absolute pulse timing is through the use of the Hilbert envelope, which is used in this pipeline. A detailed comparison of the methods is given in Sect. 3.3.5.

A sensible definition of the pulse arrival time is the measured arrival time of the maximum of the electric field strength. In practice, however, using directly max(| x(t)|²) is highly dependent on the filter characteristics of the receiving system and the sampling used. Therefore, the arrival time is defined as the position of the maximum in the amplitude envelope of the analytic signal, also called the Hilbert envelope (11)where is the Hilbert transform, or imaginary propagation, of the signal x(t) defined by (12)where ℱ denotes the Fourier transform.

In order to find the pulse maximum with subsample precision, the signal is first up-sampled by a factor 16, such that the maximum search will not be the limiting factor in the timing resolution. Subsequently, a simple maximum search is performed on the envelope. In addition, the signal-to-noise ratio is calculated where the signal is defined as the maximum and the noise as the root mean square of the envelope. 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 pulse with a signal to noise greater than three is marked as a possible cosmic-ray signal to be used for direction fitting. Because the pulse is expected to be intrinsically stronger in one of the two polarizations, depending on the angle between the shower axis and the geomagnetic field, the polarization with the highest average signal to noise (over all antennas) is first identified and only its maximum positions are used for the subsequent direction fit.

3.3.5. Arrival direction fitting

As described above, every station is processed separately, meaning that the data do not provide a large lever arm for direction fitting. However, it also means that the actual shape of the shower front is an insignificant factor in the direction fitting. For a measurement with a single station, which has a maximum baseline of 80m, the shower front can be approximated by a plane wave. Thus, to determine the arrival direction of the cosmic ray a planar wavefront is fit to the arrival times of the pulses.

This method assumes that essentially all antennas are on a single plane, which certainly holds for all LOFAR stations 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 solution for a plane wave: (13)can be found using a standard least squares approach. From A and B the Cartesian directions φ,θ can be extracted as: The plane wave fit itself is done in several iterations. After a fit is performed the residual delays are investigated and those antennas that have residual delays larger than 3 times the standard deviation on the residual delays, are removed from the set and the data are refitted. The fit is terminated when there are less than four antennas left in the set or if no further antennas need to be removed. For this best direction all residual delays, including those of removed antennas, are calculated again and used for quality cuts later.

There are several quality criteria in the pipeline related to the plane wave fit. If the fit fails, a station is not considered further. In addition, a cut is made on the remaining average residual delays with respect to the expectation of the best fit. This cut can be derived from the distribution of all occurring plane wave residual delays, as shown in the top panel of Fig. 17. The first peak with events of an average residual delay of less than 10ns corresponds to excellent events, in which a clear cosmic-ray pulse can be identified in all antennas. The largest peak corresponds to all those events in which random noise fluctuations are identified as a pulse. This can be illustrated by a small Monte Carlo simulation. A random sample is picked from the range of the search window and its residual to the middle of the search window (corresponding to a vertical shower) is calculated. This results in the distribution in the bottom panel of Fig. 17. The peak in the distribution obtained from data and the Monte Carlo distribution are centered around the same value and can therefore be identified with each other. Second order effects, being the directions of the air showers and non-infinite sampling, can influence the shape of the peak. The longer tail of the first peak (up to about 50ns) corresponds to events that have some antennas with correctly identified pulses and varying numbers of outliers, i.e. antennas where a random pulse is identified.

Therefore one can safely choose the value 90ns as a first cut for good cosmic-ray events. Further cuts for higher quality events or stations can be applied in later analyses.

The plane wave fit results now also allow for a justification of the choice of the Hilbert envelope as the method for pulse timing, as opposed to cross correlation. Figure 18 shows the ratio of the number of antennas in which a pulse has been identified by either method with respect to the remaining residuals on a test-set of randomly chosen events that contain a cosmic-ray signal. The distribution clearly shows that the Hilbert envelope finds significantly more signals in the first bin, i.e. in the correct bin with small residuals. In general, cross correlating is expected to be better for pulses with lower signal-to-noise ratio. For pulses with a high signal-to-noise, however, the Hilbert transform performs more accurately. When using the Hilbert envelope, the position of the maximum is only determined by the recorded individual pulse, whereas the peak of the cross correlation is determined by the degree to which two signals correlate. This degree of correlation can be influenced by correlations in the noise (for instance residual RFI) or lacking similarity of the pulse shape between antennas, thereby making the cross-correlation less accurate for timing of pulses with a high signal-to-noise ratio.

3.4. Coordinate transformation

After the antenna pattern unfolding cycle completes with a successful direction fit for a given station, the electric field components in on-sky polarizations E_θ(t), E_φ(t), and the shower arrival direction n are known. However, to compare measured data to air-shower simulations the three-dimensional electric field at ground level (16)is needed, where ê_x, ê_y and ê_z form the right handed coordinate system pointing east, north and up, respectively. This geometry can also be seen in Fig. 13.

Assuming the signal has no electric field component in the propagation direction −ê_n, this follows from a simple rotation (E_x,E_y,E_z)^T = R·(E_θ,E_φ,0)^T, with the rotation matrix (17)Note that this assumption is only an approximation, since the signal is measured in the near field of the shower and the source is moving. However, these are second order effects and the E_θ(t) and E_φ(t) components are expected to dominate over the E_n(t) component (Huege 2013). Moreover, since LOFAR uses a dual polarization set-up it is not possible to extract the E_n(t) component of a linearly polarized signal.

The pipeline concludes by storing pulse parameters for each antenna in the projected directions.

3.5. Extracting pulse parameters

In addition to the shower arrival direction, obtained from pulse timing, two more parameters are extracted: for each antenna the peak amplitude and integrated power of the pulse are calculated.

Without multiplicative unit conversion factors, ignored for current lack of absolute calibration, the integrated pulse power is defined through the integration of the instantaneous Poynting vector and the electric field strength as: (18)where k = (x,y,z) are the polarization components of the electric field and Δt is taken as a symmetric window around the pulse maximum.

This is calculated in discrete sampling x_i as (19)where f = 200MHz is the sampling frequency and N_signal and N_noise are the number of samples in the signal and noise windows respectively. The noise window consists of the full 327 680 ns block excluding the pulse window.

3.6. A measured air shower

The result of the reconstruction pipeline is a full three-dimensional electric field vector per antenna position as a function of time. There are various ways in which this result can be visualized. The shower footprint, Fig. 19, shows the signal strength (peak amplitude of the radio signal) at the measured antenna locations as well as the time of arrival. Here, one can see that both the radio signal strength and the arrival times are consistent with the air-shower direction and core position as determined by the scintillator array. Both effects are distinctive properties of radio emission from air showers and are not produced by RFI.

Another common way to visualize the result is in the form of the lateral distribution, shown in Fig. 20. Here the radio signal strength, in all three polarization components, is shown as a function of projected distance to the shower axis. This projection retains the spatial distribution of the antennas (i.e. stations can be seen as groups), but azimuthal symmetry in the shower plane is assumed. This rather complicated looking distribution can be explained using detailed models of the radio emission, which also include non-rotational symmetrical effects. Further details of event by event characteristics will be reported in forthcoming publications.

4. Properties of reconstructed air showers

In order to verify the data quality and the method of reconstruction a short overview of the first data taken with LOFAR is given. The data set used here (June 2011 until April 2013) contains 3341 recorded triggers, of which 1597 pass the strict quality cut for a good data reconstruction of the particle measurement. Of all triggers, 405 events contain signals of cosmic rays as identified by the pipeline, with a threshold energy of 5 × 10¹⁵ eV.

4.1. Triggers from the array of particle detectors

On the reconstruction of air showers from the particle data quality cuts are applied. The reconstruction is considered reliable, when the reconstructed shower core is contained within the array, the shower is not too horizontal (θ < 50°) and the reconstructed Molière radius⁸ falls in the range of 20−100 m. After cuts, the lowest energy of a shower that triggered a read-out of the LOFAR buffers is 1.8 × 10¹⁵ eV and the highest is 1.9 × 10¹⁸ eV. The LORA scintillator array becomes fully efficient above 2 × 10¹⁶ eV.

All triggers sent by the scintillator array follow a nearly uniform distribution in azimuth and a sin(θ)cos(θ)-distribution in zenith angle as it is expected from the geometry for a horizontal array with flat detectors.

The number of events with a detectable radio signal increases with the number of triggered particle detectors, as can be seen in Fig. 21, where the fraction of triggered events, with and without a detected radio signal, is plotted against the number of particle detectors per event. The fraction is clearly increasing with the number of triggered detectors, as shown by a fitted straight line. According to this fit, at a threshold of 13 detectors about 10% of the events contain a cosmic-ray signal.

4.2. Event rates and sensitivity

For a first estimate all reconstructed triggers are considered valid events which show radio pulses coming from a direction that agrees to 10° angular distance with the direction that was reconstructed from the arrival times measured with the particle detectors. This choice is based on the results shown in Fig. 22. This figure shows the angular difference between the two reconstructed axes for all events. A steep fall-off in number of events with an increasing angular difference can be seen. Any event that deviates more than 10σ_LORA certainly lies outside the correct distribution. The shower axis is on average reconstructed with an uncertainty σ_LORA ~ 1° from the data of the particle detectors. Thus, a quality cut of 10° is chosen.

Figure 23 shows all 405 cosmic-ray events successfully detected with the LBAs as distributed on the local sky. Visible is a clear north-south asymmetry, where 276 events arrive from the northern hemisphere. This corresponds to a probability p = 0.68 ± 0.02 for a detected event to arrive from the north. As the magnetic field at LOFAR is parallel to the north-south axis this is expected, if the main contribution 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 the fraction of detected air showers as a function of azimuth angle for the events with radio signal, as well as for all LORA triggers sent. While the events registered with the LORA detectors are uniformly distributed in azimuth, the radio events show a clear deficit from the south. Due to the orientation of the LOFAR antennas and thereby the reduced sensitivity for purely east-west polarized signals, events arriving directly form the north are not necessarily preferred, as their signal is expected to be mainly polarized in the east-west direction (Huege 2013). The detection efficiency as a function of direction follows from a deconvolution of the expected emission strength with the antenna pattern and will not be discussed in detail here.

The energies of the air showers with a detectable radio signal are shown in Fig. 25. The depicted energy is the one reconstructed from the corresponding particle data. This reconstruction has an overall systematic uncertainty of 27% and varying event by event uncertainties (Thoudam et al., in prep.). One clearly sees that below ~10¹⁷ eV the detection of air showers through their radio signal is not fully efficient, as the strength of the radio signal scales with the energy of the shower. Higher energies in this distribution are constrained by the steeply falling cosmic-ray energy spectrum and limited size of the detector array, which leads to limited event statistics at the highest energies. There are significant hints that showers of higher energies have been measured with LOFAR (especially when including the stations outside the Superterp), but these events are not well enough constrained by the data from the particle detectors in order to have a reference energy of the necessary accuracy. After a calibration of the energy of the radio measurements, those events will be used in a radio-stand-alone reconstruction.

5. Conclusions

At LOFAR cosmic-ray induced air showers are regularly measured with an array of particle detectors, LORA, and a large array of radio antennas. The cosmic-ray pipeline is routinely finding their distinctive radio signatures in the measurements and a full three-dimensional electric field vector is reconstructed for every antenna position.

A large dataset has been gathered with hundreds of identified cosmic-ray events in data from the LBAs. With up to a thousand antennas per events, these are the first highly detailed measurement of the radio signal of air showers. These measurements will be used for a detailed characterization of the shower shape and will be the benchmark data for comparison with models of radio emission in air showers.

---

Acknowledgments

The LOFAR cosmic ray key science project very much acknowledges the scientific and technical support from ASTRON, especially in constructing the LORA particle detectors. We thank the KASCADE Collaboration for providing the scintillator detectors. Furthermore, we acknowledge financial support from the Netherlands Research School for Astronomy (NOVA), the Samenwerkingsverband Noord-Nederland (SNN) and the Foundation for Fundamental Research on Matter (FOM) as well as support from the Netherlands Organization for Scientific Research (NWO), VENI grant 639-041-130. We acknowledge funding from an Advanced Grant of the European Research Council under the European Union’s Seventh Framework Program (FP/2007-2013) / ERC Grant Agreement n. 227610. LOFAR, the Low Frequency Array designed and constructed by ASTRON, has facilities in several countries, that are owned by various parties (each with their own funding sources), and that are collectively operated by the International LOFAR Telescope (ILT) foundation under a joint scientific policy. The authors would like to thank both the internal and external referees for carefully reading the manuscript. Chiara Ferrari and Giulia Macario acknowledge financial support by the “Agence Nationale de la Recherche” through grant ANR-09-JCJC-0001-01.

References

All Figures