Cosmic Ray Mass Measurements with LOFAR

In the dense core of LOFAR individual air showers are detected by hundreds of dipole antennas simultaneously. We reconstruct X max by using a hybrid technique that combines a two-dimensional fit of the radio profile to CoREAS simulations and a one-dimensional fit of the particle density distribution. For high-quality detections, the statistical uncertainty on X max is smaller than 20 g/cm 2 . We present results of cosmic-ray mass analysis in the energy regime of 10 17 - 10 17.5 eV. This range is of particular interest as it may harbor the transition from a Galactic to an extragalactic origin of cosmic rays.


Introduction
Radio detection of extended air showers is a powerful and rapidly developing technique [1] that is suitable for measuring the atmospheric depth of the shower maximum X max [2,3]. Arrays with a high density of radio antennas, like LOFAR [4] and the future SKA [5], are capable of accurately measuring the mass composition of cosmic rays below the ankle, where a transition of a Galactic to an extragalactic component is expected to occur [6,7].
LOFAR is a digital radio telescope constructed in the North of the Netherlands with satellite stations across Europe. It consists of thousands of dipole antennas, sensitive in the frequency range of 10−240 MHz. It produces extremely detailed air shower radio data, by using the dense core region, or superterp, where 384 antennas are located within a circle of 320 m diameter. A particle array, LORA [8], has been installed in the core and is used for triggering and reconstruction. Each antenna contains a ring buffer which is read-out in case of a trigger. The full waveform is stored for offline analysis [9].
The vast amount of data per shower and the complicated radiation pattern require new reconstruction techniques. By fitting two-dimensional radiation profiles to the data, LOFAR has achieved a e-mail: Stijn.Buitink@vub.ac.be resolution on X max below 20 g/cm 2 [3] which is comparable to the resolution of fluorescence detection. An analysis of the first set of high-quality LOFAR events has yielded mean X max values that are in agreement with existing data based on other techniques and is indicative of a strong contribution of light nuclei in the energy range of 10 17 -10 17.5 eV [10].

Simulations
At LOFAR, showers are reconstructed using dedicated sets of CORSIKA [11] simulations, i.e. for each measured shower we generate enough proton showers and iron showers to span the range of possible values for X max . Such simulations are computationally expensive, so we aim to minimise the number of simulated showers per events, while still ensuring that (a) the complete range of naturally occurring X max values is covered, and (b) the reconstruction is accurate enough.
We achieve this by using CONEX to quickly generate hundreds of showers [12]. For these showers only the first few interactions are simulated in CORSIKA, after which the rest of the shower development is based on parametrizations. Since the value of X max critically depends on these first interactions the set of CONEX showers contains a realistic distribution of X max values. Moreover, it is possible to select a subset of these showers for subsequent full CORSIKA simulation. We choose this subset by first selecting the two proton showers with the highest and lowest value for X max . Then, we select nine more proton showers with intermediate values, (nearly) equally spaced between the extreme values. Since the range of X max values for iron showers is much narrower, we select a total of five iron showers using the same algorithm.
Finally we add additional shower simulations close to the estimated X max of the measured shower. This estimation is given by fitting a parametrization of the radio signal [13] to the data, that has an uncertainty of ∼ 40 g/cm 2 . We select 11 additional showers around this expected value.
The subset acquired in this way is then simulated in CORSIKA including the CoREAS plugin for calculation of the radio emission [14]. The radio emission is simulated for antennas that form a star-shaped pattern in the shower plane. An antenna model is applied to the simulated electric fields to calculate the measured signals in the two dipoles of the LOFAR low band antennas. The measured pulse power is defined as the total power in a window of 55 ns centered around the pulse maximum, summed over the two polarisations. Finally, a two-dimensional map of the pulse power is created by interpolation of the star-shaped pattern. A detailed description of this procedure can be found in [3].
A GEANT4 simulation of the LORA scintillators is used to calculate the lateral distribution of the deposited energy. So, for each measured shower a dedicated set of simulations is produced, and from each simulation in a set a two-dimensional radio profile and a one-dimensional particle profile are produced.

Reconstruction
The radio and particle profiles are fitted to the data simultaneously. The fit has four free parameters: two for the shower core position and scaling parameters for both the radio and the particle profiles. Inclusion of the particle data in the fit does not have much effect on the fitted core position in practice. However, it serves two important purposes: as an additional consistency check and to determine the energy scale. The analysis presented here does not yet include the absolute calibration of the antennas [15] so the scaling parameter for the radio power has an arbitrary value. In effect, the ratio of two scaling parameters also has an arbitrary value. It should however be the same constant for Their color reflects the measured power of the radio pulse. The background is a two-dimensional radio map from interpolated CoREAS simulations. The fit quality is high when the colors inside the circles blend into the background. Right: reconstruction of X max . Every data point corresponds to a single simulated shower. Together they form the dedicated simulation set produced for this particular shower. The reduced χ 2 of the fit is plotted against X max and a parabola is fitted to find the reconstructed X max . The inlay is a zoomed in version of the minimum of the curve demonstrating the sharpness of the minimum. Purple squares correspond to iron shower and blue circles to proton shower. all measured showers. Therefore, the width of the distribution of this parameter is a measure for the energy resolution of the reconstruction 1 .
The reduced χ 2 of the fit is very sensitive to X max as is seen in Fig. 1. The minimum of the χ 2 -curve is found by fitting a parabola to the points near the lowest value. The curve is not smooth because shower-to-shower fluctuations are not just variations in X max . Other parameters that will vary include the width of the longitudinal profile, the shape of the tail of the shower, and the lateral distribution of the shower. These variations cause the jitter in the data points. If an infinite amount of simulations could be produced, the data points are expected to form a band rather than a line.
The reconstructed value for X max found in this method is the atmospheric depth corresponding to the atmosphere used in the simulation. For this analysis we used the standard US atmosphere which is also a good average for northern Europe. For each measurement we construct a more accurate atmosphere using the GDAS database [16]. A first order correction is made by assuming that the reconstructed X max corresponds to the correct altitude h above sea level. The corrected value for X max is then found by evaluating the atmospheric depth of the GDAS profile at altitude h. However, it was realised that the change in refractive index (which also depends on the humidity profile) also affects the reconstructed X max . This effect is included as a systematic effect in the present analysis, although it can be corrected for in the future [17].

