A TARGETED SEARCH FOR POINT SOURCES OF EeV NEUTRONS

The Pierre Auger Observatory measures cosmic-ray air showers near 1 EeV and higher energies (1 EeV = 10¹⁸ eV). An air shower produced by a neutron is indistinguishable from an air shower produced by a proton. Unlike protons and other nuclei, neutrons are not deflected by magnetic fields in the Galaxy, so their arrival directions point back to their sources. A flux of neutrons from a single direction can be detected as an excess of air showers arriving from that direction within the angular resolution of the Observatory. The mean decay path length for a neutron of energy E (measured in EeV) is 9.2E kpc. Above 1 EeV, the Galactic center is within the mean decay length, and above 2 EeV most of the Galactic disk is within range for neutron astronomy.

In a previous paper (Abreu et al. 2012a), the Pierre Auger Collaboration published a blind search for a neutron flux from any point of the sky with declination less than +15°, and celestial maps of flux upper limits were presented. No point stood out as statistically significant among the large number of trial source locations. In this paper, the search is limited to a small number of trials, each being a kind of "stacked analysis" for a set of candidate sources from an astrophysical catalog. The hypothesis is that many (or all) of the candidate sources of a given class are indeed emitting neutrons, so the combined signal should be more significant than that of a single target by itself. Performing the analysis on only a small number of target sets avoids a large statistical penalty.

The "ankle" of the energy spectrum (at about 5 EeV; Abraham et al. 2010a; Schulz et al. 2013) may represent a transition from Galactic to extragalactic sources of cosmic rays (Linsley 1963; Hillas 1984), although an alternative scenario holds that the "dip" of the ankle is formed by e^± pair production in a spectrum of extragalactic protons (Berezinsky et al. 2006). At energies in the range 1–5 EeV, Auger, HiRes, and Telescope Array have all found their data to be consistent with a significant component of protons among the cosmic rays (Abreu et al. 2012b; Abraham et al. 2010c; Abbasi et al. 2010; Jui et al. 2012). It is of particular interest to determine whether or not EeV protons are being produced at discrete sources in the Galaxy. A sensitive search for large scale anisotropy (Abreu et al. 2013) has found that the dipole anisotropy near 1 EeV is lower than what would be expected for protons emitted from the disk of the Galaxy. However, those expectations necessarily rely on general properties of the magnetic field in the Galaxy and its halo that have some uncertainties. The search for sources of EeV protons in the Galaxy is still highly motivated.

Any proton source is expected to produce some neutrons due to pion photo-production and nuclear interactions of the protons near the source. The detection of hadronic production of TeV gamma rays in some Galactic sources (Ackermann et al. 2013) provides direct evidence that the conditions for neutron production are favorable at least in those sources. TeV gamma-ray sources in the Galaxy are among the candidate sources targeted in this study. A 1/E² differential energy spectrum of protons from TeV to EeV in some of the H.E.S.S. sources¹⁰² would produce a neutron flux that is readily detectable by Auger. The Galactic center is a special candidate source that is well exposed to the Auger Observatory. The results here also update an earlier search for neutral particles from the Galactic center (Abraham et al. 2007).

The production of neutrons with the hadronic production of π⁺ mesons is necessarily accompanied by photons from decay of similarly produced neutral pions. The Auger Observatory can search for the existence of EeV photon fluxes using the hybrid data set that includes air fluorescence measurements. The search for point sources of photons near 1 EeV will be reported separately.

The Pierre Auger Observatory (Abraham et al. 2004) is centered at latitude 353 S, longitude 693 W, near Malargüe, Argentina, with mean altitude 1400 m above sea level (870 g cm⁻² atmospheric depth). The completed 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 2013 October 31. 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 (Abraham et al. 2010b). Periods of array instability have been omitted from the data set. The total exposure of the array with these cuts is 35,967 km² sr yr for the period of time analyzed here, yielding 854,270 events with E ⩾ 1 EeV.

The search is performed on nine target sets of astrophysically interesting directions. These classes of candidate sources are the millisecond pulsars (Manchester et al. 2005), γ-ray pulsars (Abdo et al. 2013), low-mass X-ray binaries (LMXBs; Liu et al. 2007a), high-mass X-ray binaries (HMXBs; Liu et al. 2007b), H.E.S.S. Pulsar Wind Nebulae, the other H.E.S.S. identified sources, the H.E.S.S. unidentified sources, microquasars,¹⁰³ and magnetars.¹⁰⁴ In addition to these target sets, the Galactic plane and the Galactic center are considered as 2 additional single-element sets for a total of 11 target sets. In order to have independent target sets, a source that appears in two or more sets is retained only in the most exclusive set (smallest set) while removed from the others. (The large Galactic plane target is allowed to contain targets from other sets.) Target sets are tested with and without statistical weights. A weight is assigned to each candidate source in proportion to its electromagnetic flux (recorded in the catalog), in proportion to its exposure to the Auger Observatory, and in proportion to the expected flux attenuation factor due to neutron decay. The flux attenuation factor is the fraction of emitted neutrons expected to survive decay from the distance of the candidate source assuming a 1/E² energy spectrum for emitted neutrons.¹⁰⁵ The target weights are normalized so that their sum is 1 in each target set.

