# A Search for Cosmic-ray Proton Anisotropy with the Fermi Large Area Telescope

M. Ajello, L. Baldini, G. Barbiellini, D. Bastieri, K. Bechtol, R. Bellazzini, E. Bissaldi, R.D. Blandford, R. Bonino, E. Bottacini, et al.

► To cite this version:

M. Ajello, L. Baldini, G. Barbiellini, D. Bastieri, K. Bechtol, et al. A Search for Cosmic-ray Proton Anisotropy with the Fermi Large Area Telescope. The Astrophysical Journal, 2019, 883 (1), pp.33. 10.3847/1538-4357/ab3a2e. hal-02080724

HAL Id: hal-02080724

https://hal.science/hal-02080724v1

Submitted on 12 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

CrossMark

# A Search for Cosmic-Ray Proton Anisotropy with the Fermi Large Area Telescope

M. Ajello¹, L. Baldini², G. Barbiellini³,⁴, D. Bastieri⁵,⁶, K. Bechtol⁷, R. Bellazzini⁸, E. Bissaldi⁹,¹⁰, R. D. Blandford¹¹, R. Bonino¹²,¹³, E. Bottacini¹¹,¹⁴, T. J. Brandt¹⁵, P. Bruel¹⁶, S. Buson¹⁷, R. A. Cameron¹¹, R. Caputo¹⁸, E. Cavazzuti¹⁹, S. Chen⁵,¹⁴, G. Chiaro²⁰, S. Ciprin¹¹,²², J. Cohen-Tanugi²³, D. Costantin²⁴, A. Cuoco¹²,²⁵, S. Cutini²⁶, F. D'Ammando²⁷,²⁸, P. de la Torre Luque⁹, F. de Palma¹², A. Desai¹, S. W. Digel¹¹, N. Di Lalla², L. Di Venere⁹,¹⁰, A. Domínguez²⁹, S. J. Fegan¹⁶, Y. Fukazawa³⁰, S. Funk³¹, P. Fusco⁹,¹⁰, F. Gargano¹⁰, D. Gasparrini²²,²⁶, N. Giglietto⁹,¹⁰, F. Giordano⁹,¹⁰, M. Giroletti²⁷, D. Green³², I. A. Grenier³³, S. Guiriec¹⁵,³⁴, K. Hayashi³⁵, E. Hays¹⁵, J. W. Hewitt³⁶, D. Horan¹⁶, G. Jóhannesson³⁷,³⁸, M. Kuss⁸, L. Latronico¹², J. Li³⁹, I. Liodakis¹¹, F. Longo³,⁴, F. Loparco⁹,¹⁰, P. Lubrano²⁶, S. Maldera¹², A. Manfreda², G. Martí-Devesa⁴⁰, M. N. Mazziotta¹⁰, M. Meehan⁷, I. Mereu⁴¹, M. Meyer¹¹, P. F. Michelson¹¹, N. Mirabal¹⁵,⁴², W. Mitthumsiri⁴³, T. Mizuno⁴⁴, A. Morselli²¹, M. Negro¹²,¹³, E. Nuss²³, N. Omodei¹¹, M. Orienti²⁷, E. Orlando¹¹, V. S. Paliya³⁹, D. Paneque³², M. Persic³,⁴⁵, M. Pesce-Rollins⁸, F. Piron²³, T. A. Porter¹¹, G. Principe³¹, S. Rainò⁹,¹⁰, R. Rando⁵,⁶, M. Razzano⁸,⁵⁴, S. Razzaque⁴⁶, A. Reimer¹¹,⁴⁰, O. Reimer¹¹,⁴⁰, D. Serini⁹, C. Sgrò⁸, E. J. Siskind⁴⁷, G. Spandre⁸, P. Spinelli⁹,¹⁰, D. J. Suson⁴⁸, H. Tajima¹¹,⁴⁹, J. B. Thayer¹¹, D. F. Torres⁵⁰,⁵¹, E. Troja¹⁵,⁵², J. Vandenbroucke⁷, M. Yassine³,⁴, and S. Zimmer⁴⁰,⁵³

## The Fermi-LAT Collaboration

¹ Department of Physics and Astronomy, Clemson University, Kinard Lab of Physics, Clemson, SC 29634-0978, USA

² Università di Pisa and Istituto Nazionale di Fisica Nucleare, Sezione di Pisa I-56127 Pisa, Italy

³ Istituto Nazionale di Fisica Nucleare, Sezione di Trieste, I-34127 Trieste, Italy

⁴ Dipartimento di Fisica, Università di Trieste, I-34127 Trieste, Italy

⁵ Istituto Nazionale di Fisica Nucleare, Sezione di Padova, I-35131 Padova, Italy

⁶ Dipartimento di Fisica e Astronomia "G. Galilei," Università di Padova, I-35131 Padova, Italy

⁷ Department of Physics, University of Wisconsin-Madison, Madison, WI 53706, USA; mrmeehan@wisc.edu, justin.vandenbroucke@wisc.edu

⁸ Istituto Nazionale di Fisica Nucleare, Sezione di Pisa, I-56127 Pisa, Italy

⁹ Dipartimento di Fisica "M. Merlin" dell'Università e del Politecnico di Bari, I-70126 Bari, Italy

¹⁰ Istituto Nazionale di Fisica Nucleare, Sezione di Bari, I-70126 Bari, Italy

¹¹ W.W. Hansen Experimental Physics Laboratory, Kavli Institute for Particle Astrophysics and Cosmology, Department of Physics and SLAC National Accelerator Laboratory, Stanford University, Stanford, CA 94305, USA

¹² Istituto Nazionale di Fisica Nucleare, Sezione di Torino, I-10125 Torino, Italy

¹³ Dipartimento di Fisica, Università degli Studi di Torino, I-10125 Torino, Italy

¹⁴ Department of Physics and Astronomy, University of Padova, Vicolo Osservatorio 3, I-35122 Padova, Italy

¹⁵ NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA

¹⁶ Laboratoire Leprince-Ringuet, École polytechnique, CNRS/IN2P3, F-91128 Palaiseau, France

¹⁷ Institut für Theoretische Physik und Astrophysik, Universität Würzburg, D-97074 Würzburg, Germany

¹⁸ Center for Research and Exploration in Space Science and Technology (CRESST) and NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA

¹⁹ Italian Space Agency, Via del Politecnico snc, I-00133 Roma, Italy

²⁰ INAF-Istituto di Astrofisica Spaziale e Fisica Cosmica Milano, via E. Bassini 15, I-20133 Milano, Italy

²¹ Istituto Nazionale di Fisica Nucleare, Sezione di Roma "Tor Vergata," I-00133 Roma, Italy

²² Space Science Data Center—Agenzia Spaziale Italiana, Via del Politecnico, snc, I-00133, Roma, Italy

²³ Laboratoire Univers et Particules de Montpellier, Université Montpellier, CNRS/IN2P3, F-34095 Montpellier, France

²⁴ University of Padua, Department of Statistical Science, Via 8 Febbraio, 2, I-35122 Padova, Italy

²⁵ RWTH Aachen University, Institute for Theoretical Particle Physics and Cosmology, (TTK), D-52056 Aachen, Germany

²⁶ Istituto Nazionale di Fisica Nucleare, Sezione di Perugia, I-06123 Perugia, Italy

²⁷ INAF Istituto di Radioastronomia, I-40129 Bologna, Italy

²⁸ Dipartimento di Astronomia, Università di Bologna, I-40127 Bologna, Italy

²⁹ Grupo de Altas Energías, Universidad Complutense de Madrid, E-28040 Madrid, Spain

³⁰ Department of Physical Sciences, Hiroshima University, Higashi-Hiroshima, Hiroshima 739-8526, Japan

³¹ Friedrich-Alexander-Universität Erlangen-Nürnberg, Erlangen Centre for Astroparticle Physics, Erwin-Rommel-Str. 1, D-91058 Erlangen, Germany

³² Max-Planck-Institut für Physik, D-80805 München, Germany

³³ AIM, CEA, CNRS, Université Paris-Saclay, Université Paris Diderot, Sorbonne Paris Cité, F-91191 Gif-sur-Yvette, France

³⁴ The George Washington University, Department of Physics, 725 21st St, NW, Washington, DC 20052, USA

