A SEARCH FOR POINT SOURCES OF EeV NEUTRONS

Neutrons travel on straight lines, undeflected by magnetic fields, and they produce air showers that are indistinguishable from air showers produced by protons. A flux of neutrons from a discrete source would cause an excess of cosmic-ray events around the direction to the source, clustered within the angular resolution of the observatory. Since free neutrons undergo beta decay with a mean lifetime of about 886 s at rest (Particle Data Group 2010), the mean travel distance for relativistic neutrons is 9.2 × E kpc, where E is the energy of the neutron in EeV (1 EeV = 10¹⁸ eV). The distance from Earth to the Galactic center is about 8.3 kpc (Gillessen 2009), and the radius of the Galaxy is approximately 15 kpc. Sources in part of the Galactic disk, including the Galactic center, should be detectable via neutrons above 1 EeV. Above 2 EeV, the volume for detectable neutron emitters includes most of the Galaxy.

An important unresolved issue about high-energy cosmic rays is the transition at some energy from cosmic rays produced in the Galaxy to cosmic rays from extragalactic sources. The ankle of the energy spectrum can be explained as a dip caused by e^± production in collisions of predominantly extragalactic protons with cosmic microwave background photons (Blumenthal 1970; Berezinsky et al. 2006). A time-honored alternative view is that the ankle of the spectrum near 4 EeV (Pierre Auger Collaboration 2010a) is the transition from a Galactic power-law spectrum to a harder extragalactic power-law spectrum (Hillas 1972). If sources in the Galaxy are emitting protons up to the ankle of the spectrum, they could show themselves through a flux of EeV neutrons or perhaps a flux of EeV photons.

The signature of a neutron flux is a simple excess of proton-like air showers from a single celestial direction. In contrast, special discrimination techniques should be used to optimize sensitivity to EeV photons. A search for EeV photon fluxes will be reported separately. The method used in this paper is not sufficiently sensitive to photon fluxes because muon-poor photon showers produce less signal in water-Cherenkov detectors than proton showers of the same energy. This search is optimized for neutron fluxes, and the upper limits do not apply to photon fluxes.

Concerning the production of neutrons and photons by energetic protons, in both cases the dominant process is pion-producing interactions with ambient photons, protons, or nuclei. Two photons are produced by the decay of each π⁰. Neutrons are produced by charge-exchange interactions in which a π⁺ takes the positive charge of the proton and a leading neutron emerges with most of the energy that the proton had. The production of neutrons and photons has been studied extensively (Medina Tanco & Watson 2001; Bossa et al. 2003; Aharonian & Neronov 2005; Crocker et al. 2005), especially in relation to evidence from AGASA (Hayashida et al. 1999) and SUGAR (Bellido et al. 2001) for possible fluxes from directions close to the Galactic center, which were not confirmed using Auger data (Pierre Auger Collaboration 2007). Because photons acquire only a small fraction of the proton energy, the production of neutrons exceeds the hadronic production of photons of the same energy provided the accelerated proton spectrum falls approximately like 1/E² or more steeply with energy.

Based on the energy flux of gamma rays measured at TeV energies, known Galactic sources could plausibly be producing neutron fluxes which would be detectable by the Pierre Auger Observatory. The energy flux in TeV gamma rays exceeds 1 eV cm⁻² s⁻¹ at Earth for some Galactic sources (Hinton & Hofmann 2009; H.E.S.S. Collaboration 2011). Suppose those gamma rays arise from the decay of π⁰ mesons produced by interactions of protons accelerated at the source. A source with a 1/E² differential energy spectrum puts equal energy in each decade, and most models for cosmic-ray sources favor a power-law spectrum of approximately 1/E². If there are such sources in the Galaxy producing EeV photons as well as TeV gamma rays, then the energy flux of EeV photons should also exceed 1 eV cm⁻² s⁻¹ at Earth. For sources closer than the neutron attenuation length (9.2 × E kpc), the energy flux of neutrons could be even higher than the flux of photons above energy E because the production rate of neutrons at the source should exceed the production rate of photons of the same energy, as noted in the previous paragraph. Using three different energy thresholds, the results reported here show that there are no EeV sources of neutrons that bright in the southern sky that is exposed to the Auger Observatory.

