Search for Very High-Energy Emission from the millisecond pulsar PSR J0218+4232

V. A. ACCIARI,<fnref n="1" /> S. ANSOLDI,<fnref n="2" /> L. A. ANTONELLI,<fnref n="3" /> A. ARBET ENGELS,<fnref n="4" /> M. ARTERO,<fnref n="5" /> K. ASANO,⁶

D. BAACK,⁷ A. BABIĆ,⁸ A. BAQUERO,⁹ U. BARRES DE ALMEIDA,¹⁰ J. A. BARRIO,⁹ I. BATKOVIĆ,¹¹

J. BECERRA GONZÁLEZ,¹² W. BEDNAREK,¹³ L. BELLIZZI,¹⁴ E. BERNARDINI,¹⁵ M. BERNARDOS,¹¹ A. BERTI,¹⁶

J. BESENRIEDER,¹⁷ W. BHATTACHARYYA,¹⁵ C. BIGONGIARI,<fnref n="3" /> A. BILAND,<fnref n="4" /> O. BLANCH,¹⁸ G. BONNOLI,¹⁴

Ž. BOŠNJAK,⁸ G. BUSETTO,¹¹ R. CAROSI,¹⁹ G. CERIBELLA,¹⁷ M. CERRUTI,²⁰ Y. CHAI,¹⁷

A. CHILINGARIAN,²¹ S. CIKOTA,⁸ S. M. COLAK,¹⁸ E. COLOMBO,¹² J. L. CONTRERAS,⁹ J. CORTINA,²²

S. COVINO,<fnref n="3" /> G. D'AMICO,¹⁷ V. D'ELIA,<fnref n="3" /> P. DA VELA,¹⁹,²³ F. DAZZI,<fnref n="3" /> A. DE ANGELIS,¹¹

B. DE LOTTO,<fnref n="2" /> M. DELFINO,¹⁸,²⁴ J. DELGADO,¹⁸,²⁴ C. DELGADO MENDEZ,²² D. DEPAOLI,¹⁶

F. DI PIERRO,¹⁶ L. DI VENERE,²⁵ E. DO SOUTO ESPINIIRA,¹⁸ D. DOMINIS PRESTER,²⁶ A. DONINI,<fnref n="2" />

D. DORNER,²⁷ M. DORO,¹¹ D. ELSAESSER,⁷ V. FALLAH RAMAZANI,²⁸,²⁹ A. FATTORINI,<fnref n="2" />,⁷ G. FERRARA,<fnref n="3" />

M. V. FonSECA,⁹ L. Font,³⁰ C. FRUCK,¹⁷ S. FUKAMI,⁶ R. J. GARCÍA LÓPEZ,¹² M. GARCZARCYK,¹⁵

S. GASPARYAN,³¹ M. GAUG,³⁰ N. GIGLIETTO,²⁵ F. GIORDANO,²⁵ P. GLIWNY,¹³ N. GODINOVIĆ,³²

J. G. GREEN,<fnref n="3" /> D. GREEN,¹⁷ D. HADASCH,⁶ A. HAHN,¹⁷ L. HECKMANN,¹⁷ J. HERRERA,¹² J. HOANG,⁹

D. HRUPEC,³³ M. HÜTTEN,¹⁷ T. INADA,⁶ S. INOUE,³⁴ K. ISHIO,¹⁷ Y. IWAMURA,⁶ I. JIMÉNEZ,²²

J. JORMANAINEN,²⁸ L. JOUVIN,¹⁸ Y. KAJIWARA,³⁵ M. KARJALAINEN,¹² D. KERSZBERG,¹⁸ Y. KOBAYASHI,⁶

H. KUBO,³⁵ J. KUSHIDA,³⁶ A. LAMASTRA,<fnref n="3" /> D. LELAS,³² F. LEONE,<fnref n="3" /> E. LINDFORS,²⁸ S. LOMBARDI,<fnref n="3" />

F. LONGO,<fnref n="2" />,³⁷ R. LÓPEZ-COTO,¹¹ M. LÓPEZ-MOYA,⁹ A. LÓPEZ-ORAMAS,¹² S. LOPORCHIO,²⁵

B. MACHADO DE OLIVEIRA FRAGA,¹⁰ C. MAGGIO,³⁰ P. MAJUMBAR,³⁸ M. MAKARIEV,³⁹ M. MALLAMACI,¹¹

G. MANEVA,³⁹ M. MANGANARO,²⁶ K. MANNHEIM,²⁷ L. MARASCHI,<fnref n="3" /> M. MARIOTTI,¹¹ M. MARTÍNEZ,¹⁸

D. MAZIN,⁶,¹⁷ S. MENCHIARI,¹⁴ S. MENDER,⁷ S. MIĆANOVIĆ,²⁶ D. MICELI,<fnref n="2" /> T. MIENER,⁹ M. MINEV,³⁹

J. M. MIRANDA,¹⁴ R. MIRZOYAN,¹⁷ E. MOLINA,²⁰ A. MORALEJO,¹⁸ D. MORCUENDE,⁹ V. MORENO,³⁰

E. MORETTI,¹⁸ V. NEUSTROEV,⁴⁰ C. NIGRO,¹⁸ K. NILSSON,²⁸ K. NISHIJIMA,³⁶ K. NODA,⁶

S. NOZAKI,³⁵ Y. OHTANI,⁶ T. OKA,³⁵ J. OTERO-SANTOS,¹² S. PAIANO,<fnref n="3" /> M. PALATIELLO,<fnref n="2" />

D. PANEQUE,¹⁷ R. PAOLETTI,¹⁴ J. M. PAREDES,²⁰ L. PAVLETIĆ,²⁶ P. PEÑIL,⁹ C. PERENNES,¹¹

M. PERSIC,<fnref n="2" />,⁴¹ P. G. PRADA MORONI,¹⁹ E. PRANDINI,¹¹ C. PRIYADARSHI,¹⁸ I. PULJAK,³² W. RHODE,⁷

M. RIBÓ,²⁰ J. RICO,¹⁸ C. RIGHI,<fnref n="3" /> A. RUGLIANCICH,¹⁹ L. SAHA,⁹ N. SAHAKYAN,³¹ T. SAITO,⁶

S. SAKURAI,⁶ K. SATALECKA,¹⁵ F. G. SATURNI,<fnref n="3" /> B. SCHLEICHER,²⁷ K. SCHMIDT,⁷ T. SCHWEIZER,¹⁷

J. SITAREK,¹³ I. ŠNIDARIĆ,⁴² D. SOBCZYNSKA,¹³ A. SPOLON,¹¹ A. STAMERRA,<fnref n="3" /> D. STROM,¹⁷ M. STRZYS,⁶

Y. SUDA,¹⁷ T. SURIĆ,⁴² M. TAKAHASHI,⁶ F. TAVECCHIO,<fnref n="3" /> P. TEMNIKOV,³⁹ T. TERZIĆ,²⁶ M. TESHIMA,¹⁷,⁶

L. TOSTI,⁴³ S. TRUZZI,¹⁴ A. TUTONE,<fnref n="3" /> S. UBACH,³⁰ J. VAN SCHERPENBERG,¹⁷ G. VANZO,¹²

M. VAZQUEZ ACOSTA,¹² S. VENTURA,¹⁴ V. VERGUILOV,³⁹ C. F. VIGORITO,¹⁶ V. VITALE,⁴⁴ I. VOVK,⁶

M. WILL,¹⁷ C. WUNDERLICH,¹⁴ D. ZARIĆ,³²

(MAGIC COLLABORATION)

P. A. CARAVEO,⁴⁵ I. COGNARD,⁴⁶,⁴⁷ L. GUILLEMOT,⁴⁶,⁴⁷ A. K. HARDING,⁴⁸ J. LI,⁴⁹,⁵⁰ B. LIMYANSKY,⁵¹

C. Y. Ng,⁵²,⁵³ D. F. TORRES,⁵⁴,⁵⁵ AND P. M. SAZ PARKINSON,⁵¹,⁵²,⁵³

<footnote n="1">¹Inst. de Astrofísica de Canarias, E-38200 La Laguna, and Universidad de La Laguna, Dpto. Astrofísica, E-38206 La Laguna, Tenerife, Spain</footnote>

<footnote n="2">²Università di Udine and INFN Trieste, I-33100 Udine, Italy</footnote>

<footnote n="3">³National Institute for Astrophysics (INAF), I-00136 Rome, Italy</footnote>

<footnote n="4">⁴ETH Zürich, CH-8093 Zürich, Switzerland</footnote>

<footnote n="5">⁵Institut de Física d'Altes Energies (IFAE), The Barcelona Institute of Science and Technology (BIST), E-08193 Bellaterra (Barcelona), Spain</footnote>
⁶Japanese MAGIC Group: Institute for Cosmic Ray Research (ICRR), The University of Tokyo, Kashiwa, 277-8582 Chiba, Japan

⁷Technische Universität Dortmund, D-44221 Dortmund, Germany

⁸Croatian MAGIC Group: University of Zagreb, Faculty of Electrical Engineering and Computing (FER), 10000 Zagreb, Croatia

⁹IPARCOS Institute and EMFTEL Department, Universidad Complutense de Madrid, E-28040 Madrid, Spain

¹⁰Centro Brasileiro de Pesquisas Físicas (CBPF), 22290-180 URCA, Rio de Janeiro (RJ), Brazil

¹¹Università di Padova and INFN, I-35131 Padova, Italy

