Testing a Reported Correlation between Arrival Directions of Ultra-high-energy Cosmic Rays and a Flux Pattern from nearby Starburst Galaxies using Telescope Array Data

The origins of ultra-high-energy cosmic rays (UHECRs) are still unknown. Anisotropies in the angular distribution of their arrival directions are rather small, requiring the detection of a large number of events to observe them. Furthermore, deflections of UHECRs by Galactic and intergalactic magnetic fields complicate the interpretation of anisotropies in terms of possible sources; this effect is reduced for the highest-energy cosmic rays, but the available statistics are significantly limited due to the steeply falling spectrum of UHECRs.

The two largest UHECR observatories in operation are the Telescope Array (TA; Abu-Zayyad et al. 2013a), located in Utah, USA, with approximately 700 km² effective area, and the Pierre Auger Observatory (Auger; Aab et al. 2015), located in Argentina with 3000 km² effective area. Their exposures peak in the Northern and Southern hemispheres, respectively.

Auger recently reported (Aab et al. 2018) a correlation between UHECR events with reconstructed energies above 39 EeV and a flux pattern of nearby starburst galaxies (SBGs). A model where 90.3% of the flux is isotropic and 9.7% originates from SBGs (with UHECR luminosities assumed to be proportional to their radio luminosities) and undergoes Gaussian random deflections with standard deviation 129 in each transverse dimension is favored over the purely isotropic model with a post-trial significance of 4.0σ, and over a model based on the overall galaxy distribution beyond 1 Mpc with a 3.0σ significance. In the Auger analysis it was found that different selections of candidate sources yield very similar results, as in any case over 90% of the anisotropic part of the flux weighed by the Auger directional exposure originates from four bright objects—NGC 4945, NGC 253, M83, and NGC 1068.

In this Letter, we follow up on this finding by testing UHECRs detected by TA in the Northern hemisphere against the same flux model and the best-fit values reported by Auger, and discuss possible interpretations of our result.

2.1. Cosmic-Ray Data Set

The TA is located at 393N, 1129W, in Millard County, Utah, USA, about 200 km southwest of Salt Lake City, about 1400 m above sea level (Abu-Zayyad et al. 2013a). The TA surface detector (SD) array consists of 507 plastic scintillation detectors on a square grid with 1.2 km spacing, covering an area of 700 km², and is surrounded by three fluorescence detector (FD) stations (Tokuno et al. 2012) with telescopes overlooking the SD array. It has been collecting data since 2008 May. The SD has ≈100% duty cycle, against ≈10% for the FD, so with a similar collection area the SD has about 10 times the statistics. The events detected in coincidence by both detectors are used to calibrate energy scale of the SD: SD reconstructed energies (determined by comparison to Monte Carlo simulations) are rescaled by a factor of 1/1.27 to match the FD energy scale (determined calorimetrically; Abu-Zayyad et al. 2013b; Tsunesada et al. 2017). The systematic uncertainty on the TA energy scale is 21% (Abbasi et al. 2016) and its energy and angular resolutions are 15%–20% and 10–15, respectively, depending on the event geometry and energy (Abbasi et al. 2014).

In this Letter we use data collected by the TA SD array in a 9 yr period from 2008 May to 2017 May with reconstructed energies above 43 EeV, zenith angles less than 55°, and declinations δ > −10° using the same quality cuts as in Abbasi et al. (2014). This data set comprises 284 events. We neglect the finite angular and energy resolution of TA events, and consider the detector fully efficient, i.e. with a flat response for all showers with energies and zenith angles in the considered range, so that its directional exposure ω_TA equals the geometrical one for δ > −10°, which varies with declination but not with right ascension (Sommers 2001):

