# Multiwavelength Observations of the Blazar PKS 0735+178 in Spatial and Temporal Coincidence with an Astrophysical Neutrino Candidate IceCube-211208A

A Acharyya, C.B Adams, A Archer, P Bangale, J.T Bartkoske, P Batista, W

Benbow, A Brill, J.H Buckley, J.L Christiansen, et al.

► To cite this version:

A Acharyya, C.B Adams, A Archer, P Bangale, J.T Bartkoske, et al.. Multiwavelength Observations of the Blazar PKS 0735+178 in Spatial and Temporal Coincidence with an Astrophysical Neutrino Candidate IceCube-211208A. Astrophys.J., 2023, 954 (1), pp.70. 10.3847/1538-4357/ace327. hal-04160420

HAL Id: hal-04160420

https://hal.science/hal-04160420

Submitted on 9 Nov 2023

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
Multiwavelength Observations of the Blazar PKS 0735+178 in Spatial and Temporal Coincidence with an Astrophysical Neutrino Candidate IceCube-211208A

A. Acharyya¹, C. B. Adams², A. Archer³, P. Bangale⁴, J. T. Bartkoske⁵, P. Batista⁶, W. Benbow⁷, A. Brill⁸, J. H. Buckley⁹, J. L. Christiansen¹⁰, A. J. Chromey⁷, M. Errando⁹, A. Falcone¹¹, Q. Feng⁷, G. M. Foote⁴, L. Fortson¹², A. Furniss¹³, G. Gallagher¹⁴, W. Hanlon⁷, D. Hanna¹⁵, O. Hervet¹⁶, C. E. Hinrichs¹⁷, J. Hoang¹⁶, J. Holder⁴, T. B. Humensky¹⁸, W. Jin¹, P. Kaaret¹⁹, M. Kertzman³, M. Kherlakian⁶, D. Kieda⁵, T. K. Kleiner⁶, N. Korzoun⁴, S. Kumar²⁰, M. J. Lang²¹, M. Lundy¹⁵, G. Maier⁶, C. E McGrath²², M. J. Millard¹⁹, J. Millis¹⁴, C. L. Mooney⁴, P. Moriarty²¹, R. Mukherjee²³, S. O'Brien¹⁵, R. A. Ong²⁴, M. Pohl²⁵, E. Pueschel⁶, J. Quinn²², K. Ragan¹⁵, P. T. Reynolds²⁶, D. Ribeiro¹², E. Roache⁷, I. Sadeh⁶, A. C. Sadun²⁷, L. Saha⁷, M. Santander¹, G. H. Sembroski²⁸, R. Shang²³, M. Splettstoesser¹⁶, A. Kaushik Talluri¹², J. V. Tucci²⁹, V. V. Vassiliev²⁴, A. Weinstein³⁰, D. A. Williams¹⁶, S. L. Wong¹⁵, J. Woo³¹

The VERITAS Collaboration,

and

F. Aharonian³²,³³, J. Aschersleben³⁴, M. Backes³⁵,³⁶, V. Barbosa Martins³⁷, R. Batzofin³⁸, Y. Becherini³⁹,⁴⁰, D. Berge³⁷,⁴¹, K. Bernlöhr³³, B. Bi⁴², M. Böttcher³⁶, C. Boisson⁴³, J. Bolmont⁴⁴, M. de Bony de Lavergne⁴⁵, J. Borowska⁴¹, M. Bouyahiaoui³³, F. Bradascio⁴⁶, M. Breuhaus³³, R. Brose³², F. Brun⁴⁶, B. Bruno⁴⁷, T. Bulik⁴⁸, C. Burger-Scheidlin³², S. Caroff⁴⁵, S. Casanova⁴⁹, R. Cecil⁵⁰, J. Celic⁴⁷, M. Cerruti³⁹, T. Chand³⁶, S. Chandra³⁶, A. Chen⁵¹, J. Chibueze³⁶, O. Chibueze³⁶, G. Cotter⁵², S. Dai⁵³, J. Damascene Mbarubucyeye³⁷, A. Djannati-Atai³⁹, A. Dmytriiev³⁶, V. Doroshenko⁴², S. Einecke⁵⁴, J.-P. Ernenwein⁵⁵, G. Fichet de Clairfontaine⁴³, M. Filipovic³³, G. Fontaine⁵⁶, M. Füßling³⁷, S. Funk⁴⁷, S. Gabici³⁹, S. Ghafourizadeh⁵⁷, G. Giavitto³⁷, D. Glawion⁴⁷, J. F. Glicenstein⁴⁶, P. Goswami³⁶, G. Grolleron⁴⁴, L. Haerer³³, J. A. Hinton³³, T. L. Holch³⁷, M. Holler⁵⁸, D. Horns⁵⁰, M. Jamrozy⁵⁹, F. Jankowsky⁵⁷, V. Joshi⁴⁷, I. Jung-Richardt⁴⁷, E. Kasai³⁵, K. Katarzyński⁶⁰, R. Khatoon³⁶, B. Khélifi³⁹, S. Klepser³⁷, W. Kluźniak⁶¹, K. Kosack⁴⁶, D. Kostunin³⁷, R. G. Lang⁴⁷, S. Le Stum⁵⁵, A. Lemière³⁹, J.-P. Lenain⁴⁴, F. Leuschner⁴², T. Lohse⁴¹, A. Luashvili⁴³, I. Lypova⁵⁷, J. Mackey³², D. Malyshev⁴², V. Marandon³³, P. Marchegiani⁵¹, A. Marcowith⁶², G. Martí-Devesa⁵⁸, R. Marx⁵⁷, A. Mitchell⁴⁷, R. Moderski⁶¹, L. Mohrmann³³, A. Montanari⁵⁷, E. Moulin⁴⁶, T. Murach³⁷, K. Nakashima⁴⁷, J. Niemiec⁴⁹, A. Priyana Noel⁵⁹, P. O'Brien⁶³, L. Olivera-Nieto³³, E. de Ona Wilhelmi³⁷, M. Ostrowski⁵⁹, S. Panny⁵⁸, M. Panter³³, G. Peron³⁹, D. A. Prokhorov⁶⁴, G. Pühlhofer⁴², M. Punch³⁹, A. Quirrenbach⁵⁷, P. Reichherzer⁴⁶, A. Reimer⁵⁸, O. Reimer⁵⁸, H. Ren³³, M. Renaud⁶², F. Rieger³³, B. Rudak⁶¹, E. Ruiz-Velasco³³, V. Sahakian⁶⁵, A. Santangelo⁴², M. Sasaki⁴⁷, J. Schäfer⁴⁷, F. Schüssler⁴⁶, H. M. Schutte³⁶, U. Schwanke⁴¹, J. N. S. Shapopi³⁵, A. Specovius⁴⁷, S. Spencer⁴⁷, Ł. Stawarz⁵⁹, R. Steenkamp³⁵, S. Steinmassl³³, I. Sushch³⁶, H. Suzuki⁶⁶, T. Takahashi⁶⁷, T. Tanaka⁶⁶, R. Terrier³⁹, C. van Eldik⁴⁷, M. Vecchi³⁴, J. Veh⁴⁷, C. Venter³⁶, J. Vink⁶⁴, R. White³³, A. Wierzcholska⁴⁹, Yu Wun Wong⁴⁷, M. Zacharias³⁶,⁵⁷, D. Zargaryan³², A. A. Zdziarski⁶¹, A. Zech⁴³, S. Zouari³⁹, N. Żywucka³⁶

The H.E.S.S. Collaboration, and

K. Mori²

¹Department of Physics and Astronomy, University of Alabama, Tuscaloosa, AL 35487, USA; aacharyya1@ua.edu, wjin4@crimson.ua.edu

²Physics Department, Columbia University, New York, NY 10027, USA

³Department of Physics and Astronomy, DePauw University, Greencastle, IN 46135-0037, USA

⁴Department of Physics and Astronomy and the Bartol Research Institute, University of Delaware, Newark, DE 19716, USA

⁵Department of Physics and Astronomy, University of Utah, Salt Lake City, UT 84112, USA

⁶DESY, Platanenallee 6, D-15738 Zeuthen, Germany

⁷Center for Astrophysics | Harvard & Smithsonian, Cambridge, MA 02138, USA; qi.feng@cfa.harvard.edu

⁸N.A.S.A./Goddard Space-Flight Center, Code 661, Greenbelt, MD 20771, USA

⁹Department of Physics, Washington University, St. Louis, MO 63130, USA

¹⁰Physics Department, California Polytechnic State University, San Luis Obispo, CA 94307, USA

¹¹Department of Astronomy and Astrophysics, 525 Davey Lab, Pennsylvania State University, University Park, PA 16802, USA

¹²School of Physics and Astronomy, University of Minnesota, Minneapolis, MN 55455, USA

¹³Department of Physics, California State University - East Bay, Hayward, CA 94542, USA

¹⁴Department of Physics and Astronomy, Ball State University, Muncie, IN 47306, USA

¹⁵Physics Department, McGill University, Montreal, QC H3A 2T8, Canada

¹⁶Santa Cruz Institute for Particle Physics and Department of Physics, University of California, Santa Cruz, CA 95064, USA

¹⁷Center for Astrophysics | Harvard & Smithsonian, Cambridge, MA 02138, USA and Department of Physics and Astronomy, Dartmouth College, 6127 Wilder Laboratory, Hanover, NH 03755 USA

¹⁸Department of Physics, University of Maryland, College Park, MD, USA and NASA GSFC, Greenbelt, MD 20771, USA

¹⁹Department of Physics and Astronomy, University of Iowa, Van Allen Hall, Iowa City, IA 52242, USA

²⁰Department of Physics, University of Maryland, College Park, MD, USA

²¹School of Natural Sciences, University of Galway, University Road, Galway H91 TK33, Ireland

²²School of Physics, University College Dublin, Belfield, Dublin 4, Ireland

²³Department of Physics and Astronomy, Barnard College, Columbia University, NY 10027, USA

²⁴Department of Physics and Astronomy, University of California, Los Angeles, CA 90095, USA
²⁵ Institute of Physics and Astronomy, University of Potsdam, 14476 Potsdam-Golm, Germany and DESY, Platanenallee 6, D-15738 Zeuthen, Germany;
martin.pohl@desy.de

²⁶ Department of Physical Sciences, Munster Technological University, Bishopstown, Cork T12 P928, Ireland
²⁷ Campus Box 157, P.O.Box 173364, Denver, CO 80217, USA

²⁸ Department of Physics and Astronomy, Purdue University, West Lafayette, IN 47907, USA

²⁹ Department of Physics, Indiana University-Purdue University Indianapolis, Indianapolis, IN 46202, USA
³⁰ Department of Physics and Astronomy, Iowa State University, Ames, IA 50011, USA

³¹ Columbia Astrophysics Laboratory, Columbia University, New York, NY 10027, USA; jw3855@columbia.edu
³² Dublin Institute for Advanced Studies, 31 Fitzwilliam Place, Dublin 2, Ireland

³³ Max-Planck-Institut für Kernphysik, P.O. Box 103980, D-69029 Heidelberg, Germany

³⁴ Kapteyn Astronomical Institute, University of Groningen, Landleven 12, 9747 AD Groningen, The Netherlands
³⁵ University of Namibia, Department of Physics, Private Bag 13301, Windhoek 10005, Namibia

³⁶ Centre for Space Research, North-West University, Potchefstroom 2520, South Africa

³⁷ DESY, D-15738 Zeuthen, Germany; contact.hess@hess-experiment.eu

³⁸ Institut für Physik und Astronomie, Universität Potsdam, Karl-Liebknecht-Strasse 24/25, D-14476 Potsdam, Germany

³⁹ Université de Paris, CNRS, Astroparticule et Cosmologie, F-75013 Paris, France; contact.hess@hess-experiment.eu

⁴⁰ Department of Physics and Electrical Engineering, Linnaeus University, SE-351 95 Växjö, Sweden