¹²Inst. de Astrofísica de Canarias, E-38200 La Laguna, and Universidad de La Laguna, Dpto. Astrofísica, E-38206 La Laguna, Tenerife, Spain

¹³University of Lodz, Faculty of Physics and Applied Informatics, Department of Astrophysics, 90-236 Lodz, Poland

¹⁴Università di Siena and INFN Pisa, I-53100 Siena, Italy

¹⁵Deutsches Elektronen-Synchrotron (DESY), D-15738 Zeuthen, Germany

¹⁶INFN MAGIC Group: INFN Sezione di Torino and Università degli Studi di Torino, 10125 Torino, Italy

¹⁷Max-Planck-Institut für Physik, D-80805 München, Germany

¹⁸Institut de Física d'Altes Energies (IFAE), The Barcelona Institute of Science and Technology (BIST), E-08193 Bellaterra (Barcelona), Spain

¹⁹Università di Pisa and INFN Pisa, I-56126 Pisa, Italy

²⁰Universitat de Barcelona, ICCUB, IEEC-UB, E-08028 Barcelona, Spain

²¹Armenian MAGIC Group: A. Alikhanyan National Science Laboratory

²²Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas, E-28040 Madrid, Spain

²³now at University of Innsbruck

²⁴also at Port d'Informació Científica (PIC) E-08193 Bellaterra (Barcelona) Spain

²⁵INFN MAGIC Group: INFN Sezione di Bari and Dipartimento Interateneo di Fisica dell'Università e del Politecnico di Bari, 70125 Bari, Italy

²⁶Croatian MAGIC Group: University of Rijeka, Department of Physics, 51000 Rijeka, Croatia

²⁷Universität Würzburg, D-97074 Würzburg, Germany

²⁸Finnish MAGIC Group: Finnish Centre for Astronomy with ESO, University of Turku, FI-20014 Turku, Finland

²⁹now at Ruhr-Universität Bochum, Fakultät für Physik und Astronomie, Astronomisches Institut (AIRUB), 44801 Bochum, Germany

³⁰Departament de Física, and CERES-IEEC, Universitat Autònoma de Barcelona, E-08193 Bellaterra, Spain

³¹Armenian MAGIC Group: ICRANet-Armenia at NAS RA

³²Croatian MAGIC Group: University of Split, Faculty of Electrical Engineering, Mechanical Engineering and Naval Architecture (FESB), 21000 Split, Croatia

³³Croatian MAGIC Group: Josip Juraj Strossmayer University of Osijek, Department of Physics, 31000 Osijek, Croatia

³⁴Japanese MAGIC Group: RIKEN, Wako, Saitama 351-0198, Japan

³⁵Japanese MAGIC Group: Department of Physics, Kyoto University, 606-8502 Kyoto, Japan

³⁶Japanese MAGIC Group: Department of Physics, Tokai University, Hiratsuka, 259-1292 Kanagawa, Japan

³⁷also at Dipartimento di Fisica, Università di Trieste, I-34127 Trieste, Italy

³⁸Saha Institute of Nuclear Physics, HBNI, 1/AF Bidhannagar, Salt Lake, Sector-1, Kolkata 700064, India

³⁹Inst. for Nucl. Research and Nucl. Energy, Bulgarian Academy of Sciences, BG-1784 Sofia, Bulgaria

⁴⁰Finnish MAGIC Group: Astronomy Research Unit, University of Oulu, FI-90014 Oulu, Finland

⁴¹also at INAF Trieste and Dept. of Physics and Astronomy, University of Bologna

⁴²Croatian MAGIC Group: Ruđer Bošković Institute, 10000 Zagreb, Croatia

⁴³INFN MAGIC Group: INFN Sezione di Perugia, 06123 Perugia, Italy

⁴⁴INFN MAGIC Group: INFN Roma Tor Vergata, 00133 Roma, Italy

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

⁴⁶Laboratoire de Physique et Chimie de l'Environnement et de l'Espace – Université d'Orléans / CNRS, F-45071 Orléans Cedex 02, France

⁴⁷Station de radioastronomie de Nançay, Observatoire de Paris, CNRS/INSU, F-18330 Nançay, France

⁴⁸NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA

⁴⁹CAS Key Laboratory for Research in Galaxies and Cosmology, Department of Astronomy, University of Science and Technology of China, Hefei 230026, People's Republic of China

⁵⁰School of Astronomy and Space Science, University of Science and Technology of China, Hefei 230026, People's Republic of China

⁵¹Santa Cruz Institute for Particle Physics, Department of Physics, University of California at Santa Cruz, Santa Cruz, CA 95064, USA

⁵²Department of Physics, The University of Hong Kong, Pokfulam Road, Hong Kong, China

⁵³Laboratory for Space Research, The University of Hong Kong, Hong Kong, China

⁵⁴Institute of Space Sciences (ICE, CSIC), Campus UAB, Carrer de Magrans s/n, E-08193 Barcelona, Spain; and Institut d'Estudis Espacials de Catalunya (IEEC), E-08034 Barcelona, Spain

⁵⁵Institució Catalana de Recerca i Estudis Avançats (ICREA), E-08010 Barcelona, Spain
(Accepted to ApJ, 23 August 2021)

Abstract

PSR J0218+4232 is one of the most energetic millisecond pulsars known and has long been considered
as one of the best candidates for very high-energy (VHE; >100 GeV) γ-ray emission. Using 11.5 years
of Fermi Large Area Telescope (LAT) data between 100 MeV and 870 GeV, and ~90 hours of MAGIC
observations in the 20 GeV to 20 TeV range, we have searched for the highest energy γ-ray emission
from PSR J0218+4232. Based on the analysis of the LAT data, we find evidence for pulsed emission
above 25 GeV, but see no evidence for emission above 100 GeV (VHE) with MAGIC. We present
the results of searches for γ-ray emission, along with theoretical modeling, to interpret the lack of
VHE emission. We conclude that, based on the experimental observations and theoretical modeling,
it will remain extremely challenging to detect VHE emission from PSR J0218+4232 with the current
generation of Imaging Atmospheric Cherenkov Telescopes (IACTs), and maybe even with future ones,
such as the Cherenkov Telescope Array (CTA).

Keywords: Pulsars (1306) — Millisecond pulsars (1062) — Binary pulsars (153) — Gamma-ray astron-
omy (628) — Gamma-ray sources (633)

# 1. INTRODUCTION

PSR J0218+4232 (hereafter J0218) is a millisecond pulsar (MSP) with a period of 2.3 ms in a 2-day orbit with a ~0.2*M<sub>☉</sub> white dwarf companion (Bassa et al. 2003). It was serendipitously discovered as a steep spectrum, highly polarized source in a low-frequency radio study of an unrelated supernova (Navarro et al. 1995). Its broad radio peak, with a large unpulsed component (~50%) makes it unusual, suggesting that it may be an *aligned rotator*, a pulsar in which the magnetic field is aligned with the axis of rotation. This was later supported by subsequent polarimetry studies (Stairs et al. 1999). With a characteristic age τ <0.5 Gyr and a spindown power of 2.4×10<sup>35</sup> erg s<sup>-1</sup>, it is one of the youngest and most energetic MSPs known. It has an extremely strong magnetic field at the light cylinder (B<sub>LC</sub> ~ 3.2 × 10<sup>5</sup> G), only slightly weaker than young Crab-like pulsars, but several orders of magnitude weaker at the neutron star surface (Saito et al. 1997). Like the Crab Pulsar, J0218 also displays giant pulses (Joshi et al. 2004; Knight et al. 2006), something very uncommon among MSPs. Its distance, previously estimated to be greater than 6 kpc, potentially making it the most luminous MSP known (Du et al. 2014), has since been revised downwards to ~3 kpc (Verbiest & Lorimer 2014). Table 1 provides a summary of the key properties of J0218.

J0218 was first detected as a steady source of X-rays and γ-rays using ROSAT and EGRET, respectively (Verbunt et al. 1996). It was later shown to display non-thermal pulsed X-ray emission (Kuiper et al. 1998). Like in radio, J0218 has a large unpulsed component in X-rays. It has been detected by most of the X-ray missions, including BeppoSAX (Mineo et al. 2000), Chandra (Kuiper et al. 2002), XMM-Newton (Webb et al. 2004), RXTE (Kuiper et al. 2004), and most recently NICER (Deneva et al. 2019). Non-thermal hard X-ray emission was detected with NuSTAR up to ~70 keV (Gotthelf & Bogdanov 2017), with no evidence for a spectral break above these energies, although it must turn over somewhere between 100 keV and 100 MeV, to be consistent with the γ-ray emission detected at GeV energies (Gotthelf & Bogdanov 2017).

The γ-ray emission of J0218 has often been confused with that of the blazar 3C 66A because of their close proximity (0.97° separation, see Figure 1), the poor angular resolution of many γ-ray instruments, and the varying intensity of the blazar. The Second EGRET (2EG) Catalog (Thompson et al. 1995) contained a source (2EG J0220+4228) that was more often associated with 3C 66A. In the third EGRET (3EG) catalog (Hartman et al. 1999), the source 3EG J0222+4253 was identified with 3C 66A based on its >1 GeV emission, however its low-energy flux (100-300 MeV) was found to be dominated by the pulsar, rather than the blazar, leading to the conclusion that both were likely contributing to the EGRET source (Kuiper et al. 2000; Guillemot et al. 2007). Indeed, Kuiper et al. (2000) reported marginal evidence (~3.5σ) for the detection of pulsed γ-ray emission from J0218, making it potentially the first MSP detected at these energies.