Event selection
For this analysis we have selected showers that were detected by at least four LOFAR antenna groups (each containing 48 active antennas). At trigger level, we require that 16 out of the 20 scintillator detectors find a signal. Both conditions introduce a possible bias towards light or heavy primaries in the sample. Heavy nuclei interact higher in the atmosphere, so the particle density at the ground becomes smaller, especially for more inclined showers. This leads to reduced trigger probability. On the other hand, protons and light nuclei penetrate deeper into the atmosphere and produce showers closer to the ground. Therefore, the size of the radio footprint is reduced and the probability of having a detectable signal in at least four antenna groups is reduced.
To ensure that our final event sample is bias-free, we require that all of simulated showers in the dedicated simulation set of a particular event (a) produce a trigger in the LORA array and (b) have a detectable radio signal in at least four antenna groups. This requirement removes only a few events above 10 17 eV, but below that energy the amount of removed events rapidly increases, mainly because of the trigger efficiency. Composition measurements below 10 17 eV become possible by changing the triggering strategy in the future, for example by triggering on radio-only [18] or a combination of radio and particles.
A quality cut is imposed in order to only include showers with relatively small uncertainty in the reconstruction on the core, the energy, and X max . In fact, those three uncertainties are strongly correlated and a single cut on the precision of the core suffices. The uncertainties are calculated using the dedicated set of simulations. They therefore depend on the arrival direction, core position and energy of the measured shower, but not on its value for X max . Thus, no additional bias on composition is introduced by this cut.
The final event sample counts 118 showers recorded before January 2015. The mean energy resolution of the sample is 32% and the resolution on X max is 17 g/cm 2 . The reduced χ 2 values of the fits range from 0.9 to 2.9.
Several systematic effects have been investigated. As mentioned above, humidity in the atmosphere increases the index of refraction and shifts the reconstructed X max . Since the simulations were performed for dry air this introduces an asymmetric systematic uncertainty of 10 g/cm 2 . The magnitude of the effect is expected to depend on zenith angle. Indeed, separating the event sample in two groups above or below 32 degrees zenith angle, yields a ±8 g/cm 2 difference in mean X max . The choice of hadronic interaction model also introduces another uncertainty of 5 g/cm 2 in the reconstruction method 2 . The total systematic uncertainty on X max is +14/-10 g/cm 2 and on the energy is 27%. More details are given in [3] and [10].

Results
The mean value of X max as a function of energy is shown in the left panel of Fig. 2. The error bars indicate 1σ statistical uncertainties and the shaded region shows the systematic uncertainty. The values are in agreement with other experiments using different techniques (fluorescence detection and non-imaging cherenkov). There is some tension with the data from the low-energy extension of the Pierre Auger observatory [19] towards 10 17.5 eV although the results still agree within systematic uncertainties. Further information on the composition can be extracted by studying the shape of the distribution of X max . For each shower, we calculate: where X shower is the reconstructed X max , and X proton and X iron are mean values predicted by the hadronic interaction code QGSJETII.04. The cumulative probability density function (CDF) is shown in the right panel of Fig. 2. The yellow curves at the extremes represent pure proton (centered at a=0, wide distribution) and pure iron (centered at a=1, narrow distribution) composition.
We fit theoretical CDFs based on composition with two or four mass components to the data. The test statistic in the fit is the maximum deviation between the data and the model CDFs. The p-value is given by the probability of observing this deviation, or a larger one, assuming the fitted composition model. The two-component model consists of proton and iron nuclei, and has only one free parameter: the mixing ratio. The best fit is found for a proton fraction of 62%, but it describes the data rather poorly with a p-value of 1.1 × 10 −6 .
The fit improves when using a four-component model (p+He+N+Fe), yielding a p-value of 0.17. While the best fit is found for a Helium fraction of 80%, the fit quality deteriorates only slowly when replacing helium by protons. Although there is a large uncertainty in the contribution of individual elements, the total fraction of light element (p+He) is well constrained. It has a best fit value of 0.8 and lies within the range [0.38,0.98] at a 99% confidence level based on QGSJETII.
It has been shown that the radio technique can be used to measure mass composition of cosmic rays and is competitive with other experiments. The first composition study with LOFAR yields a relatively high fraction of light particles just above 10 17 eV. The energy range can still be extended to lower and higher energies by expanding the triggering array and improving trigger strategies. Accurate composition measurements in the 10 16 -10 18 eV range are of crucial importance for understanding the transition of galactic to extragalactic origin of cosmic rays [7]. In the future, SKA will offer the opportunity of exploiting this technique with increased area, antenna density and frequency bandwidth, leading to extremely accurate measurements [5].