⁴¹ Institut für Physik, Humboldt-Universität zu Berlin, Newtonstr. 15, D-12489 Berlin, Germany

⁴² Institut für Astronomie und Astrophysik, Universität Tübingen, Sand 1, D-72076 Tübingen, Germany

⁴³ Laboratoire Univers et Théories, Observatoire de Paris, Université PSL, CNRS, Université de Paris, F-92190 Meudon, France

⁴⁴ Sorbonne Université, Université Paris Diderot, Sorbonne Paris Cité, CNRS/IN2P3, Laboratoire de Physique Nucléaire et de Hautes Energies, LPNHE, 4 Place Jussieu, F-75252 Paris, France

⁴⁵ Université Savoie Mont Blanc, CNRS, Laboratoire d'Annecy de Physique des Particules - IN2P3, F-74000 Annecy, France

⁴⁶ IRFU, CEA, Université Paris-Saclay, F-91191 Gif-sur-Yvette, France; contact.hess@hess-experiment.eu

⁴⁷ Friedrich-Alexander-Universität Erlangen-Nürnberg, Erlangen Centre for Astroparticle Physics, Erwin-Rommel-Str. 1, D-91058 Erlangen, Germany
⁴⁸ Astronomical Observatory, The University of Warsaw, Al. Ujazdowskie 4, 00-478 Warsaw, Poland

⁴⁹ Institut Fizyki Jądrowej PAN, ul. Radzikowskiego 152, 31-342 Kraków, Poland

⁵⁰ Universität Hamburg, Institut für Experimentalphysik, Luruper Chaussee 149, D-22761 Hamburg, Germany

⁵¹ School of Physics, University of the Witwatersrand, 1 Jan Smuts Avenue, Braamfontein, Johannesburg, 2050 South Africa

⁵² University of Oxford, Department of Physics, Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, UK

⁵³ School of Science, Western Sydney University, Locked Bag 1797, Penrith South DC, NSW 2751, Australia

⁵⁴ School of Physical Sciences, University of Adelaide, Adelaide 5005, Australia

⁵⁵ Aix Marseille Université, CNRS/IN2P3, CPPM, Marseille, France

⁵⁶ Laboratoire Leprince-Ringuet, École Polytechnique, CNRS, Institut Polytechnique de Paris, F-91128 Palaiseau, France
⁵⁷ Landes sternwarte, Universität Heidelberg, Königstuhl, D-69117 Heidelberg, Germany

⁵⁸ Leopold-Franzens-Universität Innsbruck, Institut für Astro- und Teilchenphysik, A-6020 Innsbruck, Austria

⁵⁹ Obserwatorium Astronomiczne, Uniwersytet Jagielloński, ul. Orla 171, 30-244 Kraków, Poland

⁶⁰ Institute of Astronomy, Faculty of Physics, Astronomy and Informatics, Nicolaus Copernicus University, Grudziadzka 5, 87-100 Torun, Poland

⁶¹ Nicolaus Copernicus Astronomical Center, Polish Academy of Sciences, ul. Bartycka 18, 00-716 Warsaw, Poland

⁶² Laboratoire Univers et Particules de Montpellier, Université Montpellier, CNRS/IN2P3, CC 72, Place Eugène Bataillon, F-34095 Montpellier Cedex 5, France
⁶³ Department of Physics and Astronomy, The University of Leicester, University Road, Leicester LE1 7RH, UK

⁶⁴ GRAPPA, Anton Pannekoeck Institute for Astronomy, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands
⁶⁵ Yerevan Physics Institute, 2 Alikhanian Brothers St., 375036 Yerevan, Armenia

⁶⁶ Department of Physics, Konan University, 8-9-1 Okamoto, Higashinada, Kobe, Hyogo 658-8501, Japan

⁶⁷ Kavli Institute for the Physics and Mathematics of the Universe (WPI), The University of Tokyo Institutes for Advanced Study (UTIAS), The University of Tokyo, 5-1-5 Kashiwa-no-Ha, Kashiwa, Chiba 277-8583, Japan

Received 2023 April 6; revised 2023 June 5; accepted 2023 June 28; published 2023 August 23

## Abstract

We report on multiwavelength target-of-opportunity observations of the blazar PKS 0735+178, located 2.2 away from the best-fit position of the IceCube neutrino event IceCube-211208A detected on 2021 December 8. The source was in a high-flux state in the optical, ultraviolet, X-ray, and GeV γ-ray bands around the time of the neutrino event, exhibiting daily variability in the soft X-ray flux. The X-ray data from Swift-XRT and NuSTAR characterize the transition between the low-energy and high-energy components of the broadband spectral energy distribution (SED), and the γ-ray data from Fermi-LAT, VERITAS, and H.E.S.S. require a spectral cutoff near 100 GeV. Both the X-ray and γ-ray measurements provide strong constraints on the leptonic and hadronic models. We analytically explore a synchrotron self-Compton model, an external Compton model, and a lepto-hadronic model. Models that are entirely based on internal photon fields face serious difficulties in matching the observed SED. The existence of an external photon field in the source would instead explain the observed γ-ray spectral cutoff in both the leptonic and lepto-hadronic models and allow a proton jet power that marginally agrees with the Eddington limit in the lepto-hadronic model. We show a numerical lepto-hadronic model with external target photons that reproduces the observed SED and is reasonably consistent with the neutrino event despite requiring a high jet power.

*Unified Astronomy Thesaurus concepts: Active galactic nuclei (16); Blazars (164); Neutrino astronomy (1100); Gamma-ray astronomy (628)*

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

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

## 1. Introduction

The IceCube Neutrino Observatory has detected a diffuse flux of astrophysical neutrinos (Aartsen et al. 2013), whose
isotropic distribution of arrival directions suggests an extra-
galactic origin, but has not firmly identified any neutrino point
sources to date despite strong evidence of TeV neutrino
emission from the nearby active galaxy NGC 1068 (Aartsen
et al. 2020; Abbasi et al. 2022). A number of extragalactic
sources are proposed as candidates for high-energy neutrino
emitters, including clusters of galaxies, active galactic nuclei
(including blazars), starburst galaxies, γ-ray bursts, super-
novae, and tidal disruption events (e.g., Kurahashi et al. 2022).
Any detection of a TeV–PeV neutrino-emitting source will
directly constrain the century-old puzzle of the origin of cosmic
rays (e.g., Mészáros 2017), as these neutrinos must be
produced by hadronic cosmic-ray interactions.

Hadronic cosmic-ray interactions produce not only neutrinos but also γ-rays (from neutral pion decays) and X-rays (from synchrotron radiation and cascading of the secondary electrons/positrons). Observations at γ-ray and X-ray bands are therefore critical for studies of neutrinos and cosmic rays. An important type of source that exhibits strong and highly variable γ-ray and X-ray emission is blazars, a subclass of active galactic nuclei with relativistic jets pointing toward Earth. High-energy emission from blazars can be produced by either leptonic (via inverse-Compton (IC) scattering of relativistic electrons) or hadronic interactions (Böttcher et al. 2013), rendering electromagnetic observations, which, by themselves, are generally insufficient to probe the origin of cosmic rays.

A more efficient approach to identify astrophysical neutrino and ultra-high-energy cosmic-ray sources is through multi-messenger observations of correlated neutrino and γ-ray events. High-power γ-ray blazars are disfavored as the dominant origin of the observed IceCube neutrinos (e.g., Aartsen et al. 2017), but individual flaring blazars still provide promising opportunities for the identification of neutrino emitters (e.g., Murase et al. 2018). Because of their highly variable nature, blazars can undergo strong flaring episodes during which an accompanying IceCube neutrino signal would stand out from the background. Besides, preselecting the time window of a blazar flare also mitigates the look-elsewhere effect associated with a blind point-source search in the entire neutrino data set (see, e.g., Abbasi et al. 2021), further improving the sensitivity of these searches.

The first evidence for a candidate TeV-PeV extragalactic neutrino source involved the γ-ray blazar TXS 0506+056. The coincident detection of the neutrino event, IceCube-170922A, and the temporally correlated γ-ray flaring activity from TXS 0506+056 in 2017 (Aartsen et al. 2018a), combined with an excess of 13 ± 5 muon-neutrino events observed between 2014 and 2015 by IceCube (Aartsen et al. 2018b), suggested a possible association of the neutrino emission with the blazar. The lack of γ-ray counterpart during the neutrino flare in 2014–2015 and the low γ-ray flux seen by Fermi-LAT from NGC 1068 compared to the neutrino flux (Abbasi et al. 2022), however, suggest that γ-rays may be absorbed in a dense radiation field that acts as an efficient target for neutrino production through photohadronic processes. Therefore, it is still unclear whether or not astrophysical neutrinos are associated with blazars and blazar flares. Further multi-messenger data of other candidate sources are important to understand the conditions and mechanisms for potential neutrino emission from blazars.