With the launch of the Large Area Telescope (LAT, Atwood et al. 2009) on board the *Fermi* satellite, the picture became significantly clearer. Next to 3C 66A, a strong source (> 19σ), 1FGL J0218.1+4232, was reported in the First *Fermi* LAT Catalog, based on 11 months of data (Abdo et al. 2010a), and a firm detection of GeV γ-ray pulsations
<figure><img src="image_1.png" /><figcaption>**Figure 1.** LAT All-sky map showing the location of PSR J0218+4232 (aka 4FGL J0218.1+4232). The background map shows the 12-year (August 4, 2008 - August 4, 2020) all-sky intensity map, generated using γ-ray data above 1 GeV, in Galactic coordinates (Credit: Seth Digel). The square inset region shows the 15° × 15° counts map centered on PSR J0218+4232, generated with the 11.5 years of data analyzed in this paper (from 2008 August 4 to 2020 February 10), using all events above 1 GeV. Note the bright γ-ray blazar 3C 66A (aka 4FGL J0222.6+4302) located less than one degree away from PSR J0218+4232.</figcaption></figure>

confirmed the identification with J0218 (Abdo et al. 2009, 2010b). The LAT detection, however, showed a broad single-peaked γ-ray light curve, bearing little similarity to the EGRET one reported in Kuiper et al. (2000), or indeed the double-peaked X-ray pulse profile (Kuiper et al. 2002; Webb et al. 2004). Despite the fact that the First LAT Catalog of sources above 10 GeV (Ackermann et al. 2013, 1FHL) contained no source associated with J0218, this pulsar was identified as having hints of pulsed emission above 10 GeV (Ackermann et al. 2013), and the Third Catalog of Hard Fermi-LAT Sources (Ajello et al. 2017a, 3FHL) contained a source associated with J0218 (3FHL J0218.3+4230), which was shown to have >10 GeV (and even >25 GeV) pulsations (Saz Parkinson et al. 2017), despite the limited statistics due to the small effective area of the LAT at such high energies, compared to ground-based γ-ray telescopes<fnref n="1" />. These preliminary LAT results provided a strong motivation for observing this pulsar at even higher energies, using ground-based γ-ray telescopes.

At very high energies (VHE >100 GeV), J0218 has been a target for the Major Atmospheric Gamma Imaging Cherenkov (MAGIC) telescope, starting in 2004, when it was observed for 13 hours during the commissioning phase (single-telescope MAGIC-1 observations), in large part due to it being in the same field of view as 3C 66A (Oña-Wilhemi et al. 2005). It was subsequently observed for 20 hours, between October 2006 and January 2007, yielding no significant detection, and a $3\sigma$ flux upper limit of < 9.4 × 10⁻¹² cm⁻² s⁻<fnref n="1" />, above 140 GeV (Anderhub et al. 2010). Since then, the performance of the MAGIC telescopes has significantly improved (Aleksić et al. 2016a).

In this paper, we report results from an analysis of 11.5 years of Fermi-LAT data, together with ~90 hours of data from new MAGIC stereoscopic observations of J0218, collected from November 2018 to November 2019, using the low-energy threshold Sum-Trigger-II system (Dazzi et al. 2021). The structure of the paper is as follows: Section 2.1

<footnote n="1">¹ Note that the sensitivity of ground-based γ-ray telescopes like MAGIC depends not only on their large effective areas, but also on their ability to reject the cosmic-ray background. Given the challenges of performing background rejection with MAGIC in the 10–100 GeV range, it is perhaps not surprising that MAGIC is less sensitive than Fermi-LAT at these energies, despite its much larger effective area.</footnote>
describes the Fermi-LAT data analysis. Section 2.2 describes our MAGIC observations and analysis. Section 3 presents our results, both based on LAT data (Section 3.1) and the MAGIC observations (Section 3.2). Section 4 discusses our theoretical modeling. Finally, in Section 5 we discuss the main conclusions of our work.

## 2. DATA ANALYSIS

### 2.1. Fermi-LAT Data Analysis

For the Fermi-LAT analysis, we used 11.5 years of Pass 8 data (Atwood et al. 2013; Bruel et al. 2018) (specifically, P8R3_SOURCE_V2), from 2008 August 4 (MJD 54682.7) to 2020 February 10 (MJD 58890).<fnref n="2" /> We used Fermipy (Wood et al. 2017) to select Source class (evclass 128), Front and Back converting events (evtype 3) with an energy range from 100 MeV to 870 GeV, and from a square region of 15° × 15°, centered on the position of 4FGL J0218.1+4232 (RA=34.5344°, DEC=42.5459°). A maximum Earth zenith angle of 90° is imposed, helping to eliminate contamination from the Earth's limb. We further ensured that the selection only included events at times where the LAT was in normal science configuration and taking good data.

Figure 1 shows the Fermi-LAT >1 GeV all-sky γ-ray intensity map, highlighting the region around J0218. As discussed in Section 1, the blazar 3C 66A, located 0.97° away, complicates the analysis of J0218. The bright γ-ray counterpart of 3C 66A (4FGL J0222.6+4302, see Figure 1) has a > 1 GeV flux of $1.6 \times 10^{-8} \text{ cm}^{-2} \text{ s}^{-1}$ (> 160σ) in 4FGL, compared to $5.2 \times 10^{-9} \text{ cm}^{-2} \text{ s}^{-1}$ (73.5σ) for J0218. Because blazars are typically variable sources, we considered excluding certain *flaring* periods of 3C 66A. Unfortunately, the Fermi-LAT light curves of J0218 and 3C 66A, over the entire time interval (11.5 years) and energy (100 MeV to 870 GeV) range (Figure 2) reveal that such a strategy would not be possible as the blazar has been quite active throughout the entire 11.5-year period of our observations. Figure 3 shows the zoomed-in 1-year period covered by our MAGIC observations, illustrating how 3C 66A is brighter than J0218 at GeV energies. We note that there is a report of quasi-periodic variability in the optical light curve of 3C 66A with a period of ~3 years, although these have not been confirmed in the γ-ray data (Otero-Santos et al. 2020).

A binned likelihood analysis was performed, utilizing spatial bins of 0.1° × 0.1° and 8 logarithmically spaced bins per decade of energy. The initial region model was seeded from the 8-year Fermi-LAT Fourth Source Catalog (4FGL, Abdollahi et al. 2020) by including sources lying within a square area of 40° × 40° centered on 4FGL J0218.1+4232. Because the LAT only measures the energy of photons with finite resolution (<10% in the 1–100 GeV range, but ~20% at 100 MeV, and ~28% at 30 MeV)<fnref n="3" /> a correction for this energy dispersion was enabled for all sources except the isotropic diffuse emission, as recommended by the Fermi Science Support Center (FSSC)<fnref n="4" />. The fitting of source spectra within the region model was carried out in an iterative process, with sources being removed from the model if they fell below a test statistic (TS) of 10. Here, TS is defined as $TS = -2\ln(\mathcal{L}_0/\mathcal{L}_1)$, where $\mathcal{L}_1$ is the probability that a model which includes the source of interest matches the data, and $\mathcal{L}_0$ is the probability that the same model without the source of interest matches the data<fnref n="5" />. The final iteration of the fit was performed with Minuit (James 1994). In addition to the spectral parameters of 4FGL J0218.1+4232, the normalization and index of the Galactic diffuse emission (modeled with a power-law spectrum), the normalization of the isotropic diffuse emission, and the normalizations of background sources with a TS of at least 100 were allowed to vary.

We modeled the J0218 spectrum using a power law with an exponential cutoff <fnref n="6" />, $\frac{dN}{dE} = N_0(\frac{E}{E_0})\gamma \exp(-aE^b)$. We set the index (b) to a fixed “sub-exponential” value of 2/3, as this source is too faint for it to be determined by maximum likelihood estimation, and 2/3 approximates the values of other, brighter pulsars (Abdollahi et al. 2020).

After obtaining our best region model, we examined events within 5° of J0218, and used the “gtsrcprob” Fermitool to assign them a probability (*weight*) of originating from either J0218 or 3C 66A relative to other sources in the model. Finally, we used TEMPO2 (Hobbs et al. 2006) with the fermi plug-in (Ray et al. 2011) to assign the pulsar rotational phases $\phi_i$, according to our pulsar ephemeris obtained with the Nançay radio telescope, given in Table 1.

We calculated the source spectrum and flux points for the theoretical modeling described in Section 4 by utilizing the aforementioned region model. Three energy bins spanning 12.38 - 28.99 GeV were combined in order to produce

<footnote n="2">² We considered the possibility of including Calorimeter-Only (Cal-Only) data (Takahashi et al. 2015; Takahashi 2019) in our analysis. However, after a preliminary look at 8 years of such data, covering the period 2008-2016, using three different Cal-Only event classes, we found no evidence for pulsed emission, most likely due to the large PSF and corresponding large cosmic-ray background level. Thus, we opted to limit our analysis to standard LAT data.</footnote>

<footnote n="3">³ https://www.slac.stanford.edu/exp/glast/groups/canda/lat_Performance.htm</footnote>

<footnote n="4">⁴ This correction occurs during the fitting of the model. One additional energy bin is added beyond the existing upper and lower energy limits, and filled with a number of photons that is extrapolated from the existing model. The number of photons in each energy bin is then scaled by a factor relating to the instrument response functions, in order to account for the possibility photons with sufficiently large energy uncertainty may have a true energy that lies in a neighboring energy bin. See https://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Pass8_edisp_usage.html</footnote>