$$\begin{eqnarray}\begin{array}{rcl}{\omega }_{\mathrm{TA}}(\delta ) & \propto & \cos {\phi }_{\mathrm{TA}}\cos \delta \sin {\alpha }_{{\rm{m}}}+{\alpha }_{{\rm{m}}}\sin {\phi }_{\mathrm{TA}}\sin \delta ,\\ {\alpha }_{{\rm{m}}} & = & \left\{\begin{array}{ll}\pi , & \xi \lt -1;\\ \arccos \xi , & -1\leqslant \xi \leqslant 1;\\ 0, & \xi \gt 1;\end{array}\right.\\ \xi & = & \displaystyle \frac{\cos {\theta }_{{\rm{m}}}-\sin {\phi }_{\mathrm{TA}}\sin \delta }{\cos {\phi }_{\mathrm{TA}}\cos \delta },\end{array}\end{eqnarray} \tag{ 1 }$$

where ϕ_TA = +393 is the detector latitude and θ_m = 55° is the maximum zenith angle accepted.

The energy threshold of E_min = 43 EeV used in this analysis corresponds to the Auger energy threshold of 39 EeV, at which the most significant correlation with SBG was found. Here we took into account the 10.4% difference between the energy scales of the two experiments as estimated by a comparison of energy spectra around 5 EeV (Verzi et al. 2017; Abu-Zayyad et al. 2018).

2.2. Source Catalog

Following the Auger analysis (Aab et al. 2018), we select the candidate sources from a sample of 63 SBGs outside the Local Group compiled by the Fermi-Large Area Telescope (LAT) Collaboration (Ackermann et al. 2012) for the gamma-ray emission search.³⁶ Imposing the cut of flux greater than 0.3 Jy at 1.4 GHz leaves 23 objects in the catalog of candidate sources. Their UHECR fluxes were assumed to be proportional to their radio fluxes at 1.4 GHz. These objects are listed in Table 1.

Table 1.  Selected Source Candidates from the SBG Catalog used in this Analysis (the same as in Aab et al. 2018)

Name	Gal. (l, b)		Distance	Flux ϕ	ϕωTA
NGC 253	974	−880	2.7 Mpc	13.6%	1.6%
M82	1414	406	3.6 Mpc	18.6%	35.7%
NGC 4945	3053	133	4.0 Mpc	16.0%	0.0%
M83	3146	320	4.0 Mpc	6.3%	0.4%
IC 342	1382	106	4.0 Mpc	5.5%	10.5%
NGC 6946	957	117	5.9 Mpc	3.4%	6.2%
NGC 2903	2087	445	6.6 Mpc	1.1%	1.4%
NGC 5055	1060	743	7.8 Mpc	0.9%	1.5%
NGC 3628	2409	648	8.1 Mpc	1.3%	1.5%
NGC 3627	2420	644	8.1 Mpc	1.1%	1.2%
NGC 4631	1428	842	8.7 Mpc	2.9%	4.4%
M51	1049	686	10.3 Mpc	3.6%	6.2%
NGC 891	1404	−174	11.0 Mpc	1.7%	2.8%
NGC 3556	1483	563	11.4 Mpc	0.7%	1.3%
NGC 660	1416	−474	15.0 Mpc	0.9%	1.0%
NGC 2146	1357	249	16.3 Mpc	2.6%	5.2%
NGC 3079	1578	484	17.4 Mpc	2.1%	3.8%
NGC 1068	1721	−519	17.9 Mpc	12.1%	9.1%
NGC 1365	2380	−546	22.3 Mpc	1.3%	0.0%
Arp 299	1419	554	46.0 Mpc	1.6%	2.9%
Arp 220	366	530	80.0 Mpc	0.8%	1.1%
NGC 6240	207	273	105.0 Mpc	1.0%	0.8%
Mkn 231	1216	602	183.0 Mpc	0.8%	1.4%

Note. The last column shows the relative source contribution weighted with the TA directional exposure ω_TA.

Download table as:  ASCIITypeset image

In the Auger analysis, the effect of energy losses by UHECRs during their propagation was found to be negligible in the SBG model, as most of the anisotropic flux originates from sources within a few Mpc; in this Letter, we neglected the losses for simplicity.

2.3. Test Statistic and Flux Model

Let $\hat{{\boldsymbol{n}}}$ be the unit vector representing a direction in the sky, pointing away from the observer. Given two flux models ${{\rm{\Phi }}}_{1}(\hat{{\boldsymbol{n}}}),{{\rm{\Phi }}}_{2}(\hat{{\boldsymbol{n}}})$ describing a null hypothesis and an alternative hypothesis, respectively, and the directional exposure $\omega (\hat{{\boldsymbol{n}}})$ of an experiment, the test statistic (hereinafter $\mathrm{TS}$) is defined as twice the log-likelihood ratio

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{TS} & = & 2\mathrm{ln}(L({{\rm{\Phi }}}_{2})/L({{\rm{\Phi }}}_{1})),\\ \mathrm{where}\,L({{\rm{\Phi }}}_{j}) & = & \displaystyle \prod _{i}\displaystyle \frac{{{\rm{\Phi }}}_{j}({\hat{{\boldsymbol{n}}}}_{i})\omega ({\hat{{\boldsymbol{n}}}}_{i})}{{\int }_{4\pi }{{\rm{\Phi }}}_{j}(\hat{{\boldsymbol{n}}})\omega (\hat{{\boldsymbol{n}}})\,{\rm{d}}{\rm{\Omega }}},\end{array}\end{eqnarray} \tag{ 2 }$$

and ${\hat{{\boldsymbol{n}}}}_{i}$ being the reconstructed arrival direction of the i-th observed event. A positive (negative) $\mathrm{TS}$ indicates that the data set is more (less) likely if the real flux is described by ${{\rm{\Phi }}}_{2}(\hat{{\boldsymbol{n}}})$ than by ${{\rm{\Phi }}}_{1}(\hat{{\boldsymbol{n}}})$.

In this analysis, the null hypothesis is an isotropic flux, ${{\rm{\Phi }}}_{1}(\hat{{\boldsymbol{n}}})={{\rm{\Phi }}}_{\mathrm{iso}}=1/4\pi$, whereas the alternative hypothesis is ${{\rm{\Phi }}}_{2}(\hat{{\boldsymbol{n}}})$ =

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\mathrm{mod}}(\hat{{\boldsymbol{n}}})={f}_{\mathrm{SBG}}{{\rm{\Phi }}}_{\mathrm{SBG}}(\hat{{\boldsymbol{n}}})+(1-{f}_{\mathrm{SBG}}){{\rm{\Phi }}}_{\mathrm{iso}},\end{eqnarray} \tag{ 3 }$$

