Polarized radio emission from extensive air showers measured with LOFAR

We present LOFAR measurements of radio emission from extensive air showers. We find that this emission is strongly polarized, with a median degree of polarization of nearly $99\%$, and that the angle between the polarization direction of the electric field and the Lorentz force acting on the particles, depends on the observer location in the shower plane. This can be understood as a superposition of the radially polarized charge-excess emission mechanism, first proposed by Askaryan and the geomagnetic emission mechanism proposed by Kahn and Lerche. We calculate the relative strengths of both contributions, as quantified by the charge-excess fraction, for $163$ individual air showers. We find that the measured charge-excess fraction is higher for air showers arriving from closer to the zenith. Furthermore, the measured charge-excess fraction also increases with increasing observer distance from the air shower symmetry axis. The measured values range from $(3.3\pm 1.0)\%$ for very inclined air showers at $25\, \mathrm{m}$ to $(20.3\pm 1.3)\%$ for almost vertical showers at $225\, \mathrm{m}$. Both dependencies are in qualitative agreement with theoretical predictions.


Introduction
Cosmic rays impinging on the Earth's atmosphere induce showers of secondary particles. The motions of electrons and positrons in the electromagnetic part of the shower produce radio emission. Such radio emission, in the form of strong coherent nanosecond timescale pulses, has been detected [1][2][3] and studied by several experiments [4][5][6]. It is now generally assumed that two mechanisms are principally responsible for the emission measured at MHz frequencies.
The dominant contribution to the radio emission is geomagnetic in origin [3,[7][8][9][10]. Electrons and positrons in the shower are accelerated in opposite directions by the Lorentz force exerted by the Earths magnetic field. The radio emission generated in this manner is linearly polarized in the direction of the Lorentz forceê v× B . Here v is the propagation velocity vector of the shower and B represents the Earths magnetic field.
A secondary contribution to the radio emission results from the net excess of electrons at the front of the shower [11]. This excess is built up from electrons that are knocked out of atmospheric molecules by interactions with shower particles and by a net depletion of positrons due to annihilation. This charge-excess contribution to the emission is also linearly polarized, but now radially with respect to the shower axis.
The resulting emission measured at ground level is the coherent sum of both components. Interference between these components may be constructive or destructive, depending on the position of the observer relative to the shower. Furthermore, the emission is strongly beamed in the forward direction due to the relativistic velocities of the particles. Additionally, the emission propagates through the atmosphere, which has a non-unity index of refraction that changes with height. This gives rise to relativistic time-compression effects most prominently resulting in a ring of amplified emission around the Cherenkov angle [12]. These time compression effects do not influence the polarization of the emission [13], but do contribute to the highly asymmetric total intensity pattern measured on the ground [14].
Recently, advances in simulations of air shower radio emission have converged to similar results [15]. Moreover, it has been demonstrated that such simulations accurately predict the complex emission pattern and can be used to derive air shower parameters [16]. This reflects the fact that the emission mechanisms are now well incorporated in them [17]. However, most models describe the emission at a microscopic level by calculating the contributions of each particle, in a full CORSIKA [18] or AIRES [19] shower simulation [14,20,21]. This makes it difficult to disentangle the precise contributions of the geomagnetic and charge-excess components.
It is important to stress here that a separation into macroscopic emission mechanisms is not necessary for a correct modeling of the emission (the charge-excess and geomagnetic effects are both captured by the distributions and motions of the simulated charged particles). However, a description in the form of macroscopic emission mechanisms can contribute to a better physical understanding of the conditions within an air shower.
In [13], it is predicted that the relative contributions of the charge-excess and geomagnetic components fall off differently with increasing opening angle from the emission maximum. This is perceived as an increase of the charge-excess fraction (the ratio between the strengths of both contributions to the total electric field) with increasing radial distance from the shower axis. Moreover, for air showers arriving at larger zenith angles the bulk of the emission is generated further away from the observer and thus the same distance to the shower axis corresponds to a smaller opening angle (where the charge-excess fraction is lower). Therefore the charge-excess fraction should also decrease with increasing zenith angle.
An experimental indication of the presence of the charge-excess component was found by the CODALEMA experiment as a shift of the radio signal maximum with respect to the particle core [22]. A measurement of a radially polarized emission component, consistent with that produced by the charge-excess mechanism, was obtained by the Auger Engineering Radio Array (AERA). This contribution was further quantified for their specific sample of air showers, giving an average strength of (14 ± 2)% in amplitude when compared to the geomagnetic contribution [23]. However, these measurements have a large antenna spacing and thus sample the radio emission of each shower only at a limited number of locations. This prohibits testing of the predicted dependence of the charge-excess fraction on the distance from the shower axis and the arrival direction of the shower.
Here we present high antenna density measurements of polarized radio emission from air showers with the LOFAR radio telescope. The instrumental setup is described in section 2. The measurements and data reduction techniques are described in sections 3 and 4. Section 5 describes the analysis procedure. Results are presented in section 6, followed by a discussion of systematic uncertainties in section 7. The article concludes in section 8.

Instrumental setup
The Low Frequency Array (LOFAR) [24], is a digital radio telescope, observing in the frequency range 10 − 240 MHz. It was designed as a flexible instrument, capable of carrying out multiple modes of observation simultaneously. One of these modes is the measurement of radio emission from extensive air showers [25]. For this purpose LOFAR is equipped with ring buffers, that store the raw voltage traces of each individual antenna in the array for up to 5 s. When a trigger is received the ring buffers are frozen and their contents copied over the network to a central storage location. This trigger can either be generated by inspecting the signal with the on-board FPGA 1 or received over the network from an external source.
For the purpose of cosmic-ray measurements a trigger is currently generated by a dedicated particle detector array. The LOFAR Radboud Airshower Array (LORA) [26], consists of 20 scintillator detectors. A trigger is sent to LOFAR when an air shower is detected with an energy exceeding E 2 · 10 16 eV. The corresponding radio emission for showers of this energy is at the lower limit of detectability by LOFAR for favorable shower geometries.
The particle detector array is located in the core of LOFAR. Here, the highest density of radio antennas is found, which makes the setup ideal for cosmic-ray detections. LOFAR consists of two types of antennas, the Low Band Antennas (LBAs) covering the frequency range 10 − 90 MHz and the High Band Antennas (HBAs) for 110 − 240 MHz (excluding the highly polluted 90 − 110 MHz commercial FM-radio band). Radio emission from cosmic rays has been measured in both frequency bands [25,27]. The analysis presented here focusses on the LBA measurements.
Each LBA consists of two inverted V-shaped dipoles labeled X and Y. These are aligned along the southwest to northeast and southeast to northwest direction respectively. They are grouped into stations, which in the core consist of 96 antennas each. Due to signal path limitations, only half of these antennas can be active during any given observation and, usually either the inner circle or the outer ring is selected. The ∼ 2 km diameter LOFAR core consists of 24 such stations with the highest density offered by the six stations located in a 350 m diameter region. An overview of the layout of all stations used in this analysis is given in figure 1.

Measurements
For this analysis, data collected with the outer rings of low-band antennas between June 2011 and January 2014 are used. The inner rings of low-band antennas are excluded for this analysis as crosstalk effects are present in some of the more closely spaced antennas, which are currently not included in the antenna simulations needed for polarization analysis. All data corresponding to a LORA trigger are stored together for offline analysis and constitute an event.
In the offline analysis these data are processed per station as described in [25]. The most important steps, for polarization analysis, are briefly summarized here.
First radio frequency interference (RFI) is identified and removed. Thereafter the received power in each dipole, outside of the cosmic-ray signal, is dominated by Galactic emission. This is used to perform a relative gain calibration between all dipoles, independently for the X and Y dipoles. An absolute calibration is currently not yet available at LOFAR. An initial estimate for the arrival direction of the air shower is given by the particle detector array. The predicted geometric signal arrival time delays for this direction are used to coherently add the signals from all dipoles in one station, thus forming a beam in this direction. This increases the signal-to-noise ratio for a cosmic-ray signal from that direction by approximately a factor of seven. When a pulse is found with a signal-to-noise ratio exceeding three the station is marked for further processing.
In the next processing step the arrival direction of the air shower is reconstructed from a plane wave fit to the arrival times of the pulse maxima (defined as the maximum of its Hilbert envelope). In this analysis only air showers where four or more stations have a successful direction reconstruction are included.
There are 230 air showers that meet this criterium. Of these a further 30 are excluded from this analysis, because they are likely influenced by thunderstorm conditions. These are defined by lightning strikes being recorded by the Royal Netherlands Meteorological Institute (KNMI) within two hours of the event. Five events do not coincide with a thunderstorm but show a similar polarization pattern and are also excluded. An additional 15 events were too weak to get a reliable estimate of the shower core position (see section 4.1). One additional event was excluded because the polarization reconstruction failed for a single station. This leaves a total of 179 air showers used for this analysis, for which the emission reaching the ground is sampled at between 192 and 528 distinct locations.

Reconstructing polarized radio emission
Using the reconstructed direction of the air shower, the full time traces measured by the X and Y dipoles are combined, while correcting for their frequency-dependent complex gain [25]. These gains are obtained from antenna simulations and are calculated for a plane wave, arriving from direction −ê v , where v is the propagation velocity vector of the air shower, and . Also depicted is the (north, east, zenith) coordinate frame corresponding to the x, y and z-axis, respectively. Furthermore the dipole antennas X and Y are shown (projected onto the ground plane). polarized alongê θ orê φ as defined in figure 2. The resulting combined signals are thus the electric field components alongê θ andê φ .
For the present analysis these are subsequently projected, following figure 3, onto thê e v× B andê v× v× B directions. Here B again represents the geomagnetic field.
Note that for this procedure it is assumed that the emission has no component along the propagation directionê v [25].

Observer positions in the shower plane
Given a position for the shower core and arrival direction, the antenna positions can be projected onto the shower plane as given by figure 4. Here each antenna i is represented by the polar coordinates φ i measured from theê v× B axis (positive in the direction of thê e v× v× B axis), and r i the distance from the shower axis. Figure 4. Geometry of the shower plane for a shower arriving with azimuth and zenith angles φ and θ respectively. The direction of the shower plane is defined by its normal vectorê v . All antenna positions (x i , y i ) are projected onto this plane giving r i , the distance to the shower axis and φ the angle with theê v× B direction. Note that the direction ofê v× B in this figure is chosen for clarity in display and does not accurately reflect the direction of the magnetic field at LOFAR.
While the particle detector array offers an initial estimate of the core position this is not reliable when the shower core is not contained within the particle detector array, as is often the case for the measured air showers. Therefore, the shower core position is determined by fitting a two-dimensional lateral distribution function to the integrated pulse power of the radio signal, following the procedure described in [28]. This provides the core position with an estimated statistical uncertainty of ∼ 15 m.
The arrival direction is the average of those obtained per station. This has an estimated statistical uncertainty of ∼ 1 • . A better angular resolution of ∼ 0.1 • can be obtained by fitting a hyperbolic wavefront to the arrival times at the antennas for all stations [29], however this is not needed for the current analysis.

Stokes parameters
Due to the expected elliptically polarized nature of the received signal [30] it is difficult to directly compare the electric field amplitudes in both polarization components. A better approach is to use time integrated quantities such as the Stokes parameters [31]: The integration is performed over n = 5 samples, of 5 ns each, centered around the pulse maximum. Here E i,j is sample i of the j−th electric field component andÊ i,j its Hilbert transform.
For an elliptically polarized signal one can calculate from the Stokes parameters the angle that the semi-major axis of the polarization ellipse makes with theê v× B axis Additionally the degree of polarization is calculated which is defined to be the fraction of the power in the polarized component of the wave

Determining the contributions of the emission mechanisms
Following [23] the expected electric field at any given time t can be written as Here E G (t) is the electric field produced by the geomagnetic contribution that is purely polarized along the direction of the Lorentz force, and E C (t) the radial electric field, produced by charge-excess. The charge-excess fraction is then defined as where |E C | is the amplitude of the electric field produced by charge-excess when the total electric field vector amplitude reaches its maximum value. The sin α factor, with α the angle between the magnetic field B and the propagation direction of the shower v, reflects the known dependence of the geomagnetic contribution [32].
Thus, the electric field at the time of the pulse maximum can be written as Therefore, the expected angle that the electric field vector at the time of the pulse maximum makes with theê v× B axis is given by This angle is equal to the angle of the semi major axis ψ = ψ of the polarization ellipse. Thus, the charge-excess fraction, a, can be determined by fitting eq. (5.4) to ψ as a function of azimuthal angle φ in the shower plane.

Results
Since we expect the air shower signal to be completely polarized, the degree of polarization should be close to unity. The distribution of p for all measured showers can be seen in the top panel of figure 5. With a degree of polarization of 98.0%, 99.2% and 99.7% in the 25%, 50% (median) and 75% quartiles respectively, the air shower signal is strongly polarized. For those cases where the ratio is less than unity the difference is expected to be caused by the unpolarized background noise. This can be seen in the bottom panel of figure 5 where the degree of polarization approaches unity for antennas with a high measured signal-to-noise ratio.
An instructive way to represent the polarization information is in the form of the polarization footprint of the shower as given in figure 6 for an example event. Here, the angle and degree of polarization are depicted as arrows for each antenna position projected in the shower plane. As expected for geomagnetically dominated emission, the arrows are roughly aligned along theê v× B direction. However, a small deviation from this direction can be seen. This deviation is interpreted as being due to the charge-excess component and causes the arrows to point upward or downward, above or below the projected shower core.

Measurement of the radially polarized emission component
In the presence of a radially polarized emission component, with strength a/ sin α relative to the geomagnetic component, we expect that the angle of polarization depends on the observer location in the shower plane according to eq. (5.4). In figure 7 this dependence can clearly be seen for two measured air showers.   Figure 6. Polarization footprint of a single air shower, as recorded with the LOFAR low-band antennas, projected onto the shower plane. Each arrow represents the electric field measured by one antenna. The direction of the arrow is defined by the polarization angle ψ with theê v× B axis and its length is proportional to the degree of polarization p. The shower axis is located at the origin (indicated by the black dot). In order to minimize systematic uncertainties, as discussed in section 7, only antennas with a measured signal-to-noise ratio exceeding three (in electric field amplitude) are included. Furthermore, a minimum of 48 antennas (one full station) are required to ensure a stable fit result. The distribution of the best fitting values for the charge-excess fractions of all 137 events surviving these cuts can be seen in figure 8. These fits give an average chargeexcess fraction of (11 ± 4)% for our sample, where the uncertainty reflects the spread of the distribution. This measured presence of a radially polarized emission component is consistent with that produced by the charge-excess mechanism. The uncertainty on the individual values of a are determined as described in appendix B and their distribution is plotted in figure 9. The fit quality, as parameterised by χ 2 r = χ 2 /dof, is given in figure 10. The overall fit quality is reasonably good χ 2 r = 2.6 ± 0.9. However, simulations predict additional dependencies on the distance to the shower axis as well as the shower arrival direction [13], these are not taken into account at this stage, which will necessarily lead to suboptimal fit results. For the same reason it is important to note that any quoted average charge-excess value will only apply to the specific set of showers for which it was measured. So while this analysis confirms that there is a radially polarized emission component, the average value itself has no meaning outside this sample. In the next section these dependencies are further investigated.

Checking for additional dependencies on the geomagnetic angle
It is important to note that eq. (5.2) assumes that the charge-excess fraction a only depends on the angle α, that the propagation axis of the shower makes with the geomagnetic field, through the strength of the geomagnetic contribution which is proportional to sin α. This assumption can now be checked by looking for an additional dependence of a to α in figure 11, where the charge-excess fraction is plotted as a function of α. No trend is seen, therefore  the geomagnetic angle and that sin α is the proper way of normalizing the geomagnetic component. Note that the scatter of the points is greater than their uncertainties suggest. This indicates an additional dependence which does not scale with the geomagnetic angle. This is almost certainly due to a dependence of the charge-excess fraction on the shower arrival direction and the distance to the shower axis, which cannot be determined on a single air shower basis due to the limited number of data points available after radial and angular binning (see section 6.3).

Dependence on shower arrival direction and radial distance to the shower axis
As is explained in [13], the relative contributions of the charge-excess and geomagnetic components fall off differently with increasing opening angle from the emission maximum. This is perceived as a fall-off of the charge-excess fraction with increasing radial distance from the shower axis. Moreover, at larger zenith angles the bulk of the emission is generated further away from the observer and thus the same distance to the shower axis corresponds to a smaller opening angle (where the charge-excess fraction is lower). Therefore, the chargeexcess fraction also decreases with increasing zenith angle.
To verify this, the antenna positions are projected onto the shower plane and subsequently their corresponding signals are grouped into radial distance bins of 50 m. For each group, a was reconstructed following the procedure from section 5. This was done separately for three different zenith angle bins. The result can be seen in table 1 and in figure 12. The uncertainty on a is determined as described in appendix B.
Both the predicted increase with increasing radial distance, as well as the decrease with Table 1. Charge-excess fraction as a function of the distance from the shower axis for three different zenith angle bins. increasing zenith angle, can clearly be seen in figure 12. Note however that the specific values obtained, and listed in table 1, still depend on the event set used due to shower-to-shower fluctuations.

Systematic uncertainties
While the addition of background noise results in an additional statistical uncertainty on the polarization angle and thus the charge-excess fraction, which is accounted for in the Monte Carlo procedures described in appendices A and B, it also introduces a systematic bias on the angle of polarization [30] which worsens with decreasing signal-to-noise ratio. While this in principle can be corrected for by subtracting the Stokes parameters calculated on background noise alone before calculating the angle of polarization, this has the downside of increasing the statistical uncertainty. For this reason it was opted to not correct for the systematic effect and instead minimize its influence by selecting only antennas with a measured signal-to-noise ratio (in electric field amplitude) exceeding three. A second possible systematic effect may be the result of using an inaccurate model to correct for the antenna response. It is however difficult to assess what effect this may have. Given that the antenna model gives an excellent agreement between data and simulation for total power (see [16]) we expect that this effect will be small.

Discussion & conclusions
Polarized radio emission from a sample of 179 air showers, measured with the LOFAR radio telescope, has been analyzed.
In total 35 air showers where excluded from this analysis due to coinciding thunderstorm activity. The strong atmospheric electric fields present during thunderstorms are expected to significantly affect the charged particle distributions and therefore the polarization of the emission [33,34]. This effect will be investigated in a future publication.
The measured emission is strongly polarized, with a median degree of polarization exceeding 99%. Because the geomagnetic and charge-excess emission are linearly polarized in different directions, their relative contributions can be determined by polarization analysis.
In all measured air showers the geomagnetic emission mechanism clearly dominates the polarization pattern. However, a sub-dominant charge-excess component can also be seen, varying in strength between showers. All dependence of the total electric field strength on the angle α, between the propagation direction of the shower and the geomagnetic field, is captured by the sin α dependence of the geomagnetic emission mechanism. No additional dependence of the charge-excess contribution on α is found.
The relative strength of both contributions is quantified by the charge-excess fraction. We find that the measured charge-excess fraction is higher for air showers arriving from closer to zenith. Furthermore, the measured charge-excess fraction also increases with increasing observer distance from the air shower symmetry axis. The measured values range from (3.47± 0.79)% for very inclined air showers at 25 m to (20.80 ± 0.98)% for almost vertical showers at 225 m. Both dependencies are in qualitative agreement with air shower simulations.
These measurements provide further confirmation that the mechanisms that produce radio emission in air showers are now well understood and incorporated in simulations which are used to reconstruct properties of the primary particle.

A Uncertainty on the angle of polarization
Both direction dependent effects and background noise contribute to the uncertainty on the angle of polarization ψ.
Any uncertainty on the direction of the shower axis translates into an uncertainty on the Stokes parameters through the combination of the measured signals in both dipoles in the antenna model. In order to estimate the uncertainty on ψ due to this effect a Monte Carlo procedure is employed. In this procedure the Stokes parameters and ψ are determined 200 times for each event while picking a random direction and core position from within the 1 • and 15 m uncertainties respectively. The spread (standard deviation) of the thus determined distribution of ψ gives the uncertainty due to the directional dependence σ ψ,d .
A second contribution to the total uncertainty on ψ results from the influence of background noise. An uncertainty σ ψ,n is assigned to each measurement based on the relation where S/N is the signal-to-noise ratio of the electric field amplitude. To derive this relation a similar Monte Carlo procedure is used. For 200 trials per event the Stokes parameters and ψ are determined after adding background noise to the signal. The spread on the thus determined distribution of ψ gives the uncertainty resulting from noise σ ψ,n .
However, because this Monte Carlo procedure lowers the signal-to-noise ratio before calculating the Stokes parameters, the resulting uncertainty contribution is overestimated (it is the contribution for S/( √ 2N ) instead of S/N ) and cannot be assigned to the measurement directly. Therefore, as a first step σ ψ,n is calculated for all measurements and plotted as a function of the actual S/N (corrected by a factor √ 2 to account for the addition of noise) in figure 13. In order to minimize the influence of outliers the median is calculated for each S/N bin. The bin size is 0.2 below S/N = 5 and increases to 1 and 5 for higher signal-to-noise ratios to ensure a sufficient number of measurements per bin. As can be seen, the uncertainty is inversely proportional to the signal-to-noise ratio and fitting for the proportionality constant gives relation (A.1).
The total uncertainty on ψ is now given by B Uncertainty on the charge-excess fraction As discussed in appendix A the uncertainty on the polarization angle ψ is caused by two effects. The influence of noise, and the propagation of the uncertainties on the arrival direction and core position of the shower through the antenna model. The uncertainty on ψ propagates into the uncertainty on the charge-excess fraction a, determined by fitting eq. (5.4), as σ a,ψ which is readily available from the covariance matrix of the fit. However there is a second effect. Any uncertainty on the core position and arrival direction of the shower also introduces an uncertainty on the azimuth of the observer position in the shower plane φ . This effect can be particularly strong for antennas close to the shower axis and needs to be taken into account in the uncertainty on a.
In order to estimate the uncertainty on the charge-excess fraction determined from a fit of eq. [rad] Figure 13. Uncertainty contribution to the polarization angle due to background noise as a function of signal-to-noise ratio in amplitude.
step, j, a core position and arrival direction were selected randomly within their respective uncertainties around the measured value giving φ i,j . Subsequently a was determined by fitting eq. (5.4) with fixed weights 1/σ ψ from eq. (A.2). This procedure was repeated 200 times. The standard deviation of the thus determined distribution of a is taken to be the uncertainty contribution σ a,φ of the uncertainties on core position and arrival direction through φ to a. The total uncertainty on a is found by adding the two contributions in quadrature σ 2 a = σ 2 a,ψ + σ 2 a,φ . (B.1)