The results without the weights are also reported since they are independent of the assumption that neutron emissions are proportional to the low-energy photon emissions and independent of the choice to make the weights proportional to each of the three factors.

Four energy ranges are used: 1 EeV ⩽E < 2 EeV (621,375 events), 2 EeV ⩽E < 3 EeV (135,444 events), E ⩾ 3 EeV (97,451 events), as well as E ⩾ 1 EeV. The first three are independent data sets, while the final cumulative data set should give maximum sensitivity to a flux that extends over the entire energy range. The Auger energy scale has a systematic uncertainty of 14% (Verzi et al. 2013).

The solid angle size for each target is optimized based on the average angular resolution for its declination and the energy range as explained in (Abreu et al. 2012a). The solid angle is a target circle of radius 1.05 times the angle within which 68% of neutron arrival directions from the candidate source should be included after the event reconstruction. For the special case of the Galactic plane, the signal-to-noise ratio is optimized by a strip centered on Galactic latitude b = 0° with half thickness |b| < 0.93ψ, where ψ is the mean angular resolution along the Galactic plane (and within the exposure of the Observatory) for a given energy range. For the four energy ranges above, the half thicknesses are, respectively, 128, 102, 072, and 117.

To recognize the existence of an excess of events in any target circle and any of the energy ranges, it is necessary to know the number that is expected in that circle without the neutron flux. The expected number of events in a given target circle is taken to be the average number found in 10,000 simulated data sets, each having the same number of events as in the actual data set. As in Abreu et al. (2012a), the arrival directions in a simulation event set are produced by sampling independently from the measured distributions for zenith angle, azimuth angle, and sidereal time. The average of many simulation data sets has no structure on small angular scales, providing a robust measure of the expected cosmic-ray background in each target.

The result for any target i is summarized by a p-value p_i. This p-value is here defined by p_i ≡ 1/2[Poisson(n, b) + Poisson(n + 1, b)], where Poisson(n, b) is the probability of getting n or more arrival directions in the target when the observed value is n, and the expected number from the background is b, as determined using simulated data sets. Averaging the values for n and n + 1 avoids the bias toward high p-values that occurs with Poisson(n, b) and the bias toward low p-values that occurs with Poisson(n + 1, b) for pure background fluctuations. When combining probabilities from a large target set with the Fisher formula (Fisher 1925), it is important that the individual p-values p_i have uniform expected distributions in the absence of any signal. Then for N targets with probabilities p_i, i = 1, 2, ..., N, the chance probability for their product (Π) not to be greater than their actual product (Π₀) is:

$$\mathbb {P} (\Pi \le \Pi _{0}) = \Pi _{0} \sum \limits _{j=0}^{N-1} \frac{(-\log \Pi _{0})^{j}}{j!} = 1-{\rm Poisson}(N,-\log \Pi _{0}).$$

For a weighted set of N targets with weights w_i, the combined p-value P_w is given by Good's formula (Good 1955). It is the chance probability for the weighted product of p-values (but sampled from uniform distributions) not to be greater than the actual weighted product Π_w. The weighted product is the product of factors $p_i^{w_i}$. Each p-value p_i is raised to the power w_i in the product of p-values, so the weight w_i can be regarded as the "number of times" the result for target i is counted relative to other targets of the set. In practice, the Good combined P_w is evaluated numerically using an ensemble of sets p_i (i = 1, 2, ..., N) with every p_i sampled randomly between 0 and 1.

Results are shown in Tables 1 and 2. The first table gives the unweighted combined p-value P for each of the 11 target sets and for each of the four energy ranges. The weighted combined p-value P_w is also given for each of the nine target sets that have multiple targets.

Table 2. Results for the Most Significant Target from Each Target Set

Class	R.A. [°]	Decl. [°]	Obs	Exp	Flux U.L.	E-Flux U.L.	p-value	p-value
					(km−2 yr−1)	(eV cm−2 s−1)		(penalized)