where ${f}_{\mathrm{SBG}}=9.7 \%$ is the fraction of the flux assumed to originate from the SBGs in the catalog (the rest being assumed to be isotropic), and

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\mathrm{SBG}}(\hat{{\boldsymbol{n}}})=\displaystyle \frac{{\displaystyle \sum }_{k}{\phi }_{k}\exp ({\hat{{\boldsymbol{n}}}}_{k}\cdot \hat{{\boldsymbol{n}}}/{\theta }^{2})}{{\int }_{4\pi }{\displaystyle \sum }_{k}{\phi }_{k}\exp ({\hat{{\boldsymbol{n}}}}_{k}\cdot \hat{{\boldsymbol{n}}}/{\theta }^{2})\,{\rm{d}}{\rm{\Omega }}}\end{eqnarray} \tag{ 4 }$$

is a weighed sum of von Mises–Fisher distributions (the spherical analog of the Gaussian distribution), where ${\phi }_{k}$ and ${\hat{{\boldsymbol{n}}}}_{k}$ are the flux and position of the k-th source from Table 1 and θ = 129 is the rms deviation in each transverse dimension, the total rms deviation being $\sqrt{2}\theta$. The exposure is assumed to be geometrical, $\omega (\hat{{\boldsymbol{n}}})={\omega }_{\mathrm{TA}}(\hat{{\boldsymbol{n}}})$ from Equation (1). In this Letter we do not optimize the parameter values but keep them fixed to the Auger best-fit values, in order not to include any freedom in the model which would require a statistical penalty. The resulting model flux is shown in Figure 1, along with the events in the TA data set.

Zoom In Zoom Out Reset image size

Substituting the coordinates of the TA events $\{{\hat{{\boldsymbol{n}}}}_{i}\}$ into Equation (2), the test statistic that we obtained was $\mathrm{TS}\,=-1.00$. In order to assess the significance of this result, we computed TS for 10⁶ Monte Carlo (MC) data sets generated assuming an isotropic flux, and found $\mathrm{TS}\geqslant -1.00$ in p =14.3% of the 10⁶ cases, corresponding to a 1.1σ significance.³⁷

We also computed test statistics for 10⁶ MC sets generated under an assumption of the Auger best-fit SBG flux model to know the range of $\mathrm{TS}$ values that could be expected in that case. The results are shown in Figure 2. We found that 92.5% of realizations in the latter case have a higher $\mathrm{TS}$ value than the TA data (corresponding to a −1.4σ significance). We also verified that, as should be by design, the ratio between the two $\mathrm{TS}$ distributions is $\exp (\mathrm{TS}/2)$. A negative $\mathrm{TS}$ means that the angular distribution in a data set resembles isotropy more than the SBG model, and a positive $\mathrm{TS}$ means the reverse, so most isotropic realizations have $\mathrm{TS}\lt 0$ and most SBG-like realizations have $\mathrm{TS}\gt 0$. $\mathrm{TS}\approx 0$ would mean that the angular distribution in a data set is about equally different from the two models considered.

Zoom In Zoom Out Reset image size