³⁵ Department of Physics and Astrophysics, Nagoya University, Chikusa-ku Nagoya 464-8602, Japan

³⁶ University of North Florida, Department of Physics, 1 UNF Drive, Jacksonville, FL 32224, USA

³⁷ Science Institute, University of Iceland, IS-107 Reykjavik, Iceland

³⁸ Nordita, Royal Institute of Technology and Stockholm University, Roslagstullsbacken 23, SE-106 91 Stockholm, Sweden

³⁹ Deutsches Elektronen Synchrotron DESY, D-15738 Zeuthen, Germany

⁴⁰ Institut für Astro- und Teilchenphysik and Institut für Theoretische Physik, Leopold-Franzens-Universität Innsbruck, A-6020 Innsbruck, Austria

⁴¹ Dipartimento di Fisica, Università degli Studi di Perugia, I-06123 Perugia, Italy

⁴² Department of Physics and Center for Space Sciences and Technology, University of Maryland Baltimore County, Baltimore, MD 21250, USA

⁴³ Department of Physics, Faculty of Science, Mahidol University, Bangkok 10400, Thailand

⁴⁴ Hiroshima Astrophysical Science Center, Hiroshima University, Higashi-Hiroshima, Hiroshima 739-8526, Japan

⁴⁵ Osservatorio Astronomico di Trieste, Istituto Nazionale di Astrofisica, I-34143 Trieste, Italy

⁴⁶ Department of Physics, University of Johannesburg, PO Box 524, Auckland Park 2006, South Africa

⁴⁷ NYCB Real-Time Computing Inc., Lattingtown, NY 11560-1025, USA

⁴⁸ Purdue University Northwest, Hammond, IN 46323, USA

⁴⁹ Institute of Space Sciences (CSIC-IEEC), Campus UAB, Carrer de Magrans s/n, E-08193 Barcelona, Spain
Abstract

The Fermi Large Area Telescope (LAT) has amassed a large data set of primary cosmic-ray protons throughout its mission. In fact, it is the largest set of identified cosmic-ray protons ever collected at this energy. The LAT's wide field of view and full-sky survey capabilities make it an excellent instrument for studying cosmic-ray anisotropy. As a space-based survey instrument, the LAT is sensitive to anisotropy in both R.A. and decl., while ground-based observations only measure the anisotropy in R.A. We present the results of the first-ever proton anisotropy search using Fermi LAT. The data set was collected over eight years and consists of approximately 179 million protons above 78 GeV, enabling it to probe dipole anisotropy below an amplitude of $10^{-3}$, resulting in the most stringent limits on the decl. dependence of the dipole to date. We measure a dipole amplitude $\delta = 3.9 \pm 1.5 \times 10^{-4}$ with a p-value of 0.01 (pretrials) for protons with energy greater than 78 GeV. We discuss various systematic effects that could give rise to a dipole excess and calculate upper limits on the dipole amplitude as a function of minimum energy. The 95% confidence level upper limit on the dipole amplitude is $\delta_{UL} = 1.3 \times 10^{-3}$ for protons with energy greater than 78 GeV and $\delta_{UL} = 1.2 \times 10^{-3}$ for protons with energy greater than 251 GeV.

Unified Astronomy Thesaurus concepts: Galactic cosmic rays (567); High-energy astrophysics (739)

1. Introduction

Cosmic rays with energy $\lesssim 10^{18}$ eV diffuse through interstellar magnetic fields toward Earth, where they arrive with a high degree of isotropy. However, a small anisotropy in the arrival directions of cosmic rays of $O(10^{-4} - 10^{-3})$ has been consistently observed over the past several decades. A large ($\sim$6%) anisotropy was recently measured at the $10^{18}$ eV scale by the Pierre Auger Observatory (The Pierre Auger Collaboration et al. 2017). The cosmic-ray anisotropy landscape has recently grown more complex as large experiments with long duty cycles have measured the anisotropy over nine decades in energy with unprecedented precision (Amenomori et al. 2005, 2006, 2017; Abdo et al. 2008, 2009; Aglietta et al. 2009; Abbasi et al. 2010, 2011, 2012; Aartsen et al. 2013, 2016; Bartoli et al. 2013, 2015; Abeysekara et al. 2014, 2018, 2019). Broadly speaking, the anisotropy in the TeV-PeV range can be described by a large-scale (dominated by dipole and quadrupole terms) feature with an energy-dependent amplitude and phase. Anisotropy at medium and small scales has also been measured down to angular scales of $\sim 10^{\circ}$, though with amplitudes an order of magnitude lower than the large-scale anisotropy (Abbasi et al. 2010; Abeysekara et al. 2014, 2018). A variety of physical mechanisms could explain the large-scale anisotropy, though there is no consensus on the exact causes of the energy dependence. Standard diffusion theory predicts a dipole in the direction of the density gradient of cosmic rays, but the predicted amplitude is up to two orders of magnitude larger than the observed anisotropy (Hillas 2005; Ptuskin et al. 2006; Blasi & Amato 2012). The observed large-scale anisotropy at Earth could be due to the particular distribution of sources nearby as well as their directions relative to the local interstellar magnetic field (Mertsch & Funk 2015; Schwadron et al. 2015). Observational effects of measuring the anisotropy with Earth-fixed observatories also obscure the true two-dimensional phase of the anisotropy. For example, the analysis techniques used to reach the appropriate sensitivity are incapable of measuring anisotropy along the decl. axis, resulting in measurements along R.A. only. Partial sky coverage also biases the measurement of large-scale

features. The all-sky anisotropy was recently measured using a combined data set from the IceCube and HAWC detectors, which demonstrated that previous measurements by either detector alone underestimate the large-scale amplitude (Abeysekara et al. 2019). It has been shown that the observed large-scale anisotropy can be explained by a combination of the astrophysical and instrumental effects described above (Ahlers 2016). Many of the systematic effects introduced by the analysis techniques used by ground-based experiments can be mitigated by studying anisotropy with a full-sky, space-based observatory.

The Fermi Large Area Telescope (LAT) scans the entire celestial sky and detects cosmic rays in the GeV-TeV energy range (in addition to its primary purpose as a gamma-ray telescope). As a space-based survey instrument, it is sensitive to cosmic-ray anisotropy in both R.A. and decl. Additionally, the subsystems of the instrument can measure the charge of the cosmic rays, enabling a proton-only measurement of the anisotropy.<fnref n="55" /> The study of cosmic-ray anisotropy in this energy range offers complementary information to that at higher energies—as well as constraints on the decl. dependence of the anisotropy, which has never been measured. Furthermore, the study of anisotropy in this energy range is complementary to the study of the energy spectrum of protons and could shed light on the unexpected spectral break at a few hundred GeV (Adriani et al. 2011; Aguilar et al. 2015).

2. Fermi Large Area Telescope

The LAT is a pair-conversion gamma-ray telescope on board the Fermi Gamma-ray Space Telescope (Fermi) mission. Its wide instantaneous field of view (2.4 sr) and full-sky survey capabilities make it an excellent instrument for studying cosmic-ray anisotropy. The Fermi spacecraft is in an equatorial orbit with an inclination of 25°.6. It rocks north and south from zenith toward the celestial poles on successive orbits, enabling the LAT to scan the entire sky every two orbits ($\sim$3 hr).<fnref n="56" /> The Fermi-LAT Collaboration has published two studies of the

<footnote n="51">⁵¹ Institució Catalana de Recerca i Estudis Avançats (ICREA), E-08010 Barcelona, Spain</footnote>

<footnote n="52">⁵² Department of Astronomy, University of Maryland, College Park, MD 20742, USA</footnote>

<footnote n="53">⁵³ University of Geneva, Département de physique nucléaire et corpusculaire (DPNC), 24 quai Ernest-Ansermet, CH-1211 Genève 4, Switzerland</footnote>

<footnote>Received 2019 March 6; revised 2019 August 8; accepted 2019 August 9; published 2019 September 18</footnote>

<footnote n="54">⁵⁴ Funded by contract FIRB-2012-RBFR12PM1F from the Italian Ministry of Education, University, and Research (MIUR).</footnote>

<footnote n="55">⁵⁵ This is in contrast to ground-based detectors, which have poor composition resolution.</footnote>

