Investigating an angular correlation between nearby starburst galaxies and UHECRs with the Telescope Array experiment

The arrival directions of cosmic rays detected by the Pierre Auger Observatory (Auger) with energies above 39~EeV were recently reported to correlate with the positions of 23 nearby starburst galaxies (SBGs): in their best-fit model, 9.7\% of the cosmic-ray flux originates from these objects and undergoes angular diffusion on a $12.9^\circ$~scale. On the other hand, some of the SBGs on their list, including the brightest one (M82), are at northern declinations outside the Auger field of view. Data from detectors in the northern hemisphere would be needed to look for cosmic-ray excesses near these objects. In this work, we tested the Auger best-fit model against data collected by the Telescope Array (TA) in a 9-year period, without trying to re-optimize the model parameters for our dataset in order not to introduce statistical penalties. The resulting test statistic (double log-likelihood ratio) was $-1.00$, corresponding to $1.1\sigma$ significance among isotropically generated random datasets, and to $-1.4\sigma$ significance among ones generated assuming the Auger best-fit model. In other words, our data is still insufficient to conclusively rule out either hypothesis. The ongoing fourfold expansion of TA will collect northern hemisphere data with much more statistics, improving our ability to discriminate between different flux models.


The Auger analysis
We will briefly illustrate the main features of the analysis by the Pierre Auger Collaboration [1], which we tried to replicate using Telescope Array data.

Assumed sources (23 starburst galaxies)
In Ref. [1], the candidate sources were selected from a list of 64 starburst galaxies (SBGs) outside the Local Group searched by Fermi-LAT [2] for gamma-ray emission (which was only significantly detected from 4 of them in that work). Among these, Ref. [1] selected the 23 objects whose radio flux at 1.4 GHz is at least 0.3 Jy (figure 1). Their UHECR luminosity was assumed to be proportional to their radio luminosity.

The flux model
In Ref. [1], the directional flux of UHECRs from SBGs is assumed to be a weighed sum of von Mises-Fisher distributions (the spherical analog of a Gaussian): where φ source is taken proportional to the radio flux at 1.4 GHz and Ψ is the RMS deviation in each transverse * e-mail: armando.di.matteo@ulb.ac.be * * e-mail: fujii@icrr.u-tokyo.ac.jp * * * e-mail: kawata@icrr.u-tokyo.ac.jp Figure 1. Distances from Earth and intrinsic luminosities of objects listed in Ref. [2], with the ones selected in the Ref. [1] analysis highlighted dimension (the total RMS deviation being √ 2Ψ). The total flux is assumed to be that from SBGs plus an isotropic background, where Φ iso = 1/4π.

The log-likelihood ratio test
Given two flux models Φ 1 (n), Φ 2 (n) and the directional exposure ω(n) of the experiment, the log-likelihood ratio test statistic is given by A positive TS indicates that the data favour the model Φ 2 over Φ 1 , and a negative TS vice versa. Ref. [1] used

Results
The best-fit parameters found by Ref. [1] were Ψ = 12.9 • , f SBG = 9.7%, and E min = 39 EeV. The corresponding flux map is shown in figure 2. This model is favoured at the 4.0σ (post-trial) level over an isotropic null hypothesis and at the 3.0σ level over the hypothesis that the UHECR emissivity is proportional to the overall matter density in the nearby extragalactic Universe.

The Telescope Array follow-up
Since several of the objects in the Auger analysis (including the brightest one, M82) are in the Northern Hemisphere outside the Auger field of view, we decided to test the same hypothesis using Telescope Array data [3].

The analysis
We computed TS for our data using the same values of Ψ and f SBG as in the best-fit of Ref. [1] (12.9 • , 9.7%), without scanning them in order to not introduce any statistical penalty. We neglected the attenuation of UHECRs in intergalactic space, which was found to be negligible by Ref. [1] because most of the flux from SBGs in this model originates from a few Mpc.

The dataset (284 events)
We used events collected by the Telescope Array (39.3 • N, 112.9 • W, 700 km 2 area) in the 9-year period from May 2008 to May 2017. We used the same quality cuts as in Ref. [4], among which zenith angle θ ≤ 55 • and declination δ ≥ 10 • . The energy threshold we used is E min = 43 EeV, corresponding to Ref. [1]'s 39 EeV when taking into account the 10% mismatch in energy scales between the two experiments [5]. We neglected the measurement resolution ( 20% on energy, 1.5 • on arrival directions) and assumed the detector to be fully efficient in the considered energy and zenith-angle range, calculating its exposure from purely geometrical considerations. This dataset includes 284 events, whose arrival directions are plotted in figure 3, along with the detector exposure multiplied by the model flux.

Result and discussion
The test statistic value we obtained was TS = −1.00, less than in 14.3% of the Monte Carlo simulations we made assuming an isotropic flux but more than in 7.5% of the simulations assuming the best-fit model from Ref. [1] is correct ( figure 4). This means our data are still insufficient to rule out either scenario. The upcoming expansion of Telescope Array, TA×4 [6], will increase its effective area by a factor of 4. If the data-taking of the TAx4 is started in 2019 as planned, by 2024 it will have collected the equivalent of 30 years of data with the current exposure, improving the statistical sensitivity proportionally to the square root of the number of events 1 , hence by a factor ≈ 1.8 compared to the present work.
The present study does not include the Local Group objects SMC, LMC, M33 and M31, which were also listed in Ref. [2] but in a separate table. Their intrinsic radio luminosities are very low compared to those in figure 1 (ranging from 0.19 to 6.2 Jy Mpc 2 ), but due to their proximities (0.05 to 0.85 Mpc) their radio fluxes are very large (3.3 to 444 Jy), and would dominate the sky if the assumed proportionality between UHECR and radio luminosities also applied to them. A possible justification for omitting them is hypothesizing that UHECRs are only accelerated during transient events, the expected number of which is proportional to the star-formation rate with a coefficient of order 0.01-0.1/(Jy Mpc 2 ), such that Poissonian fluctuations in the brightest few selected objects in figure 1 are minor but most of the time there are zero events in the Local Group. One example of such a class of events is gamma-ray bursts from supernovae (one ev-ery ∼ 2 500 years in M82, each producing UHECRs for ∼ 3 000 years) [7].