A recent opportunity to explore the connection between neutrinos and high-energy blazars came with the announcement of the neutrino event IceCube-211208A (IceCube Collaboration 2021) detected by IceCube as a track-like event with energy $E_ν \approx 171$ TeV and a 50.2% probability of having an astrophysical origin<fnref n="68" /> on 2021 December 8. The γ-ray blazar PKS 0735+178 (redshift $z=0.45$) is located immediately outside of the 90% error region (2°13 statistical error only) for the neutrino event, 2°2 away from the best-fit position. Additionally, the Baikal-GVD experiment detected a high-energy neutrino candidate event with energy $E_ν \approx 43$ TeV (Dzhilkibaev et al. 2021) approximately four hours after the IceCube event, about 4°7 (with an estimated 68% containment point-spread function of 8°1) from the position of PKS 0735+170. KM3NeT neutrino detectors found one up-going muon-neutrino candidate (~18 TeV) on December 15 in spatial coincidence with PKS 0735+178 with a *p*-value of 0.14 (Filippini et al. 2022). The Baksan Underground Scintillation Telescope reported the observation of a GeV neutrino candidate event four days before IceCube-211208A (Petkov et al. 2021). Multiwavelength observations of PKS 0735+178 revealed flaring states in the radio band (Kadler et al. 2021), optical band (Zhirkov et al. 2021), X-ray band (D'Ammando 2021; Feng et al. 2021; Santander & Buson 2021), and GeV γ-ray band (Garrappa et al. 2021), which are potentially associated with the neutrino event IceCube-211208A. To follow up, we triggered NuSTAR observations and provided precise measurements on 2021 December 11 and 13 of the X-ray spectrum that was found to be harder than that seen with Swift-XRT (Feng et al. 2021). In addition, the imaging atmospheric Cherenkov telescopes (IACTs) the Very Energetic Radiation Imaging Telescope Array System (VERITAS) and the High Energy Stereoscopic System (H.E.S.S.) performed target-of-opportunity (ToO) observations that yielded upper limits above 100 GeV.

Recently, Sahakyan et al. (2022) explored the connection between PKS 0735+178 and IceCube-211208A and reported on Fermi-LAT, Swift, and optical observations of the blazar. They compared a proton-synchrotron model and two lepto-hadronic models to explain the spectral energy distribution (SED) of the source and found that a lepto-hadronic model with PeV protons (proton luminosity about three times the Eddington luminosity) interacting with an external target photon field yielded the highest neutrino rate (0.067 neutrinos in three weeks). A similar study including NuSTAR data was carried out by Prince et al. (2023), who argued through numerical modeling that neutrino production is insufficient without external photons, assuming the neutrino event and the photons are associated.

In this work, we report on multiwavelength observations of the blazar PKS 0735+178, including those from NuSTAR, VERITAS, and H.E.S.S., contemporaneous with the IceCube astrophysical neutrino candidate IceCube-211208A. Since the blazar is outside of the 90% error region of the IceCube event by ~0°1, and there are a large number of γ-ray blazars that exhibit strong flares, it is uncertain whether or not there is an association between the flaring events across the electromagnetic spectrum and the neutrino event. We interpret the broadband SED of the source in the context of both leptonic

<footnote n="68">⁶⁸ https://gcn.gsfc.nasa.gov/notices_amon_g_b/136015_21306805.amon</footnote>
and lepto-hadronic models, and discuss whether the neutrino event could originate from the blazar.

## 2. Observations and Data Analysis

### 2.1. VERITAS

VERITAS is an array of four IACTs located in southern Arizona (31°40′30″N, 110°57′07″W, 1.3 km a.s.l.; Park 2016). It is capable of detecting γ-rays with energies from 85 GeV to >30 TeV, with an energy resolution of ~15% (at 1 TeV) and a point-spread function of ~0.1 (68% containment at 1 TeV). A point source of 1% the Crab Nebula flux is detectable by VERITAS at a statistical significance of five standard deviations ($5\sigma$) in ~25 hr.

Independently of the IceCube trigger, the blazar PKS 0735+178 was previously observed by VERITAS between 2010 December and 2011 February for 5.2 hr. The source was not detected and a differential flux upper limit at 99% confidence level (CL) of $9 \times 10^{-12} \text{cm}^{-2}\text{s}^{-1}\text{TeV}^{-1}$ at 260 GeV was derived (Archambault et al. 2016).

In response to the IceCube and multiwavelength alerts, VERITAS collected about 20 hr of quality-selected data on PKS 0735+178 between 2021 December 9 and 2022 January 8, at an average zenith angle of 20°2. Observations were performed using the standard “wobble” observation mode (Fomin et al. 1994) with a 0°5 offset. The source was not significantly detected. A statistical significance of $3.4\sigma$ on the excess was derived using events within a 0°1 region around the source. The integral flux upper limit above 220 GeV, which is the highest energy threshold among these observations, is $2.55 \times 10^{-12} \text{cm}^{-2}\text{s}^{-1}$ at 99% CL assuming a power-law spectrum with a photon index of 3 (following Rolke et al. 2005a).

The new VERITAS data were analyzed using the software described in Cogan (2008), with a shower-image template maximum-likelihood reconstruction method and cuts optimized for lower-energy γ-ray events (see, e.g., Christiansen 2017), and independently confirmed with another analysis software described in Maier & Holder (2017). Two nearby Fermi-LAT sources, 4FGL J0738.4+1539 and 4FGL J0743.1+1713, were excluded from the background estimation with an exclusion radius of 0°3.

### 2.2. H.E.S.S.

Located in the Khomas Highland of Namibia (23°16′18″ south, 16°30′00″ east), at an elevation of 1800 m above sea level, H.E.S.S. is the only IACT array in the Southern Hemisphere (Aharonian et al. 2006). It consists of four 12 m telescopes (CT1-4) placed in a square with a side length of 120 m, with a field of view (FoV) of 5°. A fifth 28 m telescope (CT5; Bolmont et al. 2014) was added in 2012 at the center, but it is not used in this analysis.

Since 2012, H.E.S.S. has conducted a neutrino ToO program searching for spatial and temporal correlations between neutrinos and very-high-energy γ-ray emission. Triggered by the ToO alert from IceCube on 2021 December 8, H.E.S.S. observed in the direction of IceCube-211208A for a total of 16 hr from 2021 December 8 to 15, at an average zenith angle of 42°2. Only 3.8 hr of data were selected based on strict criteria on weather conditions and instrumental status. Observations were performed in wobble mode at an offset from the center of the camera of 0°5 (Aharonian et al. 2006).

The data were analyzed using the method described in de Naurois & Rolland (2009) with standard γ-hadron separation and event selection cuts. A circular region of interest (RoI) of 0°1 centered on the position of PKS 0735+178 was defined, and two regions of 0°25 radius around the nearby sources 4FGL J0738.4+1539 and 4FGL J0743.1+1713 were excluded from the background estimation. The background was determined using the standard “reflected background” technique (Berge et al. 2007). The results were validated by a second analysis which uses an independent event calibration and reconstruction (Parsons & Hinton 2014).

No significant γ-ray excess above the expected background was detected from the direction of PKS 0735+178. The integral flux upper limit at 99% CL in the energy range between 0.1 TeV and 10 TeV is $1.82 \times 10^{-11} \text{cm}^{-2}\text{s}^{-1}$ assuming a power-law spectrum with a photon index of 3 (following Rolke et al. 2005). For the upper limit calculations, the minimum energy was chosen as the energy where the effective area reaches 10% of its maximum value, while the maximum energy was chosen such that the number of background events $N_{\text{OFF}} \ge 10$.

### 2.3. Fermi-LAT

The Large Area Telescope (LAT) on board the Fermi satellite is a pair-conversion γ-ray telescope sensitive to γ-rays with energies from ~20 MeV to >300 GeV from a >2 sr FoV (Atwood et al. 2009).

A binned likelihood analysis was performed for the Fermi-LAT data using the Fermi Science Tools version 11 May 3$^{69}$ and FermiPy version 1.0.1$^{70}$ (Wood et al. 2017) in conjunction with the latest PASS 8 instrument response functions (Atwood et al. 2013). Photon events were selected from an energy range between 100 MeV and 300 GeV and an RoI of radius 15° centered on the location of PKS 0735+178, within a maximum zenith angle of 90°. A spatial binning of 0°1 pixel$^{-1}$ and four logarithmic energy bins per decade were used. The initial model consisted of all sources within 20° of the center of the RoI, based on the spatial positions obtained from the 4FGL Data Release 2 (4FGL-DR2; Abdollahi et al. 2020) catalog, and the templates for isotropic and Galactic diffuse emission, iso_P8R3_SOURCE_V2_v1.txt and gll_iem_v07.fits, respectively. An iterative optimization of the model was performed, removing weak sources with $\text{TS} < 10$ and searching for any additional point sources ($\text{TS} \ge 10$) not accounted for in the 4FGL-DR2 catalog at each iteration.

PKS 0735+178 is included in the 4FGL catalog as 4FGL J0738.1+1742, with a peak energy $E_{\text{peak}} \approx (2.1 \pm 0.3)$ GeV in the $\nu F_{\nu}$ representation estimated from the preferred log-parabola model $dN/dE = N_0(E/E_{\text{b}})^{-\alpha-\beta \log(E/E_{\text{b}})}$, where $E_{\text{b}}$ is the energy scale, $E$ is the energy, $\alpha + \beta \log(E/E_{\text{b}})$ is the energy-dependent photon index, and $N_0$ is the normalization. In this work, the Fermi-LAT spectrum of the source between 2021 December 1 and 28 was fit with a log-parabola model with $E_{\text{b}}$ fixed at 1.54 GeV. The best-fit parameters are $N_0 = (1.95 \pm 0.16) \times 10^{-11} \text{MeV}^{-1} \text{cm}^{-2} \text{s}^{-1}$, $\alpha = 2.03 \pm 0.05$, and $\beta = 0.04 \pm 0.03$, yielding the best-fit peak energy as $E_{\text{peak}} \approx 0.7$ GeV and its $1\sigma$ CL interval as roughly 200 MeV–2.5 GeV. The uncertainties are statistical only. Monthly and daily binned light curves between 100 MeV and 300 GeV were

<footnote>$^{69}$ http://fermi.gsfc.nasa.gov/ssc/data/analysis/software</footnote>

<footnote>$^{70}$ http://fermipy.readthedocs.io</footnote>
made from likelihood analyses, freezing the spectral curvature parameters $\alpha$ and $\beta$ while leaving the normalization free.

To check for GeV spectral variability, we conducted a likelihood analysis with daily bins between 2021 December 1 and 28, leaving $\alpha$ and $\beta$ free. The best-fit values for $\alpha$ and $\beta$ are consistent with being constant with probabilities (p-values) of 0.55 and 0.88, respectively. The curvature parameter, $\beta$, is poorly constrained in the daily spectra due to limited statistics and is consistent with 0. On 2021 December 10 and 13, the best-fit $\alpha$ values are $2.30 \pm 0.32$ and $2.07 \pm 0.19$, respectively, in agreement with the average value of $2.03 \pm 0.05$ used for modeling.

## 2.4. NuSTAR

NuSTAR is a hard X-ray space telescope consisting of two coaligned optics and two focal plane detectors (FPMA and FPMB) covering a $13' \times 13'$ FoV. NuSTAR is sensitive to 3-79 keV photons with an on-axis point-spread function of 18" (FWHM; Harrison et al. 2013).

Two NuSTAR ToO observations of PKS 0735+178 were triggered by the IceCube-211208A alert, in combination with the flux increase from PKS 0735+178 detected by Swift-XRT and Fermi-LAT. Observation 1 (ObsID = 80701621002, exposure = 22 ks) was performed on 2021 December 11, and observation 2 (ObsID=80701621004, exposure = 22 ks) was performed on 2021 December 13. The two observations were separated by two days to test the flux and spectral variability of PKS 0735+178.

The NuSTAR data were processed using NuSTAR data analysis software (NuSTARDAS) version 2.1.1 contained within HEASOFT version 6.29 along with the NuSTAR calibration database (CALDB) version 20211202. The source and background spectra were extracted from a circular region of radius 30" and a box region of $2' \times 2'$, respectively. The source spectra were binned such that each bin has a minimum of 50 counts. Using XSPEC (Arnaud 1996), the binned spectra from FPMA and FPMB of each observation were simultaneously fit between 3 and 40 keV beyond which the background begins to dominate. Observation 1 yielded 493/517 (net/total) counts for FPMA and 469/521 counts for FPMB. Observation 2 yielded 386/412 counts for FPMA and 395/433 counts for FPMB.

The spectra were well fit by a single power law with a constant factor (cross-normalization) to account for the difference between the two detectors (cons\*po). Adding an absorption component (cons\*tbabs\*po) does not improve the fit quality or constrain the hydrogen column density. The total neutral hydrogen column density is low in the direction of the source ($4.48 \times 10^{20}$ atoms cm$^{-2}$; Willingale et al. 2013),<fnref n="71" /> thus the absorption in the NuSTAR hard X-ray band is negligible. The best-fit photon index $\Gamma$ of observation 1 is $1.85 \pm 0.06$ (1$\sigma$ statistical error) with a reduced $\chi^2 = 1.08$ for 31 degrees of freedom (dof). Observation 2 exhibits slightly harder spectrum than observation 1, yielding a best-fit $\Gamma = 1.70 \pm 0.07$ with a reduced $\chi^2 = 1.06$ for 26 dof.

## 2.5. Swift-XRT

The X-Ray Telescope (XRT) on the Neil Gehrels Swift Observatory, sensitive to energies from ~0.2 keV to 10 keV

(Gehrels et al. 2004; Burrows et al. 2005), observed PKS 0735+178 nine times in photon counting mode for a total exposure of ~15.2 ks from 2021 December 10 to 2022 January 6 in response to the IceCube alert. Prior to this event, 12 Swift observations were taken between 2009 December and 2011 October, with a total XRT exposure of ~21.3 ks. The Swift-XRT light curve was retrieved from the public online tool "the Swift-XRT data products generator"<fnref n="72" /> (Evans et al. 2007, 2009) and is shown in Figure 1. The Swift-XRT spectral analysis was performed on data taken on 2021 December 10 and 13 (ObsID 00036372014 and 00036372016, respectively) using HEASOFT version 6.29 and CALDB version 20210915. The source and background spectra were extracted from a circular region of radius 20 pixels and an annulus region with inner and outer radii of 70 and 120 pixels, respectively. The source spectra for the observation on December 10 and 13 were grouped requiring a minimum of 20 and 10 counts per bin, respectively. An absorbed power law with the neutral hydrogen column density frozen at $4.48 \times 10^{20}$ atoms cm$^{-2}$ was used to fit the spectra, yielding photon indices of $2.71 \pm 0.09$ ($\chi^2$/dof = 11.8/18) and $2.68 \pm 0.16$ ($\chi^2$/dof = 18.3/15), respectively. Note in Figure 2 that the Swift-XRT spectrum on December 13 captures the low-energy end of the high-energy SED component above a few keV, consistent with the NuSTAR measurements on the same night.

## 2.6. Swift-UVOT

The Swift Ultraviolet/Optical Telescope (UVOT; Roming et al. 2005b) is capable of detecting optical to UV photons with six filters with central wavelengths of 5468 Å (V), 4392 Å (B), 3465 Å (U), 2600 Å (UVW1), 2246 Å (UVM2), and 1928 Å (UVW2). The Swift-UVOT observations of PKS 0735+178 were analyzed using HEASOFT version 6.29 and CALDB version 20211108. A source region with a radius of 5.0 and a background region of the same size were used to extract signal and background counts. The magnitude of the source was then computed using uvotsource and converted to flux using the zero-point for each of the UVOT filters from Poole et al. (2008). An extinction correction was applied following Roming et al. (2009), using the color excess $E(B-V) = 0.0292$ (Schlafly & Finkbeiner 2011).

## 2.7. Other Multiwavelength Facilities

PKS 0735+178 was monitored by many optical facilities. For this study, we used the optical magnitudes of PKS 0735+178 in the g band from the publicly available aperture photometry results from the All-Sky Automated Survey for Supernovae (ASAS-SN; Shappee et al. 2014; Kochanek et al. 2017),<fnref n="73" /> in the R band from the Asteroid Terrestrial-impact Last Alert System (ATLAS; Heinze et al. 2018; Tonry et al. 2018), in the B and R bands from the Automatic Telescope for Optical Monitoring<fnref n="74" /> (ATOM; Hauser et al. 2004), and in the Ks band reported in an Astronomer's Telegram from the Nordic Optical Telescope (Lindfors et al. 2021). We also used the 37 GHz results from the Owens Valley Radio Observatory (Kadler et al. 2021). The ASAS-SN, ATLAS, and ATOM light curves are shown in Figure 1, and the fluxes in all bands mentioned above measured between one and eight days after

<footnote n="71">⁷¹ https://www.swift.ac.uk/analysis/nhtot/index.php</footnote>

<footnote n="72">⁷² https://www.swift.ac.uk/user_objects/index.php</footnote>

<footnote n="73">⁷³ https://asas-sn.osu.edu/</footnote>

<footnote n="74">⁷⁴ https://www.lsw.uni-heidelberg.de/projects/hess/ATOM/</footnote>
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 1.** The light curves of PKS 0735+178 measured by Fermi-LAT, NuSTAR, Swift-XRT, ASAS-SN, ATLAS, and ATOM. (a) The long-term variability of PKS 0735+178. The flux above 100 MeV measured by Fermi-LAT starting in 2008 (top) is binned by month, the count rate between 0.3 and 10 keV measured by Swift-XRT (middle) and the optical *B*-, *g*-, *V*-, and *R*-band magnitudes from ASAS-SN, ATLAS, and ATOM (bottom) are binned by observations. The dashed lines show the time of the neutrino event IceCube-211208A. (b) The light curves focusing on a ~2 month period after the neutrino event. The Fermi-LAT flux is binned nightly. The vertical dashed lines show the time of the neutrino event, the dotted lines show the times of the two NuSTAR observations presented in this work and the green hatch-filled regions show the period during which the VERITAS and H.E.S.S. data were taken. The dashed curve in the middle panel shows the best-fit model with three exponential components used to describe the Swift-XRT data.</figcaption></figure>

IceCube-211208A are shown in the SEDs in Figures 2 and 3. Archival data from the radio to UV bands were taken from the ASDC SED Builder Tool of Italian Space Agency (Stratta et al. 2011), including data from various catalogs and databases (Kuehr et al. 1981; Ochsenbein et al. 2000; Wright et al. 2009, 2010; Planck Collaboration et al. 2011; Giommi et al. 2012).

## 3. Results

### 3.1. Temporal Variability

The monthly γ-ray and the observation-wise X-ray and optical light curves are shown in Figure 1(a), which illustrate that around the time of the neutrino event IceCube-211208A, the blazar PKS 0735+178 exhibited the highest fluxes in all these data sets. The flux above 100 MeV reached a peak of (3.1 ± 0.3) × 10⁻⁷ photon cm⁻² s⁻¹ in the monthly bin centered at 2021 December 24, roughly seven times higher than the average flux in the 4FGL catalog and twice as bright as the previously highest state in 2014.

Figure 1(b) shows the daily γ-ray and observation-wise X-ray and optical light curves during a ~2 month period after the IceCube event, focusing on the highest-flux state of the source. The soft X-ray flux exhibited daily variability, starting at the highest count rate of 0.28 ± 0.01 s⁻¹ on 2021 December 10, followed by a rapid decay down to about 0.1 s⁻¹ on December 12 and 13 (with a flux halving time of 2.2 ± 0.4 days), and a rapid rise to the second highest count rate of 0.26 ± 0.01 s⁻¹ on December 17 (with a flux doubling time of 0.8 ± 0.4 days, which was used to constrain the size of

the emitting region in the 4). A model with three exponential components (two decaying and one increasing) was used to characterize the soft X-ray flux variability (see Figure 1(b)). The peak times were fixed to the two XRT observations with the highest count rate, and five parameters, two peak times and three variability timescales, were left free. The optical observations also revealed a decay in the brightness of the blazar by ~1 magnitude within three weeks after IceCube-211208A, while exhibiting daily variability with a smaller amplitude.

Daily variability above 100 MeV was present during this highest-flux period, with a constant-flux model yielding a poor fit with $\chi^2 = 79$ for 25 dof, which is consistent with the variability above 300 MeV reported by Sahakyan et al. (2022). The strongest evidence for fast GeV γ-ray variability comes from the low flux measured on MJD 59557 and the immediately following high fluxes measured on MJD 59558 and 59559. The variability timescale, however, is model dependent and not well constrained given the large uncertainty in the flux measurement. A fit to a model with three flux peaks yielded the shortest best-fit variability time of 0.38 ± 0.46 days, consistent with the X-ray variability mentioned above and the value of 0.35 days reported by Prince et al. (2023). No evidence for short-term spectral variability was found in the GeV spectrum (see Section 2.3).

No hard X-ray variability was observed between the two NuSTAR observations on December 11 and 13, the 3–40 keV fluxes of which are (3.2 ± 0.2) and (3.0 ± 0.2) in units of 10⁻¹² erg cm⁻² s⁻¹, respectively.
<figure><img src="image_5.png" /><figcaption>**Figure 2.** The broadband SED of PKS 0735+178 in 2021 December along with archival data from radio to UV bands (shown as gray filled circles). The VERITAS and H.E.S.S. spectra were averaged over the entire period of the ToO observations (see Section 2), with the observed SED and those corrected for EBL absorption shown as the unfilled and filled symbols, respectively. The Fermi-LAT spectrum was averaged over a four week period from 2021 December 1 to 28. The magenta shaded region shows the best-fit log-parabola model to the Fermi-LAT spectrum and its $1\sigma$ uncertainty.</figcaption></figure>

## 3.2. Broadband SED

The broadband SED of PKS 0735+178 contemporaneous with the IceCube neutrino event is shown in Figure 2. The SED includes Fermi-LAT, VERITAS, and H.E.S.S. data averaged over ~1 month after IceCube-211208A, and data from X-rays down to optical bands for each observation on a given night between 2021 December 9 and 16. Among these observations, an approximately simultaneous SED can be constructed using data from NuSTAR, Swift, and ATLAS taken on 2021 December 13. Given the daily variability, as discussed in the previous section, the SED on December 13 should be the focus of the numerical modeling (discussed in Section 4.4).

The synchrotron and the high-energy peak frequencies are loosely constrained by the rather flat optical/UV and GeV γ-ray spectra. Nominally, the peak energy flux is observed at a few times $10^{14}$ Hz and at ~$10^{24}$ Hz, respectively, but the spectral cutoff occurs at a somewhat higher frequency, in particular for the γ-ray component. The measured soft and hard X-ray spectra fully constrain the tail of the synchrotron emission and the beginning of the high-energy component of the SED, respectively. The transition occurs at a few keV. To be noted from the SED is the steep decline of the flux near 100 GeV. The VERITAS upper limit at 330 GeV is particularly relevant. After

correction for absorption by extragalactic background light (EBL; shown as the darker red squares in Figure 2; Domínguez et al. 2011), the flux limit is about a factor of 10 below the log-parabola extrapolation of the Fermi-LAT spectrum. There are at least three ways to explain this finding:

1. There could be an intrinsic cutoff at about 100 GeV in the observed photon energy.

2. The redshift $z = 0.45 \pm 0.06$ (Nilsson et al. 2012) could be an underestimate, and the true value might be $z \approx 0.8$. This would bring the upper limit derived by VERITAS, after correcting for EBL absorption, in line with the power-law extrapolation by Fermi-LAT. The redshift of the source was derived assuming that BL Lacertae host galaxies can be used as standard candles. A lower limit of 0.424 was derived from a Mg II absorption-line doublet (Carswell et al. 1974). The identity of the absorber toward PKS 0735+178 is yet to be determined.

3. There could be significant γ-ray absorption within the broad-line or narrow-line region (BLR/NLR), whichever exists, which would require γ-ray production within a parsec or so from the central engine.

In Section 4 we shall discuss the observed SED in the context of the neutrino association.
<figure><img src="image_6.png" /><figcaption>**Figure 3.** The lepto-hadronic model with an external photon field, as described in Section 4.4. The solid blue, orange, and green lines are the synchrotron, SSC, and external Compton components of the leptonic emission from the primary electrons, respectively. The dashed red line is the combined synchrotron and IC emission from Bethe-Heitler pairs. The dotted and dashed-dotted lines are the neutral and charged pion cascade components, respectively. The solid gray line is the sum of all electromagnetic components, and the solid pink line above a few TeV is the neutrino emission. The yellow star marks the nominal flux of 1.5, 170 TeV neutrinos per year using an effective area of 30 m² (see Sections 4.2.1 and 4.4) to guide the eye.</figcaption></figure>

# 4. Discussion

## 4.1. Model Basics and SSC

We shall first discuss the implications of an intrinsic γ-ray cutoff at 100 GeV, assuming a one-zone synchrotron self-Compton (SSC) scenario. The synchrotron spectrum seems to peak at a few $10^{14}$ Hz in frequency, or $\epsilon_{peak} = \epsilon_{peak,eV}$ eV in photon energy, where $\epsilon_{peak,eV}$ is close to 1. The synchrotron peak is quite broad and loosely constrained by the rather flat optical to UV spectrum, with UVOT data suggesting an optical/near-UV flux that is a slowly declining function of frequency. Beyond the near-UV band, there are no measurements until the soft X-ray band. The Swift-XRT flux at 300 eV is a factor of 30 below the UV flux at 3 eV, indicating that the cutoff energy is likely below 100 eV, or $\epsilon_{peak,eV} \lesssim 100$, and likely far lower based on the slowly declining UV spectrum.

There has to be an SSC contribution to the SED. If the Thomson limit applies, the ratio between the IC and the synchrotron peak frequencies (100 GeV/$\epsilon_{peak}$) gives the square of the peak Lorentz factor of electrons, leading to

$$ \gamma_{peak} \approx \frac{3 \times 10^5}{\sqrt{\epsilon_{peak,eV}}}. \qquad (1) $$

The Thomson limit at $\epsilon_{\gamma} = 100$ GeV requires $\gamma \gg \epsilon_{\gamma}/(Dm_e c^2) \approx (8 \times 10^3)/D_{25}$, where $D = 25 D_{25}$ is the Doppler factor, which renders the Thomson limit applicable for a

synchrotron peak up to the hard-UV band. If the true synchrotron peak energy were close to or beyond this limit, which is unlikely given the UVOT and XRT observations, the Klein-Nishina transition would add to the cutoff at 100 GeV.

In the following, any quantity measured in the jet frame is denoted with a prime (e.g., $R'$). For the nominal redshift (0.45), the synchrotron peak frequency and the peak Lorentz factor of electrons constrain the magnetic-field strength in the emission zone

$$ B'D \approx (10^{-3} \text{ G}) \epsilon_{peak,eV}^2. \qquad (2) $$

Equation (2) indicates that the magnetic field would have to be very weak, unless the true synchrotron peak frequency were large, $\epsilon_{peak,eV} \gg 1$. Alternatively, the γ-ray cutoff may be at a few GeV, but that would leave the observed flux at 10 GeV–100 GeV unexplained, and so a weak magnetic field is needed in the one-zone SSC scenario. At the same time, the jet-frame photon energy density at the synchrotron peak would have to be a factor of two or three larger than that in the magnetic field to reproduce the $\nu F_{\nu}$ flux ratio of the IC and the synchrotron peak. We shall now discuss whether or not that is possible.

Assuming isotropic emission in the jet frame, at the distance $d \approx 8 \times 10^{27}$ cm (see Equation (28) in Planck Collaboration et al. 2020) from the observer, the intrinsic photon density per logarithmic energy interval and the energy density in the
radiation field are

$$U' = \frac{3 d^2}{c R'^2 D^4} \nu F_{\nu} = \epsilon' n'_{\ln \epsilon}. \quad (3)$$

Because of the flat observed optical spectrum in the $\nu F_{\mu}$ representation, the energy density of the synchrotron peak should be roughly similar to that at $\epsilon \approx 1$ eV, corresponding to $\epsilon' \approx D^{-1}$ eV. We find

$$U' \approx \frac{5 \times 10^{-3} \text{ erg cm}^{-3}}{D_{25}^{4} R_{16}^{2}}, \quad (4)$$

where we denote the radius of the emission zone as $R' = R_{16}(10^{16}$ cm).

The energy density in the magnetic field (see Equation (2)) required for the SSC peak to be close to 100 GeV and the synchrotron turnover at $\epsilon_{peak,eV}$ is

$$U'_{B} \approx (6 \times 10^{-11} \text{ erg cm}^{-3}) \frac{\epsilon_{\text{peak,eV}}^{4}}{D_{25}^{2}}, \quad (5)$$

which agrees with the energy density of the photon field at the synchrotron peak only if the synchrotron cutoff energy were high, $\epsilon_{peak,eV} \gtrsim 10$, and likewise the size of the emission zone, $R_{16} \gtrsim 10$. More precisely, our new TeV-band data require that

$$\epsilon_{\text{peak,eV}} \approx \frac{100}{\sqrt{D_{25} R_{16}}}, \quad (6)$$

otherwise the high-energy component of the SSC scenario will roll over far below 100 GeV. As the UVOT data suggest $\epsilon_{peak,eV} \approx 1$, and the source radius cannot be made arbitrarily large while keeping SSC dominant, the SSC model has difficulty reproducing a 100 GeV cutoff in the $\gamma$-ray spectrum.

The timescale of the flux variability at $\epsilon_{peak,eV}$ should be at least as large as the total electron cooling time in the observer frame, here about a third of the synchrotron cooling time, and be commensurate with the observed daily flux variations (see Section 3.1). We find in the observer frame

$$t_{\text{loss}} \gtrsim \frac{1}{3} t_{\text{syn}} \approx (0.7 \text{ hr}) D_{25}^{11/4} R_{16}^{7/4}, \quad (7)$$

where in the last step we used Equations (1), (2), and (6). Evidently, a very large emission zone will not allow rapid energy losses. We should add that the flat $\nu F_{\nu}$ spectrum between the optical and the peak at $\epsilon_{peak,eV}$ is likely explained by cooling, which requires an active state that is at least as long as the loss time of electrons radiating in the optical band. This time is about ten times longer than that at the best-fitting cutoff frequency (see Equation (6)), but still shorter than one day. The fastest observed variability timescale is $\sim 1$ day at soft X-ray energies, at which the emitting electrons have an even shorter synchrotron cooling time than that at $\epsilon_{peak,eV}$. The causality argument limits the radius to $R_{16} \lesssim 5 D_{25}$ for daily variations. Diffusive escape primarily affects the IC component in an SSC model and here is sufficiently slow, even if the diffusion coefficient is a few hundred times the Bohm rate.

We conclude that for a synchrotron peak that extends far into the UV band, the SED and in particular the cutoff near 100 GeV could be reasonably well explained with a one-zone SSC model. But if the synchrotron peak indeed lies in the optical, as the data appear to show, the high-energy $\gamma$-ray emission would

fall off at significantly lower energy than is observed, $\epsilon_{IC,peak} \ll 100$ GeV, and the SSC model has difficulty reproducing the observed $\gamma$-ray spectrum.

## 4.2. Viability of Neutrino Emission

Since a simple SSC model cannot explain the SED of the radiation from PKS 0735+178 unless the true synchrotron peak is hidden in the unobserved far-UV band, we shall now explore the viability of associated neutrino emission, which cannot be accounted for by the SSC scenario alone. A significant number of accelerated ions or protons may also exist in the jet. Neutrino emission through the pp channel requires a very large particle density in the jet and hence a large jet power (e.g., Pohl & Schlickeiser 2000). More likely, neutrino production occurs through $p\gamma$ interactions (Gaisser et al. 1995). The implications in the case of PKS 0735+178 depend on the origin of the photon field with which the energetic protons would interact (Sahakyan et al. 2022).

### 4.2.1. Cos spatially Produced Target Photons

If there were no absorption of $\gamma$-rays with a few hundred GeV by pair production with ambient photons from the BLR/NLR, the low observed flux upper limit at 350 GeV would strongly constrain the possibility of neutrino production. This is due to moderate EBL absorption at this energy (for the nominal redshift of $z=0.45$), as well as negligible internal absorption, unless $D_{25}^5 R_{16} \lesssim 0.03$. The effective area of IceCube for $\gamma$-ray follow-up Bronze alerts with neutrino energy $E_{\nu} \approx 170$ TeV is $A_{\text{eff}} \approx 30 \text{ m}^2$ (Blaufuss et al. 2019), and one 170 TeV neutrino event within a month would correspond to a neutrino flux of

$$EF_{\nu}(E) \approx 1.5 \times 10^{-10} \text{ erg cm}^{-2} \text{ s}^{-1}, \quad (8)$$

which extends over a decade in neutrino energy. The true neutrino flux, however, should be much lower than that estimate, otherwise the hard X-ray flux from secondary electrons ($p\gamma$ pair production) in PKS 0735+178 would exceed the observed level (e.g., Gao et al. 2019).

The cascade emission of charged and neutral pions in the 100 TeV band provides a $\gamma$-ray flux at a few hundred GeV that is similar to the neutrino flux at 170 TeV (Gao et al. 2017). The VERITAS upper limits are about a factor 30 lower than that. If there is no strong $\gamma$-ray absorption above about 100 GeV, then the expected rate of neutrinos would be only 0.03 per month and the association would likely be a chance event.

The proton energy is about 20 times the energy of the emergent neutrino, assuming interaction through $\Delta$ resonance (e.g., Hümmer et al. 2010). This is justified by (a) the steeply declining number spectrum of target photons between the optical and hard X-rays and (b) as protons of higher energies would also produce neutrinos at higher energy. In fact, at higher energies the IceCube effective area is even larger, thus leading to a large expected number of observed neutrinos at an energy higher than that of the event in question. The proton Lorentz factor in the jet frame is

$$\gamma'_{p} \approx \frac{1+z}{D} \frac{20 E_{\nu}}{m_{p} c^{2}} \approx 2 \times 10^{5} D_{25}^{-1}. \quad (9)$$

The typical energy of the target photons in $p\gamma$ interactions is

$$\epsilon'_{\text{target}} \approx \frac{200 \text{ MeV}}{\gamma'_{p}} \approx D_{25}(1 \text{ keV}). \quad (10)$$
In the observer frame, that would be 25 keV, i.e., in the energy band covered by NuSTAR. Assuming the hard X-rays are produced cospatially with the neutrinos and isotropically in the jet frame, we can estimate the jet-frame photon density as a function of the radius of the emission zone, $R'$. In analogy to the derivation of Equation (4), but for the hard X-rays with flux $\nu F_\nu \approx 2 \times 10^{-12}$ erg cm$^2$ s$^{-1}$, we find for the jet-frame photon density