msec PSRs	260.27	−24.95	237	214	0.019	0.14	0.058	0.98
γ-ray PSRs	8.59	−5.58	176	149	0.024	0.18	0.016	0.70
LMXB	264.57	−26.99	265	219	0.028	0.20	0.0012	0.10
HMXB	152.45	−58.29	283	248	0.019	0.14	0.014	0.49
H.E.S.S. PWN	128.75	−45.60	275	248	0.018	0.13	0.043	0.53
H.E.S.S. other	269.72	−24.05	235	211	0.019	0.14	0.054	0.59
H.E.S.S. UNID	266.26	−30.37	251	227	0.018	0.13	0.055	0.57
Microquasars	262.75	−26.00	247	216	0.022	0.16	0.020	0.23
Magnetars	81.50	−66.08	268	241	0.016	0.11	0.040	0.48
Gal. center	266.42	−29.01	234	223	0.014	0.10	0.24	...
Gal. plane	|Gal. lat.| < 117		16965	17197	0.077	0.56	0.96	...

Note. The upper limits are computed at 95% confidence level.

Download table as:  ASCIITypeset image

The second table presents specific information for each target set about the candidate source that had the smallest individual p-value p_i for the full energy range E ⩾ 1 EeV. The direction of the source is given together with the observed number of events in the target, the expected number, the neutron flux upper limit, energy flux upper limit (assuming a 1/E² spectrum), and the p-value. The final column gives the penalized p-value p* = 1 − (1 − p)^N. This is the chance probability that one or more of the N candidate sources in the target set would have a p-value less than p if each p-value were randomly sampled from the uniform probability distribution.

The method to compute the limits is the same as the one explained in (Abreu et al. 2012a), where the definition of the upper limit in the number of neutrons is that of Zech (Zech 1989), using a 95% confidence level.

None of the candidate source classes tested in this study reveals compelling evidence for fluxes of EeV neutrons. Moreover, the minimum p-value for each target set is not statistically significant when penalized for the number of targets in the set. Neither collectively nor individually is there evidence for EeV neutron fluxes from the candidate sources that were examined.

The upper limits on the energy flux from these candidate sources, as shown in Table 2, are below the energy fluxes detected from TeV gamma-ray sources in the Galaxy. If such a source were accelerating protons in the same environment to EeV energies with the 1/E² dependence expected for Fermi acceleration, then the energy flux of neutrons in the EeV energy decade would exceed the energy flux in TeV gamma rays, since neutrons are produced more efficiently than gamma rays of equal energy.

The upper limit on the neutron flux from the Galactic plane provides a stringent constraint on models for continuous production of EeV protons in the Galaxy. The emission rate in the disk of the Galaxy must be sufficient to replace those protons as fast as they escape, providing an estimate on the required emission rate. The concomitant neutron emission rate is model dependent. It could exceed the proton emission rate if protons are magnetically bound to the sources and only the produced neutrons escape, yielding EeV protons by their later decays. More likely, however, the neutron luminosity at any fixed energy is less than the proton luminosity. Based on an estimate of the proton emission rate, the limits here on neutron flux from the Galactic plane can be used to derive a qualified upper limit on this ratio η = neutron luminosity/proton luminosity. The ratio η is closely related to the average optical depth of the source regions to escaping EeV protons. A very low value for η would imply that protons escape remarkably freely without interacting with photons or nuclei.

The expected neutron flux from the Galactic plane can be estimated from the luminosity density of protons ω in the Galactic plane together with the unknown ratio η. The luminosity density ω is the rate of proton production per unit area of the plane. It is given in terms of the proton density ρ, the thickness H of the cosmic-ray disk, and the proton residence time τ, by ω = ρH/τ. The density of protons ρ can be written as (4π/c)f_pI, where I is the cosmic-ray intensity and f_p denotes the fraction of cosmic rays that are protons. Substituting this expression for ρ yields:

$$\omega = 4\pi f_p I \frac{H}{c\tau }$$

for any given energy range. From any infinitesimal area element of the Galactic plane dA at distance r, the expected neutron flux is dF = η(ωdA/4πr²)exp (− r/D), where D is the mean neutron decay distance. (Note that the Auger Observatory has little exposure to the anti-center direction but good exposure for the Galactic plane in directions looking toward inner parts and the far side of the Galaxy, so the plane is treated crudely as emitting uniformly from the exposed annulus out to 20 kpc.) Integrating over an annulus of the Galactic plane centered on Earth gives an expected flux of neutrons:

$$\begin{eqnarray*} F_{{\rm GP}} &=& \eta 4\pi f_p I \frac{H}{c\tau } \int _{R_{{\rm min}}}^{R_{{\rm max}}}\frac{\exp (-r/D)}{4\pi r^2}2\pi r {\rm dr} \nonumber\\ &=& \eta 2\pi f_p I \frac{H}{c\tau } \int _{R_{{\rm min}}}^{R_{{\rm max}}}\frac{\exp (-r/D)}{r}{\rm dr}. \end{eqnarray*}$$