<footnote n="5">⁵ https://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Cicerone/Ciceroe_Likelihood/Likelihood_overview.html</footnote>

<footnote n="6">⁶ See https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/source_models.html</footnote>
<table><caption>**Table 1.** Timing ephemeris for PSR J0218+4232, obtained with the Nançay radio telescope. We used the DE436 Solar System ephemeris, with time units in barycentric dynamic time (TDB) and the ELL1 binary model for low eccentricity orbits, where EPS1 and EPS2 represent the first and second Laplace-Lagrange parameters (Lange et al. 2001). We refer the reader to the TEMPO2 manual (Hobbs et al. 2006) for the detailed definition of all parameters included in our timing model.</caption><thead><tr><th colspan="2">Timing and binary parameters</th></tr></thead><tbody><tr><td>R.A., α (J2000.0)</td><td>02<sup>h</sup> 18<sup>m</sup> 06.35863(1)<sup>s</sup></td></tr><tr><td>Decl., δ (J2000.0)</td><td>+42°32′17.3722(2)"</td></tr><tr><td>Frequency, F0 (Hz)</td><td>430.46105998103612106(6)</td></tr><tr><td>1st frequency derivative, F1, (Hz s<sup>-1</sup>)</td><td>-1.434128(1) × 10<sup>-14</sup></td></tr><tr><td>PMRA (α/cosδ, mas yr<sup>-1</sup>)</td><td>5.32(3)</td></tr><tr><td>PMDEC (δ, mas yr<sup>-1</sup>)</td><td>-3.68(6)</td></tr><tr><td>PEPOCH (MJD)</td><td>56000</td></tr><tr><td>POSEPOCH (MJD)</td><td>56000</td></tr><tr><td>DMEPOCH (MJD)</td><td>56000</td></tr><tr><td>DM (cm<sup>-3</sup> pc)</td><td>61.2374(7)</td></tr><tr><td>DM1 (cm<sup>-3</sup> pc s<sup>-1</sup>)</td><td>-0.0004(2)</td></tr><tr><td>BINARY MODEL</td><td>ELL1</td></tr><tr><td>PB (d)</td><td>2.0288460845(6)</td></tr><tr><td>A1 (lt-s)</td><td>1.9844348(2)</td></tr><tr><td>TASC (MJD)</td><td>49148.5799767(2)</td></tr><tr><td>EPS1</td><td>5.0(2)×10<sup>-6</sup></td></tr><tr><td>EPS2</td><td>4.9(2)×10<sup>-6</sup></td></tr><tr><td>START (MJD)</td><td>53579.2</td></tr><tr><td>FINISH (MJD)</td><td>58960.5</td></tr><tr><td>UNITS</td><td>TDB (Barycentric Dynamical Time)</td></tr><tr><td>EPHEM</td><td>DE436</td></tr></tbody><thead><tr><th colspan="2">Derived parameters</th></tr></thead><tbody><tr><td>Period, P (ms)</td><td>2.32309053</td></tr><tr><td>1st period derivative, Ṗ (s s<sup>-1</sup>)</td><td>7.739 × 10<sup>-20</sup></td></tr><tr><td>Characteristic age, τ<sub>c</sub> (yr)</td><td>4.8 × 10<sup>8</sup></td></tr><tr><td>Spin-down power, Ė (erg s<sup>-1</sup>)</td><td>2.4 × 10<sup>35</sup></td></tr><tr><td>Surface B-field strength, B<sub>S</sub> (G)</td><td>4.3 × 10<sup>8</sup></td></tr><tr><td>Light-cylinder B-field, B<sub>LC</sub> (G)</td><td>3.1 × 10<sup>5</sup></td></tr><tr><td>Distance, d (kpc)</td><td>3.15<sup>+0.85</sup><sub>-0.60</sub></td></tr><tr><td>ON pulse region</td><td>(0.34–0.98)</td></tr><tr><td>OFF pulse region</td><td>[0,0.34)∪(0.98,1]</td></tr></tbody></table>

<table><caption>**Table 2.** Gamma-ray spectral parameters for the total emission from PSR J0218+4232. Photon and energy flux cover the entire 100 MeV - 870 GeV energy range.</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>N<sub>0</sub> (ph cm<sup>-2</sup> s<sup>-1</sup> MeV<sup>-1</sup>)</td><td>(2.07 ± 0.03) × 10<sup>-11</sup></td></tr><tr><td>γ</td><td>-1.76 ± 0.01</td></tr><tr><td>E<sub>0</sub> (MeV)</td><td>821.6 (fixed)</td></tr><tr><td>a</td><td>(6.19751 ± 0.00007) × 10<sup>-3</sup></td></tr><tr><td>b</td><td>0.6667 (fixed)</td></tr><tr><td>Photon flux (photons cm<sup>-2</sup> s<sup>-1</sup>)</td><td>(7.67 ± 0.15) × 10<sup>-8</sup></td></tr><tr><td>Energy flux (MeV cm<sup>-2</sup> s<sup>-1</sup>)</td><td>(3.05 ± 0.04) × 10<sup>-5</sup></td></tr></tbody></table>
<figure><img src="image_2.png" /></figure>

**Figure 2.** Light curve showing PSR J0218+4232 (green circles) and 3C 66A (blue squares). The LAT data ranges from 2008 August 4 (MJD 54682.7) to 2020 February 10 (MJD 58890), and covers the 100 MeV - 870 GeV energy bands. The time period of MAGIC observations (MJD 58424 – 58791) is shown in cyan. Note the larger variability and γ-ray flux of 3C 66A. To generate this plot, background sources were fixed to the value in the region model, and the normalizations of 3C 66A and J0218 were allowed to vary.

a flux point instead of an upper limit. To extract the overall spectrum of J0218, the normalization of the isotropic and Galactic diffuse emission components were allowed to vary while Minuit fit the spectral parameters (summarized in Table 2). To generate flux points, the index of the Galactic diffuse emission was also allowed to vary, along with the normalizations of background sources with a TS of at least 500 or which lie within 5° of 4FGL J0218.1+4232. The spectrum of 4FGL J0218.1+4232 is replaced by a power law with an index of -2, and Minuit is used to fit the normalization of this modified spectrum within each energy bin. The result is interpreted as either a flux point or an upper limit, depending on the significance with which the power law source was detected.

## 2.2. *MAGIC Observations and Data Analysis*

We used the MAGIC telescopes to search for the VHE emission component of J0218. MAGIC consists of two imaging atmospheric Cherenkov telescopes (IACTs) of 17m diameter located at the Roque de Los Muchachos Observatory in La Palma, Canary Islands, Spain (Aleksić et al. 2016a). We collected the data in stereoscopic mode with the Sum-Trigger-II system. This system is designed to improve the performance of the telescopes in the sub-100 GeV energy range (Dazzi et al. 2021).

We observed the source from 2018 November 2 to 2019 November 4 (MJD 58424 – 58791) with a zenith angle range from 13 to 30 degrees, for maximum sensitivity at low-energies. Wobble mode was used for robust flux and
<figure><img src="image_3.png" /><figcaption>**Figure 3.** Fermi-LAT light curve showing PSR J0218+4232 (green circles) and 3C 66A (blue squares), zoomed in on the time period of MAGIC observations (MJD 58424 – 58791). Note that 3C 66A has a significantly larger flux than J0218 most of the time.</figcaption></figure>

background estimation by pointing the telescopes 0.4° away from the source (Fomin et al. 1994). Weather conditions were monitored simultaneously by measuring the atmospheric transmission with the LIDAR system operating together with the MAGIC Telescopes (Fruck et al. 2014). The Cherenkov radiation produced by the sub-100 GeV particle showers is fainter; consequently, they are more affected by the lower atmospheric transmissions. Therefore, we set a strict requirement of excellent atmospheric conditions, a minimum of 0.85 atmospheric transmission at an altitude of 9km. After discarding around 20% of the total data, 87 hours of good quality data remained.

The data were analyzed using the Magic Standard Analysis Software, MARS (Zanin et al. 2013). We applied the Sum-Trigger-II dedicated algorithm for calibration and image cleaning, which enabled us to improve the performance and achieve an energy threshold of 20 GeV. We performed the higher-level analysis following the standard pipeline (Aleksić et al. 2016b).

## 3. RESULTS

### 3.1. *Fermi-LAT results*

Figure 4 shows the LAT spectrum obtained in our analysis. Table 2 reports the best-fit spectral parameters and Table 3 gives the spectral values, including upper limits. Note that the spectrum falls steeply at energies above 10 GeV with only upper limits reported for energies above ~29 GeV. These results are consistent with previous LAT results reported in 4FGL (Abdollahi et al. 2020) and 3FHL (Ajello et al. 2017a), the latter of which reported an index of $\Gamma = -4.5$, when fitting the >10 GeV data with a simple power law.
3.1.1. *Search for high-energy pulsation in LAT data*