$$n'_{\ln \epsilon} \approx \frac{2 \times 10^5 \text{ cm}^{-3}}{D_{25}^4 R_{16}^2}. \qquad (11)$$

The cross section of $p\gamma$ interaction is $\sigma_{p\gamma} \approx 5 \times 10^{-28}$ cm$^2$. With an inelasticity of $K=0.2$ (Mücke et al. 2000) we find for the energy-loss time of the radiating protons

$$t'_{p\gamma} \approx \frac{1}{K c n'_{\ln \epsilon} \sigma_{p\gamma}} \approx D_{25}^4 R_{16}^2 (1.7 \times 10^{12} \text{ s}). \qquad (12)$$

In the observer frame, that would correspond to a few thousand years. Protons will diffusively escape from the emission zone in a much shorter time, even for Bohm diffusion, and hence $p\gamma$ interactions are an inefficient loss process. To reach any kind of reasonable output, we need a source radius $R_{16} \lesssim 1$. That would also be called for, if the emitter were located inside the BLR/NLR.

The most important unknown input in an estimate of the proton injection luminosity is the observed duration of the active state, $T_{\text{act}}$. Building on Equation (8) we find the neutrino luminosity that gives one neutrino as

$$L'_{\nu} \approx \frac{6 \times 10^{41} \text{ erg s}^{-1}}{D_{25}^4} \left( \frac{30 \text{ days}}{T_{\text{act}}} \right). \qquad (13)$$