Following a description of the data set (Section 2), the method used (Section 3), and the uncertainties (Section 4), the results of this blind search for a point-like neutron flux are reported in Section 5. Differential and integral plots of the Li–Ma significances (Li & Ma 1983) are presented, and upper limits are plotted on maps of the exposed sky. These results are then summarized and discussed in Section 6. A preliminary version of this study, using a slightly smaller data set, has been reported (Rouillé-d'Orfeuil 2011).

The Pierre Auger Observatory (Pierre Auger Collaboration 2004) is located in Malargüe, Argentina, at latitude 35.2 S, longitude 69.5 W, and mean altitude 1400 m above sea level (870 g cm⁻² atmospheric depth). The surface detector array consists of 1660 water-Cherenkov stations covering an area of about 3000 km² on a triangular grid with 1.5 km spacing, allowing secondary muons, electrons, and photons to be sampled at ground level with a duty cycle of nearly 100%.

The data set analyzed here consists of events recorded by the surface detector (SD) from 2004 January 1 to 2011 September 30. During this time, the size of the Observatory increased from 154 to 1660 SD stations. Events used in this analysis have zenith angles less than 60°. Moreover, an event is accepted only if all six nearest neighbors of the station with the highest signal were operational at the time of the event. This is the standard geometrical aperture cut that ensures good event reconstruction (Pierre Auger Collaboration 2010b). Periods of array instability have been omitted from the data set. The total exposure of the array with these cuts is 24,880 km² sr yr for the period of time analyzed here, yielding 429,138 events with E ⩾ 1 EeV.

The arrival direction of a cosmic ray is determined from a fit to the arrival times of the shower front at the SD stations. The precision achieved in the arrival direction depends on the clock resolution of each detector and on the fluctuations in the time of arrival of the first particle (Bonifazi 2009). The angular resolution is defined as the radius of the circular solid angle that would include 68% of the reconstructed events that arrive from a fixed direction. The angular resolution depends on energy as stated in Section 3.2.

The energy of a given air shower is measured by fitting for the ground signal S(1000) that a station would have measured at 1000 m from the core. This is converted to the energy parameter S₃₈, which is independent of zenith angle. The energy parameter S₃₈ has been calibrated using the quasi-calorimetric air fluorescence detector. See Pierre Auger Collaboration (2008) and Pesce (2011) for details about the SD energy determination. There is a systematic uncertainty of 22% in the absolute energy calibration. Statistical uncertainty in the SD energy determination is approximately 15%.

3.1. Energy Cuts

3.2. Target Sizes and Angular Resolutions

3.3. Simulation Data Sets

To recognize the existence of an excess of events in any solid angle "target," it is necessary to know the number that is expected in that target without the neutral flux. Simulation data sets are used for this. The expected number of events in a given target is taken to be the average number found in 10,000 simulated data sets.

The simulation data sets are obtained from the actual arrival directions, for each energy range, by a scrambling procedure that thoroughly smoothes out any small-scale anisotropy. Each simulation data set has the same number of arrival directions as the actual data set. An arrival direction is produced by randomly sampling a sidereal time from the set of measured sidereal times, a zenith angle from the set of measured zenith angles, and an azimuthal angle from a uniform distribution over 2π radians. Each simulation data set should be equivalent to the actual data aside from statistical fluctuations, unless astrophysical fluxes have imprinted small-scale anisotropy in the actual data.

Figure 1 shows the expected number of events per target, averaging over the targets with centers in 3° bands of declination. For each energy, the expected number depends on declination partly because of the declination dependence of the target size, but primarily because the directional exposure varies with declination.

Zoom In Zoom Out Reset image size

3.4. Li–Ma Significance

The statistical significance S of an excess (or deficit) in a given target is based on the number of events n observed in the target, the number b expected in the target from background cosmic rays, and the Li–Ma parameter α:

$$\begin{equation} S = \frac{n-b}{|n-b|} \sqrt{2} \left\lbrace n \ln \left(\frac{n+\alpha n}{b+ \alpha n}\right) + \frac{b}{\alpha } \ln \left(\frac{b + \alpha b}{b + \alpha n } \right) \right\rbrace ^{1/2}. \end{equation} \tag{ 1 }$$

This formula is Equation (17) of Li & Ma (1983) using N_on ≡ n and N_off ≡ b/α. In gamma-ray astronomy, α is the ratio of time spent observing on-source to the time spent observing an equivalent off-source solid angle. For the analysis here, all off-source regions are used in estimating the background, so α is taken to be the expected number in the target region divided by the expected number in the remainder of the sky.

3.5. Upper Limit Calculation

There are alternative ways to define the upper limit s_UL of confidence level CL for the expected signal s when an observation results in a count n in the presence of a Poisson background distribution with mean value b. The definition for s_UL adopted here is that of Zech (Zech 1989):

$$\begin{equation} P(\le n | b + s_{{\rm UL}}) = (1 - {\rm CL}) \times P(\le n | b), \end{equation} \tag{ 2 }$$

where CL is the fractional confidence level (e.g., 95% confidence level ⇒CL = 0.95 and 1 − CL = 0.05). This definition avoids unphysical negative upper limits that occur in the classical definition when the observed number is a strong downward fluctuation of the background. The frequentist interpretation (Zech 1989) of the above equation is as follows: "For an infinitely large number of repeated experiments looking for a signal with expectation s_UL and background with mean b, where the background contribution is restricted to a value less than or equal to n, the frequency of observing n or fewer events is α." This definition of the upper limit agrees, in this case of a Poisson process, with the Bayesian upper limit with flat prior:

$$\begin{equation} \int _0^{s_{{\rm UL}}} P(n|b+s) ds = ({\rm CL})\int _0^{\infty } P(n|b+s)ds. \end{equation} \tag{ 3 }$$

This upper limit s_UL is for the expected number of events from the source that would be contained within the top-hat target region. For the top-hat radius and assumed point-spread function described in Section 3.2, that top-hat region is expected to encompass 71.8% of the total signal. The upper limit for the total expected number of events from a source in the direction of the target center is obtained by scaling s_UL by the factor 1/0.718 = 1.39.

3.6. Flux Upper Limit

The flux upper limit is the upper limit on the number of events (as described above) divided by the directional exposure at the target center. The exposure depends on the trigger efficiency: this is 100% for energies above 3 EeV (Pierre Auger Collaboration 2004). Below that energy, the efficiency can depend on energy, zenith angle, and primary mass of the cosmic ray. The directional exposure for any celestial direction is given by b/ωI. Here b is the expected number (obtained empirically from the average of simulation data sets) in the target of solid angle ω, and I is the cosmic-ray intensity with units (km² sr yr)⁻¹ calculated by integrating the known energy spectrum (Salamida 2011) over the relevant energy range. The directional exposure is measured in units of km² yr. The dependence of the directional exposure on declination is shown in Figure 2 for the four different energy ranges. (Note that celestial points within 5° of the south pole are constantly exposed to the array at zenith angles between 50° and 60°.) For energy ranges that include energies below 3 EeV, the empirically derived b includes an implicit efficiency factor for triggering (generally less than unity) which depends on energy. The trigger efficiency depends also on the zenith angle of arrival, causing a slight dependence of the trigger efficiency on declination since the distribution of zenith angles varies with declination.

Zoom In Zoom Out Reset image size

This empirical determination of the directional exposure for an energy interval implicitly uses the cosmic-ray energy spectrum in weighting the average energy dependence of the directional exposure over that interval, whereas the (unknown) energy spectrum of a possible neutron flux would provide the ideal set of weights. This is not an issue above 3 EeV or for any narrow energy interval. Results are reported for two separate energy bins below 3 EeV partly to reduce this uncertainty in deriving a flux upper limit from the upper limit on the number of particles.

The flux upper limit is an upper limit on the time-averaged flux based on the cumulative Auger exposure. Periodic, episodic, or transient fluxes may exceed these time-averaged limits.

3.7. Pixelation and Target Spacing

Statistical uncertainty comes from Poisson fluctuations in the number of events in each target. The Li–Ma significance and the flux upper limit at a fixed confidence level are designed to account appropriately for the Poisson statistical fluctuations. As such, statistical error is not a concern.

A systematic error in the angular resolution can affect the overall normalization of the upper limits. The mean particle upper limit scales approximately with the square root of the expected number b in a target, and is therefore proportional to the target radius or the assumed angular resolution. The systematic uncertainty in the angular resolution for each energy bin is approximately 10%. Also, a genuine signal would be underestimated (overestimated) if measured within a fixed target radius proportional to an assumed angular resolution that is too small (large). The fractional error in the top-hat measurement is 0.714 times the fractional error in angular resolution for small fractional errors, which means about 7% uncertainty due to a 10% uncertainty in the angular resolution. The specific shape of the point-spread function assumed in Section 3.2 (a Rayleigh distribution) has little impact on the results. Any similar distribution having 68% containment within angular radius ψ should contain a fraction of the signal flux within the top-hat radius 1.05ψ that is not much different from the 71.8% Rayleigh expectation.

In principle, there could be error in the expected background counts due, for example, to imperfection in the smoothing procedure based on simulation data sets. However, those uncertainties are minuscule compared to Poisson fluctuations in all parts of the sky and for all energy ranges. Alternative methods for obtaining the expected counts give almost identical results. The uncertainty in expected counts is small and has a negligible impact on the results. Similarly, there could be some systematic uncertainty in the Li–Ma significances stemming from how the Li–Ma α parameter has been identified in this context. Small uncertainty in the background expectations corresponds to small α, and the Li–Ma significances are stable against changes in α as long as it remains small. The results here would be the same for alternative identifications of α that similarly imply insignificant uncertainty in the background counts.

The results are presented for fixed energy intervals. If there is a systematic error in the energy normalization, then the Li–Ma significance and the flux upper limit for each target pertain to a different true energy range. The Auger energy scale presently has a systematic uncertainty of 22% (Salamida 2011). It should be noted that the SD is not fully efficient below 3 EeV, and the trigger inefficiency can introduce a systematic energy error by favoring upward fluctuating signals in the surface stations. Events measured in hybrid mode by air fluorescence profiles as well as by the SD alone indicate that SD energy assignments are systematically high by about 2% for 2–3 EeV and close to 7% for the 1–2 EeV range.

In addition, fluctuations in energy measurements can cause unequal migration of signal and background events into (and out of) an energy range, thereby affecting apparent signals and upper limits. The effect depends on the exact shape of the arriving neutron signal spectrum, including its suppression at low energies due to in-flight neutron decays. A non-negligible underestimation in an upper limit could exist if there were no arriving neutrons to spill upward into an energy bin which does acquire background events by upward fluctuation in energy measurements. The errors are not large, however, since the energy measurement fluctuations are small compared to the size of the energy ranges, and the background contamination from outside an energy range does not exceed a few percent.

As explained in Section 3.6 and seen in Figure 2, the directional exposure has some dependence on energy. Ideally one would use the weighted average for the spectrum of a hypothetical neutron flux, but the method here uses a weighted average using the cosmic ray energy spectrum. For example, a neutron flux from a distant source might be fully attenuated below 3 EeV, so the trigger efficiency is 100% for that flux. For a measurement above 1 EeV, however, the directional exposure using the cosmic ray spectrum might be 20% lower (see Figure 2). The upper limit would therefore be (conservatively) too high by about 20%. This systematic uncertainty can be reduced by considering separately the three differential energy ranges where the range of energies is small for the intervals below 3 EeV.

Some systematic uncertainty exists due to the uncertain cosmic ray composition. There is good evidence for a mixed composition which includes protons throughout the EeV energy decade (Pierre Auger Collaboration 2010c), and neutron showers are reconstructed the same as proton showers. If the composition were purely heavy nuclei, the limits here would pertain to somewhat higher neutron energies than stated. For a mixed composition with a substantial proton component, however, energies reconstructed for neutron showers based on SD data do not differ systematically more than about 5% from the background cosmic ray energies which are calibrated using air fluorescence measurements.

5.1. Li–Ma Significances

Statistical results for the ensemble of celestial targets are shown in Figure 3. Each row of plots represents one of the four energy ranges. Red lines show the distribution of Li–Ma significance obtained from the data, while blue lines show the expectation obtained by averaging over simulation data sets. (Each simulation data set is analyzed exactly as the real data, using all of the other simulation data sets to determine the background for every target.) Also shown in each figure is the Gaussian function that the Li–Ma distribution is expected to approximate if deviations from expected values are due only to statistical fluctuations.

Zoom In Zoom Out Reset image size

For each energy range there are two plots. On the left is the differential histogram, binned in increments of Li–Ma significance. On the right are two (unbinned) integral distributions of the same Li–Ma significances. One focuses on the tail of high significance by plotting (for each Li–Ma significance) the total number of targets of equal or greater significance. The other focuses on the tail of low significances by plotting the total number of targets that had equal or lower Li–Ma significance. The shaded bands are 95% containment bands for simulation data sets. For any number of targets (plotted vertically), the shaded band extends horizontally over 95% of the simulation data sets; 2.5% of the simulation integral curves were to the left of the band at that vertical level, and 2.5% of the simulation integral curves were to the right of the band.

The fact that the red curve does not lie to the right of the right-hand shaded region means that this search has not identified obviously significant hot spots. The deviation from the Gaussian curve for negative significances in the case E ⩾ 3 EeV is caused by the very low statistics in many targets.

5.2. Upper Limits

Flux upper limits (95% CL) for each target direction are displayed in the color sky plots of Figure 4. Each limit is calculated according to the method explained in Sections 3.5 and 3.6, and it is the upper limit on the time-averaged neutron flux from that celestial direction.

Zoom In Zoom Out Reset image size

The mean flux upper limit is shown as a function of declination in Figure 5 for each of the energy ranges. The upper limits tend to be greater (weaker) for the northern declinations where the directional exposure (shown in Figure 2) is reduced. The limits are lowest (strongest) near the South Pole (−90° declination) where the directional exposure is maximum, but the mean value is less accurately determined in that region because there are relatively few targets in a declination band.

Zoom In Zoom Out Reset image size

The blind search for a flux of neutral particles using the Auger SD data set finds no candidate point on the sky that stands out among the large number of trial targets. Upper limits have been calculated for all parts of the sky using four different energy ranges. Three of those ranges are independent data sets and the fourth is the combination of the other three. These upper limits pertain to neutrons, with systematic uncertainties as discussed in Section 4. (The methods used in this paper are less sensitive to photons.)

The upper limits are generally more stringent where the directional exposure is relatively high, but they are strong enough to be of considerable astrophysical interest in all parts of the exposed sky. Above 1 EeV, the typical (median) flux upper limit is 0.0114 neutron km⁻² yr⁻¹. That corresponds to an energy flux limit of 0.083 eV cm⁻² s⁻¹ (or 0.026 EeV km⁻² yr⁻¹) in the EeV energy decade if the differential neutron spectrum is proportional to 1/E². Even for the regions of minimum sensitivity, the flux upper limit does not exceed 0.046 particles km⁻² yr⁻¹, corresponding to 0.34 eV cm⁻² s⁻¹ (or 0.106 EeV km⁻² yr⁻¹) for a 1/E² spectrum.

As noted in the introduction, this energy flux limit is well below what is observed from some Galactic TeV gamma-ray sources, and hadronic production of photons by protons with a 1/E² spectrum should have equal power in each energy decade. The luminosity emitted in neutrons should be at least as great as the luminosity emitted in hadronically produced photons. The upper limits on neutron fluxes at EeV energies indicate that TeV gamma-ray emission from those sources might be of electromagnetic origin or else their proton spectra are not as hard as 1/E² up to EeV energies.

At EeV energies, there is evidence that the cosmic ray composition includes a strong proton component (Pierre Auger Collaboration 2010c). With reasonable assumptions about the containment time of such protons in the Galaxy, it can be shown that a neutron flux should be detectable if the sources of those protons are in the Galaxy and continuously emitting protons in all directions, assuming the neutron luminosity of a source is not negligible compared to its proton luminosity and there are not more than a few such sources per cubic kiloparsec. The absence of any detectable neutron flux might suggest that the sources are extragalactic, or transient, or emitting in jets, or optically thin to escaping protons, or individually weak but densely distributed. The flux limits on neutrons couples with the absence of detectable large-scale anisotropy at EeV energies (Pierre Auger Collaboration 2011) to constrain models in which EeV protons are produced by a low density of strong sources in the Galaxy.

Although no flux of neutrons has been detected in this blind search of the exposed sky, it is possible that there is a measurable flux of neutrons from some candidate source of cosmic rays. There are some targets with Li–Ma significance greater than 4σ in this study, but their number is consistent with what is expected from normal statistical fluctuations. The blind search performed here necessarily involves a very large number of trials. It is sensible to look carefully at a small number of astrophysically motivated candidate source locations. Results of a targeted search will be reported separately.

The successful installation, commissioning, and operation of the Pierre Auger Observatory would not have been possible without the strong commitment and effort from the technical and administrative staff in Malargüe.

We are very grateful to the following agencies and organizations for financial support: Comisión Nacional de Energía Atómica, Fundación Antorchas, Gobierno De La Provincia de Mendoza, Municipalidad de Malargüe, NDM Holdings and Valle Las Leñas, in gratitude for their continuing cooperation over land access, Argentina; the Australian Research Council; Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), Financiadora de Estudos e Projetos (FINEP), Fundação de Amparo à Pesquisa do Estado de Rio de Janeiro (FAPERJ), Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP), Ministério de Ciência e Tecnologia (MCT), Brazil; AVCR AV0Z10100502 and AV0Z10100522, GAAV KJB100100904, MSMT-CR LA08016, LG11044, LC527, 1M06002, MSM0021620859 and RCPTM—CZ.1.05/2.1.00/03.0058, Czech Republic; Centre de Calcul IN2P3/CNRS, Centre National de la Recherche Scientifique (CNRS), Conseil Régional Ile-de-France, Département Physique Nucléaire et Corpusculaire (PNC-IN2P3/CNRS), Département Sciences de l'Univers (SDU-INSU/CNRS), France; Bundesministerium für Bildung und Forschung (BMBF), Deutsche Forschungsgemeinschaft (DFG), Finanzministerium Baden-Württemberg, Helmholtz-Gemeinschaft Deutscher Forschungszentren (HGF), Ministerium für Wissenschaft und Forschung, Nordrhein-Westfalen, Ministerium für Wissenschaft, Forschung und Kunst, Baden-Württemberg, Germany; Istituto Nazionale di Fisica Nucleare (INFN), Ministero dell'Istruzione, dell'Università e della Ricerca (MIUR), Italy; Consejo Nacional de Ciencia y Tecnología (CONACYT), Mexico; Ministerie van Onderwijs, Cultuur en Wetenschap, Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO), Stichting voor Fundamenteel Onderzoek der Materie (FOM), The Netherlands; Ministry of Science and Higher Education, grant Nos. N N202 200239 and N N202 207238, Poland; Fundação para a Ciência e a Tecnologia, Portugal; Ministry for Higher Education, Science, and Technology, Slovenian Research Agency, Slovenia; Comunidad de Madrid, Consejería de Educación de la Comunidad de Castilla La Mancha, FEDER funds, Ministerio de Ciencia e Innovación and Consolider-Ingenio 2010 (CPAN), Xunta de Galicia, Spain; Science and Technology Facilities Council, United Kingdom; Department of Energy, Contract Nos. DE-AC02-07CH11359, DE-FR02-04ER41300, National Science Foundation, grant No. 0450696, The Grainger Foundation, USA; ALFA-EC/HELEN, European Union 6th Framework Program, grant No. MEIF-CT-2005-025057, European Union 7th Framework Program, grant No. PIEF-GA-2008-220240, and UNESCO.