<footnote n="56">⁵⁶ The rocking angle of the instrument (angle between zenith and instrument axis) was 35° for the first year of data taking and increased to 50° thereafter.</footnote>
anisotropy of cosmic-ray electrons and positrons (CREs), the latter of which provides the most stringent constraints to date on the dipole anisotropy (Ackermann et al. 2010; Abdollahi et al. 2017a). The LAT has also amassed a large sample of hadronic cosmic rays that can be studied using similar techniques to those in previous analyses.

The LAT has three subsystems (Atwood et al. 2009; Ackermann et al. 2012): an anticoincidence detector (ACD) to reject the charged particle background, a tracker (TKR) to promote conversion of gamma-rays to $e^+/e^-$ and measure their incident direction, and an electromagnetic calorimeter (CAL) to measure the energy of the resulting particle shower. The anticoincidence detector consists of 89 segmented plastic scintillator panels covering the top and sides of the LAT. In cosmic-ray analyses, the ACD can be used to measure the charge of incident particles via their ionization losses through the scintillator tiles. The tracker is composed of 18 layers of x-y silicon strip detectors (SSDs) with interwoven tungsten foils to promote conversion of gamma-rays into $e^+/e^-$ pairs. We used the tracker for direction reconstruction as well as a second, independent measurement of cosmic-ray charge. The calorimeter sits at the bottom of the LAT and consists of 1536 CsI (Tl) crystal bricks in a hodoscopic arrangement, which allows for 3D reconstruction of the particle shower and is crucial for lepton-hadron separation. The imaging capabilities of the calorimeter also provide an independent, complementary direction measurement, which we used to improve the angular resolution of the data set.

While the LAT's reconstruction algorithms are optimized for gamma-rays, the same basic principles for reconstructing particle direction and energy also apply to leptons and hadrons, though some care must be taken in understanding biases when measuring hadrons. The largest difference between reconstructing electromagnetic and hadronic showers is in the energy estimate, because the calorimeter is relatively shallow and does not fully contain hadronic showers. The calorimeter is 8.6 radiation lengths deep on axis, but only ~0.5 hadronic interaction lengths (Ackermann et al. 2012). In general, the calorimeter contains the electromagnetic portion of the particle shower and underestimates the energy of hadronic showers. We apply a scaling relation developed with a Geant4 (Agostinelli et al. 2003) Monte Carlo (MC) simulation of protons interacting with the detector to account for the missing energy in the reconstruction and remove the bias when estimating proton energy.

### 3. Data Set

In this section, we describe the event selection and resulting data set used in the analysis. We selected eight years of data from 2008 December 26 to 2016 December 26 for this analysis. To ensure at least the minimum quality necessary to reconstruct events, a set of basic cuts was applied: we required that the events pass the onboard gamma filter,<fnref n="57" /> a track must be found, events must traverse a minimum of $4X_0$ of calorimeter integrated along their path length, and events must deposit a minimum of 21 GeV in the calorimeter.<fnref n="58" /> We also cut on a

classification tree-based variable trained to identify accurately reconstructed events. These quality cuts are based on the event selection developed for the Fermi-LAT proton spectral measurement (Green 2016; Green & Hays 2017). Additionally, we required that the LAT was in standard survey mode<fnref n="59" /> and that the rocking angle of the spacecraft was <52°.<fnref n="60" />

To separate protons from helium nuclei and other heavy cosmic rays, we applied cuts on two independent measurements of the cosmic-ray charge, Z, using the tracker and the ACD. The signals in each subsystem are correlated with the charge due to the $Z^2$-scaling of ionization loss. In the ACD, we measured the energy deposited in the single scintillator tile that is intersected by the best-fit track. A geometric correction was applied to account for the path length of each track. The average pulse height in the TKR provides a second, independent measurement of the charge. A two-dimensional cut on both variables was used to separate $Z = 1$ particles from $Z > 1$ particles. An additional cut on the energy deposited in ACD tiles within a 15° cone of the best-fit track is also applied to remove residual heavy nuclei, most of which have large incidence angles or enter the bottom of the detector. Residual contamination from helium and other $Z > 1$ nuclei is estimated to be well below 1% (Green & Hays 2017).

The charge-based selection described above yields a data set of protons and electrons. We used the electron classifier developed for the Fermi-LAT Pass 8 cosmic-ray electron and positron spectral and anisotropy analyses to separate protons from the remaining electrons (Abdollahi et al. 2017a, 2017b). The multivariate classifier uses the differences in the morphology of electromagnetic and hadronic showers to separate the two event types. For example, the transverse width of hadronic showers is wider on average than that of electromagnetic showers. We refer to the most recent LAT CRE spectral measurement for more details on the classifier (Abdollahi et al. 2017b).

To mitigate the effects of the geomagnetic field, we imposed energy-dependent off-axis angle ($\theta$) cuts to reduce the LAT's field of view. This removes cosmic rays coming from the horizon, which have larger deflection angles than those with incident directions closer to zenith and are subjected to charge-dependent Earth shadowing. The maximum allowed off-axis angle is 45° for protons with energies between 78 and 139 GeV, and 50° for all other events. Details of this selection are described in Section 5.1.1.

A detailed understanding of the instrument's point-spread function (PSF) is critical for a measurement of anisotropy. Compared to gamma-rays, the LAT's angular resolution for protons is excellent: the 68% containment angle is ~0.01° above 78 GeV. However, MC studies of the effect of the tails of the PSF on the measured anisotropy showed that a small fraction of events with large reconstruction errors can generate a false-positive signal. We used a second, independent measurement of the event direction from the 3D imaging calorimeter to reject poorly reconstructed events, primarily those that entered the bottom of the detector. We required that the angle between the best-fit track in the TKR and the best-fit direction in the CAL be <0.2 radians (11°5). The procedure for determining this value without biasing the measurement is described in detail in Section 5.1.2.

<footnote n="57">⁵⁷ All showering cosmic rays should pass at these energies, because any event with >20 GeV of raw energy (i.e., uncalibrated) deposited in the CAL passes the gamma filter.</footnote>

<footnote n="58">⁵⁸ The raw 20 GeV threshold used by the onboard filter is uncalibrated and affected by decreased light yields in the CAL crystals over time. The 21 GeV threshold on a calibrated quantity mitigates the time dependence of the onboard filter. Details of this effect are discussed in Section 5.1.4.</footnote>

<footnote n="59">⁵⁹ LAT_MODE==5, LAT_CONFIG==1, DATA QUAL==1, which are publicly available in the LAT pointing history (FT2) files.</footnote>

<footnote n="60">⁶⁰ This is to ensure that the Earth limb was well outside the field of view.</footnote>
# 4. Anisotropy Search Method

## 4.1. Reference Map

To achieve the sensitivity necessary to measure anisotropy at these energies, the instrument's exposure cannot be calculated using simulation, because the resulting uncertainties are significantly larger than the expected signal $O(10^{-3})$. The observed counts map is instead compared to a reference map, which is an estimate of the instrument's response to an isotropic cosmic-ray flux. Many data-driven methods of creating reference maps have been developed over the years to avoid relying on MC simulation. The method we adopted uses the time-averaged event rate, $R_{\text{avg}}$, and distribution of detected event directions in instrument coordinates, $P(\theta, \phi)$, as empirical estimates of the detector's efficiency. In order to generate the reference map, the pointing history of the instrument was divided into 1 s time bins. For each 1 s bin with live time $lt_{\text{bin}}$, the expected number of events was drawn from a Poisson distribution with mean $R_{\text{avg}} \times lt_{\text{bin}}$ and event directions were randomly drawn from $P(\theta, \phi)$. The sky directions for each event were then calculated for the position and orientation of the instrument at that time. To ensure that any anisotropy signal in the rate and $P(\theta, \phi)$ was adequately averaged out, we averaged over long time intervals of one year, i.e., $R_{\text{avg}}$ and $P(\theta, \phi)$ were calculated for each of the eight years in the data set. The choice of an integer number of years also minimizes contamination from the Compton–Getting dipole created by Earth's motion around the Sun, which cancels out in each complete year (Compton & Getting 1935). Additionally, the use of year-long time bins, rather than a single bin for the entire data set, mitigates the effects of decreasing light yields in the calorimeter, which cause a small, monotonic decrease in the total event rate over the course of the LAT mission (Bregeon et al. 2013). One hundred independent reference maps were created and averaged to decrease the statistical uncertainty in the resulting eight-year map.