The required proton injection luminosity is at least a factor $t'_{p\gamma}/T'_{\text{act}}$ larger than the neutrino luminosity

$$L'_{p} \approx \frac{R_{16}^2}{D_{25}} \left( \frac{30 \text{ days}}{T_{\text{act}}} \right)^2 (1.5 \times 10^{46} \text{ erg s}^{-1}). \qquad (14)$$

In the AGN frame, the kinetic power of the jet is approximately a factor $D^2$ larger than that and would exceed the Eddington limit (see, e.g., Ghisellini et al. 2014), which is $L_{\text{edd}} \approx 10^{47}$ erg s$^{-1}$ for $M_{\text{BH}} \approx 8 \times 10^8 M_{\odot}$ (Ghisellini et al. 2010), unless the activity state lasts much longer than 30 days, in fact at least half a year. The same problem arises for TXS 0506+056. Given the high X-ray flux (which could largely be produced by secondary electrons), it is possible to reproduce the SED of PKS 0735+178 with a lepto-hadronic scenario, if the power requirement could be reduced (see also Sahakyan et al. 2022).

### 4.2.2. Target Photons from the BLR/NLR

External photon fields may potentially reduce the power requirement. For example, spine-sheath models (e.g., Tavecchio & Ghisellini 2015) require less power but imply a coherent bright sheath that is at least 10 pc long, otherwise the activity cannot last more than a few weeks.

The featureless optical spectrum of PKS 0735+178 suggests that it may be a BL Lacertae object. It has been noted in BL Lacertae objects that the jet may pass through a region harboring a significant jet-external photon field (Foffano et al. 2022) that for simplicity we henceforth refer to as the BLR.

The BLR is an external source of photons, with which high-energy protons could interact and produce neutrinos. If BLR emission would be the dominant source of target photons, then the particle acceleration site cannot be comoving with the jet. If it were, it would leave the BLR within a few days of observed time, as in that time span the jet travels about a parsec for a Doppler factor of $D=25$. This decreases the Doppler boosting of the emission, but relaxes the factor $D^2$ that would otherwise govern the relation between the jet-frame proton injection luminosity and the kinetic power of the jet, and likewise for leptonic models (e.g., Sokolov et al. 2004).