To test for possible pulsed emission above 10 GeV, an analysis analogous to what was carried out in the First Fermi-LAT Catalog of >10 GeV sources (1FHL) (Ackermann et al. 2013) and in Saz Parkinson et al. (2017) was performed. We defined a low-energy probability density function (PDF), which we refer to as PDF<sub>LE</sub>, based on the best estimate fit of the 1-10 GeV events (see Figure 5, top panel). For the high-energy PDF, we considered the family of distributions given by PDF<sub>HE</sub>(φ) = (1-x)+x·PDF<sub>LE</sub>(φ), with 0 ≤ x ≤ 1. We maximized the unbinned likelihood function derived from this PDF, with respect to x, obtaining L(ˆx), and comparing it to the null hypothesis, for x=0, that there is no pulsation (i.e. PDF<sub>HE</sub>(φ)=1). By construction, L(0) = 1, so the test statistic (TS=-2ln(L(0)/L(ˆx))) can be simplified to TS = 2 lnL(ˆx). We converted the measured TS value into a tail probability (or p-value), by assuming, following Wilks' theorem (Wilks 1938), that for the null hypothesis the TS follows a χ<sup>2</sup> distribution with 1 degree of freedom.

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

**Figure 4.** Spectrum of the total emission of millisecond pulsar J0218+4232 measured by Fermi-LAT (green points and upper limits) and the MAGIC telescopes (light blue upper limits). The green dashes represent the fit of the Fermi-LAT data with an exponentially cutoff power-law model. Note that the width of the error region is narrower than the dashes showing the best-fit model. Although included, it is difficult to distinguish in this plot. For the MAGIC analysis we assumed a spectral index Γ = -4.5 obtained from the spectral index of the power-law fit to the high-energy (>10 GeV) part of the Fermi-LAT spectrum.

Using the LAT data set described in Section 2.1, we first selected events in the 1–10 GeV energy range (with a probability >50% of coming from the pulsar, as obtained via the Fermi Science Tool<sup>7</sup> gtsrcprob) and used this histogram to generate a ‘low-energy template’. We used the non-parametric SOPIE (Sequential Off-Pulse Interval Estimation) R package (Schutte & Swanepoel 2016), to obtain a *smooth* kernel density estimator of our histogram (which we defined as our “low-energy template”), and also derived an estimate of the off-pulse interval of the light curve, using the median value of the results obtained from four different goodness-of-fit tests: Kolmogorov–Smirnov, Cramér–von Mises, Anderson–Darling, and Rayleigh test statistics. Figure 5 (top panel) shows our results, including the 1–10 GeV histogram, along with the resulting *low-energy* template, and the estimated off-pulse interval, all calculated using SOPIE.

To test for emission at higher energies, we looked at the >10 GeV events arriving within the 95% containment radius of the point-spread function (0.5/0.8 degrees for front/back converting events) and performed a likelihood test to

<footnote n="7">7 https://fermi.gsfc.nasa.gov/ssc/data/analysis/software/</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 5.** Top panel – Histogram of the 1–10 GeV events for PSR J0218+4232, along with the smooth circular kernel density estimator (red line) fitted to the data, which we define as the *low-energy* template in our subsequent searches for pulsed emission above 10 GeV. Two rotation cycles are shown, with 100 bins per cycle. The blue brackets indicate the estimated off-pulse interval, $[0-0.34) \cup (0.98,1]$, obtained using SOPIE (Schutte & Swanepoel 2016). Bottom panel – Search for high-energy pulsations using LAT standard events above 10 GeV. The blue histogram are the 58 events above 10 GeV in energy, while the pink histogram are the 17 events above 25 GeV. Two rotation cycles are shown, with 65 bins per cycle.</figcaption></figure>

determine whether they are likely to come from a similar distribution function, as represented by the lower energy template. We set a threshold p-value of 0.05 to claim evidence for emission at a specific energy. We carried out the same test with events of energies greater than 25 GeV. Figure 5 (bottom panel) shows the distribution of 58 (17) events above 10 (25) GeV, leading to a p-value of 1e-4 (0.01), thus showing evidence for emission above 10 and, marginally, above 25 GeV. We also tested for possible emission above 30 GeV but found that, despite the presence of 10 events above this energy, their distribution in phase yielded a p-value that was not significant ($p > 0.05$).

## 3.2. MAGIC Results

We analyzed our MAGIC data to search for possible pulsed and un-pulsed γ-ray emission above 20 GeV. The skymap is shown in Figure 6, where no emission is observed from J0218. The high emission spot observed in the image is the blazar 3C 66A, which is significantly detected as a by-product of the observations centered on J0218.
We used the TEMPO2 package (Hobbs et al. 2006) to assign the rotational phase to each event using the same ephemeris as in our LAT analysis described in Section 3.1 given in Table 1.

Given the broad pulse shape in the high-energy band, the use of the off-pulse region to estimate the background would lead to large uncertainties due to its smaller phase extent of 0.36 compared to the on-pulse width of 0.64. Therefore, we chose three source-free reflected-region backgrounds located at the same distance from the FoV center, which we expect to have the same acceptance as the region containing the source (Berge et al. 2007). The upper limits (ULs) to the differential flux were obtained by following the Rolke & López (2001) method under the assumption of a
<figure><img src="image_6.png" /></figure>

**Figure 6.** MAGIC skymap of the region around PSR J0218+4232 (indicated by a purple square) above 20 GeV. The relative flux (in arbitrary units) is calculated by the number of smeared excess events divided by the residual background flux within 0.1 degrees (Zanin et al. 2013)). Although no VHE emission is detected from J0218, the blazar 3C 66A (green cross), a well-known VHE source (Acciari et al. 2009; Aliu et al. 2009), is seen with high significance.

Gaussian systematic uncertainty in the detection efficiency, with a standard deviation of 30% systematic uncertainty in the flux level. Hereafter, the ULs will be given at 95% confidence level (CL). We assumed a spectral index Γ=-4.5 obtained from the power-law fit to the high-energy (>10 GeV) Fermi-LAT data, as reported in the Third Hard Source Catalog (Ajello et al. 2017b, 3FHL).

Figure 4 shows our MAGIC upper limits, indicated with blue arrows, along with the green points and upper limits from the Fermi-LAT analysis. The numerical values are reported in Table 3.