Here R_min is the distance at which most of the source disk fits within the angle (based on the angular resolution) used in collecting the signal from the Galactic plane, and R_max is the distance to the edge of the Galaxy. Roughly, R_min ∼ 1 kpc and R_max ∼ 20 kpc. Using D = 9.2E kpc (with E measured in EeV units), the final integral on the right has the following values: 1.97, 2.32, 2.53, and 2.27 respectively for the four energy ranges 1–2, 2–3, ⩾3, and ⩾1 EeV. Each is somewhat less than ln (20/1) = 2.996, which would be obtained in the limit D ≫ 20 kpc, and it should be noted that the expected flux has only an (approximately) logarithmic dependence on R_max/R_min, so the flux estimates are not terribly sensitive to the particular values adopted here. The proton fraction f_p just below the ankle of the cosmic-ray energy spectrum is believed to be at least 30% (Abreu et al. 2012b; Abraham et al. 2010c). A rough estimate for H/cτ is 0.1, based on proton escape times and density scale heights using model magnetic fields in the Galaxy. Using these estimates, together with I ≃ 29 (km² sr yr)⁻¹ for E ⩾ 1 EeV, the upper limit F_UL on a neutron flux from the Galactic plane yields a qualified upper limit on η given by η_UL = 0.08 F_UL. The upper limit F_UL ≃ 0.077 (km² yr)⁻¹ in Table 2 then implies η_UL ≃ 0.006, which is a significant constraint on models for continuous production of EeV protons in the Galaxy.

While the scope of this paper is a search for Galactic point sources, it should be noted that a flux of neutrons may be produced diffusely at a small level by interactions of EeV cosmic-ray protons in the Galactic disk. Our bound on η might therefore be slightly conservative.

A positive detection of neutron sources would have identified sources of EeV protons in the Galaxy. The null results here leave open the question of whether the observed EeV protons are produced in the Galaxy or whether they fill the space also between galaxies. They might be produced in transient events within the Galaxy, like supernova explosions or rare gamma ray bursts. The Auger Observatory could detect a strong flux of neutrons only at times of the bursts. The non-detection of any neutron flux averaged over the Auger exposure time does not constrain models for infrequent transient sources of EeV protons in the Galaxy. Also, protons emitted in jets would produce neutron jets with possibly none of them pointing toward Earth. Models with transient sources or such jet sources are better constrained by the absence of any strong anisotropy of the protons themselves (Abreu et al. 2011, 2012c).

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), São Paulo Research Foundation (FAPESP) grant Nos. 2010/07359-6, 1999/05404-3, Ministério de Ciência e Tecnologia (MCT), Brazil; MSMT-CR LG13007, 7AMB14AR005, CZ.1.05/2.1.00/03.0058 and the Czech Science Foundation grant 14-17501S, 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), Gran Sasso Center for Astroparticle Physics (CFA), CETEMPS Center of Excellence, 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; National Centre for Research and Development, grant Nos. ERA-NET-ASPERA/01/11 and ERA-NET-ASPERA/02/11, National Science Centre, Grant Nos. 2013/08/M/ST9/00322, 2013/08/M/ST9/00728, and 2013/10/M/ST9/00062, Poland; Portuguese national funds and FEDER funds within COMPETE-Programa Operacional Factores de Competitividade through Fundação para a Ciência e a Tecnologia, Portugal; Romanian Authority for Scientific Research ANCS, CNDI-UEFISCDI partnership projects nr.20/2012 and nr.194/2012, project nr.1/ASPERA2/2012 ERA-NET, PN-II-RU-PD-2011-3-0145-17, and PN-II-RU-PD-2011-3-0062, the Minister of National Education, Programme for research—Space Technology and Advanced Research—STAR, project number 83/2013, Romania; Slovenian Research Agency, Slovenia; Comunidad de Madrid, FEDER funds, Ministerio de Educación y Ciencia, Xunta de Galicia, Spain; The Leverhulme Foundation, Science and Technology Facilities Council, United Kingdom; Department of Energy, Contract No. DE-AC02-07CH11359, DE-FR02-04ER41300, and DE-FG02-99ER41107, National Science Foundation, grant No. 0450696, The Grainger Foundation, USA; NAFOSTED, Vietnam; Marie Curie-IRSES/EPLANET, European Particle Physics Latin American Network, European Union 7th Framework Program, grant No. PIRSES-2009-GA-246806; and UNESCO.