To estimate the photon density inside the BLR, we assume that the accretion disk produces 10% of the observed UV flux, about 10% of which is rescattered. For a BLR of radius $R_{\text{BLR}}=2$ pc and mono-energetic photons of energy $\epsilon_{\text{BLR}} \approx 15$ eV, the photon density within the BLR in the AGN frame is

$$n_{\ln \epsilon} \approx 7 \times 10^5 \text{ cm}^{-3}. \qquad (15)$$

The pair-production opacity within the BLR would be around unity for photons of about 200 GeV (Poutanen & Stern 2010), consistent with a scenario in which the sharp drop of observed flux near 100 GeV is at least partially caused by absorption in the BLR.

In the jet frame, the photon density would increase by the jet Lorentz factor, $\Gamma$, and the mean energy of BLR photons would be

$$n'_{\ln \epsilon} \approx \Gamma_{25}(2 \times 10^7 \text{ cm}^{-3}), \quad \epsilon'_{\text{BLR}} \approx \Gamma_{25}(500 \text{ eV}), \qquad (16)$$

which is reasonably close to that estimated in Equation (10) and hence permits a rescaling of the estimate for internal photon fields.

The energy-loss time of the radiating protons would be

$$t'_{p\gamma} \approx \frac{1.7 \times 10^{10} \text{ s}}{\Gamma_{25}}, \qquad (17)$$

and the proton injection luminosity would be

$$L'_{p} \approx \frac{1.5 \times 10^{44} \text{ erg s}^{-1}}{\Gamma_{25} D_{25}^{5}} \left( \frac{30 \text{ days}}{T_{\text{act}}} \right)^{2}. \qquad (18)$$

The jet luminosity is about a factor $D^2$ larger (about $10^{47}$ erg s$^{-1}$) and in marginal agreement with the Eddington luminosity (Xie et al. 2002).

### 4.3. External Compton Scattering

The emission zone radius is not constrained in an external hadronic model and can be determined by the leptonic emission processes. IC scattering of BLR photons would be strongly Klein–Nishina-suppressed for electrons with a Lorentz factor exceeding a few hundred. Radiation modeling of blazars often requires a hard electron spectrum, $N(\gamma) \propto \gamma^{-5}$ or similar, with Lorentz factors up to a few thousand. The IC output would then be a narrow spectral component at the Klein–Nishina transition, $E_{\text{KN}} \approx 10$ GeV, with a hard spectrum at lower energies, $\nu F_\nu \propto \nu^{0.5}$, that asymptotically transitions to a steep decline well above a few tens of GeV (e.g., Sikora 2009), which is well in line with the new flux limits in the TeV band that we present. The same pair absorption in the BLR as discussed in Section 4.2 also introduces a spectral cutoff around 200 GeV. The Klein–Nishina effect and the pair-absorption process will
both contribute to the steep decline of the observed flux around 100 GeV. The peak flux of the upscattered BLR photons at about 10 GeV can be scaled to the observed flux at 1 MeV, to which the same electrons would up-scatter the optical photons at the synchrotron peak in the Thomson regime.

For our fiducial parameters $D_{25}$ and $R_{16}$ the $\nu F_{\nu}$ flux at 10 GeV from up-scattering BLR/NLR photons may be a factor of 10 higher than the SSC flux at 1 MeV, and even higher for a larger Doppler factor. It then appears viable to produce the observed radiation at 10 to 100 GeV with external Compton scattering, which permits the SSC component to fall off at a few GeV. The corresponding reduction in $\gamma_{\text{peak}}$ (see Equation (1)) permits a stronger magnetic field than estimated in Equation (2). Reproducing a comparable $\nu F_{\nu}$ flux of the synchrotron and the $\gamma$-ray component in the SSC scenario would no longer require a high value of $\epsilon_{\text{peak},eV}$. This alleviates the requirement of a strong unseen UV synchrotron component that in an SSC model we found necessary to explain all the $\gamma$-ray data up to 100 GeV (see Section 4.1). It also implies that the cutoff near 100 GeV is partially caused by pair production with BLR photons and partially by the maximum energy of the external Compton emission component.

### 4.4. Numerical SED Modeling

Starting from the analytical estimates discussed in the previous section, we performed a numerical simulation of the photon and neutrino emission from PKS 0735+178 during the 2021 December flaring event. We used the code described in Cerruti et al. (2015) that computes the emission at equilibrium by a spherical plasmoid in a jet, parameterized by its radius $R'$ and Doppler factor $D$, filled with a homogeneous magnetic field $B'$. The code assumes that primary populations of electrons and protons are present in the emitting region. Hadronic interactions spawn secondary particles (photons from $\pi^0$ decay and electrons/positrons from $\pi^{\pm}$ decay and Bethe-Heitler pair production) in the emitting region; these trigger electron-positron pair cascades that radiate via both synchrotron and IC scattering. Neutrinos instead leave the emitting region with a negligible probability of interaction. Protons interact with both the internal photon fields (synchrotron and IC emission by primary electrons) and with external ones that are parameterized by a graybody spectrum with temperature $T' = 2 \times 10^6$ K, commensurate with the mean photon energy quoted in Equation (16).

Primary protons are assumed to be described by a simple power-law energy distribution with an exponential cutoff, while primary electrons are described by a broken power law with an exponential cutoff. This choice is motivated by the fact that primary electrons are certainly cooled due to synchrotron radiation, while primary protons are supposed here to be not cooled; this hypothesis is confirmed a posteriori by comparing all relevant timescales. The primary electron distribution is thus defined by six free parameters (the two power-law indices $\alpha_{e,1}$ and $\alpha_{e,2}$, the minimum/break/maximum Lorentz factors $\gamma'_{\text{min},e}$, $\gamma'_{\text{br},e}$, and $\gamma'_{\text{max},e}$, respectively, and the normalization number density $n'_e$), and the primary proton distribution by an additional four (the index $\alpha_p$, minimum/maximum Lorentz factors $\gamma'_{\text{min},p}$ and $\gamma'_{\text{max},p}$, respectively, and the normalization number density $n'_p$).

Together with the three parameters ($B'$, $D$, and $R'$) describing the emitting region, and the normalization of the external

<table><caption>**Table 1**

Parameters for the Lepto-hadronic Modeling of the SED</caption><thead><tr><th>Parameter<sup>a</sup></th><th>Value</th></tr></thead><tbody><tr><td>z</td><td>0.45</td></tr><tr><td>D</td><td>30</td></tr><tr><td>B'<sub>e</sub> [G]</td><td>0.3</td></tr><tr><td>R'<sub>e</sub> [cm]</td><td>2.1 &times; 10<sup>16</sup></td></tr><tr><td>n'<sub>e</sub> [1/cm<sup>3</sup>]</td><td>4 &times; 10<sup>3</sup></td></tr><tr><td>n'<sub>e</sub>/n'<sub>p</sub></td><td>1.0</td></tr><tr><td>&alpha;<sub>e,1</sub></td><td>2</td></tr><tr><td>&alpha;<sub>e,2</sub></td><td>3</td></tr><tr><td>&gamma;'<sub>min,e</sub></td><td>100</td></tr><tr><td>&gamma;'<sub>br,e</sub></td><td>3.5 &times; 10<sup>3</sup></td></tr><tr><td>&gamma;'<sub>max,e</sub></td><td>1.3 &times; 10<sup>4</sup></td></tr><tr><td>&alpha;<sub>p</sub></td><td>2.0</td></tr><tr><td>&gamma;'<sub>min,p</sub></td><td>1</td></tr><tr><td>&gamma;'<sub>max,p</sub></td><td>1 &times; 10<sup>7</sup></td></tr><tr><td>T'<sub>bbody</sub> [K]</td><td>2 &times; 10<sup>6</sup></td></tr><tr><td>U'<sub>bbody</sub> [erg cm<sup>-3</sup>]</td><td>0.16</td></tr><tr><td>*U'<sub>e</sub> [erg cm<sup>-3</sup>]</td><td>1.2 &times; 10<sup>-2</sup></td></tr><tr><td>*U'<sub>p</sub> [erg cm<sup>-3</sup>]</td><td>3.6 &times; 10<sup>-3</sup></td></tr><tr><td>*U'<sub>p</sub> [erg cm<sup>-3</sup>]</td><td>93</td></tr><tr><td>*U'<sub>e</sub>/U'<sub>p</sub></td><td>3.5</td></tr><tr><td>*U'<sub>p</sub>/U'<sub>p</sub></td><td>2.6 &times; 10<sup>4</sup></td></tr><tr><td>*L'<sub>jet</sub> [erg s<sup>-1</sup>]</td><td>2.0 &times; 10<sup>48</sup></td></tr><tr><td>*&nu;<sub>BRONZE</sub> [year<sup>-1</sup>]</td><td>1.49</td></tr></tbody></table>

**Note.**

<sup>a</sup> The quantities flagged with an asterisk are derived quantities and not model parameters. The luminosity of the emitting region has been calculated as $L = 2\pi R'c(D/2)^2(U'_b + U'_p + U'_e)$, where $U'_b$, $U'_p$, and $U'_e$, are the energy densities of the magnetic field, the electrons, and the protons, respectively.

photon density $U'_{bbody}$, the total amount of free parameters of the model is 14. Instead of fitting a model with many free parameters, we aimed at a solution that describes the photon SED and maximizes the neutrino output. The indices of the injected primary particles are fixed to $\alpha_{e,1} = \alpha_p = 2.0$, while that of the cooled electrons is fixed to 3.0 as expected from simple synchrotron cooling. The Doppler factor of the emitting region is fixed at 30, a typical value for blazars. The other free parameters are manually optimized to reproduce the SED via dominant SSC and external-IC components in the $\gamma$-ray band, and a subdominant hadronic component. The latter cannot be arbitrarily luminous; it is mainly constrained by the Bethe-Heitler component that emerges in the X-ray band. Neutrino rates are computed using the effective area $A_{\text{eff}} \approx 30 \text{ m}^2$. In Figure 3 we show the SED resulting from the lepto-hadronic modeling, with both the photon and neutrino components.

A detailed list of model parameters is provided in Table 1. The expected neutrino rate is ~1.5 events per year, or 0.125 events in a 30 day activity period, marginally consistent with the detection of a single event by IceCube. The power of the jet needed to produce the amount of photons and neutrinos is $2 \times 10^{48} \text{ erg s}^{-1}$—roughly 20 times higher than the Eddington luminosity of a supermassive black hole with a mass of $10^9 \text{ M}_\odot$, but it is very dependent on the choice of $\alpha_p$ and $\gamma_{p,\text{min}}$. The values assumed here, 2.0 and 1, respectively, are conservative with respect to the total jet power. A harder proton injection or a higher minimum Lorentz factor of protons significantly reduces the energetic requirement and would bring
the particle energy density closer to that of the magnetic field. As an example, an equivalent SED modeling, with similar neutrino rates, is possible assuming $\alpha_p = 1.9$ and $L_{\text{jet}} = 8.5 \times 10^{47}$ erg s$^{-1}$.

## 5. Summary

A few days after the neutrino event IceCube-211208A, detected by IceCube with an energy $E_ν \approx 171$ TeV on 2021 December 8, the blazar PKS 0735+178 was observed at an elevated flux state from the optical to γ-ray bands, reaching the highest soft X-ray and GeV γ-ray fluxes among all measured values since 2007. PKS 0735+178 is located 2°2 away from the best-fit position of the IceCube event, immediately outside the 90% error region (2°13). Its soft X-ray flux exhibited fast variability on a one day timescale. The active flux state of the source lasted from roughly two weeks to a month.

The broadband SED of PKS 0735+178 near the time of the IceCube event was well measured by extensive follow-up observations across the electromagnetic spectrum. In particular, the X-ray data from Swift-XRT and NuSTAR characterize the tail of the synchrotron SED peak and the onset of the high-energy SED peak. The γ-ray data from Fermi-LAT, VERITAS, and H.E.S.S. require a spectral cutoff near 100 GeV, after taking into account EBL absorption.