Given the soft spectrum observed by Fermi-LAT, we concentrated our search for pulsed signal in the lowest energy decade starting at threshold energy for the MAGIC observations. Hence, the phase-folded light curve of J0218 (see Figure 7) is computed in the energy range from 20 GeV to 200 GeV. We performed the same unbinned likelihood test described in Section 3.1, to determine whether they are likely to come from the distribution function represented by the lower energy template (see Figure 5, top panel), obtaining no evidence for pulsation (p-value ≫ 0.05). In addition, we chose to carry out a standard pulsation search, looking at ON and OFF events. The on-pulse region was selected as the phase interval between 0.34-0.98 and shown as the gold area in Figure 7, as defined by our Fermi-LAT analysis (see the top panel of Figure 5 and Table 1). The same source-free reflected-region backgrounds, shown with the grey horizontal band in Figure 7, were used for calculating the significance of the excess events using Eq.17 of (Li & Ma 1983)), and no significant (0.057σ) pulsation was found. Moreover, we applied region-independent signal tests (de Jager et al. 1989) (χ², and H-test), also with null results (5.54, for 11 degrees of freedom, and 0.05σ, respectively).
<table><caption>**Table 3.** Fermi-LAT and MAGIC spectral points and Upper Limits. Centers of energy bins are reported. Fermi-LAT data utilizes 32 logarithmically spaced bins between 100 MeV and 870 GeV. Three bins, spanning 12.38 - 28.99 GeV, were combined in order to produce a flux point instead of an upper limit. As such, a total of 30 bins are reported for the Fermi-LAT. MAGIC utilizes 14 logarithmically spaced bins between 20 GeV and 63 TeV. Note that we did not obtain upper limits for the last five MAGIC bins (i.e. E> 3.56 TeV) because they have zero counts and such limits would be considered too unreliable.</caption><thead><tr><th colspan="3">Fermi-LAT</th><th colspan="3">MAGIC</th></tr><tr><th>E<br>[GeV]</th><th>E<sup>2</sup>dN/(dEdAdt)</th><th>E<sup>2</sup>dN/(dEdAdt)</th><th>E<br>[GeV]</th><th>E<sup>2</sup>dN/(dEdAdt)</th><th>E<sup>2</sup>dN/(dEdAdt)</th></tr><tr><th></th><th>[TeV cm<sup>-2</sup> s<sup>-1</sup>]</th><th>[TeV cm<sup>-2</sup> s<sup>-1</sup>]</th><th></th><th>[TeV cm<sup>-2</sup> s<sup>-1</sup>]</th><th>[TeV cm<sup>-2</sup> s<sup>-1</sup>]</th></tr></thead><tbody><tr><td>0.12</td><td>(7.24 ± 1.58) × 10<sup>-12</sup></td><td>...</td><td>33.40</td><td>&lt; 3.79 × 10<sup>-13</sup></td><td>...</td></tr><tr><td>0.15</td><td>(7.69 ± 1.02) × 10<sup>-12</sup></td><td>...</td><td>44.35</td><td>&lt; 4.71 × 10<sup>-13</sup></td><td>...</td></tr><tr><td>0.20</td><td>(8.11 ± 0.73) × 10<sup>-12</sup></td><td>...</td><td>45.93</td><td>...</td><td>&lt; 8.32 × 10<sup>-12</sup></td></tr><tr><td>0.27</td><td>(7.77 ± 0.57) × 10<sup>-12</sup></td><td>...</td><td>58.88</td><td>&lt; 5.82 × 10<sup>-13</sup></td><td>...</td></tr><tr><td>0.36</td><td>(8.25 ± 0.48) × 10<sup>-12</sup></td><td>...</td><td>78.18</td><td>&lt; 7.69 × 10<sup>-13</sup></td><td>...</td></tr><tr><td>0.48</td><td>(8.15 ± 0.42) × 10<sup>-12</sup></td><td>...</td><td>81.68</td><td>...</td><td>&lt; 2.36 × 10<sup>-12</sup></td></tr><tr><td>0.63</td><td>(8.34 ± 0.38) × 10<sup>-12</sup></td><td>...</td><td>103.80</td><td>&lt; 1.33 × 10<sup>-12</sup></td><td>...</td></tr><tr><td>0.84</td><td>(7.50 ± 0.36) × 10<sup>-12</sup></td><td>...</td><td>137.80</td><td>&lt; 1.38 × 10<sup>-12</sup></td><td>...</td></tr><tr><td>1.11</td><td>(8.22 ± 0.37) × 10<sup>-12</sup></td><td>...</td><td>145.25</td><td>...</td><td>&lt; 6.57 × 10<sup>-13</sup></td></tr><tr><td>1.48</td><td>(6.79 ± 0.36) × 10<sup>-12</sup></td><td>...</td><td>183.00</td><td>&lt; 1.80 × 10<sup>-12</sup></td><td>...</td></tr><tr><td>1.96</td><td>(5.34 ± 0.35) × 10<sup>-12</sup></td><td>...</td><td>243.00</td><td>&lt; 2.40 × 10<sup>-12</sup></td><td>...</td></tr><tr><td>2.61</td><td>(5.37 ± 0.38) × 10<sup>-12</sup></td><td>...</td><td>258.30</td><td>...</td><td>&lt; 3.50 × 10<sup>-13</sup></td></tr><tr><td>3.46</td><td>(4.05 ± 0.38) × 10<sup>-12</sup></td><td>...</td><td>322.60</td><td>&lt; 3.19 × 10<sup>-12</sup></td><td>...</td></tr><tr><td>4.59</td><td>(3.07 ± 0.37) × 10<sup>-12</sup></td><td>...</td><td>428.30</td><td>&lt; 4.29 × 10<sup>-12</sup></td><td>...</td></tr><tr><td>6.10</td><td>(2.90 ± 0.40) × 10<sup>-12</sup></td><td>...</td><td>459.34</td><td>...</td><td>&lt; 1.94 × 10<sup>-13</sup></td></tr><tr><td>8.10</td><td>(1.28 ± 0.30) × 10<sup>-12</sup></td><td>...</td><td>568.70</td><td>&lt; 5.83 × 10<sup>-12</sup></td><td>...</td></tr><tr><td>10.75</td><td>(9.24 ± 3.01) × 10<sup>-13</sup></td><td>...</td><td>755.00</td><td>&lt; 8.10 × 10<sup>-12</sup></td><td>...</td></tr><tr><td>18.95</td><td>(4.66 ± 1.71) × 10<sup>-13</sup></td><td>...</td><td>816.84</td><td>...</td><td>&lt; 3.30 × 10<sup>-13</sup></td></tr><tr><td>25.83</td><td>...</td><td>&lt; 9.12 × 10<sup>-11</sup></td><td>1452.58</td><td>...</td><td>&lt; 1.06 × 10<sup>-13</sup></td></tr><tr><td></td><td></td><td></td><td>2583.09</td><td>...</td><td>&lt; 1.05 × 10<sup>-13</sup></td></tr></tbody></table>
<figure><img src="image_7.png" /><figcaption>**Figure 7.** Search for VHE pulsations using MAGIC events between 20 - 200 GeV, shown with the pink histogram. We used the same on-pulse interval as LAT analysis, [0.34-0.98], presented with the gold area. The grey horizontal dashed line within the one sigma uncertainty band indicates the average number of OFF events collected from three reflected-region backgrounds in the FoV. No significant pulsation is detected.</figcaption></figure>
## 4. THEORETICAL MODELING

We modeled the broadband spectrum of J0218 from UV to VHE γ-rays (14 orders of magnitude in energy) using a numerical force-free magnetosphere model for the global magnetic field, computing the individual trajectories of particles injected at the neutron star surface. Two populations of particles are injected: primary electrons/positrons along field lines that connect to the current sheet and are accelerated by an assumed parallel electric field distribution, and secondary electrons/positrons from polar cap pair cascades along field lines where there is no accelerating electric field. The dynamics and radiation of the particles are followed from the neutron star surface to a distance of 2 light cylinder radii ($2R_{lc}$) and radiated photons are stored in energy-dependent sky maps of observer angle vs. rotation phase (Harding & Kalapotharakos 2015; Harding et al. 2018a).

All particles radiate by synchro-curvature (SC) and inverse Compton (IC) emission. The pitch angles for SC are maintained through cyclotron resonant absorption of radio photons emitted above the polar cap (PC) (Harding et al. 2008). The SC also assumes the radius of curvature of the particle trajectory in the inertial observers' frame. The IC requires that trajectories be followed twice, once to store the SC radiation emissivity and another to compute the local photon densities from the stored emissivity and radiate IC (Harding & Kalapotharakos 2015; Harding et al. 2018a). The main assumptions of the model are the parallel electric field ($E_‖$) distribution, the source of pairs, pair multiplicity and their injection distribution on the PC, and the mechanism for generating pitch angle. Apart from this, the model requires the observed parameters of the pulsar ($P$ and $\dot{P}$). The magnetic and electric field distribution assumptions are based on results of Particle-In-Cell (PIC) simulations showing that pulsars producing high pair multiplicity have near-force-free magnetospheres and that the highest parallel electric fields are in the current sheet. Increasing the pair multiplicity increases the pair SC (mostly synchrotron radiation (SR)) and the pair IC (mostly synchrotron self-Compton, SSC). Increasing the E parallel increases the SC of primaries, increases the high-energy (GeV) cutoff, and increases the IC (at 10 TeV). For J0218, we assumed a magnetic inclination angle $\alpha = 45^\circ$, a viewing angle $\zeta = 65^\circ$, and a pair multiplicity of $M_+ = 1 \times 10^5$. Figure 8 shows the model predictions, including the various individual emission components.

We have also used a synchro-curvature model where all unknowns are reduced to just a few parameters that represent the observed spectrum (see Torres 2018; Torres et al. 2019, for details). The model follows particle trajectories in a generic region of a pulsar magnetosphere threaded by an accelerating parallel electric field, ($E_‖$). The region is located around the light cylinder, and particles are assumed to enter it at $x_{in}$ with a (sizeable) pitch angle $\alpha$. The model parameterizes the magnetic field by a power law $B(x) = B_s(R_s/x)^b$ (see the discussion in Viganò et al. (2015a)), where $x$ is the distance along the field line, $b$ is referred to as the magnetic gradient, $B_s$ is the surface magnetic field, and $R_s$ is the pulsar radius. Given ($E_‖, b$) as free parameters, and the period and period derivative ($P, \dot{P}$), the model solves the equations of motion that balance acceleration and losses by SC radiation (see Cheng & Zhang (1996); Viganò et al. (2015b)), computing the total emission. The model assumes that the distribution of particles emitting towards us can be parameterized as $dN_e/dx \propto e^{-(x-x_{in})/x_0}$ where the inverse of $x_0/R_{lc}$ ($R_{lc}$ is the light cylinder radius) is referred to as the contrast.

Figure 9 shows the results of the model with best fit parameters, $log(E_‖/V m^{-1})=10.92$, $log(x_0/R_{lc})=-4.20$ and $b=3.70$. The agreement between the model description and the broad-band data is acceptable (the fractional residual errors are of the order ~10%), despite the significant increase in the precision of each spectral measurement and the number of data points.

## 5. DISCUSSION AND CONCLUSIONS

The detection by ground-based Cherenkov telescopes of pulsed emission from the Crab pulsar (Aliu et al. 2008; VERITAS Collaboration et al. 2011), most recently detected up to TeV energies (Ansoldi et al. 2016), followed by the detection of pulsations from Vela (H. E. S. S. Collaboration et al. 2018), also up to TeV energies (Harding et al. 2018b), has led to a flurry of activity by pulsar experts to develop self-consistent models able to explain the detected emission over such a broad range of energies.

In addition, the development of the Sum-Trigger-II system in MAGIC has significantly improved the sensitivity of the telescopes below 100 GeV, something which has enabled the detection of Geminga between 15 GeV and 75 GeV (MAGIC
<figure><img src="image_8.png" /></figure>

**Figure 8.** Model predictions for the spectrum of phase-averaged emission from accelerated particles and pairs in PSR J0218+4232, for an assumed magnetic inclination angle $\alpha = 45^\circ$ and viewing angle $\zeta = 65^\circ$. The solid orange line represents the predicted ICS component due to accelerated SC-emitting primaries scattering the pair SR component (blue dotted line). The thick blue line identifies the overall emission model. Data points show the soft (XMM) and hard (NuSTAR) X-ray emission (from Gotthelf & Bogdanov 2017), as well as the LAT spectral points and MAGIC upper limits obtained in this work. Note that the LAT and MAGIC spectral points and upper limits represent the total (pulsed plus unpulsed) emission, however, given the broad peak and large pulsed fraction (see Figure 5), the differences between this and the *pulsed* spectrum would be marginal, and would not significantly affect the model fit parameters. Note that the $\gamma$-ray data are identical to those in Figure 4.

Collaboration et al. 2020), making this the third $\gamma$-ray pulsar (and first *middle-aged* one) detected with ground-based telescopes $^8$.