A limitation in this analysis is the exclusion of Local Group objects (Small Magellanic Cloud (SMC), Large Magellanic Cloud (LMC), M33, and M31), which were listed in Ackermann et al. (2012) but in a separate table. These objects are not particularly intrinsically luminous (several times less than the dimmest objects in Table 1), but due to their proximity (D = 0.06, 0.05, 0.85, and 0.78 Mpc, respectively) they appear very bright. If the assumed proportionality between the UHECR luminosity, the star formation rate and the radio luminosity also applied to them, then the LMC and SMC would outshine all other objects combined in the Auger sky, and M33 and M31 would be the second- and third-brightest objects in the TA sky; but no excess of events is apparent in the vicinity of either pair of objects in our data or in Aab et al. (2018). A discussion about possible theoretical astrophysical motivations for not including these objects in the sample is outside the scope of this Letter.

Aab et al. (2018) also tested their data for correlations with gamma-ray loud active galactic nuclei from the second catalog of hard Fermi-LAT sources (Ackermann et al. 2016). The best fit (${E}_{\min }=60$ EeV, f_γAGN = 6.7%, θ = 69) is favored over isotropy at the 2.7σ level. Unlike with SBGs, UHECR energy losses in propagation are not negligible in this case because the unattenuated flux is not dominated by nearby objects. Testing TA data for correlations with this catalog would not be very useful, because the attenuated flux at Earth is dominated by Cen A, way outside of the TA field of view (at δ = −43°), leaving the flux in the northern hemisphere very nearly isotropic, and therefore requiring a very large number of events for an experiment in the northern hemisphere to detect the correlation; also, the Auger best-fit energy threshold found with this catalog (${E}_{\min }=60$ EeV) was higher than with the SBGs, further reducing the available statistics.

This Letter presents the result of a search for a correlation between arrival directions of UHECRs observed by TA and the flux pattern of SBGs. The SBG sample, anisotropic fraction, and angular scale were fixed to be the best-fit values as in the Auger study. The energy threshold of 43 EeV was determined by taking into account of the energy-scale difference between two experiments (Abu-Zayyad et al. 2018), corresponding to 39 EeV, at which the most significant correlation was reported in Auger. The result of this test was inconclusive, being compatible both with isotropy to within 1.1σ and with the Auger result to within 1.4σ. This means that the current TA data is not capable of discriminating between these two hypotheses. The ongoing expansion of TA (Kido 2018) will increase its effective area by a factor of 4, allowing us to reduce the statistical uncertainties and possibly to discriminate between different hypothesis about the UHECR origin.

We thank Pierre Auger collaboration members Jonathan Biteau and Olivier Deligny for useful discussions about their analysis.

The Telescope Array experiment is supported by the Japan Society for the Promotion of Science (JSPS) through Grants-in-Aid for Priority Area 431, for Specially Promoted Research JP21000002, for Scientific Research (S) JP19104006, for Specially Promoted Research JP15H05693, for Scientific Research (S) JP15H05741 and for Young Scientists (A) JPH26707011; by the joint research program of the Institute for Cosmic Ray Research (ICRR), The University of Tokyo; by the U.S. National Science Foundation awards PHY-0601915, PHY-1404495, PHY-1404502, and PHY-1607727; by the National Research Foundation of Korea (2016R1A2B4014967, 2016R1A5A1013277, 2017K1A4A3015188, 2017R1A2A1A05071429); by the Russian Academy of Sciences, RFBR grant 16-02-00962a (INR), IISN project No. 4.4502.13, and Belgian Science Policy under IUAP VII/37 (ULB). The foundations of Ezekiel R. and Edna Wattis Dumke, Willard L. Eccles, and George S. and Dolores Doré Eccles all helped with generous donations. The State of Utah supported the project through its Economic Development Board, and the University of Utah through the Office of the Vice President for Research. The experimental site became available through the cooperation of the Utah School and Institutional Trust Lands Administration (SITLA), U.S. Bureau of Land Management (BLM), and the U.S. Air Force. We appreciate the assistance of the State of Utah and Fillmore offices of the BLM in crafting the Plan of Development for the site. Patrick Shea assisted the collaboration with valuable advice on a variety of topics. The people and the officials of Millard County, Utah have been a source of steadfast and warm support for our work which we greatly appreciate. We are indebted to the Millard County Road Department for their efforts to maintain and clear the roads that get us to our sites. We gratefully acknowledge the contribution from the technical staffs of our home institutions. An allocation of computer time from the Center for High Performance Computing at the University of Utah is gratefully acknowledged.