We analytically demonstrated that the observed SED of the blazar, especially the γ-ray cutoff in the TeV data, constitutes a challenge to a simple one-zone SSC model, unless, in the unlikely scenario, the synchrotron peak extends into the far-UV band that the observations do not cover. Alternatively, it can be explained by an SSC/external Compton scenario, which naturally provides the observed 100 GeV cutoff through the Klein–Nishina effect and γ-γ pair absorption. It could also be explained by a lepto-hadronic mechanism with an external photon field, which marginally agrees with the Eddington limit on the jet luminosity and the observed IceCube neutrino rate.

We presented a numerical lepto-hadronic model with external target photons that is consistent with the observed SED and marginally consistent with the neutrino event. In this model, the electromagnetic emission is dominated by the leptonic components, while the subdominant hadronic components are constrained by the observed X-ray spectrum. The jet power is significantly higher than the Eddington luminosity, unless the spectrum of the parent protons is harder than 2.0. The expected total neutrino rate is 1.5 per year, or 0.125 per month.

Some challenges in searching for neutrino-emitting blazars are common in the cases of PKS 0735+178 and TXS 0506+056; the limited localization precision of the IceCube Observatory and the large number of γ-ray blazars as potential counterparts make the association between neutrino events and γ-ray blazars difficult. The electromagnetic emission can often be explained by leptonic models alone without the need for a hadronic component; the jet power and the proton luminosity required in lepto-hadronic models are often too high, exceeding the Eddington limit, for a short period of activity. These challenges can only be addressed with additional data and continued searches for astrophysical neutrino emitters, including follow-up observations of flaring blazars in temporal and spatial coincidence with astrophysical neutrino events. Such follow-up efforts remain a focus in the multimessenger community (see, e.g., Acciari et al. 2022; Garrappa & Buson 2022; McBride et al. 2022).

## Acknowledgments

This research is supported by grants from the U.S. Department of Energy Office of Science, the U.S. National Science Foundation, and the Smithsonian Institution, by NSERC in Canada, and by the Helmholtz Association in Germany. This research used resources provided by the Open Science Grid, which is supported by the National Science Foundation and the U.S. Department of Energy's Office of Science, and resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility operated under Contract No. DE-AC02-05CH11231. We acknowledge the excellent work of the technical support staff at the Fred Lawrence Whipple Observatory and at the collaborating institutions in the construction and operation of the instrument.

This work was supported by NASA grants 80NSSC22K0573, 80NSSC22K1515, 80NSSC22K0950, 80NSSC20K1587, and 80NSSC20K1494, and NSF grant PHY-1806554.

The support of the Namibian authorities and of the University of Namibia in facilitating the construction and operation of H.E.S.S. is gratefully acknowledged, as is the support by the German Ministry for Education and Research (BMBF), the Max Planck Society, the German Research Foundation (DFG), the Helmholtz Association, the Alexander von Humboldt Foundation, the French Ministry of Higher Education, Research and Innovation, the Centre National de la Recherche Scientifique (CNRS/IN2P3 and CNRS/INSU), the Commissariat à l'énergie atomique et aux énergies alternatives (CEA), the U.K. Science and Technology Facilities Council (STFC), the Irish Research Council (IRC) and the Science Foundation Ireland (SFI), the Knut and Alice Wallenberg Foundation, the Polish Ministry of Education and Science, agreement no. 2021/WK/06, the South African Department of Science and Technology and National Research Foundation, the University of Namibia, the National Commission on Research, Science & Technology of Namibia (NCRST), the Austrian Federal Ministry of Education, Science and Research and the Austrian Science Fund (FWF), the Australian Research Council (ARC), the Japan Society for the Promotion of Science, the University of Amsterdam, and the Science Committee of Armenia grant 21AG-1C085. We appreciate the excellent work of the technical support staff in Berlin, Zeuthen, Heidelberg, Palaiseau, Paris, Saclay, Tübingen, and in Namibia in the construction and operation of the equipment. This work benefited from services provided by the H.E.S.S. Virtual Organisation, supported by the national resource providers of the EGI Federation.

This work made use of data supplied by the UK Swift Science Data Centre at the University of Leicester.

This work has made use of data from the Asteroid Terrestrial-impact Last Alert System (ATLAS) project. The Asteroid Terrestrial-impact Last Alert System (ATLAS) project is primarily funded to search for near earth asteroids through NASA grants NN12AR55G, 80NSSC18K0284, and 80NSSC18K1575; by-products of the NEO search include images and catalogs from the survey area. This work was partially funded by Kepler/K2 grant J1944/80NSSC19K0112 and HST GO-15889, and STFC grants ST/T000198/1 and ST/S006109/1. The ATLAS science products have been made possible through the contributions of the University of Hawaii Institute for Astronomy, the Queens University Belfast, the
Space Telescope Science Institute, the South African Astronomical Observatory, and The Millennium Institute of Astrophysics (MAS), Chile.

**Facilities:** VERITAS, H.E.S.S., Fermi-LAT, NuSTAR, Swift-XRT, Swift-UVOT, ASAS-SN, ATLAS, NOT, and OVRO.

**Software:** Astropy (Astropy Collaboration et al. 2013, 2018), NumPy (van der Walt et al. 2011), Matplotlib (Hunter 2007), SciPy (Jones et al. 2001), VEGAS (Cogan 2008), Eventdisplay (Maier & Holder 2017), Model Analysis (de Naurois & Rolland 2009), Fermitools (Fermi Science Support Development Team 2019), and HEAsoft (Nasa High Energy Astrophysics Science Archive Research Center (Heasarc) 2014).

## ORCID iDs

A. Acharyya https://orcid.org/0000-0002-2028-9230

C. B. Adams https://orcid.org/0000-0002-9021-6192

P. Bangale https://orcid.org/0000-0002-3886-3739

J. T. Bartkoske https://orcid.org/0000-0002-9675-7328

P. Batista https://orcid.org/0000-0001-8138-1391

W. Benbow https://orcid.org/0000-0003-2098-170X

A. Brill https://orcid.org/0000-0002-6208-5244

J. H. Buckley https://orcid.org/0000-0001-6391-9661

M. Errando https://orcid.org/0000-0002-1853-863X

A. Falcone https://orcid.org/0000-0002-5068-7344

Q. Feng https://orcid.org/0000-0001-6674-4238

G. M. Foote https://orcid.org/0000-0002-2944-6060

L. Fortson https://orcid.org/0000-0002-1067-8558

A. Furniss https://orcid.org/0000-0003-1614-1273

W. Hanlon https://orcid.org/0000-0002-0109-4737

D. Hanna https://orcid.org/0000-0002-8513-5603

O. Hervet https://orcid.org/0000-0003-3878-1677

C. E. Hinrichs https://orcid.org/0000-0001-6951-2299

J. Holder https://orcid.org/0000-0002-6833-0474

T. B. Humensky https://orcid.org/0000-0002-1432-7771

W. Jin https://orcid.org/0000-0002-1089-1754

P. Kaaret https://orcid.org/0000-0002-3638-0637

D. Kieda https://orcid.org/0000-0003-4785-0101

T. K. Kleiner https://orcid.org/0000-0002-4260-9186

N. Korzoun https://orcid.org/0000-0002-4289-7106

S. Kumar https://orcid.org/0000-0001-6243-7599

M. J. Lang https://orcid.org/0000-0003-4641-4201

M. Lundy https://orcid.org/0000-0003-3802-1619

G. Maier https://orcid.org/0000-0001-9868-4700

M. J. Millard https://orcid.org/0000-0001-7106-8502

C. L. Mooney https://orcid.org/0000-0001-5937-446X

P. Moriarty https://orcid.org/0000-0002-1499-2667

R. Mukherjee https://orcid.org/0000-0002-3223-0754

S. O'Brien https://orcid.org/0000-0002-9296-2981

R. A. Ong https://orcid.org/0000-0002-4837-5253

M. Pohl https://orcid.org/0000-0001-7861-1707

E. Pueschel https://orcid.org/0000-0002-0529-1973

J. Quinn https://orcid.org/0000-0002-4855-2694

K. Ragan https://orcid.org/0000-0002-5351-3323

D. Ribeiro https://orcid.org/0000-0002-7523-7366

I. Sadeh https://orcid.org/0000-0003-1387-8915

A. C. Sadun https://orcid.org/0000-0001-8086-7242

L. Saha https://orcid.org/0000-0002-3171-5039

M. Santander https://orcid.org/0000-0001-7297-8217

R. Shang https://orcid.org/0000-0002-9856-989X

D. A. Williams https://orcid.org/0000-0003-2740-9714

S. L. Wong https://orcid.org/0000-0002-2730-2733

J. Woo https://orcid.org/0009-0001-6471-1405

M. Backes https://orcid.org/0000-0002-9326-6400

V. Barbosa Martins https://orcid.org/0000-0002-5085-8828

R. Batzofin https://orcid.org/0000-0002-5797-3386

Y. Becherini https://orcid.org/0000-0002-2115-2930

D. Berge https://orcid.org/0000-0002-2918-1824

K. Bernlöhr https://orcid.org/0000-0001-8065-3252

M. Böttcher https://orcid.org/0000-0002-8434-5692

C. Boisson https://orcid.org/0000-0001-5893-1797

J. Bolmont https://orcid.org/0000-0003-4739-8389

F. Bradascio https://orcid.org/0000-0002-7750-5256

M. Breuhaus https://orcid.org/0000-0003-0268-5122

R. Brose https://orcid.org/0000-0002-8312-6930

F. Brun https://orcid.org/0000-0003-0770-9007

T. Bulik https://orcid.org/0000-0003-2045-4803

S. Caroff https://orcid.org/0000-0002-1103-130X

S. Casanova https://orcid.org/0000-0002-6144-9122

M. Cerruti https://orcid.org/0000-0001-7891-699X

A. Chen https://orcid.org/0000-0001-6425-5692

J. Chibueze https://orcid.org/0000-0002-9875-7436

G. Cotter https://orcid.org/0000-0002-9975-1829

S. Dai https://orcid.org/0000-0002-9618-2499

J. Damascene Mbarubucyeye https://orcid.org/0000-0002-4991-6576

A. Djannati-Ataï https://orcid.org/0000-0002-4924-1708

V. Doroshenko https://orcid.org/0000-0001-8162-1105

G. Fichet de Clairfontaine https://orcid.org/0000-0003-1143-3883

M. Filipovic https://orcid.org/0000-0002-4990-9288

G. Fontaine https://orcid.org/0000-0002-6443-5025

S. Funk https://orcid.org/0000-0002-2012-0080

G. Giavitto https://orcid.org/0000-0002-7629-6499

D. Glawion https://orcid.org/0000-0003-4865-7696

J. F. Glicenstein https://orcid.org/0000-0003-2581-1742

T. L. Holch https://orcid.org/0000-0001-5161-1168

M. Jamrozy https://orcid.org/0000-0002-0870-7778

V. Joshi https://orcid.org/0000-0003-4467-3621

B. Khélifi https://orcid.org/0000-0001-6876-5577

S. Klepsen https://orcid.org/0000-0002-8949-4275

K. Kosack https://orcid.org/0000-0001-8424-3621

D. Kostunin https://orcid.org/0000-0002-0487-0076

J.-P. Lenain https://orcid.org/0000-0001-7284-9220

F. Leuschner https://orcid.org/0000-0001-9037-0272

A. Luashvili https://orcid.org/0000-0003-4384-1638

J. Mackey https://orcid.org/0000-0002-5449-6131

D. Malyshev https://orcid.org/0000-0001-9689-2194

V. Marandon https://orcid.org/0000-0001-9077-4058

M. Marchegiani https://orcid.org/0000-0001-7487-8287

G. Martí-Devesa https://orcid.org/0000-0003-0766-6473

R. Marx https://orcid.org/0000-0002-6557-4924

A. Mitchell https://orcid.org/0000-0003-3631-5648