We compared the rate-based reference map method to the time-scrambling (Abbasi et al. 2010, 2011, 2012; Aartsen et al. 2013, 2016) or event-shuffling (Ackermann et al. 2010; Abdollahi et al. 2017a) method commonly used for cosmic-ray anisotropy measurements. Using the time-scrambling method, the reference map is created by randomly shuffling the times of events in the data set and calculating new sky directions for each event with the shuffled time. We performed an MC study to compare the performance of our rate-based method to the time-scrambling method when measuring a dipole using an ideal, conical detector sensitive at off-axis angles up to 60°. One million events were injected with a dipolar angular distribution with an amplitude of 0.01 and its maximum oriented at an angle $\alpha \in [0^\circ, 45^\circ, 90^\circ]$ relative to the North Celestial Pole. A null hypothesis data set was also produced by generating events isotropically. Reference maps were created using both the rate-based and time-scrambling methods and the dipole amplitude was reconstructed for each method. The results of one thousand realizations of this study are summarized in Figure 1. The time-scrambling method consistently underestimates the true dipole amplitude by a factor of ~2, while the rate-based method is unbiased. This is consistent with what was seen in the most recent Fermi-LAT anisotropy study of CREs, which compared four separate reference map techniques (Abdollahi et al. 2017a). A similar effect emerges with ground-based observatories at middle

latitudes, and iterative likelihood methods have been developed to un bias the measurement (Ahlers et al. 2016; Abeysekara et al. 2018). This bias occurs for any anisotropy with an angular scale larger than, or comparable to, the detector's field of view if the reference map method preserves the event rate within the detector's field of view. This results in a reference map that preserves some of the anisotropy, and therefore attenuates the anisotropy signal in the difference between the observed map and reference map. In both panels of Figure 1, the distribution of measured dipole amplitudes under the null hypothesis is centered around a non-zero value. This is the expected dipole amplitude due to Poisson noise in the data set, and it determines the sensitivity of the analysis. Furthermore, neither method shows bias with respect to the direction of the dipole, i.e., the angle $\alpha$, indicating that both methods are sensitive to the R.A. and decl. components of a dipole anisotropy. The rate-based reference map method is therefore sensitive to the two-dimensional direction of dipole anisotropy.

## 4.2. Spherical Harmonic Analysis

After we generated the reference map, we compared it to the data map and performed a spherical harmonic analysis. Spherical harmonics offer a convenient way to characterize anisotropy at multiple angular scales. First, we calculated the relative intensity between the data map and the reference map:

$$ \delta I_i(\theta, \phi) = \frac{N_i}{\langle N \rangle_i} - 1, \quad (1) $$

where $N_i$ and $\langle N \rangle_i$ are the counts in the $i$th pixel of the data map and reference map, respectively.

The relative intensity was then decomposed into spherical harmonics:

$$ \delta I(\theta, \phi) = \sum_{\ell=0}^{\infty} \sum_{m=-\ell}^{\ell} a_{\ell m} Y_{\ell m}(\theta, \phi). \quad (2) $$

The coefficients of the expansion were computed using the anafast algorithm in the HEALPix$^{61}$ library (Gorski et al. 2005). The coefficients were then converted to coefficients of the real-valued, or tesseral, spherical harmonic functions for a more natural geometric interpretation. All of the sky maps were binned using the HEALPix scheme with an Nside parameter of 16 (3072 pixels), which corresponds to an angular extent of approximately 4°.

The coefficients of the spherical harmonic decomposition can be used to characterize any anisotropy at each angular scale. The angular power at each multipole is calculated directly from the coefficients of the multipole expansion:

$$ C_{\ell} = \frac{1}{2\ell + 1} \sum_{m=-\ell}^{\ell} |a_{\ell m}|^2. \quad (3) $$

The measured angular power $C_{\ell}$ contains contributions from two terms: the true anisotropy power in the map, $C_{\ell}^{\text{true}}$, which we are interested in measuring, and a noise term from randomly correlated statistical fluctuations in the map, $C_N$. For an isotropic sky, the variance of the $a_{\ell m}$ is the angular power spectrum $C_{\ell}$, which has an expected value of $C_N$ for all $\ell$. The value of $C_N$ can be calculated by propagating errors in the relative intensity map to the final quantities (Knox 1995). The

<footnote>$^{61}$ http://healpix.sf.net</footnote>
following expression accounts for pixel-to-pixel variation introduced by nonuniform exposure in the maps (Fornasa et al. 2016; Abdollahi et al. 2017a):

$$C_N = \frac{4\pi}{N_{\text{pix}}^2} \sum_{i=1}^{N_{\text{pix}}} \left[ \frac{N_i}{\langle N \rangle_i^2} + \alpha \frac{N_i^2}{\langle N \rangle_i^3} \right], \quad (4)$$

where $N_{\text{pix}}$ is the number of pixels in the map, $\alpha = 1/n$, and $n$ is the number of reference maps created and averaged to estimate the isotropic expectation, which is 100 in this analysis. Note that, in the case of uniform exposure and very large $n$, the formula simplifies to $C_N = \frac{4\pi}{N}$, where $N$ is the total number of events in the map.<fnref n="62" /> To estimate the true anisotropy power in the map, we subtract the noise contribution, which results in the maximum likelihood estimator for the true anisotropy:

$$\hat{C}_{\ell}^{\text{true}} = C_{\ell} - C_{N}. \quad (5)$$

The angular power spectrum characterizes the total anisotropy at each angular scale, with angular features at each multipole, $\ell$, ~180°/ℓ. Any excess in the angular power spectrum compared to the isotropic expectation indicates anisotropy at that angular scale. The expected distribution of angular power at each multipole under the null hypothesis can be calculated by assuming that $C_\ell$ follows a $\chi^2_{2\ell+1}$ distribution with mean $C_N$ (Knox 1995).

As described in the introduction, the dipole anisotropy is especially interesting scientifically. The amplitude of the dipole is defined as:

$$\delta = 3 \sqrt{\frac{C_1}{4\pi}}, \quad (6)$$

where $C_1$ is the angular power at multipole $\ell = 1$. The estimator for the true dipole amplitude, $\hat{\delta}$, is calculated by inputting $\hat{C}_1^{\text{true}}$ into the above equation and imposing the requirement that $\hat{C}_1^{\text{true}} > 0$.

In addition to the total dipole amplitude, we can also calculate the full two-dimensional direction from the spherical harmonic coefficients. The R.A. of the maximum of the dipole is given by:

$$\text{R.A.} = \arctan\left(\frac{a_1 - 1}{a_{11}}\right), \quad (7)$$

while the decl. of the maximum is given by:

$$\text{Decl.} = \frac{\pi}{2} - \arccos\left(\sqrt{\frac{3}{4\pi}} \frac{a_{10}}{\delta}\right). \quad (8)$$

The accompanying uncertainty in each quantity above is calculated by propagating the statistical uncertainty in the relative intensity map to the derived quantities. Note that, under the null hypothesis, the $a_{\ell m s}$s should be normally distributed with a mean of zero.

## 5. Results

In this section, we present the results of the analysis performed using the methods described in the previous section. The data were divided into eight logarithmically spaced energy

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Results from an MC study comparing two reference map methods using an ideal detector that has a conical field of view with a maximum off-axis angle of 60°. Each realization consists of one million events injected with a 1% dipole oriented at an angle $\alpha$ relative to the North Celestial Pole. Reference maps were created using the time-scrambling method (top) and the rate-based method (bottom). The null hypothesis refers to realizations with events injected isotropically.</figcaption></figure>

bins in estimated gamma-ray energy, which correspond to approximately 78 GeV–10 TeV in estimated proton energy. Reference maps were generated for each energy bin using an averaging period, or time bin, of one year as described in Section 4.1. The data maps and reference maps were summed over all eight years of the data set, and summed cumulatively in energy to maximize the sensitivity of the analysis. The data map and reference map for the minimum cumulative energy bin, i.e., the energy bin spanning the full data set, can be seen in Figure 2. The structure seen in both the data and reference maps is a result of the LAT’s exposure, which is biased toward the northern and southern celestial poles due to the instrument’s rocking profile. There is also left-right asymmetry in the exposure due to a modified survey strategy to study the Galactic center that began in 2013 December and lasted one year. Note that the large pixel size results in 58,000 events per pixel, on average. The relative intensity map, which is used for the spherical harmonic analysis, and a significance map, are shown in Figure 3.