J0218 is part of a small but diverse population of MSPs with a well-characterized broad-band non-thermal energy distribution. Several of these MSPs have an ‘inverted’ spectrum in X-rays (where $E^2dF/dE$ decreases with energy), quite different from that of J0218 as we have reported here, see Figure 3 of Coti Zelati et al. (2020). We also note that the fitted magnetic gradient $b$ for J0218 and other MSPs within the synchro-curvature model (Torres 2018) is larger than for normal pulsars. This is perhaps the result of the larger $B_{lc}$ of MSPs compared to typical pulsars, due to the smaller size of $R_{lc}$. This needs to be taken into account when making predictions for their observability at lower energies based only on the $\gamma$-ray data. Fits to the $\gamma$-ray data alone are mostly insensitive to the value of the magnetic gradient, and assuming a lower $b$ could lead to incorrect predictions that an MSP is undetectable in the X-ray band. In accordance with previous studies, we also find here that the relevant scales for the production of the pulsar’s spectrum (given by $x_0$) is small in comparison with the light cylinder radius.

This is true in general for MSPs, for which the light cylinder is already orders of magnitude smaller than in normal pulsars, i.e., the $x_0/R_{lc}$-values imply a relevant region of emission $\ll 1$ km.

<footnote>$^8$ The H.E.S.S. Collaboration has also reported at the 36th International Cosmic Ray Conference (ICRC2019) the detection of $\gamma$-ray pulsations up to $\sim 70$ GeV from PSR B1706-44 (Spir-Jacob et al. 2019), which, if confirmed, would bring the total number of ground-based detected $\gamma$-ray pulsars to four.</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 9.** Broad-band spectrum of PSR J0218+4232, from the X-ray (XMM and NuSTAR) to the γ-ray (Fermi-LAT and MAGIC) range, along with the best fit to the synchro-curvature model (Torres 2018). The model is described by these parameters: $log(E_{\parallel}/V m^{-1})=10.92$, $log(x_0/R_{lc})=-4.20$, and $b=3.70$. Note that the X-ray and γ-ray data are identical to those in Figure 8.</figcaption></figure>

Instead of subtracting the background events from the off-pulse region, we applied reflected-region background subtraction approach for the MAGIC analysis due to the large on-pulse interval of the LAT phaseogram. No evidence of emission (either pulsed or unpulsed) is apparent in the MAGIC data, and the measured MAGIC upper limits are well above our two theoretical model predictions for VHE emission. The curvature radiation component from particles accelerated mostly in the current sheet is expected to fall to flux levels too low at VHE energies for detection by MAGIC and the ICS components from both pairs (mostly SSC) and accelerated primaries are predicted to be at even lower flux levels.

Most models for γ-ray emission from pulsars do not predict high levels of ICS and SSC emission for MSPs. In the model we used here, described by Harding & Kalapotharakos (2015) and Harding et al. (2018a), for example, the pairs that come from the PC cascade and MSP surface magnetic fields are so low that the photons need to have much higher energies to produce pairs by one-photon magnetic pair production than do normal pulsars. The MSP pair spectra are thus shifted to much higher energies (typically $\gamma \sim 10^4 - 10^7$) (Harding & Muslimov 2011). This will produce higher energy SR near the light cylinder. Since VHE emission is most likely ICS or SSC, and both particles and photons have higher energies, the VHE emission will be Klein-Nishina limited and therefore suppressed. This is also a problem for outer gap (OG) models since the latest models have pairs also produced near the PC since otherwise, MSPs cannot sustain OGs (Harding 2021). Observationally, we see that the SR spectra seem to extend to higher energy in MSPs (at least the energetic ones that have non-thermal emission). So the SR photons and the particles that produce them must be at higher energy. The Cherenkov Telescope Array (CTA) is expected to have significantly better sensitivity than MAGIC in the 10–100 GeV range, and this and other pulsars will thus be prime targets for observation (Burtovoi et al. 2017). On the other hand, pulsars like J0218 are also good sources for MeV telescopes, such as AMEGO (McEnery et al. 2019), that can detect the predicted SR peaks around 1 - 10 MeV.
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), High Energy Accelerator Research Organization (KEK) and 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. This research was partially carried out using the HKU Information Technology Services research computing facilities that are supported in part by the Hong Kong UGC Special Equipment Grant (SEG HKU09). P.S.P. and C.Y.N. are supported at HKU by a grant from the Big Data Project Fund (BDPF) and a GRF grant (Project 17307618) from the Hong Kong Government. AKH acknowledges resources supporting the theoretical modeling provided by the NASA High-End Computing Program through the NASA Center for Climate Simulation. DFT acknowledges support from the Spanish grants PGC2018-095512-B-I00, SGR2017-1383; as well as from the Chinese Academy of Sciences Presidential Fellowship Initiative 2021VMA0001. The Nançay Radio Observatory is operated by the Paris Observatory, associated with the French Centre National de la Recherche Scientifique (CNRS). We acknowledge financial support from the "Programme National Hautes Energies" (PNHE) of CNRS/INSU, France.

We would like to thank the Instituto de Astrofísica de Canarias for the excellent working conditions at the Observatorio del Roque de los Muchachos in La Palma. The financial support of the German BMBF, MPG and HGF; the Italian INFN and INAF; the Swiss National Fund SNF; the ERDF under the Spanish Ministerio de Ciencia e Innovación (MICINN) (FPA2017-87859-P, FPA2017-85668-P, FPA2017-82729-C6-5-R, FPA2017-90566-REDC, PID2019-104114RB-C31, PID2019-104114RB-C32, PID2019-105510GB-C31, PID2019-107847RB-C41, PID2019-107847RB-C42, PID2019-107988GB-C22); the Indian Department of Atomic Energy; the Japanese ICRR, the University of Tokyo, JSPS, and MEXT; the Bulgarian Ministry of Education and Science, National RI Roadmap Project DO1-268/16.12.2019 and the Academy of Finland grant nr. 320045 is gratefully acknowledged. This work was also supported by the Spanish Centro de Excelencia "Severo Ochoa" SEV-2016-0588 and CEX2019-000920-S, and "María de Maeztu" CEX2019-000918-M, the Unidad de Excelencia "María de Maeztu" MDM-2015-0509-18-2 and the "la Caixa" Foundation (fellowship LCF/BQ/PI18/11630012) and by the CERCA program of the Generalitat de Catalunya; by the Croatian Science Foundation (HrZZ) Project IP-2016-06-9782 and the University of Rijeka Project 13.12.1.3.02; by the DFG Collaborative Research Centers SFB823/C4 and SFB876/C3; the Polish National Research Centre grant UMO-2016/22/M/ST9/00382; and by the Brazilian MCTIC, CNPq and FAPERJ. We thank Matthew Kerr and Philippe Bruel for their careful reading of the manuscript and useful comments and suggestions, which greatly improved the paper. We thank Slavko Bogdanov for providing the X-ray data points used in this work. Finally, we sincerely thank the journal referee for taking the time and effort to provide a detailed review of our paper, providing some very helpful feedback, which helped us improve the final version of the article.

*Software:* Fermipy (Wood et al. 2017), Minuit (James 1994), MARS (Zanin et al. 2013), R (R Core Team 2013), SOPIE (Schutte & Swanepoel 2016), Tempo2 (Hobbs et al. 2006)

*Facilities:* Fermi-LAT, MAGIC, ADS

## REFERENCES

Abdo, A. A., Ackermann, M., Ajello, M., et al. 2009, Science, 325, 848, doi: 10.1126/science.1176113

-. 2010a, ApJS, 188, 405, doi: 10.1088/0067-0049/188/2/405

-. 2010b, ApJS, 187, 460, doi: 10.1088/0067-0049/187/2/460

Abdollahi, S., Acero, F., Ackermann, M., et al. 2020, ApJS, 247, 33, doi: 10.3847/1538-4365/ab6bcb

Acciari, V. A., Aliu, E., Arlen, T., et al. 2009, ApJL, 693, L104, doi: 10.1088/0004-637X/693/2/L104

Ackermann, M., Ajello, M., Allafort, A., et al. 2013, ApJS, 209, 34, doi: 10.1088/0067-0049/209/2/34
Ajello, M., Atwood, W. B., Baldini, L., et al. 2017a, ApJS, 232, 18, doi: 10.3847/1538-4365/aa8221

-. 2017b, ApJS, 232, 18, doi: 10.3847/1538-4365/aa8221

Aleksić, J., Ansoldi, S., Antonelli, L. A., et al. 2016a, Astroparticle Physics, 72, 61, doi: 10.1016/j.astropartphys.2015.04.004

-. 2016b, Astroparticle Physics, 72, 76, doi: 10.1016/j.astropartphys.2015.02.005

Aliu, E., Anderhub, H., Antonelli, L. A., et al. 2008, Science, 322, 1221, doi: 10.1126/science.1164718

-. 2009, ApJL, 692, L29, doi: 10.1088/0004-637X/692/1/L29

Anderhub, H., Antonelli, L. A., Antoranz, P., et al. 2010, ApJ, 710, 828, doi: 10.1088/0004-637X/710/1/828

Ansoldi, S., Antonelli, L. A., Antoranz, P., et al. 2016, A&A, 585, A133, doi: 10.1051/0004-6361/201526853

Atwood, W., Albert, A., Baldini, L., et al. 2013, in Proceedings of the 4th Fermi Symposium. https://arxiv.org/abs/1303.3514