L. Mohrmann https://orcid.org/0000-0002-9667-8654

A. Montanari https://orcid.org/0000-0002-3620-0173

E. Moulin https://orcid.org/0000-0003-4007-0145

T. Murach https://orcid.org/0000-0003-1128-5008

J. Niemiec https://orcid.org/0000-0001-6036-8569

L. Olivera-Nieto https://orcid.org/0000-0002-9105-0518

M. Ostrowski https://orcid.org/0000-0002-9199-7031

S. Panny https://orcid.org/0000-0001-5770-3805

G. Pühlhofer https://orcid.org/0000-0003-4632-4644

M. Punch https://orcid.org/0000-0002-4710-2165
P. Reichherzer https://orcid.org/0000-0003-4513-8241

A. Reimer https://orcid.org/0000-0001-8604-7077

B. Rudak https://orcid.org/0000-0003-0452-3805

E. Ruiz-Velasco https://orcid.org/0000-0001-6939-7825

V. Sahakian https://orcid.org/0000-0003-1198-0043

A. Santangelo https://orcid.org/0000-0003-4187-9560

M. Sasaki https://orcid.org/0000-0001-5302-1866

F. Schüssler https://orcid.org/0000-0003-1500-6571

H. M. Schutte https://orcid.org/0000-0002-1769-5617

J. N. S. Shapopi https://orcid.org/0000-0002-7130-9270

A. Specovius https://orcid.org/0000-0002-1156-4771

S. Spencer https://orcid.org/0000-0001-5516-1205

S. Steinmassl https://orcid.org/0000-0002-2865-8563

I. Sushch https://orcid.org/0000-0002-2814-1257

T. Tanaka https://orcid.org/0000-0002-4383-0368

R. Terrier https://orcid.org/0000-0002-8219-4667

C. van Eldik https://orcid.org/0000-0001-9669-645X

J. Veh https://orcid.org/0000-0003-4736-2167

A. Wierzcholska https://orcid.org/0000-0003-4472-7204

M. Zacharias https://orcid.org/0000-0001-5801-3945

D. Zargaryan https://orcid.org/0000-0002-2876-6433

A. A. Zdziarski https://orcid.org/0000-0002-0333-2452

S. Zouari https://orcid.org/0000-0002-5333-2004

N. Żywucka https://orcid.org/0000-0003-2644-6441

## References

Aartsen, M. G., Abbasi, R., Abdou, Y., et al. 2013, *Sci*, 342, 1242856

Aartsen, M. G., Abraham, K., Ackermann, M., et al. 2017, *ApJ*, 835, 45

Aartsen, M. G., Ackermann, M., Adams, J., et al. 2018a, *Sci*, 361, eea1378

Aartsen, M. G., Ackermann, M., Adams, J., et al. 2018b, *Sci*, 361, 147

Aartsen, M. G., Ackermann, M., Adams, J., et al. 2020, *PhRvL*, 124, 051103

Abbasi, R., Ackermann, M., Adams, J., et al. 2021, *ApJL*, 920, L45

Abbasi, R., Ackermann, M., Adams, J., et al. 2022, *Sci*, 378, 538

Abdollahi, S., Acero, F., Ackermann, M., et al. 2020, *ApJS*, 247, 33

Acciari, V. A., Ansoldi, S., Antonelli, L. A., et al. 2022, *ICRC*, 37, 960

Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2006, *A&A*, 457, 899

Archambault, S., Archer, A., Benbow, W., et al. 2016, *AJ*, 151, 142

Arnaud, K. A. 1996, in ASP Conf. Ser. 101, Astronomical Data Analysis Software and Systems V, ed. G. H. Jacoby & J. Barnes (San Francisco, CA: ASP), 17

Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, *AJ*, 156, 123

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, *A&A*, 558, A33

Atwood, W., Albert, A., Baldini, L., et al. 2013, *arXiv:1303.3514*

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

Berge, D., Funk, S., & Hinton, J. 2007, *A&A*, 466, 1219

Blaufuss, E., Kintscher, T., Lu, L., & Tung, C. F. 2019, *ICRC*, 36, 1021

Bolmont, J., Corona, P., Gauron, P., et al. 2014, *NIMPA*, 761, 46

Böttcher, M., Reimer, A., Sweeney, K., & Prakash, A. 2013, *ApJ*, 768, 54

Burrows, D. N., Hill, J. E., Nousek, J. A., et al. 2005, *SSRv*, 120, 165

Carswell, R. F., Strittmatter, P. A., Williams, R. E., Kinman, T. D., & Serkowski, K. 1974, *ApJL*, 190, L101

Cerruti, M., Zech, A., Boisson, C., & Inoue, S. 2015, *MNRAS*, 448, 910

Christiansen, J. 2017, *ICRC (Busan)*, 35, 789

Cogan, P. 2008, *ICRC*, 3, 1385

D'Ammando, F. 2021, *ATel*, 15109, 1

de Naurois, M., & Rolland, L. 2009, *APh*, 32, 231

Domínguez, A., Primack, J. R., Rosario, D. J., et al. 2011, *MNRAS*, 410, 2556

Dzhilkibaev, Z. A., Suvorova, O. & Baikal-GVD Collaboration 2021, *ATel*, 15112, 1

Evans, P. A., Beardmore, A. P., Page, K. L., et al. 2007, *A&A*, 469, 379

Evans, P. A., Beardmore, A. P., Page, K. L., et al. 2009, *MNRAS*, 397, 1177

Feng, Q., Jin, W., Mori, K., et al. 2021, *ATel*, 15113, 1

Fermi Science Support Development Team 2019, *Fermitools: Fermi Science Tools, Astrophysics Source Code Library*, ascl:1905.011

Filippini, F., Illuminati, G., Heijboer, A., et al. 2022, *ATel*, 15290, 1

Foffano, L., Vittorini, V., Tavani, M., & Menegoni, E. 2022, *ApJ*, 926, 95

Fomin, V. P., Stepanian, A. A., Lamb, R. C., et al. 1994, *APh*, 2, 137

Gaisser, T. K., Halzen, F., & Stanev, T. 1995, *PhR*, 258, 173

Gao, S., Fedynitch, A., Winter, W., & Pohl, M. 2019, *NatAs*, 3, 88

Gao, S., Pohl, M., & Winter, W. 2017, *ApJ*, 843, 109

Garrappa, S., & Buson, S. 2022, *ATel*, 15570, 1

Garrappa, S., Buson, S., Sinapius, J., & Kadl, M. 2021, *ATel*, 15099, 1

Gehrels, N., Chincarini, G., Giommi, P., et al. 2004, *ApJ*, 611, 1005

Ghisellini, G., Tavecchio, F., Foschini, L., et al. 2010, *MNRAS*, 402, 497

Ghisellini, G., Tavecchio, F., Maraschi, L., Celotti, A., & Sbarrato, T. 2014, *Natur*, 515, 376

Giommi, P., Polenta, G., Lähteenmäki, A., et al. 2012, *A&A*, 541, A160

Harrison, F. A., Craig, W. W., Christensen, F. E., et al. 2013, *ApJ*, 770, 103

Hauser, M., Möllenhoff, C., Pühlhofer, G., et al. 2004, *AN*, 325, 659

Heinze, A. N., Tonry, J. L., Denneau, L., et al. 2018, *AJ*, 156, 241

Hümmer, S., Rüger, M., Spanier, F., & Winter, W. 2010, *ApJ*, 721, 630

Hunter, J. D. 2007, *CSE*, 9, 90

IceCube Collaboration 2021, *GCN*, 31191, 1

Jones, E., Oliphant, T., Peterson, P., et al. 2001, *SciPy: Open source scientific tools for Python*, http://www.scipy.org/

Kadler, M., Benke, P., Gokus, A., et al. 2021, *ATel*, 15105, 1

Kochanek, C. S., Shappee, B. J., Stanek, K. Z., et al. 2017, *PASP*, 129, 104502

Kuehr, H., Witzel, A., Pauliny-Toth, I. I. K., & Nauber, U. 1981, *A&AS*, 45, 367

Kurahashi, N., Murase, K., & Santander, M. 2022, *ARNPS*, 72, 365

Lindfors, E., Hovatta, T., Pursimo, T., et al. 2021, *ATel*, 15136, 1

Maier, G., & Holder, J. 2017, *ICRC (Busan)*, 35, 747

McBride, F., Fox, D. B., Cowen, D., et al. 2022, *GCN*, 31687, 1

Mészáros, P. 2017, *ARNPS*, 67, 45

Mücke, A., Engel, R., Rachen, J. P., Protheroe, R. J., & Stanev, T. 2000, *CoPhC*, 124, 290

Murase, K., Oikonomou, F., & Petropoulou, M. 2018, *ApJ*, 865, 124

Nasa High Energy Astrophysics Science Archive Research Center (Heasarc) 2014, HEAsoft: Unified Release of FTOOLS and XANADU, Astrophysics Source Code Library, ascl:1408.004

Nilsson, K., Pursimo, T., Villforth, C., et al. 2012, *A&A*, 547, A1

Ochsenbein, F., Bauer, P., & Marcout, J. 2000, *A&AS*, 143, 23

Park, N. 2016, *ICRC*, 34, 771

Parsons, R., & Hinton, J. 2014, *APh*, 56, 26

Petkov, V. B., Novoseltsев, Y. F., Novoseltsевse, R. V. & Baksan Underground Scintillation Telescope Group 2021, *ATel*, 15143, 1

Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2011, *A&A*, 536, A7

Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, *A&A*, 641, A6

Pohl, M., & Schlickeiser, R. 2000, *A&A*, 354, 395

Poole, T. S., Breeveld, A. A., Page, M. J., et al. 2008, *MNRAS*, 383, 627

Poutanen, J., & Stern, B. 2010, *ApJL*, 717, L118

Prince, R., Das, S., Gupta, N., Majumdar, P., & Czerny, B. 2023, *arXiv:2301.06565*

Rolke, W. A., López, A. M., & Conrad, J. 2005, *NIMPA*, 551, 493

Rolke, W. A., López, A. M., & Conrad, J. 2005a, *NIMPA*, 551, 493

Roming, P. W. A., Kennedy, T. E., Mason, K. O., et al. 2005b, *SSRv*, 120, 95

Roming, P. W. A., Koch, T. S., Oates, S. R., et al. 2009, *ApJ*, 690, 163

Sahakyan, N., Giommi, P., Padovani, P., et al. 2022, *MNRAS*, 519, 1396

Santander, M., & Buson, S. 2021, *ATel*, 15102, 1

Schlafly, E. F., & Finkbeiner, D. P. 2011, *ApJ*, 737, 103

Shappee, B. J., Prieto, J. L., Grupe, D., et al. 2014, *ApJ*, 788, 48

Sikora, M., Stawarz, Ł., Moderski, R., Nalewajko, K., & Madejski, G. M. 2009, *ApJ*, 704, 38

Sokolov, A., Marscher, A. P., & McHardy, I. M. 2004, *ApJ*, 613, 725

Stratta, G., Capalbi, M., Giommi, P., et al. 2011, *arXiv:1103.0749*

Tavecchio, F., & Ghisellini, G. 2015, *MNRAS*, 451, 1502

Tonry, J. L., Denneau, L., Heinze, A. N., et al. 2018, *PASP*, 130, 064505

van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, *CSE*, 13, 22

Willingale, R., Starling, R. L. C., Beardmore, A. P., Tanvir, N. R., & O'Brien, P. T. 2013, *MNRAS*, 431, 394

Wood, M., Caputo, R., Charles, E., et al. 2017, *ICRC*, 35, 824

Wright, E. L., Chen, X., Odegard, N., et al. 2009, *ApJS*, 180, 283

Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, *AJ*, 140, 1868

Xie, G. Z., Liang, E. W., Xie, Z. H., & Dai, B. Z. 2002, *AJ*, 123, 2352

Zhirkov, K., Lipunov, V., Gorbovskoy, E., et al. 2021, *ATel*, 15098, 1