<footnote n="62">⁶² Uniform exposure refers to the case where $\langle N_i \rangle = \langle N_j \rangle$ for all i, j. In the large statistics limit, i.e., $N_i = \langle N_i \rangle$ and $\langle N_j \rangle = N/N_{\text{pix}}$, and with uniform exposure, Equation (5) results in this expression.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Data and reference sky maps in equatorial coordinates (J2000) for all events in the data set. Sky maps are created using the HEALPix pixelization scheme with 3072 pixels and use the astronomical convention of R.A. increasing to the left.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Relative intensity and significance maps in equatorial coordinates (J2000) for all events in the data set. Statistical fluctuations are smaller toward the equatorial poles because the exposure is greater toward the poles. Sky maps are created using the HEALPix pixelization scheme with 3072 pixels and use the astronomical convention of R.A. increasing to the left.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Measured $a_{1m}$ coefficients of the spherical harmonic expansion vs. minimum energy. From left to right: $a_{1-1}$, $a_{10}$, and $a_{11}$. Note that the $a_{10}$ coefficient corresponds to the dipole component aligned with the polar axis, which is unconstrained by ground-based measurements. The error bars are the $1\sigma$ statistical uncertainty on the measured coefficients and the colored bands represent the distribution of expected results under the null hypothesis, i.e., isotropic sky.</figcaption></figure>

We performed the spherical harmonic analysis detailed in Section 4.2 on the relative intensity maps corresponding to each of the cumulative energy bins. The measured $a_{\ell m}$ coefficients, which contain the direction of the large-scale anisotropy, are plotted in Figure 4 and their measured values are given in Table 1. Note that measured coefficients are not independent, due to the cumulative energy binning, which could create the energy dependence seen in Figure 4. The angular power spectrum for the entire energy range, i.e., $E > 78$ GeV, is shown in Figure 5. We calculated the angular power spectrum up to $\ell = 30$, which corresponds to an angular scale of $\sim 6^\circ$. There is a weak excess in the dipole power, $C_1$, with a $p$-value of 0.01 (pretrials). There is an

<table><caption>**Table 1**
Spherical Harmonic Coefficients</caption><thead><tr><th>Min. Energy (GeV)</th><th>a<sub>1-1</sub> [10<sup>-3</sup>]</th><th>a<sub>10</sub> [10<sup>-3</sup>]</th><th>a<sub>11</sub> [10<sup>-3</sup>]</th><th>σ<sub>stat</sub> [10<sup>-3</sup>]</th></tr></thead><tbody><tr><td>78</td><td>-0.39</td><td>-0.61</td><td>-0.57</td><td>0.28</td></tr><tr><td>139</td><td>-0.27</td><td>-0.16</td><td>-0.32</td><td>0.36</td></tr><tr><td>251</td><td>0.37</td><td>0.72</td><td>0.26</td><td>0.56</td></tr><tr><td>455</td><td>1.29</td><td>1.34</td><td>-0.95</td><td>0.92</td></tr><tr><td>830</td><td>0.83</td><td>0.65</td><td>0.90</td><td>1.53</td></tr><tr><td>1522</td><td>2.33</td><td>-2.13</td><td>2.47</td><td>2.58</td></tr><tr><td>2810</td><td>1.95</td><td>-1.87</td><td>8.68</td><td>4.58</td></tr><tr><td>5218</td><td>-4.78</td><td>-8.40</td><td>17.66</td><td>9.29</td></tr></tbody></table>