Atwood, W. B., Abdo, A. A., Ackermann, M., et al. 2009, ApJ, 697, 1071, doi: 10.1088/0004-637X/697/2/1071

Bassa, C. G., van Kerkwijk, M. H., & Kulkarni, S. R. 2003, A&A, 403, 1067, doi: 10.1051/0004-6361:20030384

Berge, D., Funk, S., & Hinton, J. 2007, A&A, 466, 1219, doi: 10.1051/0004-6361:20066674

Bruel, P., Burnett, T. H., Digel, S. W., et al. 2018, in Proceedings of the 8th Fermi Symposium. https://arxiv.org/abs/1810.11394

Burtovoi, A., Saito, T. Y., Zampieri, L., & Hassan, T. 2017, MNRAS, 471, 431, doi: 10.1093/mnras/stx1582

Cheng, K. S., & Zhang, J. L. 1996, ApJ, 463, 271, doi: 10.1086/177239

Coti Zelati, F., Torres, D. F., Li, J., & Viganò, D. 2020, MNRAS, 492, 1025, doi: 10.1093/mnras/stz3485

Dazzi, F., Schweizer, T., Ceribella, G., et al. 2021, IEEE Transactions on Nuclear Science, 1, doi: 10.1109/TNS.2021.3079262

de Jager, O. C., Raubenheimer, B. C., & Swanepoel, J. W. H. 1989, A&A, 221, 180

Deneva, J. S., Ray, P. S., Lommen, A., et al. 2019, ApJ, 874, 160, doi: 10.3847/1538-4357/ab0966

Du, Y., Yang, J., Campbell, R. M., et al. 2014, ApJL, 782, L38, doi: 10.1088/2041-8205/782/2/L38

Fomin, V. P., Stepanian, A. A., Lamb, R. C., et al. 1994, Astroparticle Physics, 2, 137, doi: 10.1016/0927-6505(94)90036-1

Fruck, C., Gaug, M., Zanin, R., et al. 2014, 33rd International Cosmic Ray Conference (ICRC), Rio de janerio 2013, arXiv:1403.3591. https://arxiv.org/abs/1403.3591

Gotthelf, E. V., & Bogdanov, S. 2017, ApJ, 845, 159, doi: 10.3847/1538-4357/aa813c

Guillemot, L., Lonjou, V., Dumora, D., et al. 2007, in American Institute of Physics Conference Series, Vol. 921, The First GLAST Symposium, ed. S. Ritz, P. Michelson, & C. A. Meegan, 395-396, doi: 10.1063/1.2757369

H. E. S. S. Collaboration, Abdalla, H., Aharonian, F., et al. 2018, A&A, 620, A66, doi: 10.1051/0004-6361/201732153

Harding, A. K. 2021, arXiv e-prints, arXiv:2101.05751. https://arxiv.org/abs/2101.05751

Harding, A. K., & Kalapotharakos, C. 2015, ApJ, 811, 63, doi: 10.1088/0004-637X/811/1/63

Harding, A. K., Kalapotharakos, C., Barnard, M., & Venter, C. 2018a, ApJL, 869, L18, doi: 10.3847/2041-8213/aaf3b2

-. 2018b, ApJL, 869, L18, doi: 10.3847/2041-8213/aaf3b2

Harding, A. K., & Muslimov, A. G. 2011, ApJ, 743, 181, doi: 10.1088/0004-637X/743/2/181

Harding, A. K., Stern, J. V., Dyks, J., & Frackowiak, M. 2008, ApJ, 680, 1378, doi: 10.1086/588037

Hartman, R. C., Bertsch, D. L., Bloom, S. D., et al. 1999, ApJS, 123, 79, doi: 10.1086/313231

Hobbs, G. B., Edwards, R. T., & Manchester, R. N. 2006, MNRAS, 369, 655, doi: 10.1111/j.1365-2966.2006.10302.x

James, F. 1994, CERN Report, CERN-D-506

Joshi, B. C., Kramer, M., Lyne, A. G., McLaughlin, M. A., & Stairs, I. H. 2004, in IAU Symposium, Vol. 218, Young Neutron Stars and Their Environments, ed. F. Camilo & B. M. Gaensler, 319. https://arxiv.org/abs/astro-ph/0310285

Knight, H. S., Bailes, M., Manchester, R. N., Ord, S. M., & Jacoby, B. A. 2006, ApJ, 640, 941, doi: 10.1086/500292

Kuiper, L., Hermsen, W., & Stappers, B. 2004, Advances in Space Research, 33, 507, doi: 10.1016/j.asr.2003.08.019

Kuiper, L., Hermsen, W., Verbunt, F., & Belloni, T. 1998, A&A, 336, 545

Kuiper, L., Hermsen, W., Verbunt, F., et al. 2002, ApJ, 577, 917, doi: 10.1086/342220

-. 2000, A&A, 359, 615. https://arxiv.org/abs/astro-ph/0005338

Lange, C., Camilo, F., Wex, N., et al. 2001, MNRAS, 326, 274, doi: 10.1046/j.1365-8711.2001.04606.x

Li, T. P., & Ma, Y. Q. 1983, ApJ, 272, 317, doi: 10.1086/161295

MAGIC Collaboration, Acciari, V. A., Ansoldi, S., et al. 2020, A&A, 643, L14, doi: 10.1051/0004-6361/202039131

McEnergy, J., van der Horst, A., Dominguez, A., et al. 2019, in Bulletin of the American Astronomical Society, Vol. 51, 245. https://arxiv.org/abs/1907.07558
Mineo, T., Cusumano, G., Kuiper, L., et al. 2000, A&A,
355, 1053

Navarro, J., de Bruyn, A. G., Frail, D. A., Kulkarni, S. R.,
& Lyne, A. G. 1995, ApJL, 455, L55, doi: 10.1086/309816

Oña-Wilhemi, E., de los Reyes, R., Contreras, J. L., et al. 2005, in International Cosmic Ray Conference, Vol. 4, 29th International Cosmic Ray Conference (ICRC29), Volume 4, 247

Otero-Santos, J., Acosta-Pulido, J. A., Becerra González, J., et al. 2020, MNRAS, 492, 5524, doi: 10.1093/mnras/staa134

R Core Team. 2013, R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing, Vienna, Austria.
http://www.R-project.org/

Ray, P. S., Kerr, M., Parent, D., et al. 2011, ApJS, 194, 17, doi: 10.1088/0067-0049/194/2/17

Rolke, W. A., & López, A. M. 2001, Nuclear Instruments and Methods in Physics Research A, 458, 745, doi: 10.1016/S0168-9002(00)00935-9

Saito, Y., Kawai, N., Kamae, T., et al. 1997, ApJL, 477, L37, doi: 10.1086/310512

Saz Parkinson, P., Belfiore, A., Fidalgo, D., et al. 2017, in Proceedings of the 7th International Fermi Symposium, 8. https://arxiv.org/abs/1712.06808

Schutte, W. D., & Swanepoel, J. W. H. 2016, MNRAS, 461, 627, doi: 10.1093/mnras/stw1335

Spir-Jacob, M., Djannati-Ataï, A., Mohrmann, L., et al. 2019, arXiv e-prints, arXiv:1908.06464. https://arxiv.org/abs/1908.06464

Stairs, I. H., Thorsett, S. E., & Camilo, F. 1999, ApJS, 123, 627, doi: 10.1086/313245

Takahashi, M. 2019, in International Cosmic Ray Conference, Vol. 36, 36th International Cosmic Ray Conference (ICRC2019), 606

Takahashi, M., Caputo, R., Paneque, D., & Sgrò, C. 2015, 5th Fermi Symposium, Nagoya, Japan, 20-24 Oct 2014, arXiv:1503.01364. https://arxiv.org/abs/1503.01364

Thompson, D. J., Bertsch, D. L., Dingus, B. L., et al. 1995, ApJS, 101, 259, doi: 10.1086/192240

Torres, D. F. 2018, Nature Astronomy, 2, 247, doi: 10.1038/s41550-018-0384-5

Torres, D. F., Viganò, D., Coti Zelati, F., & Li, J. 2019, MNRAS, 489, 5494, doi: 10.1093/mnras/stz2403

Verbiest, J. P. W., & Lorimer, D. R. 2014, MNRAS, 444, 1859, doi: 10.1093/mnras/stu1560

Verbunt, F., Kuiper, L., Belloni, T., et al. 1996, A&A, 311, L9

VERITAS Collaboration, Aliu, E., Arlen, T., et al. 2011, Science, 334, 69, doi: 10.1126/science.1208192

Viganò, D., Torres, D. F., Hirotani, K., & Pessah, M. E. 2015a, MNRAS, 447, 2649, doi: 10.1093/mnras/stu2565 —. 2015b, MNRAS, 447, 1164, doi: 10.1093/mnras/stu2456

Webb, N. A., Olive, J. F., & Barret, D. 2004, A&A, 417, 181, doi: 10.1051/0004-6361:20040032

Wilks, S. S. 1938, Ann. Math. Statist., 9, 60, doi: 10.1214/aoms/1177732360

Wood, M., Caputo, R., Charles, E., et al. 2017, in International Cosmic Ray Conference, Vol. 301, 35th International Cosmic Ray Conference (ICRC2017), 824. https://arxiv.org/abs/1707.09551

Zanin, R., Carmona, E., Sitarek, J., et al. 2013, in International Cosmic Ray Conference, Vol. 33, International Cosmic Ray Conference, 2937