Note. There is a maximum energy of $\sim 10$ TeV. The uncertainties are the $1\sigma$ statistical errors on the measured coefficients.
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Angular power spectrum calculated for the relative intensity map with energy greater than 78 GeV. The horizontal axis is the multipole *ℓ* of the spherical harmonic expansion and the vertical axis is angular power at that *ℓ* value. The angular scale of each multipole is ~180°/*ℓ*. Here, *C*<sub>*ℓ*</sub> is the angular power in the map, which includes an anisotropy signal component and a noise component, *C*<sub>*N*</sub>. The angular power due to Poisson noise is subtracted from the measured power, which is the maximum likelihood estimator for the anisotropy signal included in the map. The error bars are the statistical uncertainty on the measured noise-subtracted angular power. The colored bands represent the distribution of expected results under the null hypothesis, i.e., isotropic sky.</figcaption></figure>

additional excess at *ℓ* = 16 with a *p*-value of 0.008 (pretrials); however, this is likely to be a statistical fluctuation. Each of the 30 angular power measurements is independent. Under the null hypothesis, the expected value for the number of measurements that should exceed the 95% interval due to random chance is 1.5. We note that either excess could be a statistical fluctuation, though it is unlikely that an astrophysical anisotropy would be measured as an excess at large *ℓ* without additional deviations from isotropy at intermediate angular scales. The known astrophysical anisotropy measured by other detectors has the most power at large angular scales, i.e., small *ℓ*, and the angular power spectrum is smoothly decreasing as *ℓ* increases (Abeysekara et al. 2019). The large-scale anisotropy measured by other detectors is near the sensitivity of our analysis, and there are systematics that could potentially create a false-positive dipole excess that will be discussed in Section 5.1. The interpretation of the *ℓ* = 1, i.e., dipole, measurement is therefore more nuanced than that at *ℓ* = 16 and will be discussed in detail below.

Dipole amplitudes for each cumulative energy bin were measured using the dipole power ($C_1$) calculated for each bin and can be seen in Figure 6. The dipole excess at $E > 78$ GeV corresponds to the excess seen in the angular power spectrum in Figure 5. The measured dipole amplitudes in the remaining energy intervals are all consistent with an isotropic sky. The exact amplitude and direction of the dipole excess are shown in Table 2.

In Section 5.1, we discuss the major sources of systematic uncertainty in the analysis that could lead to a dipole excess. We do not expect any of the systematics to create an excess at the level seen in our data. However, we cannot completely rule out the possibility that the excess is caused by an unaccounted systematic effect. We therefore computed upper

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Top: measured dipole amplitude, $\delta$, for each cumulative energy bin. The dipole amplitude is calculated directly from the measured dipole power ($C_1$). The error bars are the statistical uncertainty on the measured quantities, and the colored bands represent the distribution of expected results under the null hypothesis, i.e., isotropic sky. Note that the dipole amplitude in this plot is not noise-subtracted. Bottom: 95% confidence level upper limits on the dipole amplitude for each cumulative energy bin. The dashed line represents the expected upper limits under the null hypothesis. The colored bands show the central 68% and 95% range of upper limits expected from an ensemble of measurements under the null hypothesis.</figcaption></figure>

<table><caption>**Table 2**

Observed Amplitude and Direction of the Maximum of the Dipole Excess and 95% CL Upper Limit</caption><thead><tr><th>Min. Energy (GeV)</th><th>δ<sub>obs</sub> [10<sup>-4</sup>]</th><th>R.A. (°)</th><th>Decl. (°)</th><th>δ<sup>95%</sup><sub>UL</sub> [10<sup>-3</sup>]</th></tr></thead><tbody><tr><td>78</td><td>3.9 ± 1.5</td><td>215 ± 23</td><td>-51 ± 21</td><td>1.3</td></tr></tbody></table>

limits on the total dipole amplitude, for all cumulative energy bins, which can be seen in Figure 6. The 95% CL upper limits were calculated using the frequentist likelihood ratio approach used in Abdollahi et al. (2017a). The upper limits on the observed dipole power are calculated after enforcing that $\hat{C}_{\ell}^{\text{true}} > 0$ and then converted to upper limits on the dipole amplitude using Equation (6). The 95% CL upper limit on the dipole amplitude at a minimum energy of 78 GeV is $\delta_{\text{UL}} = 1.3 \times 10^{-3}$. This calculated upper limit is
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Comparison of a histogram of the significance in each pixel of the sky map in the right panel of Figure 3 to a Gaussian distribution with a mean of zero and variance of unity. The $\chi^2$ between the histogram and the Gaussian model was calculated and is summarized in the legend of the plot, showing that the data are consistent with a Gaussian distribution.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Sky maps in altitude-azimuth coordinates for events with a maximum off-axis angle cut of 45° (top) and 60° (bottom). Altitude is the complement of the zenith angle and azimuth is measured from $N = 0°$ and increases toward the E. The “east-west” effect is clearly visible in the data when a larger field of view is used.</figcaption></figure>

considerably larger than the median expected upper limit of $5.8 \times 10^{-4}$ because it was calculated from an observed data point with an excess.

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Results of simulations of an eight-year observation with decreasing tracker-calorimeter angle thresholds. Top: the quadrupole power ($C_2$) is plotted vs. the maximum tracker-calorimeter angle used in the event selection. A significant quadrupole excess is detected at larger thresholds due to events from the PSF tail. Bottom: the $\ell = 2$ coefficients of the spherical harmonic transform are plotted vs. the maximum tracker-calorimeter angle used in the event selection. All of the excess quadrupole power seen in the plot on the left is in the $m = 0$ moment, i.e., $a_{20}$.</figcaption></figure>

We note that the angular power spectrum only contains all of relevant information about the anisotropy if the random field on the sphere is Gaussian and statistically isotropic.<fnref n="63" /> This may not be the case for cosmic rays, where correlations in the spherical harmonic coefficients, i.e., $a_{\ell m}$, may occur due to magnetic effects (Ahlers & Mertsch 2017). As a separate, assumption-free test, we binned the significance in each pixel of the sky map in the right panel of Figure 3 and compared the resulting histogram to a Gaussian distribution with a mean of zero and variance of unity. The $\chi^2$ between the data and Gaussian model was calculated to perform a quantitative test of their consistency. The results, which are summarized in Figure 7, show that the sky map is consistent with a Gaussian distribution and that no significant features exist at the scale of the pixel size in the map ($Nside = 16$). We also performed

<footnote n="63">⁶³ Statistical isotropy means that the statistical properties of the random field on the sphere are invariant under rotations.</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Results of the scan over the tracker-calorimeter angle parameter using flight data. We scanned maximum tracker-calorimeter angle thresholds and ran the resulting data sets through the full anisotropy pipeline. To not bias the measurement of the other moments of the spherical harmonic analysis, we only used the $a_{20}$ in the scan, remaining blind to the other spherical harmonic parameters.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 11.** Third-degree polynomial fits of the relative event rate in four energy bins. The eight analysis energy bins are combined to yield similar event statistics in each combined bin. The relative rate is the total eight-year event rate in each McIlwain L bin divided by the average rate across all McIlwain L bins.</figcaption></figure>

the test with *Nside* values of 32 and 64, which yielded similar results.

## 5.1. Systematics

In the following subsections, we describe the major sources of systematic uncertainty in the anisotropy measurement and the techniques we used to mitigate or quantify them in a data-driven way. We describe the “east-west” effect seen in our data, as well as the selection employed to reduce the tails of the PSF. Finally, we discuss the stability of the event rate and its effect on the results.

### 5.1.1. East-West Effect

Although the cosmic rays in this data set are well above the vertical cutoff rigidity, geomagnetic effects are not completely

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Angular power spectra from an MC study of the McIlwain L-dependent event rate for events with energy greater than 78 GeV. Ten simulated data sets were generated using the McIlwain L-dependent rates in Figure 11 and processed with the anisotropy pipeline. The data points are the mean angular power at each multipole from the ten realizations. The error bars, which are hardly visible, are the standard error on the mean. The dashed lines represent the $1\sigma$ spread among the ten realizations. The colored bands represent the distribution of expected results under the null hypothesis, i.e., isotropic sky.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 13.** Comparison of the amplitude ($A_1$) and phase ($\Phi_1$) of the one-dimensional, i.e., R.A.-projected, dipole term of the harmonic expansion measured by Fermi-LAT above 78 GeV compared to previous measurements by underground muon telescopes (Bercovitch & Agrawal 1981; Thambyahpillai 1983; Swinson & Nagashima 1985; Ueno et al. 1990; Mori et al. 1995; Munakata et al. 1995). The amplitude and phase were calculated from the $a_{1-1}$ and $a_{11}$ coefficients of the spherical harmonic expansion using Equations (5) and (8) from Abeysekara et al. (2018). The Fermi measurement and upper limit from this work are both shown, for completeness. The upper limit was calculated by projecting the 95% C.L. upper limit on the total dipole amplitude onto the R.A. axis using the best-fit direction from Table 2.</figcaption></figure>

negligible. Positively charged cosmic rays arriving from near the horizon from the east are preferentially blocked by the Earth because their trajectories bend downward into the atmosphere. This famous “east-west” effect is visible in our data set if it is not accounted for. We impose energy-dependent off-axis angle (instrument theta) thresholds to mitigate this effect, which are similar to those used in the most recent Fermi-LAT $e^+/e^-$ anisotropy search (Abdollahi et al. 2017a). In the aforementioned study, the off-axis angle thresholds were
optimized using simulations that estimated the geomagnetic influence by back-tracing events through a model of the geomagnetic field. However, this method relies on the accuracy of the energy estimation to properly trace the particle trajectories. The relatively poor energy resolution for protons in the LAT does not allow for accurate back-tracing to quantify the effect to the desired precision. Therefore, the selection was determined by analyzing the data in altitude-azimuth coordinates where the geomagnetic effects are maximal. We created reference maps in differential energy bins for a range of maximum off-axis angle cuts in 5° increments. We then analyzed significance maps in altitude-azimuth coordinates and performed a $\chi^2$ comparison between the significance distribution and a standard normal distribution, requiring reduced $\chi^2 \sim 1$. The analysis resulted in a maximum off-axis angle of 45° in the lowest energy bin and 50° for all others. Significance maps for the final selection and a selection where the east-west effect is visible in the lowest energy bin are shown in Figure 8. The east-west effect is clearly visible in the lower panel, where there is a >6σ deficit of cosmic rays from the east measured in many pixels, and a corresponding excess from the west. After reducing the field of view (top panel), there is no longer any visible anisotropy.

### 5.1.2. Point-spread Function Tail

Preliminary results of this analysis from 2017 measured a significant quadrupole excess (Meehan et al. 2017). Further exploration indicated that the excess was caused by events from the tail of the PSF with large reconstruction errors. As stated in Section 3, the angular resolution for protons detected by the LAT is ~0.01. However, the tail of the distribution is non-negligible for an analysis at this level of sensitivity. The tail extends out to an angular error of 180°, which primarily consists of events that entered the bottom of the detector but were reconstructed as if they entered the top.<fnref n="64" /> In this final version of the analysis, we used an additional variable to remove events from the PSF tail. We compared independent direction measurements of each event from the tracker and calorimeter, and used the angle between them (hereafter called the “tracker-calorimeter angle”) as a proxy for the quality of the reconstruction. We created a detailed, Geant4-based simulation of our data set that includes realistic detector effects to accurately determine the tracker-calorimeter angle threshold necessary to reduce contamination from poorly reconstructed events. Figure 9 shows the quadrupole power versus maximum tracker-calorimeter angle when thresholds were applied to this simulated data set and the rest of the anisotropy analysis was applied. There is a significant quadrupole excess for large maximum threshold; it is reduced for smaller thresholds, yielding the expected behavior. Additionally, the lower panel of Figure 9 shows that the entirety of the excess is in the $a_{20}$ coefficient of the multipole expansion. This is the moment in which exposure-related systematics are expected to exist. Recall that the LAT’s exposure is primarily quadrupolar and aligned with the celestial poles (see Figure 2). Events that are reconstructed close to 180° off of their true direction are likely to pile up at the equatorial poles, creating an $a_{20}$ excess. These mistaken reconstructions are not accounted for in the reference map algorithm because they would need to be quantified with 0.1% precision, which is beyond the capability of the detector

simulation used to study this effect. We instead opted to remove poorly reconstructed events from the data set using the method described below.

While the simulation described above qualitatively explains the quadrupole, data/MC agreement was not precise enough to use the MC to determine an appropriate cut value. To tune the cut while staying unbiased to the parameters of interest, namely the dipole components of the analysis, we performed a parameter scan similar to that above, but only observed the $a_{20}$ moment of the spherical harmonic expansion. Figure 10 shows the results of this parameter scan on flight data. The $a_{20}$ component is very significant for large maximum tracker-calorimeter angle thresholds and decreases as the threshold decreases. The final cut of 0.2 radians is indicated by the dashed line in the plot, and the observed $a_{20}$ value is consistent with isotropy. In principle, this method could mask a real quadrupolar anisotropy, because the cut was tuned on the data. However, it biases only the measurement of the $a_{20}$ component of the anisotropy. Any real quadrupolar anisotropy in the four remaining components would be unaffected, though their contribution to the total quadrupole power depends on the orientation of the anisotropy. Most importantly, the dipole measurement remains unbiased by this method.

### 5.1.3. Event Rate Stability

The algorithm to construct the reference maps for the anisotropy search described in Section 4.1 uses the event rate averaged over each year of the data set to generate events. While the rate fluctuates on short timescales, the assumption is that it is stable over the entire year and also does not exhibit variation that is correlated with sky direction. Variation in the rate that is not properly accounted for in the reference map algorithm could over- or underestimate the exposure and result in a false-positive anisotropy. We tested both of these assumptions, and describe our findings in this section.

The geomagnetic cutoff varies continuously along the orbital path of Fermi. The increased rate of background cosmic rays in regions where the cutoff is lower could affect the efficiency of the various subsystems. To test this, we measured the event rate as a function of McIlwain L value of the LAT’s location, which is directly correlated with the geomagnetic cutoff. The observed event rate can vary by as much as 1% over the range of McIlwain L experienced by Fermi. To estimate the effect of this variation on the measured anisotropy, we simulated ten data sets by first fitting the rate versus McIlwain L curves in four energy bins and then generating events according to these curves.<fnref n="65" /> Figure 11 shows the observed event rate and third-degree polynomial fits for each of the four energy bins. As a full-circle test, we then processed the simulated data sets with the anisotropy search pipeline and measured the angular power spectrum for each realization. The results of this simulation are summarized in Figure 12. The mean dipole power among the realizations was ~1σ above the expectation under the null hypothesis. In other words, we expect a systematic 1σ dipole excess due to the McIlwain L-dependent event rate. We therefore conclude that this effect is not likely to create the observed dipole excess.

<footnote n="64">⁶⁴ The LAT has no simple capability to distinguish whether an event entered the top or bottom of the detector, from timing.</footnote>

<footnote n="65">⁶⁵ The exact curves used to seed this MC study were fit to data using the event selection from Meehan et al. (2017), which is slightly different than the one described here. The shape and magnitude of the McIlwain L dependence in the MC study is statistically consistent with that in the current data set.</footnote>
### 5.1.4. Raw Energy Threshold

As described in Section 3, the onboard filter uses a trigger that passes any event that deposits a minimum of 20 GeV of raw energy in the calorimeter to filter for high-energy particles. Radiation damage degrades the light yields of the calorimeter crystals over time. The raw energy used by the onboard filter is not calibrated for this decrease and the effective energy to pass the filter thus increases over time. This results in a time-dependent event rate for events with energies near the threshold. Without correcting for this effect, we observe a 0.5–1%/yr decrease in the event rate for events with energy in the range 78–139 GeV. To mitigate this effect, we imposed a threshold of 21 GeV of calibrated energy deposited in the calorimeter, i.e., energy calibrated for the decreasing light yields. The 21 GeV threshold is above the effective threshold of the uncalibrated onboard filter energy. This mitigates most of the time-dependent effect, resulting in a total event rate that decreases by only ~0.1%/yr. We performed simulation studies similar to those described in Section 5.1.3, to test the effect of this time dependence on the observed anisotropy. The study predicts that the systematic uncertainty in the anisotropy introduced by the monotonically decreasing rate is negligible compared to both the statistical uncertainty in the data set and the systematic effect of the McIlwain L variation, and therefore will not affect the results.

### 5.1.5. Stability of the Angular Distribution

In addition to the event rate, the observed angular distribution of events is assumed to be stable on long timescales and with changing geographic location to estimate the reference sky maps. In principle, the time-dependent raw energy threshold described in Section 5.1.4 could create a time-dependent incidence angle distribution, because the raw energy deposited in the calorimeter is a function of the event’s path length. However, we do not measure any significant time dependence in this distribution. Additionally, we searched for McIlwain L dependence of the two-dimensional angular distribution of events in detector coordinates. The distribution does not show any significant variation correlated with McIlwain L, and is not expected to introduce any systematic uncertainties into the measurement.

## 6. Conclusion

We have presented the results of the first search for cosmic-ray proton anisotropy using data from the Fermi LAT. The eight-year data set is the largest single-instrument, full-sky data set studying cosmic-ray anisotropy at these energies to date. It provides the most stringent constraints on the decl. dependence of the dipole anisotropy, which is not accessible by ground-based observatories.

Interpretation of the measured dipole excess is difficult due to its marginal statistical significance. We discussed three potential sources of systematics and our method for quantifying or mitigating them: the “east–west” effect, poorly reconstructed events from the PSF tail, and the McIlwain L-dependent event rate due to the varying rate of background cosmic rays. Of these, only the last is expected to have a measurable effect on the dipole anisotropy. Our simulation study described in Section 5.1.3 predicts a $1\sigma$ dipole excess due to this effect, but this is insufficient to explain the observed excess. Similar

geomagnetic effects were seen in cosmic-ray anisotropy searches by the Alpha Magnetic Spectrometer (AMS-02) experiment, where a method was developed to correct for the systematic shift in the measured anisotropy (Bindel et al. 2017; Gebauer et al. 2017). However, it is important to note that the orbit of AMS-02, which is onboard the International Space Station (ISS), has an inclination of 51°6 (see the LAT’s orbital inclination of 25°6) and travels through geographic locations with lower rigidity cutoffs than the LAT does. Geomagnetic effects are therefore expected to be less significant for the LAT anisotropy search.

Residual contamination from other particles could also introduce systematic uncertainties into the measurement. Estimates from the Geant4 MC simulation place the contamination from CREs in this data set at less than 0.1%, and the upper limit on the dipole anisotropy of CREs measured by Fermi-LAT is ~3 × 10⁻³ (Abdollahi et al. 2017a). The upper limit on the anisotropy created by CRE contamination in this data set is therefore ~3 × 10⁻⁶, well below the sensitivity of this analysis. Proton energies in this analysis were estimated by rescaling the estimated gamma-ray energy to account for the missing portion of the hadronic shower in the calorimeter. This will also rescale the accurately estimated energies of CREs by a factor of ~3, thereby introducing low-energy CRE contamination into the data set. Low-energy CREs could be affected by the heliospheric magnetic field, but it is difficult to quantify the effect on the anisotropy measurement. There is also residual contamination, estimated to be less than 1%, from Helium nuclei (Section 3). As a cross-check, we performed the anisotropy analysis on a selection of helium nuclei, which yielded a null result, i.e., consistent with isotropy. We therefore conclude that no systematic uncertainty should be introduced into the proton anisotropy measurement by residual, isotropic helium in the data set.

The statistical excess of the dipole amplitude is in a regime that makes it difficult to make a strong statement about its interpretation. The measured dipole can be described in terms of its amplitude and two-dimensional direction in equatorial coordinates: ($\delta$, R.A., decl.) = $3.9 \pm 1.5 \times 10^{-4}$, $215^\circ \pm 23^\circ$, $-51^\circ \pm 21^\circ$. Previous measurements of cosmic-ray anisotropy in the 100 GeV energy range by underground muon telescopes observed dipole amplitudes ~2 × 10⁻⁴ with maxima at right ascensions in the range $\phi_1 \in [45^\circ, 135^\circ]$ (Swinson & Nagashima 1985; Hall et al. 1999), where $\phi_1$ is the phase of the one-dimensional dipole fit typically performed by ground-based experiments. The phase of the TeV anisotropy described in Section 1 is similar to that measured by underground muon telescopes and is typically in the range $\phi_1 \in [30^\circ, 50^\circ]$. For the sake of comparison, we calculated the one-dimensional amplitude and phase of the first harmonic of the anisotropy using the cumulative Fermi data set. The resulting amplitude and phase are shown compared to previous measurements from underground muon telescopes in Figure 13. We note that the direction of the dipole measured in this analysis is in tension with the measurements by muon telescopes from decades ago, but stress the many differences between the analyses. These telescopes typically scanned a small patch of overhead sky and recorded the daily sidereal variation in the counting rate as a function of R.A., while the Fermi-LAT analysis measures the all-sky, two-dimensional anisotropy. Additionally, ground-based experiments have poor composition resolution compared to the LAT, and they measure the all-particle anisotropy, rather
than the measurement of protons only as presented here.<fnref n="66" /> The myriad differences between these measurements are important to understand to elucidate the origin of the anisotropy. Given the lack of a definitive signal, we set upper limits on the total dipole amplitude: the 95% CL upper limit on the dipole amplitude is $δ_{UL} = 1.3 \times 10^{-3}$ for protons with energy greater than 78 GeV and $δ_{UL} = 1.2 \times 10^{-3}$ for protons with energy greater than 251 GeV. Recently, a nearly all-sky measurement of cosmic-ray anisotropy was performed at ~10 TeV by combining data sets from the IceCube and HAWC detectors (Abeysekara et al. 2019), resulting in the least biased measurement of the anisotropy to date. However, this measurement is still insensitive to the decl. component, due to the limitations of ground-based measurements mentioned in Section 1. Our upper limits are therefore the most constraining on the decl. component of the anisotropy by any experiment.

The Fermi LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT as well as scientific data analysis. These include: the National Aeronautics and Space Administration and the Department of Energy in the United States; the Commissariat à l'Energie Atomique and the Centre National de la Recherche Scientifique/Institut National de Physique Nucléaire et de Physique des Particules in France; the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy; the Ministry of Education, Culture, Sports, Science, and Technology (MEXT), the High-energy Accelerator Research Organization (KEK), and the Japan Aerospace Exploration Agency (JAXA) in Japan; and the K. A. Wallenberg Foundation, the Swedish Research Council, and the Swedish National Space Board in Sweden.

Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d'Études Spatiales in France. This work performed in part under DOE Contract DE-AC02-76SF00515.

We would like to thank the INFN GRID Data Centers of Pisa, Trieste and CNAF-Bologna, the DOE SLAC National Accelerator Laboratory Computing Division, and the CNRS/IN2P3 Computing Center (CC-IN2P3 - Lyon/Villeurbanne), in partnership with CEA/DSM/Irfu, for their strong support in performing the massive simulations necessary for this work.

## ORCID iDs

M. Ajello https://orcid.org/0000-0002-6584-1703

T. J. Brandt https://orcid.org/0000-0003-4087-1786

J. W. Hewitt https://orcid.org/0000-0001-5254-2248

M. Kuss https://orcid.org/0000-0003-1212-9998

I. Liodakis https://orcid.org/0000-0001-9200-4006

S. Maldera https://orcid.org/0000-0002-0698-4421

N. Omodei https://orcid.org/0000-0002-5448-7577

V. S. Paliya https://orcid.org/0000-0001-7774-5308

M. Pesce-Rollins https://orcid.org/0000-0003-1790-8018

## References

Aartsen, M. G., Abbasi, R., Abdou, Y., et al. 2013, ApJ, 765, 55

Aartsen, M. G., Abraham, K., Ackermann, M., et al. 2016, ApJ, 826, 220

Abbasi, R., Abdou, Y., Abu-Zayyad, T., et al. 2010, ApJL, 718, L194

Abbasi, R., Abdou, Y., Abu-Zayyad, T., et al. 2011, ApJ, 740, 16

Abbasi, R., Abdou, Y., Abu-Zayyad, T., et al. 2012, ApJ, 746, 33

Abdo, A. A., Allen, B., Aune, T., et al. 2008, PhRvL, 101, 221101

Abdo, A. A., Allen, B. T., Aune, T., et al. 2009, ApJ, 698, 2121

Abdollahi, S., Ackermann, M., Ajello, M., et al. 2017a, PhRvL, 118, 091103

Abdollahi, S., Ackermann, M., Ajello, M., et al. 2017b, PhRvD, 95, 082007

Abeysekara, A. U., Alfaro, R., Alvarez, C., et al. 2014, ApJ, 796, 108

Abeysekara, A. U., Alfaro, R., Alvarez, C., et al. 2018, ApJ, 865, 57

Abeysekara, A. U., Alvarez, C., Álvarez, J. D., et al. 2019, ApJ, 871, 96

Ackermann, M., Ajello, M., Albert, A., et al. 2012, ApJS, 203, 4

Ackermann, M., Ajello, M., Atwood, W. B., et al. 2010, PhRvD, 82, 092003

Adriani, O., Barbarino, G. C., Bazilevskaya, G. A., et al. 2011, Sci, 332, 69

Aglietta, M., Alekseenko, V. V., Alessandro, B., et al. 2009, ApJL, 692, L130

Agostinelli, S., Allison, J., Amako, K., et al. 2003, NIMPA, 506, 250

Aguilar, M., Aisa, D., Alpat, B., et al. 2015, PhRvL, 114, 171103

Ahlers, M. 2016, PhRvL, 117, 151103

Ahlers, M., BenZvi, S. Y., Desiati, P., et al. 2016, ApJ, 823, 10

Ahlers, M., & Mertsch, P. 2017, PhPNP, 94, 184

Amenomori, M., Ayabe, S., Bi, X. J., et al. 2006, Sci, 314, 439

Amenomori, M., Ayabe, S., Cui, S. W., et al. 2005, ApJL, 626, L29

Amenomori, M., Bi, X. J., Chen, D., et al. 2017, ApJ, 836, 153

Atwood, W. B., Abdo, A. A., Ackermann, M., et al. 2009, ApJ, 697, 1071

Bartoli, B., Bernardini, P., Bi, X. J., et al. 2013, PhRvD, 88, 082001

Bartoli, B., Bernardini, P., Bi, X. J., et al. 2015, ApJ, 809, 90

Bercovitch, M., & Agrawal, S. 1981, Proc. ICRC (Essonne), 10, 246

Bindel, F. K., Gebauer, I., Zeissler, S., & Graziani, M. 2017, Proc. ICRC (Busan), 35, 202

Blasi, P., & Amato, E. 2012, JCAP, 2012, 011

Bregeon, J., Charles, E., Wood, M. & Fermi-LAT Collaboration 2013, arXiv:1304.5456

Compton, A. H., & Getting, I. A. 1935, PhRv, 47, 817

Fornasa, M., Cuoco, A., Zavala, J., et al. 2016, PhRvD, 94, 123005

Gebauer, I., AMS Collaboration, Bindel, K. F., et al. 2017, Proc. ICRC (Busan), 35, 186

Gorski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759

Green, D. M. 2016, PhD thesis, Univ. Maryland, doi:10.13016/M2ZK06

Green, D. M., & Hays, E. A. 2017, Proc. ICRC (Busan), 35, 159

Hall, D. L., Munakata, K., Yasue, S., et al. 1999, JGRA, 104, 6737

Hillas, A. M. 2005, JPhG, 31, R95

Knox, L. 1995, PhRvD, 52, 4307

Meehan, M., Vandenbroucke, J. & Fermi-LAT Collaboration 2017, Proc. ICRC (Busan), 35, 170

Mertsch, P., & Funk, S. 2015, PhRvL, 114, 021101

Mori, S., Yasue, S., Munakata, K., et al. 1995, Proc. ICRC (Rome), 4, 648

Munakata, K., Yasue, S., Mori, S., et al. 1995, Proc. ICRC (Rome), 4, 639

Ptuskin, V. S., Jones, F. C., Seo, E. S., & Sina, R. 2006, AdSpR, 37, 1909

Schwadron, N. A., Frisch, P., Adams, F. C., et al. 2015, ASTRP, 2, 9

Swinson, D. B., & Nagashima, K. 1985, P&SS, 33, 1069

Thambyahpillai, T. 1983, Proc. ICRC (Bangalore), 3, 383

The Pierre Auger Collaboration, Aab, A., Abreu, P., et al. 2017, Sci, 357, 1266

Ueno, H., Fujii, Z., & Yamada, T. 1990, Proc. ICRC (Adelaide), 6, 361

<footnote n="66">⁶⁶ Note that the chemical composition of the data sets used to measure cosmic-ray anisotropy is not necessarily representative of the true chemical abundance, due to differences in detector efficiency for each cosmic-ray species. For example, the estimated fractional composition of the HAWC data set used in Abeysekara et al. (2019) is ~1/3 He to 2/3 protons by particle number at a median energy of 10 TeV.</footnote>