# Gamma Rays from Fast Black-hole Winds

M. Ajello, L. Baldini, J. Ballet, G. Barbiellini, D. Bastieri, R. Bellazzini, A. Berretta, E. Bissaldi, R.D. Blandford, E.D. Bloom, et al.

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

## To cite this version:

M. Ajello, L. Baldini, J. Ballet, G. Barbiellini, D. Bastieri, et al., Gamma Rays from Fast Black-hole Winds. The Astrophysical Journal, 2021, 921 (2), pp.144. 10.3847/1538-4357/ac1bb2. hal-03260969

HAL Id: hal-03260969

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

Submitted on 1 Aug 2025

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

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

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

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

# Gamma Rays from Fast Black-hole Winds

M. Ajello¹, L. Baldini², J. Ballet³, G. Barbiellini⁴,⁵, D. Bastieri⁶,⁷, R. Bellazzini⁸, A. Berretta⁹, E. Bissaldi¹⁰,¹¹, R. D. Blandford¹², E. D. Bloom¹², R. Bonino¹³,¹⁴, P. Bruel¹⁵, S. Buson¹⁶, R. A. Cameron¹², D. Caprioli¹⁷, R. Caputo¹⁸, E. Cavazzuti¹⁹, G. Chartas²⁰, S. Chen⁶,²¹, C. C. Cheung²², G. Chiaro²³, D. Costantin²⁴, S. Cutini²⁵, F. D'Ammando²⁶, P. de la Torre Luque¹⁰, F. de Palma²⁷,²⁸, A. Desai²⁹, R. Diesing¹⁷, N. Di Lalla¹², F. Dirirsa³⁰, L. Di Venere¹⁰,¹¹, A. Domínguez³¹, S. J. Fegan¹⁵, A. Franckowiak³², Y. Fukazawa³³, S. Funk³⁴, P. Fusco¹⁰,¹¹, F. Gargano¹¹, D. Gasparrini³⁵,³⁶, N. Giglietto¹⁰,¹¹, F. Giordano¹⁰,¹¹, M. Giroletti²⁶, D. Green³⁷, I. A. Grenier³, S. Guiriec¹⁸,³⁸, D. Hartmann¹, D. Horan¹⁵, G. Jóhannesson³⁹,⁴⁰, C. Karwin¹, M. Kerr²², M. Kovačević²⁵, M. Kuss⁸, S. Larsson⁴¹,⁴²,⁴³, L. Latronico¹³, M. Lemoine-Goumard⁴⁴, J. Li⁴⁵, I. Liodakis⁴⁶, F. Longo⁴,⁵, F. Loparco¹⁰,¹¹, M. N. Lovellette²², P. Lubrano²⁵, S. Maldera¹³, A. Manfreda², S. Marchesi⁴⁷, L. Marcotulli¹, G. Martí-Devesa⁴⁸, M. N. Mazziotta¹¹, I. Mereu⁹,²⁵, P. F. Michelson¹², T. Mizuno⁴⁹, M. E. Monzani¹², A. Morselli³⁵, I. V. Moskalenko¹², M. Negro⁵⁰,⁵¹, N. Omodei¹², M. Orienti²⁶, E. Orlando¹²,⁵², V. Paliya⁵³,⁵⁴, D. Paneque³⁷, Z. Pei⁷, M. Persic⁴,⁵⁵, M. Pesce-Rollins⁸, T. A. Porter¹², G. Principe⁴,⁵,²⁶, J. L. Racusin¹⁸, S. Rainò¹⁰,¹¹, R. Rando⁶,²¹,⁵⁶, B. Rani¹⁸,⁵⁷,⁵⁸, M. Razzano⁸,⁷², A. Reimer¹²,⁴⁸, O. Reimer⁴⁸, P. M. Saz Parkinson⁵⁹,⁶⁰,⁶¹, D. Serini¹⁰, C. Sgrò⁸, E. J. Siskind⁶², G. Spandre⁸, P. Spinelli¹⁰,¹¹, D. J. Suson⁶³, D. Tak¹⁸,⁶⁴, D. F. Torres⁶⁵,⁶⁶,⁶⁷, E. Troja¹⁸,⁶⁸, K. Wood⁶⁹,⁷⁰, G. Zaharijas⁵²,⁷¹, and J. Krake¹

¹ Department of Physics and Astronomy, Clemson University, Kinard Lab of Physics, Clemson, SC 29634-0978, USA; majello@clemson.edu, ckarwin@clemson.edu

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

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

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

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

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

⁷ Dipartimento di Fisica e Astronomia “G. Galilei”, Università di Padova, I-35131 Padova, Italy

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

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

¹⁰ Dipartimento di Fisica “M. Merlin” dell’Università e del Politecnico di Bari, via Amendola 173, I-70126 Bari, Italy

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

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

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

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

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

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

¹⁷ Department of Astronomy and Astrophysics, University of Chicago, Chicago, IL 60637, USA; caprioli@uchicago.edu, rrdiesing@uchicago.edu

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

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

²⁰ Department of Physics and Astronomy of the College of Charleston, Charleston, SC 29424, USA; chartasg@cofc.edu

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

²² Space Science Division, Naval Research Laboratory, Washington, DC 20375-5352, USA

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

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

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

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

²⁷ Dipartimento di Matematica e Fisica “E. De Giorgi”, Università del Salento, Lecce, Italy

²⁸ Istituto Nazionale di Fisica Nucleare, Sezione di Lecce, I-73100 Lecce, Italy

²⁹ Department of Physics, University of Wisconsin-Madison, Madison, WI 53706, USA

³⁰ Laboratoire d'Annecy-le-Vieux de Physique des Particules, Université de Savoie, CNRS/IN2P3, F-74941 Annecy-le-Vieux, France

³¹ Grupo de Altas Energías, Universidad Complutense de Madrid, E-28040 Madrid, Spain

³² Ruhr University Bochum, Faculty of Physics and Astronomy, Astronomical Institute (AIRUB), D-44780 Bochum, Germany

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

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

³⁵ Istituto Nazionale di Fisica Nucleare, Sezione di Roma “Tor Vergata”, I-00133 Roma, Italy

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

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

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

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

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

⁴¹ Department of Physics, KTH Royal Institute of Technology, AlbaNova, SE-106 91 Stockholm, Sweden

⁴² The Oskar Klein Centre for Cosmoparticle Physics, AlbaNova, SE-106 91 Stockholm, Sweden

⁴³ School of Education, Health and Social Studies, Natural Science, Dalarna University, SE-791 88 Falun, Sweden

⁴⁴ Centre d'Études Nucléaires de Bordeaux Gradignan, IN2P3/CNRS, Université Bordeaux 1, BP120, F-33175 Gradignan Cedex, France

⁴⁵ Department of Astronomy, School of Physical Sciences, University of Science and Technology of China, Hefei, Anhui 230026, People's Republic of China

⁴⁶ Finnish Centre for Astronomy with ESO (FINCA), University of Turku, FI-21500 Piikiiö, Finland

⁴⁷ INAF—Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, Via Piero Gobetti, 93/3, I-40129, Bologna, Italy

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

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

⁵⁰ Center for Research and Exploration in Space Science and Technology (CRESST) and NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA
## Abstract

Massive black holes at the centers of galaxies can launch powerful wide-angle winds that, if sustained over time, can unbind the gas from the stellar bulges of galaxies. These winds may be responsible for the observed scaling relation between the masses of the central black holes and the velocity dispersion of stars in galactic bulges. Propagating through the galaxy, the wind should interact with the interstellar medium creating a strong shock, similar to those observed in supernovae explosions, which is able to accelerate charged particles to high energies. In this work we use data from the Fermi Large Area Telescope to search for the γ-ray emission from galaxies with an ultrafast outflow (UFO): a fast ($v \sim 0.1$ c), highly ionized outflow, detected in absorption at hard X-rays in several nearby active galactic nuclei (AGN). Adopting a sensitive stacking analysis we are able to detect the average γ-ray emission from these galaxies and exclude that it is due to processes other than UFOs. Moreover, our analysis shows that the γ-ray luminosity scales with the AGN bolometric luminosity and that these outflows transfer $\sim 0.04\%$ of their mechanical power to γ-rays. Interpreting the observed γ-ray emission as produced by cosmic rays (CRs) accelerated at the shock front, we find that the γ-ray emission may attest to the onset of the wind-host interaction and that these outflows can energize charged particles up to the transition region between galactic and extragalactic CRs.

*Unified Astronomy Thesaurus concepts: Active galactic nuclei (16); Galactic winds (572); Gamma-rays (637)*

## 1. Introduction

Accreting supermassive black holes (SMBHs) at the centers of galaxies, often called active galactic nuclei (AGN), have been observed to launch and power outflows, which can have a dramatic impact on the host galaxies themselves, the inter-galactic medium, and the intracluster medium (Silk & Rees 1998; McNamara & Nulsen 2007; Somerville et al. 2008; Hopkins & Elvis 2010; McCarthy et al. 2010). One spectacular, well-observed, type of outflow are relativistic jets, where particles are accelerated to near the speed of light in narrow collimated beams (often with an opening angle of $\sim 1^{\circ}$), which can extend up to megaparsec scales. These relativistic jets shine at all wavelengths, but are easily studied in radio, X-rays, and γ-rays when the jet axis is not far from our line of sight. Black-hole winds (King & Pounds 2015), on the other hand, are AGN outflows that are not collimated and are generally more difficult to detect, although no less important. Indeed, AGN winds have been proposed as the mechanism able to regulate the coevolution of the galaxy and its central SMBH, which is observed in the scaling of the black-hole mass and the bulge velocity dispersion (Gebhardt et al. 2000; Ferrarese & Ford 2005; Kormendy & Ho 2013). AGN

winds that are powerful enough can heat up and eject the gas from the galaxy, regulating the growth of both the galaxy itself and the black hole.

The most powerful AGN winds can reach velocities of $\sim 0.1-0.3c$ (Chartas et al. 2002; Pounds et al. 2003; Reeves et al. 2003; Tombesi et al. 2010b) and can carry enough energy to unbind the gas of the stellar bulge (King & Pounds 2015). Some of these winds have been identified in nearby AGN through X-ray observations of blueshifted Fe K-shell absorption lines (Reeves et al. 2003; Tombesi et al. 2010b, 2010a, 2012; Gofford et al. 2013).

These winds, which have been dubbed ultrafast outflows (UFOs), are made of highly ionized gas and are likely launched from near the SMBH (King & Pounds 2003). Their wide solid angle [$\Omega/2\pi \approx 0.4$, (Gofford et al. 2015)] and fast velocity allow UFOs to transfer a significant amount of kinetic energy from the AGN to the host galaxy. They are also believed to be common in nearby AGN (King & Pounds 2015).

UFOs, while traveling outward, interact and shock the interstellar medium (ISM; King 2010), producing a reverse shock and a forward shock. The reverse shock decelerates the wind itself while the forward shock travels through the galaxy with a velocity in the $\sim 200-1000$ km s$^{-1}$ range and leads to the formation of a bubble of hot, tenuous gas, see, e.g., Zubovas & King (2012).

<footnote n="51">51 Department of Physics and Center for Space Sciences and Technology, University of Maryland Baltimore County, Baltimore, MD 21250, USA</footnote>

<footnote n="52">52 Istituto Nazionale di Fisica Nucleare, Sezione di Trieste, and Università di Trieste, I-34127 Trieste, Italy</footnote>

<footnote n="53">53 Aryabhatta Research Institute of Observational Sciences (ARIES), Manora Peak, Nainital-263 129, Uttarakhand, India</footnote>

<footnote n="54">54 Deutsches Elektronen Synchrotron DESY, D-15738 Zeuthen, Germany</footnote>

<footnote n="55">55 Osservatorio Astronomico di Trieste, Istituto Nazionale di Astrofisica, I-34143 Trieste, Italy</footnote>

<footnote n="56">56 Center for Space Studies and Activities “G. Colombo”, University of Padova, Via Venezia 15, I-35131 Padova, Italy</footnote>

<footnote n="57">57 Korea Astronomy and Space Science Institute, 776 Daedeokdae-ro, Yuseong-gu, Daejeon 30455, Republic of Korea</footnote>

<footnote n="58">58 Department of Physics, American University, Washington, DC 20016, USA</footnote>

<footnote n="59">59 Santa Cruz Institute for Particle Physics, Department of Physics and Department of Astronomy and Astrophysics, University of California at Santa Cruz, Santa Cruz, CA 95064, USA</footnote>

<footnote n="60">60 Department of Physics, The University of Hong Kong, Pokfulam Road, Hong Kong, People's Republic of China</footnote>

<footnote n="61">61 Laboratory for Space Research, The University of Hong Kong, Hong Kong, People's Republic of China</footnote>

<footnote n="62">62 NYCB Real-Time Computing Inc., Lattingtown, NY 11560-1025, USA</footnote>

<footnote n="63">63 Purdue University Northwest, Hammond, IN 46323, USA</footnote>

<footnote n="64">64 Department of Physics, University of Maryland, College Park, MD 20742, USA</footnote>

<footnote n="65">65 Institute of Space Sciences (ICE, CSIC), Campus UAB, Carrer de Magrans s/n, E-08193 Barcelona, Spain</footnote>

<footnote n="66">66 Institut d'Estudis Espacials de Catalunya (IEEC), E-08034 Barcelona, Spain</footnote>

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

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

<footnote n="69">69 Praxis Inc., Alexandria, VA 22303, USA</footnote>

<footnote n="70">70 Naval Research Laboratory, Washington, DC 20375, USA</footnote>

<footnote n="71">71 Center for Astrophysics and Cosmology, University of Nova Gorica, Nova Gorica, Slovenia</footnote>

<footnote>Received 2021 May 6; revised 2021 July 1; accepted 2021 August 5; published 2021 November 10</footnote>

<footnote n="72">72 Funded by contract FIRB-2021-RBFR12PM1F from the Italian Ministry of Education, University and Research (MIUR).</footnote>
Because of the cooling, the phase and velocity of the outflow should change, eventually leading to the formation of low-velocity molecular outflows, commonly observed in many ultraluminous infrared (IR) galaxies (see e.g., Cicone et al. 2014; Feruglio et al. 2015). Indeed, there are a handful of objects like IRAS 17020+4544 (Longinotti et al. 2018) and Mrk 231 (Feruglio et al. 2015) where both a UFO and molecular outflow have been detected and found in agreement with the prediction of the energy-conserving outflow model, which is the basis of AGN feedback (Fabian 2012).

UFOs have velocities comparable to (or even larger than) those of the ejecta launched in supernova explosions, which are known to shock the ISM and accelerate cosmic rays (CRs). Gamma-ray emission is a signature of the interaction of relativistic charged particles with ambient gas and photon fields and has been observed in many cases in supernova remnants (Acero et al. 2016). Given the similarity, in this work we search for the γ-ray emission from UFOs using the Large Area Telescope (LAT; Atwood et al. 2009) on board the Fermi Gamma-ray Space Telescope (Atwood et al. 2009).

Models of the γ-ray emission from AGN outflows (Wang & Loeb 2016a; Lamastra et al. 2017) show them to be weak emitters, with γ-ray luminosities of ≈10⁴⁰ erg s⁻¹, which explains why UFOs have not yet been detected by the LAT.<fnref n="73" /> Here, we adopt a different strategy and search for the collective γ-ray emission from a sample of UFOs using a stacking technique.

The paper is organized as follows. In Sections 2 and 3, we describe the sample selection and the data analysis. Results are presented in Section 4, with additional tests discussed in Section 5. Section 6 reports the theoretical interpretation of the observed γ-ray emission, while a discussion is given in Section 7. Finally, Section 8 gives our conclusions.

## 2. Sample Selection

We start from a sample of 35 sources that have been identified as UFOs through X-ray observations (Reeves et al. 2003; Tombesi et al. 2010b, 2010a, 2012; Gofford et al. 2013). We have verified that none of the objects are positionally coincident with any known γ-ray sources reported in the Fourth Fermi Large Area Telescope (4FGL) source catalog (Abdollahi et al. 2020). From the initial sample we make the following cuts. First, we only keep the radio-quiet sources (as specified in the original references) to avoid contamination of the signal from the relativistic jet. Furthermore, we only select sources that are nearby ($z < 0.1$) with a mildly relativistic wind velocity ($v > 0.1c$). The former cut is motivated by the expected low luminosity of the UFO emission (Wang & Loeb 2016a), and the latter cut is motivated by the fact that the γ-ray emission is predicted to scale with the kinetic power of the outflow (Wang & Loeb 2016a; Lamastra et al. 2017). After making these cuts we are left with 11 sources, which we use as our benchmark sample. The details of these sources are reported in Table 1.

Table 2 reports additional properties of our sample of UFOs, including the bulge velocity dispersion, 1.4 GHz radio flux and total (8–1000 µm) IR luminosity. Figure 1 shows that the UFOs considered here obey the M–σ relation well (Gültekin et al. 2009; Woo et al. 2010), strengthening the evidence that these outflows operated in the energy-conserving phase in the past (King & Pounds 2015). Finally, the origin of the radio emission

in radio-quiet AGN is not very clear and it is likely due to a number of phenomena, including AGN winds, star formation, free-free emission from photoionized gas, and AGN coronal activity (Panessa et al. 2019). For these reasons, the radio fluxes reported in Table 2 are interpreted as upper limits to the synchrotron emission from accelerated electrons, as discussed in Section 6.

We note that there are alternative models explaining the absorption features as produced not by an outflowing wind, but as resonant absorption by highly ionized iron in the accretion disk (Gallo & Fabian 2011). However, this model has difficulties explaining several of the observed properties of the UFO features like the presence of P Cygni profiles (Nardini et al. 2015; Chartas et al. 2016), or the correlation between outflow velocity and the AGN bolometric luminosity (Saez & Chartas 2011; Matzeu et al. 2017).

## 3. Data Analysis

### 3.1. Data

We analyze data collected by Fermi-LAT between 2008 August 4 and 2019 September 10 (11.1 yr). The events have energies in the range of 1–800 GeV and are binned in eight bins per decade. The pixel size is 0.08°. To reduce contamination from the Earth’s limb, we use a maximum zenith angle of 105°. We define a 10° × 10° region of interest (ROI) centered at the position of each UFO source. We use the standard data filters: DATA QUAL > 0 and LAT_CONFIG==1. The analysis is performed using Fermipy (v0.18.0),<fnref n="74" /> which utilizes the underlying Fermitools (v1.2.23).

We select photons corresponding to the P8R3_SOURCE_V2 class (Atwood et al. 2013). In order to optimize the sensitivity of our stacking technique we implement a joint likelihood analysis with the four point-spread function (PSF) event types available in the Pass 8 data set.<fnref n="75" /> The data is divided into quartiles corresponding to the quality of the reconstructed direction, from the lowest quality quartile (PSF0) to the best quality quartile (PSF3). Each sub-selection has its own binned likelihood instance that is combined in a global likelihood function for the ROI. This is easily implemented in Fermipy by specifying the components section in the configuration file. Each PSF type also has its own corresponding isotropic spectrum, namely, iso_P8R3_SOURCE_V2_PSFᵢ_v1, for i ranging from 0–3. The Galactic diffuse emission is modeled using the standard component (gll_iem_v07), and the point-source emission is modeled using the 4FGL catalog (gll_psc_v20). In order to account for photon leakage from sources outside of the ROI due to the PSF of the detector, the model includes all 4FGL sources within a 15° × 15° region. The energy dispersion correction (edisp_bins=-1) is enabled for all sources except the isotropic component.

### 3.2. Analysis

In the Local Universe ($z < 0.1$) UFOs are predicted to have a γ-ray luminosity of ~10⁴⁰ erg s⁻¹ (Wang & Loeb 2016a), making them too faint to be detected individually by Fermi-LAT. Indeed, adopting the average photon index in the 4FGL

<footnote n="73">⁷³ No γ-ray source from the 4FGL catalog (Abdollahi et al. 2020) is associated with a UFO.</footnote>

<footnote n="74">⁷⁴ Available at https://fermipy.readthedocs.io/en/latest/.</footnote>

<footnote n="75">⁷⁵ For more information on the different PSF types see https://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Cicerone/Cicerone_Data/LAT_DP.html.</footnote>
<table><caption>**Table 1**
UFO Source Sample</caption><thead><tr><th>Name<br>(1)</th><th>R.A. (deg)<br>[J2000]<br>(2)</th><th>Decl. (deg)<br>[J2000]<br>(3)</th><th>Type</th><th>Redshift<br>[z]<br>(5)</th><th>Velocity<br>[v/c]<br>(6)</th><th>logM<sub>BH</sub><br>[M<sub>☉</sub>]<br>(7)</th><th>logE&#x0302;<sub>K</sub><sup>Min</sup><br>[erg s<sup>-1</sup>]<br>(8)</th><th>logE&#x0302;<sub>K</sub><sup>Max</sup><br>[erg s<sup>-1</sup>]<br>(9)</th><th>logL<sub>Bol</sub><br>[erg s<sup>-1</sup>]<br>(10)</th><th>95% UL (×10<sup>-11</sup>)<br>[ph cm<sup>-2</sup> s<sup>-1</sup>]<br>(11)</th></tr></thead><tbody><tr><td>Ark 120<sup>a,c</sup></td><td>79.05</td><td>-0.15</td><td>Sy1</td><td>0.033</td><td>0.27</td><td>8.2 ± 0.1</td><td>>43.1</td><td>46.2 ± 1.3</td><td>45.0<sup>f</sup><br>44.2<sup>h</sup><br><strong>44.6</strong></td><td>7.5</td></tr><tr><td>MCG-5-23-16<sup>a,c</sup></td><td>146.92</td><td>-30.95</td><td>Sy2</td><td>0.0084</td><td>0.12</td><td>7.6 ± 1.0</td><td>42.7 ± 1.0</td><td>44.3 ± 0.2</td><td>44.1<sup>l</sup></td><td>4.3</td></tr><tr><td>NGC 4151<sup>a,c</sup></td><td>182.64</td><td>39.41</td><td>Sy1</td><td>0.0033</td><td>0.105</td><td>7.1 ± 0.2</td><td>>41.9</td><td>43.1 ± 0.5</td><td>44.1<sup>g</sup><br>42.9<sup>h</sup><br>43.9<sup>i</sup><br>42.9<sup>j</sup><br>43.2<sup>k</sup></td><td>10.6</td></tr><tr><td>PG 1211+143<sup>a,c</sup></td><td>183.57</td><td>14.05</td><td>Sy1</td><td>0.081</td><td>0.13</td><td>8.2 ± 0.2</td><td>43.7 ± 0.2</td><td>46.9 ± 0.1</td><td>45.7<sup>f</sup><br>44.8<sup>h</sup><br>44.7<sup>j</sup><br>45.0<sup>k</sup><br><strong>45.1</strong></td><td>3.7</td></tr><tr><td>NGC 4507<sup>a,c</sup></td><td>188.90</td><td>-39.91</td><td>Sy2</td><td>0.012</td><td>0.18</td><td>6.4 ± 0.5</td><td>>41.2</td><td>44.6 ± 1.1</td><td>44.3<sup>e</sup></td><td>3.4</td></tr><tr><td>NGC 5506<sup>b,d</sup></td><td>213.31</td><td>-3.21</td><td>Sy1.9</td><td>0.006</td><td>0.25</td><td>7.3 ± 0.7</td><td>43.3 ± 0.1</td><td>44.7 ± 0.5</td><td>44.3<sup>e</sup></td><td>6.4</td></tr><tr><td>Mrk 290<sup>a,c</sup></td><td>233.97</td><td>57.90</td><td>Sy1</td><td>0.030</td><td>0.14</td><td>7.7 ± 0.5</td><td>43.4 ± 0.9</td><td>45.3 ± 1.2</td><td>44.4<sup>e</sup></td><td>4.5</td></tr><tr><td>Mrk 509<sup>a,c</sup></td><td>311.04</td><td>-10.72</td><td>Sy1</td><td>0.034</td><td>0.17</td><td>8.1 ± 0.1</td><td>>43.2</td><td>45.2 ± 1.0</td><td>45.2<sup>e</sup><br>44.3<sup>h</sup><br>45.3<sup>i</sup><br>44.3<sup>j</sup><br>44.5<sup>k</sup></td><td>9.5</td></tr><tr><td>SWIFT J2127.4<br>+5654<sup>b,d</sup></td><td>321.94</td><td>56.94</td><td>Sy1</td><td>0.014</td><td>0.23</td><td>~7.2</td><td>42.8 ± 0.1</td><td>45.6 ± 0.5</td><td>44.7<sup>d</sup></td><td>9.1</td></tr><tr><td>MR 2251-178<sup>b,d</sup></td><td>343.52</td><td>-17.58</td><td>Sy1</td><td>0.064</td><td>0.14</td><td>8.7 ± 0.1</td><td>43.3 ± 0.1</td><td>46.7 ± 0.7</td><td>45.8<sup>f</sup></td><td>7.4</td></tr><tr><td>NGC 7582<sup>a,c</sup></td><td>349.60</td><td>-42.37</td><td>Sy2</td><td>0.0052</td><td>0.26</td><td>7.1 ± 1.0</td><td>43.4 ± 1.1</td><td>44.9 ± 0.4</td><td>43.3<sup>e</sup></td><td>4.7</td></tr></tbody></table>

**Notes.** Our sample comprises 11 sources with $z < 0.1$ and $v > 0.1c$. The first superscript on the source name indicates the reference for the detection, and the second superscript indicates the reference for the UFO parameters (Columns 6–9), where $\dot{E}_K^{\text{min}}$ and $\dot{E}_K^{\text{max}}$ are the minimum and maximum kinetic powers. Values for the bolometric luminosity ($L_{\text{Bol}}$) are taken from the literature, with the reference indicated by the superscript. For sources with numerous determinations we also give the mean value in boldface text. The $\gamma$-ray flux (1–800 GeV) upper limit (UL) is calculated at the 95% confidence level, using a photon index of $-2.0$.

a Tombesi et al. (2010a).

b Gofford et al. (2013).

c Tombesi et al. (2012).

d Gofford et al. (2015).

e Vasudevan et al. (2010).

f Vasudevan & Fabian (2007).

g Vasudevan & Fabian (2009).

h Peterson et al. (2004).

i Crenshaw & Kraemer (2012).

j (Kaspi et al. 2005, 5100 Å flux density).

k (Kaspi et al. 2005, 1450 Å flux density).

l Alonso-Herrero et al. (2011).

catalog of $\Gamma = -2.2$ we derive $a > 1$ GeV flux of $3.3 \times 10^{-12}$ ph cm$^{-2}$ s$^{-1}$, for a source with a luminosity of $10^{40}$ erg s$^{-1}$ at $z = 0.014$ (the median redshift of our sample). This flux is $\sim 2.5$ times fainter than the weakest source reported in the 4FGL catalog. We therefore analyze our source sample using a stacking technique. This technique has been developed previously and has been successfully employed for multiple studies, i.e., upper limits on dark matter interactions (Ackermann et al. 2011), detection of the extragalactic background light (Abdollahi et al. 2018), extreme blazars (Paliya et al. 2019), and star-forming galaxies (Ajello et al. 2020a).

The main assumption that we make for the stacking technique is that the sample of UFOs we are considering can be characterized by average quantities like the average flux and

the average photon index (when we model their spectra with a power law). There are then two steps to the method. In the first step, the model components are optimized for each ROI using a maximum likelihood fit. We evaluate the significance of each source in the ROI using the test statistic (TS), which is defined as

$$ \text{TS} = -2 \log(L_0/L), \quad (1) $$

where $L_0$ is the likelihood for the null hypothesis, and $L$ is the likelihood for the alternative hypothesis.<fnref n="76" /> For the first iteration

<footnote n="76">⁷⁶ For a more complete explanation of the TS resulting from a likelihood fit see Mattox et al. (1996) and https://fermi.gsfc.nasa.gov/ssc/data/analysis/documen tation/Cicerone/Cicerone_Likelihood/.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Bulge stellar velocity dispersion vs. black-hole mass for our UFO sample, with values taken from the literature. Measurements were found for 8/11 sources. The error bars are statistical plus systematic, where the systematic uncertainty comes from different independent estimates. Information for the velocity dispersion measurements is provided in Table 2. To quantify the systematic uncertainty in the black-hole mass, we use minimum and maximum values from the different references provided in Table 2, as well as the values given in Table 1. The solid and dashed lines show the scaling relations for active and quiescent galaxies, from Woo et al. (2010) and Gültekin et al. (2009), respectively.</figcaption></figure>

<table><caption>**Table 2**
Additional UFO Properties</caption><thead><tr><th rowspan="2">Name<br>(1)</th><th>Velocity<br>Dispersion<br>[km s<sup>-1</sup>]</th><th colspan="2">1.4 GHz Radio</th></tr><tr><th>(2)</th><th>Flux<sup>h</sup><br>[mJy]</th><th>IR Lum.<sup>i</sup><br>[log(L<sub>☉</sub>)]</th></tr></thead><tbody><tr><td>Ark 120</td><td>184, 238<sup>b,c</sup></td><td>12.4</td><td>11.0</td></tr><tr><td>MCG-5-23-16</td><td>152, 192<sup>b,d</sup></td><td>14.3</td><td>9.6</td></tr><tr><td>NGC 4151</td><td>94, 119<sup>b,d</sup></td><td>347.6</td><td>10.2</td></tr><tr><td>PG 1211+143</td><td>...</td><td>4.3</td><td>...</td></tr><tr><td>NGC 4507<sup>a</sup></td><td>146, 156<sup>c</sup></td><td>67.4</td><td>10.5</td></tr><tr><td>NGC 5506</td><td>160, 200<sup>c</sup></td><td>355</td><td>10.5</td></tr><tr><td>Mrk 290</td><td>109, 111<sup>f</sup></td><td>5.32</td><td>&lt;10.3</td></tr><tr><td>Mrk 509</td><td>172, 196<sup>c</sup></td><td>19.2</td><td>10.5</td></tr><tr><td>SWIFT J2127.4<br>+5654</td><td>...</td><td>6.4</td><td>10.4</td></tr><tr><td>MR 2251-178</td><td>...</td><td>16</td><td>&lt;10.5</td></tr><tr><td>NGC 7582</td><td>110, 116<sup>c</sup></td><td>270</td><td>10.6</td></tr></tbody></table>

**Notes.** The second column gives velocity dispersion measurements taken from the literature, with the references indicated by the superscripts. Measurements were found for 8/11 sources, and we provide minimum and maximum values (separated by a comma). For sources with just one reference, the range is due to statistical error only, and for sources with two references, the range also includes the systematic error due to the different estimates.

<sup>a</sup> Note that most published estimates of the black-hole mass for NGC 4507 are based on velocity dispersion and [O III] line widths, and thus they are not independent measures. In quantifying the uncertainty in Figure 1, we also use black-hole mass values from Bian & Gu (2007); Beifiori et al. (2012); Nicastro et al. (2003).

<sup>b</sup> Woo et al. (2010).

<sup>c</sup> Grier et al. (2013).

<sup>d</sup> Onken et al. (2014).

<sup>e</sup> Marinucci et al. (2012).

<sup>f</sup> Bennert et al. (2015).

<sup>g</sup> Hyperleda.

<sup>h</sup> NVSS (Condon et al. 1998).

<sup>i</sup> IRAS (Kleinmann et al. 1986; Moshir et al. 1990).

of the fit, the spectral parameters of the Galactic diffuse component (index and normalization) and the isotropic component are freed. In addition, we free the normalizations of all 4FGL sources with TS ≥ 25 that are within 5° of the ROI center, as well as sources with TS ≥ 500 and within 7°. Lastly, the UFO source is fit with a power-law spectral model, and the spectral parameters (normalization and index) are also freed. In the first step, we also find new point sources using the Fermi function *find.sources*, which generates TS maps and identifies new sources based on peaks in the TS. The TS maps are generated using a power-law spectral model with an index of -2.0. The minimum separation between two point sources is set to 0.5°, and the minimum TS for including a source in the model is set to 16.

In the second step, 2D TS profiles are generated for the spectral parameters of each UFO source, where the TS is defined as in Equation (1). We scan photon indices from -1 to -3.3 with a spacing of 0.1 and total integrated photon flux (between 1 and 800 GeV) from $10^{-13}-10^{-9}$ ph cm$^{-2}$ s$^{-1}$ with 40 logarithmically spaced bins, freeing just the parameters of the diffuse components. For this step, the power-law spectra of the UFOs are defined in terms of the total flux ($F_{tot}$), integrated between the minimum energy ($E_{min}$) and the maximum energy ($E_{max}$):

$$ \frac{dN}{dE} = \frac{F_{tot}(\Gamma + 1)E^{\Gamma}}{E_{\max}^{\Gamma+1} - E_{\min}^{\Gamma+1}} \quad (2) $$

Note that the likelihood value for the null hypothesis is calculated at the end of the first step by removing the UFO source from the model. Since we perform a joint likelihood in the different PSF event types (PSF0-PSF3), the total profile for each source is obtained by adding the profiles from each of the four event types. Lastly, the TS profiles for all sources are added to obtain the stacked profile. The TS is an additive quantity, and so the stacked profile gives the statistical significance for the combined signal.

We validated the stacking method relying on a set of Monte Carlo simulations that reproduce the Fermi-LAT observations. In these tests, the simulations include the isotropic and Galactic emission, as well as an isotropic population of point sources resembling blazars, which account for the vast majority of sources detected by Fermi-LAT. Faint, below-threshold blazars are included in the synthetic sky following the models of Ajello et al. (2015). Using this setup, two different tests were performed. The stacking analysis was performed at 60 random empty positions, i.e., positions away from bright detected sources. This analysis yielded no detection, confirming that the technique does not generate spurious detections. The second set of tests was aimed at characterizing the detected signal. The stacking was performed for 60 simulated sources whose flux was extracted from a power-law distribution with index -2.5 and minimum and median flux of, respectively, $4 \times 10^{-10}$ and $6.4 \times 10^{-10}$ ph cm$^{-2}$ s$^{-1}$. The photon indices were extracted from a Gaussian distribution with average -2.21 and dispersion of 0.2. The values derived from the stacking analysis (flux = $7.0_{-0.7}^{+0.6} \times 10^{-10}$ ph cm$^{-2}$ s$^{-1}$ and index of -2.24 ± 0.05) are in agreement with the inputs, showing that our analysis successfully retrieves the average quantities of a population of sources. Moreover, the likelihood profile would
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Stacked TS profile for the sample of UFOs. The color scale indicates the TS, and the plus sign indicates the location of the maximum value, with a TS = 30.1 (5.1σ). Significance contours (for 2 degrees of freedom (dof)) are overlaid on the plot showing the 68%, 90%, and 99% confidence levels, corresponding to ΔTS = 2.30, 4.61, and 9.21, respectively.</figcaption></figure>

not show a significant peak if those average quantities were not representative of the population.

# 4. Results

## 4.1. Stacked TS Profile for the Benchmark Sample

The log-likelihoods (i.e., logL) are maximized with the optimizer MINUIT (James & Roos 1975), and we have verified that each fit converges properly, as indicated by the MINUIT outputs of quality = 3 and status = 0. The 95% flux upper limits from the preprocessing step are reported in Table 1.

The stacked profile for our UFO sample is shown in Figure 2. The maximum TS is 30.1 (5.1σ)<fnref n="77" /> corresponding to a best-fit index of −2.1 ± 0.3 and a best-fit photon flux (1–800 GeV) of 2.5<sup>+1.5</sup><sub>-0.9</sub> × 10<sup>-11</sup> ph cm<sup>-2</sup> s<sup>-1</sup>. The 68%, 90%, and 99% significance contours are overlaid on the map, and as can be seen the spectral parameters are well constrained. The source with the overall highest individual TS is NGC 4151, having a maximum value of 21.2 (4.2σ), corresponding to a best-fit index of −1.9<sup>+0.5</sup><sub>-0.3</sub> and a best-fit flux of 6.3<sup>+3.7</sup><sub>-3.8</sub> × 10<sup>-11</sup> ph cm<sup>-2</sup> s<sup>-1</sup>. The stacking analysis excluding NGC 4151 yields a maximum TS of 15.1 (3.5σ), corresponding to a best-fit index of −2.2 ± 0.4 and a best-fit flux of 2.0<sup>+2.0</sup><sub>-1.0</sub> × 10<sup>-11</sup> ph cm<sup>-2</sup> s<sup>-1</sup>.

## 4.2. Spectral Energy Distribution (SED) of UFOs

The best-fit SED for our UFO sample is shown in Figure 3. The butterfly plot is constructed by sampling the range of parameter values that are within the 68% confidence contour of the stacked profile. In addition, we calculate the SED flux in three logarithmically spaced bins between 1 and 800 GeV. In

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Best-fit UFO SED (black solid line) with 1σ uncertainty envelope (gray band). The tan data points show the UFO energy flux calculated in four different energy bins. The dashed cyan line shows our hadronic model (see Section 6), corresponding to an outflow that has propagated to ~20 pc. The effective redshift z = 0.013 was used to convert the γ-ray flux into luminosity.</figcaption></figure>

every bin, we fix the power-law index of the UFOs to −2.0 and leave all other parameters free to vary. As can be seen, these data points are in agreement with the best-fit SED model. To characterize the UFO spectrum at low energy we repeat the stacking analysis in the energy range 0.1–1 GeV, which yields a 95% flux upper limit (Δlog L = 2.71/2) of 5.7 × 10<sup>-10</sup> ph cm<sup>-2</sup> s<sup>-1</sup>. We also overlay our best-fit hadronic model presented in Section 6.

## 4.3. Bins of Bolometric Luminosity and Kinetic Power

We test whether the γ-ray emission from UFOs scales with AGN bolometric luminosity and outflow kinetic power. To properly take the distance of each source into account, we stack in the luminosity-index space. We take estimates of the bolometric luminosity from the literature, as reported in Table 1. Such estimates can be obtained by applying a correction factor to a certain flux, typically the 5100 Å optical emission, the 1450 Å UV emission, or the 2–10 keV X-ray emission. Alternatively, the bolometric luminosity can be determined by fitting an SED to the broadband emission. In any case, the absorption from the host galaxy must be corrected for, which has a large dependence on the viewing angle of the source, and can introduce a rather significant uncertainty. In addition, the contribution from the host galaxy emission also needs to be corrected for (i.e., UV/IR/optical emission from the galactic disk). Most of the AGN emission is observed in the optical/UV, while <10% is emitted in the X-ray, and thus a broadband SED fitting ensures a more accurate determination of the bolometric luminosity. We therefore search the literature for the most reliable estimates of the bolometric luminosity, and rely on the X-ray determination for only two sources (MCG-5-23-16 and SWIFT J2127.4+5654) for which no other estimates could be found. For sources with multiple estimates we take the geometric mean. The mean of the bolometric luminosity of our sample is 2.5 × 10<sup>44</sup> erg s<sup>-1</sup>, and we create two bins around this value.

The stacked profiles for the two bins are shown in Figure 4. The first bin has five sources, with a mean redshift of 0.007. The maximum TS is 28.5 (5.0σ), corresponding to a best-fit index of −1.9<sup>+0.3</sup><sub>-0.4</sub> and a best-fit luminosity of 1.6<sup>+0.9</sup><sub>-0.8</sub> × 10<sup>40</sup> erg s<sup>-1</sup>. The second bin has six sources, with a mean redshift of 0.04. The maximum TS is 9.9 (2.7σ), corresponding to a best-fit index of

<footnote n="77">⁷⁷ The conversion from TS to σ has been performed on the assumption that the TS behaves asymptotically as a χ² distribution with 2 dof (Mattox et al. 1996). Additionally, the Akaike information criterion test also shows the null hypothesis to be highly disfavored with a relative likelihood of 2 × 10<sup>-6</sup>.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Stacked profiles for bins of bolometric luminosity (the mean kinetic power bins are also the same). The left and right panels show the stacking for sources with bolometric luminosity (or kinetic power) below and above the average, respectively. The color scale indicates the TS and is set to the maximum value for each bin. The black plus sign gives the best-fit parameters. The first bin consists of five sources, with a maximum TS of 28.5 (5.0σ); and the second bin consists of six sources, with a maximum TS of 9.9 (2.7σ).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Stacked profiles for bolometric efficiency (left) and kinetic power efficiency (right). The color scale indicates the TS and is set to the maximum value. The black plus sign gives the best-fit parameters. Significance contours (for 2 dof) are overlaid on the plot showing the 68%, 90%, and 99% confidence levels, corresponding to $\Delta$TS = 2.30, 4.61, and 9.21, respectively.</figcaption></figure>

$-2.4_{-0.5}^{+0.6}$ and a best-fit luminosity of $2.5_{-1.5}^{+1.5} \times 10^{41}$ erg s$^{-1}$. The total TS (bin 1 + bin 2) for the stacking in bins is 38.4, compared to 30.1 for the full stack.

We also stack the γ-ray luminosity in bins of kinetic power. In general the kinetic power as determined from X-ray observations has a large uncertainty, as can be seen in Table 1. Minimum and maximum values are typically reported, corresponding to minimum and maximum radii of the outflow. We use the geometric mean of the minimum and maximum estimates for our calculations (also incorporating statistical

uncertainties in the range). We create two bins around the mean kinetic power, which has a value of $1.8 \times 10^{44}$ erg s$^{-1}$. The stacked profiles for the two bins turn out to be the same as those of the bolometric bins, as shown in Figure 4.

To further verify the relations found above for the stacking in bins, we perform the stacking analysis using both bolometric efficiency ($\epsilon_{\text{Bol}} = L_{\gamma}/L_{\text{Bol}}$) and kinetic power efficiency ($\epsilon_{\text{E}_K} = L_{\gamma}/L_{\text{E}_K}$). This is done by evaluating for each source the TS of a given $\epsilon_{\text{Bol}}$ (or $\epsilon_{\text{E}_K}$) and using that efficiency value, the bolometric luminosity (or kinetic power), and the distance
<figure><img src="image_9.png" /><figcaption>**Figure 6.** γ-ray luminosity vs. bolometric luminosity (left) and kinetic power (right). The black data points result from stacking in γ-ray luminosity, and the uncertainty in the x-axis corresponds to the bin widths. The gray dashed-dotted vertical lines show the value used to divide the bins. The solid green line shows the best fit resulting from stacking in efficiency, with the green band showing the 1σ confidence level. For reference, the blue lines show a range of efficiencies within roughly an order of magnitude of the best fit. The orange bar in both plots shows the average one-sided uncertainty in individual measurements of AGN bolometric luminosity (left) and kinetic power (right). In the left panel we also overlay the predicted efficiency derived from (Liu et al. 2018, dashed purple line). See the text for more details.</figcaption></figure>

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

of the source to transform to γ-ray flux (for a given photon index). Results for these fits are shown in Figure 5. The left panel shows the bolometric efficiency, with a best-fit value of $3.2^{+1.6}_{-1.5} \times 10^{-4}$, corresponding to a best-fit index of $-1.9^{+0.3}_{-0.4}$, and a maximum TS of 28.2 (5σ). The right panel of Figure 5 shows the kinetic power efficiency, with a best-fit value of $4.0^{+2.3}_{-2.0} \times 10^{-4}$, corresponding to a best-fit index of $-1.8^{+0.3}_{-0.4}$, and a maximum TS of 23.0 (4.4σ). We note that the best-fit index from the efficiency analysis is slightly harder than the one found by the flux-index stacking, but compatible within 1σ uncertainties. The small shift observed in the best-fit index value is due to how the TS profiles are weighted differently when stacking in efficiency with respect to flux.

The result for stacking in bolometric luminosity and kinetic power are summarized in Figure 6. The left panel shows the γ-ray luminosity versus bolometric luminosity, and the right panel shows the γ-ray luminosity versus UFO kinetic power. The black data points are for stacking in bins, and the corresponding best-fit efficiency, along with the 1σ confidence interval, is plotted with the green band. Also plotted are lines for different efficiencies under the assumption of a linear scaling. As can be seen, the results on the efficiencies are in very good agreement with the stacking in bins.

In the left panel of Figure 6 we also overlay the predicted scaling of $L_γ$ with $L_{Bol}$ from Liu et al. (2018).<fnref n="78" /> As can be seen, Liu et al. (2018) predict a nearly linear scaling between the logarithms of the two luminosities (over their $L_{Bol}$(erg s⁻¹) = $10^{42} – 10^{45}$ range) with an efficiency of ~8 × $10^{-4}$, which is in reasonably good agreement with the one measured here.

### 4.4. Representative Luminosity of the Sample

Because the 11 UFO galaxies are detected at fairly different distances, we adopt a weighting scheme to compute the representative luminosity of the sample. In this framework

$\overline{L}_γ = \frac{\sum_{i=1}^{11} L_{\gamma,i} \times TS_i}{TS_{tot}}$, where $L_{γ,i}$ and $TS_i$ are the luminosity and the TS for the ith galaxy at the global best-fit position (1–800 GeV flux of $2.5 \times 10^{-11}$ ph cm⁻² s⁻¹ and photon index of −2.1) and $TS_{tot} = 30.1$. The representative luminosity is found to be $\overline{L}_γ = 7.9^{+5.1}_{-2.9} \times 10^{40}$ erg s⁻¹ and would correspond to an effective redshift of $z = 0.013$ (adopting the above best-fit parameters). This luminosity is in very good agreement with the one obtained scaling the average bolometric luminosity $L_{Bol} = 2.5 \times 10^{44}$ erg s⁻¹ by the best-fit efficiency ($\epsilon_{Bol} = 3.2 \times 10^{-4}$). The effective redshift is also very close to the median redshift of the sample ($z = 0.013$ versus $z = 0.014$) making the TS-weighted luminosity compatible with the median γ-ray luminosity of the sample.

### 4.5. Simulations

The results presented here are validated using Monte Carlo simulations. We simulate the fields of the 11 UFOs considering the Galactic and isotropic emission (modeled as gll_iem_v07 and iso_P8R3_SOURCE_V2_v1, respectively), background sources from the 4FGL catalog, and our test source at the position of the UFO in each ROI. The UFO spectral parameters are set to be the same as the best-fit values from the data, i.e., $\text{index} = -2.1$ and flux = $2.5 \times 10^{-11}$ ph cm⁻² s⁻¹. For simplicity we use the standard event type (evtype = 3), i.e., we do not use the four different PSF event types. The data is simulated using the simulate_roi function from Fermipy. The simulation is created by generating an array of Poisson random numbers, where the expectation values are drawn from the model cube.<fnref n="79" /> Finally, we run our stacking pipeline on the simulated data. We recover the input values, with a best-fit index of $-2.2^{+0.4}_{-0.2}$, a best-fit flux of $3.2^{+1.8}_{-1.6} \times 10^{-11}$ ph cm⁻² s⁻¹, and a maximum TS of 21.2 (4.2σ). The stacked profile is shown in Figure 7. Overall, the results from the simulation are consistent with the real data.

<footnote n="78">⁷⁸ Our derivation is made converting the peak 1 GeV luminosities (reported in their Figure 5) to the 1–800 GeV energy range using the best-fit spectral index of −2.1.</footnote>

<footnote n="79">⁷⁹ More information on generating the simulations is available at https://fermipy.readthedocs.io/en/latest/fermipy.html.</footnote>
<figure><img src="image_11.png" /><figcaption>**Figure 7.** Stacked profile for our simulation run, in which the UFO sources are simulated with an index of −2.1 and a flux of $2.5 \times 10^{-11} \text{ ph cm}^{-2} \text{ s}^{-1}$. The color scale indicates the TS, and the plus sign indicates the location of the maximum value, with a TS = 21.2 (4.2σ). Significance contours (for 2 dof) are overlaid on the plot showing the 68%, 90%, and 99% confidence levels, corresponding to ΔTS = 2.30, 4.61, and 9.21, respectively. The maximum TS of the color scale is set to 30.1 (the maximum value from Figure 2).</figcaption></figure>

# 5. Additional Tests

## 5.1. Control Sample

We repeat the analysis with a sample of 20 low redshift ($z < 0.1$) radio-quiet AGN that do not have UFOs. The sources were selected from the samples of Tombesi et al. (2010a) and Igo et al. (2020) for which no UFO was found. The sample of Tombesi et al. (2010a) is based on absorption features, while the sample of Igo et al. (2020) uses the excess variance method. Of the 20 sources in our control sample, there are 10 sources in common between the two studies, four additional sources from Tombesi et al. (2010a), and six additional sources from Igo et al. (2020). For reference, the list of sources in the control sample is given in Table 3. Figure 8 shows that the benchmark and control samples are well matched in X-ray luminosity and redshift.

Results for the stacked profile are shown in Figure 9. No signal is detected, with a maximum TS of 1.1. Using the profile likelihood method and a photon index of −2.0, the upper limit on the flux (1–800 GeV) at the 95% confidence level is $8.8 \times 10^{-12} \text{ ph cm}^{-2} \text{ s}^{-1}$. This supports the interpretation of the γ-ray emission being due to the outflow rather than other processes in AGN.

## 5.2. Alternative UFO Samples

The fractional excess variance method was recently used in Igo et al. (2020) to search for UFOs in the samples of Tombesi et al. (2010a) and Kara et al. (2016). Overall, the results are in agreement with the past literature, finding that UFOs are a relatively widely observed phenomena in nearby AGN. However, there are differences with respect to previous studies in regards to which sources are classified as UFOs, and the corresponding UFO parameters.

<table><caption>**Table 3**
Control Sample</caption><thead><tr><th>Name<br>(1)</th><th>R.A.<br>(2)</th><th>Decl.<br>(3)</th><th>Redshift<br>(4)</th><th>IR Lumin.<br>[log (L<sub>☉</sub>)]<br>(5)</th><th>1.4 GHz flux<br>[mJy]<br>(6)</th></tr></thead><tbody><tr><td>ESO<br>198-G024</td><td>39.58</td><td>-52.19</td><td>0.046</td><td>...</td><td>...</td></tr><tr><td>Fairall 9</td><td>20.94</td><td>-58.81</td><td>0.047</td><td>...</td><td>...</td></tr><tr><td>H 0557-385</td><td>89.51</td><td>-38.33</td><td>0.034</td><td>...</td><td>...</td></tr><tr><td>MCG+8-<br>11-11</td><td>88.72</td><td>46.44</td><td>0.020</td><td>11.1</td><td>286</td></tr><tr><td>Mrk 590</td><td>33.64</td><td>-0.77</td><td>0.026</td><td>...</td><td>...</td></tr><tr><td>Mrk 704</td><td>139.61</td><td>16.31</td><td>0.029</td><td>...</td><td>...</td></tr><tr><td>NGC 526A</td><td>20.98</td><td>-35.07</td><td>0.019</td><td>10.5</td><td>13.9</td></tr><tr><td>NGC 5548</td><td>214.50</td><td>25.14</td><td>0.017</td><td>...</td><td>...</td></tr><tr><td>NGC 7172</td><td>330.51</td><td>-31.87</td><td>0.0090</td><td>10.4</td><td>37.6</td></tr><tr><td>NGC 7469</td><td>345.82</td><td>8.874</td><td>0.016</td><td>11.6</td><td>181</td></tr><tr><td>ESO<br>113-G010</td><td>16.32</td><td>-58.44</td><td>0.027</td><td>...</td><td>...</td></tr><tr><td>ESO<br>362-G18</td><td>79.90</td><td>-32.66</td><td>0.012</td><td>...</td><td>...</td></tr><tr><td>IRAS<br>17020<br>+4544</td><td>255.88</td><td>45.68</td><td>0.060</td><td>11.6</td><td>129</td></tr><tr><td>MS22549-<br>3712</td><td>344.41</td><td>-36.94</td><td>0.039</td><td>...</td><td>...</td></tr><tr><td>NGC 1365</td><td>53.40</td><td>-36.14</td><td>0.0055</td><td>10.9</td><td>534</td></tr><tr><td>NGC 4748</td><td>193.05</td><td>-13.41</td><td>0.015</td><td>10.4</td><td>14.3</td></tr><tr><td>Mrk 110</td><td>141.30</td><td>52.29</td><td>0.035</td><td>...</td><td>...</td></tr><tr><td>IRAS<br>05078<br>+1626</td><td>77.69</td><td>16.50</td><td>0.018</td><td>10.8</td><td>6.3</td></tr><tr><td>ESO<br>511-G30</td><td>214.84</td><td>-26.64</td><td>0.022</td><td>...</td><td>...</td></tr><tr><td>NGC 2110</td><td>88.05</td><td>-7.46</td><td>0.0078</td><td>10.3</td><td>300</td></tr></tbody></table>

**Note.** See Tombesi et al. (2010a) and Igo et al. (2020) for further details of the sources. The IR luminosity is reported in the 8–1000 µm range and derived from IRAS (Kleinmann et al. 1986; Moshir et al. 1990). The radio fluxes are derived from NVSS (Condon et al. 1998).

As the authors mention in Igo et al. (2020), their method relies on the variability of the strength of the emission (or absorption) features and is less sensitive in detecting cases where these features may vary in energy. The excess variance method is well suited for detecting UFOs in objects that show small changes in the energy of the UFO, but large changes of the equivalent width for the same energy. This is one reason why the excess variance method can potentially miss objects that were detected in spectral-timing analyses that model individual spectra in single epochs.

As an additional a posteriori test we perform our stacking analysis with the UFO sample determined in Igo et al. (2020), relying on sources classified as either likely outflows or possible outflows therein. Additionally, we use the same selection criterion as for our benchmark sample, i.e., $z < 0.1$ and $v > 0.1c$. This gives a sample of 18 sources. The maximum TS is 13.0 (3.2σ), corresponding to a best-fit flux of $\sim 2.0 \times 10^{-11} \text{ ph cm}^{-2} \text{ s}^{-1}$ and a best-fit index of $\sim -2.4$. These results, although less significant, are in good agreement with those from our benchmark sample and show that there is γ-ray emission associated with UFOs independently of how these sources were selected.
<figure><img src="image_12.png" /><figcaption>**Figure 8.** Comparison of redshift and X-ray luminosity (4–10 keV) for the control sample and benchmark sample, as indicated in the legend.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 9.** Stacked profile for our control sample consisting of 20 nearby (z < 0.1) radio-quiet AGN with no UFOs (i.e., a UFO has been searched for but none has been detected). No signal is detected, with a maximum TS of 1.1.</figcaption></figure>

### 5.3. Emission from Star Formation Activity

Star-forming galaxies are known γ-ray emitters because of their CR population, which is accelerated at the shock fronts of supernova remnants and pulsar wind nebulae (Ajello et al. 2020a). The ensuing γ-ray emission is known to correlate well with the total IR luminosity (8–1000 µm), which is a tracer of star formation.

We find that the average total IR luminosity is $log(L_\odot) = 10.4$ (see Table 2). According to the correlation reported in Ajello et al. (2020a), this implies an average γ-ray luminosity (>1 GeV) of $2.2 \times 10^{39}$ erg s⁻¹. This is about 40 times smaller than the observed luminosity and implies that the contamination due to star formation activity to the signal observed in the UFO sample is negligible.

As an additional test we searched for IR fluxes for the galaxies in the control sample (see Table 3). We could find data for nine galaxies with an average total IR luminosity of $log(L_\odot) = 10.8$ (compared to 10.4 for the benchmark sample).

The stacking of this subset of galaxies in the control sample yields no detection (TS = 0.04 and 95% flux UL = 1.1 × 10⁻¹¹ ph cm⁻² s⁻¹) confirming that the contamination of the signal due to star formation is negligible.

### 5.4. Emission from Potential Jets in Radio-quiet AGN

The vast majority of the γ-ray sources detected by the LAT are powered by relativistic jets closely aligned to the line of sight (Ajello et al. 2020b). Some of the sources in our sample, particularly NGC 4151, may have a jet. However, there are several reasons why the γ-ray emission that we observe is unlikely to be produced by the jets, which may be present in these radio-quiet AGN. The best-studied system<fnref n="80" /> is NGC 4151, for which an elongated series of knots, possibly associated with a jet, have been detected in radio (Johnston et al. 1982; Wilson & Ulvestad 1982). This jet has an angle of ≈40° with respect to the line of sight and a speed ≈0.04c (Williams et al. 2017). This is among the lowest speeds measured for a jet and indicates nonrelativistic motion, likely due to thermal plasma (Ulvestad et al. 2005). NGC 4151's jet lies on the opposite end of the spectrum of jets detected by the LAT, which are aligned often within 1°–2° (Pushkarev et al. 2017), highly relativistic (Lister et al. 2016), dominated by nonthermal emission, and found only in radio-loud AGN (Ajello et al. 2020b).

Moreover, the emission from jets is not expected to correlate with the bolometric luminosity of radio-quiet AGN or the outflow kinetic power. It should also be noted that the sources in our sample follow the $L_{22\,GHz}/L_{14-195\,keV} \sim 10^{-5}$ trend, indicating a contribution to the radio luminosity from the hot AGN corona (Smith et al. 2020). Finally, the analysis of winds and jets in a sample of radio-loud AGN provides evidence for a wind-jet bimodality, where winds are the strongest when jets are the weakest (as measured by the radio-loudness parameter in Mehdipour & Costantini 2019).

More importantly, the same nine galaxies in the control sample for which we could find IR data also have 1.4 GHz fluxes (see Table 3). This sample is well matched in terms of radio fluxes and redshift to our benchmark sample and as reported above yields no γ-ray detection.

## 6. SED Modeling

We assume, as in Wang & Loeb (2016a) and Lamastra et al. (2017), that the γ-ray emission is dominated by hadronic processes resulting from diffusive shock acceleration (DSA). In order to model these processes in detail, we first calculate proton distributions using the Cosmic Ray Analytical Fast Tool (CRAFT), a code that uses a semi-analytical formalism for DSA described in Blasi (2002), Amato & Blasi (2006), Caprioli et al. (2010), and references therein. CRAFT self-consistently solves the diffusion-convection equation (e.g., Skilling 1975) for the transport of nonthermal particles in a quasi-parallel, nonrelativistic shock, including the dynamical effects of both accelerated particles and the magnetic turbulence they generate Caprioli (2011, 2012). CRAFT also uses microphysical information (particle injection, diffusion, magnetic field amplification) tuned on self-consistent kinetic plasma simulations of nonrelativistic shocks

<footnote n="80">⁸⁰ Other sources like NGC 5506 and NGC 7582 do not have resolved radio jets down to 0.1", while MCG-5-23-16 has a resolved morphology suggesting the presence of a jet (Orienti & Prieto 2010).</footnote>
<figure><img src="image_14.png" /><figcaption>**Figure 10.** Ambient gas density (red dotted line) and velocity (blue solid line) profiles used in our UFO model. Both profiles come from the models calculated in Liu et al. (2018), for an AGN bolometric luminosity of $L_{\text{Bol}} = 2.5 \times 10^{44}$ erg s⁻¹.</figcaption></figure>

(Caprioli & Spitkovsky 2014a, 2014b, 2014c; Caprioli et al. 2015; Haggerty & Caprioli 2019). Thus, given basic information about UFO shock hydrodynamics (age, velocity, and ambient density), CRAFT self-consistently predicts an instantaneous proton distribution.

To model the cumulative photon distribution of a UFO, we use the hydrodynamic model for the forward shock evolution calculated in Liu et al. (2018) and shown in Figure 10. More specifically, Liu et al. (2018) calculate the forward shock evolution (as in Lamastra et al. 2017; Wang & Loeb 2016a, 2016b, and elsewhere in the literature) using the thin-shell approximation, in which a spherically symmetric shell of negligible thickness expands due to the pressure of a hot bubble inside it. Liu et al. (2018) adopt a broken power-law density profile for the ambient gas, ∝R⁻², inside the disk radius and ∝R⁻³.⁹⁵ outside the disk. However, Liu et al. (2018) also include a flat core in the inner 100 pc of the galaxy to prevent high central densities that are inconsistent with observations, as well as a constant density beyond the virial radius of the galaxy to account for the presence of the intergalactic medium. This profile reproduces well the stellar velocity dispersion in the bulge of the galaxies in our sample (see Figure 1). Both the forward shock evolution and density profile apply to the case of an AGN with a bolometric luminosity of $L_{\text{Bol}} = 2.5 \times 10^{44}$ erg s⁻¹ (consistent with our measurement) and are both shown in Figure 10. Of course, the use of a 1D model has limitations; it cannot account for a more complex ambient medium, meaning that inferred values such as the forward shock age and radius are only approximate. However, given that the model in Liu et al. (2018) yields γ-ray spectra in good agreement with observations, this calculation demonstrates that the γ-ray emission reported in this work can be explained by a population of UFOs with reasonable parameters.

After using CRAFT to calculate the proton distribution produced at each time step of the shock evolution (see Figure 10), the resulting instantaneous distributions are weighted and shifted in energy to account for adiabatic losses as in Diesing & Caprioli (2019). Energy losses due to proton–proton collisions—which are eventually responsible for the UFO’s γ-ray emission—are also taken into account by calculating the collision rate for each distribution at each time step, assuming a target proton density given by the

adiabatically expanded postshock density of a given shell. We further assume that a proton loses half its energy in a single collision (i.e., we assume an inelasticity $\kappa = 0.5$, as in Liu et al. 2018). In other words, the accelerated proton population is treated as a series of adiabatically expanding shells, with the outermost shell located at the forward shock. Each of these shells experience proton–proton collisions—and by extension —produce γ-rays at every time step. Thus, to calculate a UFO’s γ-ray spectrum at a given time, we simply take these weighted proton distributions and convert them to photon spectra using the radiative processes code *naima* (Zabalza 2015). We then add these photon spectra together to produce a cumulative SED. Note that Coulomb losses are neglected in this calculation, as they are subdominant for protons with energies ≳1 GeV (Mannheim & Schlickeiser 1994).

The result is an estimate of a UFO’s SED at every stage of its evolution, as shown in the left panel of Figure 11. We obtain γ-ray luminosities consistent with those calculated in Liu et al. (2018) and find that the observed γ-ray emission can be explained by a forward shock that has traveled a distance between 0.02 and 0.3 kpc from the SMBH (age of $t = (0.3 – 10) \times 10^5$ yr). The modeled CR and γ-ray light curves of the UFO are also shown in the right panel of Figure 11. It is worth noting that the total energy in CRs—and thus the UFO’s γ-ray luminosity—naturally cuts off after roughly 10 Myr due to the fact that the ambient density in the reference galaxy decreases substantially with radius, thereby reducing the available energy flux across the shock.

We also estimate instantaneous electron distributions from our instantaneous proton distributions by using the formalism in Zirakashvili & Aharonian (2007) and accounting for the effects of both adiabatic and synchrotron losses in our weighting (see Diesing & Caprioli 2019). To confirm that the UFO’s synchrotron emission remains below the average radio upper limit from Table 2, the relative normalization of these electron distributions is taken to be a factor of a few larger than that needed to fit observations of Tycho’s supernova remnant (Morlino & Caprioli 2012). Again using *naima*, we then calculate the leptonic emission of a typical UFO from the weighted electron distributions, adding together the contribution of each shell to produce a cumulative SED at a given time step. As shown in Figure 11, the resulting synchrotron emission always falls below the measured average radio emission of the galaxies and the inverse-Compton and bremsstrahlung emissions are a factor >25 fainter than the $\pi_0$ emission. Note that the inverse-Compton emission is estimated by assuming that electrons scatter off the cosmic microwave background and starlight approximated by a blackbody with temperature $T = 3000$ K. This emission remains a factor >10 below the $\pi_0$ one even with an artificially enhanced stellar radiation field of energy density 100 eV cm⁻³. We also model the inverse-Compton emission assuming electrons scatter off the AGN photon field described in Sazonov et al. (2004), normalized to the bolometric luminosity of the AGN sample, and find that this emission remains a factor >3 below that produced by $\pi_0$ decay (and with a much softer spectrum above 30 GeV).

Proton–proton interactions produce γ-rays with energy $E_{\gamma} \approx E_{\text{p}}/10$ (Kelner et al. 2006), and thus the observed γ-ray SED indicates a firm detection of CR protons with energies reaching at least as high as ≈10¹²–¹³ eV. Within our hadronic emission model we derive that the maximum energy of protons accelerated at the forward shock is ≈10¹⁷ eV. This makes AGN
<figure><img src="image_15.png" /><figcaption>**Figure 11.** Left: predicted multiwavelength SED of the UFO's nonthermal emission as a function of time. Synchrotron emission (dotted curves), bremsstrahlung emission (dashed curves), inverse-Compton emission (thin solid curves), and emission from $\pi^0$ decay (thick solid curves) are shown. The inverse-Compton emission remains subdominant despite assuming an artificially enhanced stellar radiation field of energy density 100 eV cm$^{-3}$. Also overlaid is the observed $\gamma$-ray flux as shown in Figure 3 and the average radio upper limit from Table 2. Note that the leptonic emission produced at early times often does not appear as it falls below the plot range. Right: light curve of a UFO-powered forward shock moving through a representative galaxy. The total energy in CRs is shown before and after proton–proton losses are included (blue dotted and dashed lines, respectively), as is the $\gamma$-ray luminosity at 1 GeV (red solid line).</figcaption></figure>

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

winds a potential source of CRs with energies beyond the *knee* of the CR spectrum (i.e., $3 \times 10^{15}$ eV) and also likely contributors to the IceCube neutrino flux (Aartsen et al. 2013).

## 7. Discussion

This work has provided evidence for the existence of a new population of $\gamma$-ray emitters produced by AGN-driven outflows, which in the interaction with the ISM can create strong shocks able to energize charged particles potentially up to the transition region between Galactic and extragalactic CRs. These charged particles produce the observed $\gamma$-rays in the interaction with the ISM. According to our and other available models (Liu et al. 2018), the observed emission can be explained by a forward shock that has traveled $\gtrsim 20$ pc from the central SMBH. As such, the $\gamma$-ray emission from UFOs may signal the onset of the wind-host interaction. We caution the reader, however, that our model relies on a one-dimensional description of a galaxy and that it does not include the complex environment in the immediate vicinity of the SMBH. Nevertheless, our results are found to be in reasonably good agreement with previous predictions (Lamastra et al. 2016; Liu et al. 2018).

Most of the outflow energy is deposited in the bubble of hot gas rather than CRs. More precisely, this energy is a factor of $\sim 10$ larger than what is transferred to accelerated CRs, which is reported in the right panel of Figure 11. For an AGN with $\log L_{\text{Bol}} (\text{erg s}^{-1}) = 44.4$ the timescale to transfer $10^{56}$ erg of energy to the bubble is about 3 Myr. This timescale would grow slightly for less powerful AGN. Sgr A*, the SMBH in the center of the Milky Way, has an estimated mass of $4 \times 10^6 M_{\odot}$ (Abuter et al. 2019), and was very likely in an active state up to a few hundred years ago (Sunyaev et al. 1993; Koyama et al. 1996). Adopting a bolometric luminosity of $\log L_{\text{Bol}} (\text{erg s}^{-1}) = 42.7$ (about 1% of its Eddington luminosity), we find that the timescale to deposit an energy of $10^{56}$ erg to the thermal gas is $\approx 10$ Myr. At that point in time, its $\gamma$-ray luminosity would be $\log L_{\gamma} (\text{erg s}^{-1}) \approx 38$ and would decline modestly in a few hundred years after the end of the AGN activity (see also Liu et al. 2018). This is in reasonable

agreement with the luminosity of the bubbles discovered by Fermi in our galaxy ($\log L_{\gamma} = 37.6 \text{ erg s}^{-1} \text{ Su et al. 2010}$) and the energetics of the thermal gas contained in the larger bubbles ($\sim 10^{56}$ erg) recently discovered by the Extended Roentgen Survey with an Imaging Telescope Array (eROSITA; Predehl et al. 2020). Thus, the Fermi and eROSITA bubbles may be the remnant of past UFO-like activity from the SMBH in the center of our galaxy.

It is important to note that the physical implications that can be inferred from the $\gamma$-ray detection are limited by the incompleteness of the current sample of UFOs, as well as the inherent uncertainty relating to the time variability of the UFOs. Indeed, detection of UFOs is limited to $\lesssim 50$ AGN (Tombesi et al. 2010b; Kara et al. 2016; Igo et al. 2020), which is by far not a complete sample. Moreover, UFOs have been found in these AGN to vary with time and energy and this has been interpreted as a series of expanding shells (see, e.g., King & Pounds 2015) rather than a continuous outflow like we have assumed here. Although we consider it unlikely, the variability may also impact the selection of a control sample as those AGN may not show a UFO precisely at the time when they were observed, but otherwise have an active UFO. A solution to these issues will be provided with the more sensitive observations that the X-ray Imaging and Spectroscopy Mission (XRISM; XRISM Science Team 2020) and the Advanced Telescope for High-Energy Astrophysics (Athena$^{81}$) will provide.

## 8. Summary and Conclusion

To search for the collective UFO emission, a stacking technique that has been used with success in the past (Abdollahi et al. 2018; Paliya et al. 2019; Ajello et al. 2020a) is adopted. Our sample consists of all radio-quiet UFOs with $z < 0.1$ and $v > 0.1c$, which gives 11 sources in total. We model the UFO ($\gamma$-ray) spectrum with a power law, and we assume that the population can be characterized by an average flux and photon index. A fit to all the regions then optimizes

<footnote>$^{81}$ https://www.the-athena-x-ray-observatory.eu</footnote>
these parameters. We find a TS of 30.1, which corresponds to a detection significance for the UFO emission of $5.1\sigma$ (2 dof). The best-fit parameters are measured to be $\Gamma = -2.1 \pm 0.3$ and flux (1–800 GeV) $= 2.5_{-0.9}^{+1.5} \times 10^{-11} \text{ ph cm}^{-2} \text{ s}^{-1}$.

We performed several tests to confirm that the γ-ray emission is truly related to the presence of UFOs in this sample of galaxies. We employed a control sample of AGN with similar properties to those of the 11 galaxies used above, but lacking UFOs. This sample yields no detectable γ-ray emission with a (1–800 GeV) flux upper limit of $8.8 \times 10^{-12} \text{ ph cm}^{-2} \text{ s}^{-1}$. We also use a sample of UFOs selected in a different way (Igo et al. 2020) than our benchmark sample. These galaxies show a γ-ray signal whose parameters are in good agreement with those reported above. Moreover, adopting a control sample matched in X-ray flux, IR luminosity, radio flux, and redshifts we can exclude that the observed γ-ray emission arises from star formation activity or the presence of a weak jet. These tests allow us to conclude that the observed emission is associated to the presence of UFOs in these galaxies.

Observations of AGN winds have shown that AGN transfer a small fraction (~1%–5%) of their bolometric luminosity to the winds. As our analysis indicates, a portion of this transferred luminosity in turn accelerates CRs and produces γ-rays. We find that AGN convert $\approx 3 \times 10^{-4}$ of their bolometric luminosity into γ rays. We also find that $\approx 4 \times 10^{-4}$ of the wind mechanical power is transferred to γ-rays. For comparison, in the Milky Way galaxy, supernova explosions transfer $\approx 2 \times 10^{-4}$ of their mechanical energy to γ-rays. This shows that AGN winds, if sustained for a few million years, can energize a large fraction of the CR population within a galaxy.

The physical model for the UFO SED is calculated by assuming that the γ-ray emission is dominated by hadronic processes resulting from DSA. For typical UFO shock velocities and densities, a leptonic origin of the γ-ray emission is disfavored, in that inverse-Compton scattering and bremsstrahlung of relativistic electrons would produce steeper γ-ray spectra with a lower normalization. The observed γ-ray SED indicates a firm detection of CR protons with energies reaching at least as high as $\approx 10^{12-13}$ eV.

Within our hadronic emission model we derive that on average the forward shock has traveled ~20–300 pc (~65–980 lt-yr) away from the SMBH and that the maximum energy of protons accelerated at the forward shock is $\approx 10^{17}$ eV. This makes AGN winds a potential source of CRs with energies beyond the knee of the CR spectrum (i.e., $3 \times 10^{15}$ eV) and also likely contributors to the IceCube neutrino flux (Aartsen et al. 2013; Padovani et al. 2018). Lastly, our results support the hypothesis that the Fermi and eROSITA bubbles may be the remnant of past UFO-like activity from the SMBH in the center of our galaxy.

M.A. and C.K. acknowledge support from NSF and NASA through grants AST-1715256, 80NSSC21K1915, and 80NSSC18K1718. R.D. and D.C. acknowledge the Eugene & Niesje Parker Fellowship Fund, NASA (grants NNX17AG30G, 80NSSC18K1218, and 80NSSC18K1726) and the NSF (grants AST-1714658, AST-1909778). G.C. acknowledges financial support from NASA grants 80NSSC20K0438 and 80NSSC19K095.

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. Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d'Études Spatiales in France. This work performed in part under DOE Contract DE-AC02-76SF00515. Work at NRL is supported by NASA.

## ORCID iDs

M. Ajello https://orcid.org/0000-0002-6584-1703
L. Baldini https://orcid.org/0000-0002-9785-7726
J. Ballet https://orcid.org/0000-0002-8784-2977
D. Bastieri https://orcid.org/0000-0002-6954-8862
R. Bellazzini https://orcid.org/0000-0002-2469-7063
E. Bissaldi https://orcid.org/0000-0001-9935-8106
R. Bonino https://orcid.org/0000-0002-4264-1215
D. Caprioli https://orcid.org/0000-0003-0939-8775
G. Chartas https://orcid.org/0000-0003-1697-6596
C. C. Cheung https://orcid.org/0000-0002-4377-0174
F. D'Ammando https://orcid.org/0000-0001-7618-7527
P. de la Torre Luque https://orcid.org/0000-0002-4150-2539
A. Desai https://orcid.org/0000-0001-7405-9994
R. Diesing https://orcid.org/0000-0002-6679-0012
N. Di Lalla https://orcid.org/0000-0002-7574-1298
L. Di Venere https://orcid.org/0000-0003-0703-824X
A. Domínguez https://orcid.org/0000-0002-3433-4610
A. Franckowiak https://orcid.org/0000-0002-5605-2219
Y. Fukazawa https://orcid.org/0000-0002-0921-8837
S. Funk https://orcid.org/0000-0002-2012-0080
P. Fusco https://orcid.org/0000-0002-9383-2425
F. Gargano https://orcid.org/0000-0002-5055-6395
D. Gasparini https://orcid.org/0000-0002-5064-9495
N. Giglietto https://orcid.org/0000-0002-9021-2888
M. Giroletti https://orcid.org/0000-0002-8657-8852
S. Guiriec https://orcid.org/0000-0001-5780-8770
D. Hartmann https://orcid.org/0000-0002-8028-0991
G. Jóhannesson https://orcid.org/0000-0003-1458-7036
C. Karwin https://orcid.org/0000-0002-6774-3111
M. Kerr https://orcid.org/0000-0002-0893-4073
M. Kovačević https://orcid.org/0000-0003-4928-4510
M. Kuss https://orcid.org/0000-0003-1212-9998
S. Larsson https://orcid.org/0000-0003-0716-107X
L. Latronico https://orcid.org/0000-0002-0984-1856
M. Lemoine-Goumand https://orcid.org/0000-0002-4462-3686
F. Longo https://orcid.org/0000-0003-2501-2270
F. Loparco https://orcid.org/0000-0002-1173-5673
P. Lubrano https://orcid.org/0000-0003-0221-4806
S. Maldera https://orcid.org/0000-0002-0698-4421
A. Manfreda https://orcid.org/0000-0002-0998-4953

S. Marchesi https://orcid.org/0000-0001-5544-0749

L. Marcotulli https://orcid.org/0000-0002-8472-3649

M. N. Mazziotta https://orcid.org/0000-0001-9325-4672

I. Mereu https://orcid.org/0000-0003-0219-4534

T. Mizuno https://orcid.org/0000-0001-7263-0296

M. E. Monzani https://orcid.org/0000-0002-8254-5308

A. Morselli https://orcid.org/0000-0002-7704-9553

I. V. Moskalenko https://orcid.org/0000-0001-6141-458X

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

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

M. Persic https://orcid.org/0000-0003-1853-4900

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

S. Rainò https://orcid.org/0000-0002-9181-0345

R. Rando https://orcid.org/0000-0001-6992-818X

B. Rani https://orcid.org/0000-0001-5711-084X

M. Razzano https://orcid.org/0000-0003-4825-1629

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

O. Reimer https://orcid.org/0000-0001-6953-1385

D. Serini https://orcid.org/0000-0002-9754-6530

C. Sgrò https://orcid.org/0000-0001-5676-6214

D. Tak https://orcid.org/0000-0002-9852-2469

D. F. Torres https://orcid.org/0000-0002-1522-9065

E. Troja https://orcid.org/0000-0002-1869-7817

## References

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

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

Abdollahi, S., Ackermann, M., Ajello, M., et al. 2018, *Sci*, **362**, 1031

Abuter, R., Amorim, A., Bauböck, M., et al. 2019, *A&A*, **625**, L10

Acero, F., Ackermann, M., Ajello, M., et al. 2016, *ApJS*, **224**, 8

Ackermann, M., Ajello, M., Albert, M., et al. 2011, *PhRvL*, **107**, 241302

Ajello, M., Angioni, R., Axelsson, M., et al. 2020b, *ApJ*, **892**, 105

Ajello, M., Di Mauro, M., Paliya, V. S., & Garrappa, S. 2020a, *ApJ*, **894**, 88

Ajello, M., Gasparrini, D., Sanchez-Conde, M., et al. 2015, *ApJL*, **800**, L27

Alonso-Herrero, A., Almeida, C. R., Mason, R., et al. 2011, *ApJ*, **736**, 82

Amato, E., & Blasi, P. 2006, *MNRAS*, **371**, 1251

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

Atwood, W. B., Atwood, W. B., Baldini, L., et al. 2013, *ApJ*, **774**, 76

Beifiori, A., Courteau, S., Corsini, E., & Zhu, Y. 2012, *MNRAS*, **419**, 2497

Bennert, V. N., Treu, T., Auger, M. W., et al. 2015, *ApJ*, **809**, 20

Bian, W., & Gu, Q. 2007, *ApJ*, **657**, 159

Blasi, P. 2002, *ApH*, **16**, 429

Caprioli, D. 2011, *JCAP*, **05**, 026

Caprioli, D. 2012, *JCAP*, **7**, 38

Caprioli, D., Amato, E., & Blasi, P. 2010, *APh*, **33**, 307

Caprioli, D., Pop, A., & Spitkovsky, A. 2015, *ApJ*, **798**, 28

Caprioli, D., & Spitkovsky, A. 2014a, *ApJ*, **783**, 91

Caprioli, D., & Spitkovsky, A. 2014b, *ApJ*, **794**, 46

Caprioli, D., & Spitkovsky, A. 2014c, *ApJ*, **794**, 47

Chartas, G., Brandt, W. N., Gallagher, S. C., & Garmire, G. P. 2002, *ApJ*, **579**, 169

Chartas, G., Cappi, M., Hamann, F., et al. 2016, *ApJ*, **824**, 53

Cicone, C., Maiolino, R., Sturm, E., et al. 2014, *A&A*, **562**, A21

Condon, J. J., Cotton, W. D., Greisen, E. W., et al. 1998, *AJ*, **115**, 1693

Crenshaw, D. M., & Kraemer, S. B. 2012, *ApJ*, **753**, 75

Diesing, R., & Caprioli, D. 2019, *PhRvL*, **123**, 071101

Fabian, A. C. 2012, *ARA&A*, **50**, 455

Ferrarese, L., & Ford, H. 2005, *SSRv*, **116**, 523

Feruglio, C., Fiore, F., Carniani, S., et al. 2015, *A&A*, **583**, A99

Gallo, L. C., & Fabian, A. C. 2011, *MNRAS*, **418**, L59

Gebhardt, K., Bender, R., Bower, G., et al. 2000, *ApJL*, **539**, L13

Gofford, J., Reeves, J., McLaughlin, D., et al. 2015, *MNRAS*, **451**, 4169

Gofford, J., Reeves, J. N., Tombesi, F., et al. 2013, *MNRAS*, **430**, 60

Grier, C., Martini, P., Watson, L., et al. 2013, *ApJ*, **773**, 90

Gültekin, K., Richstone, D. O., Gebhardt, K., et al. 2009, *ApJ*, **698**, 198

Haggerty, C. C., & Caprioli, D. 2019, *ApJ*, **887**, 165

Hopkins, P. F., & Elvis, M. 2010, *MNRAS*, **401**, 7

Igo, Z., Parker, M. L., Matzeu, G. A., et al. 2020, *MNRAS*, **493**, 1088

James, F., & Roos, M. 1975, *CoPhC*, **10**, 343

Johnston, K. J., Elvis, M., Kjer, D., & Shen, B. S. P. 1982, *ApJ*, **262**, 61

Kara, E., Alston, W. N., Fabian, A. C., et al. 2016, *MNRAS*, **462**, 511

Kaspi, S., Maoz, D., Netzer, H., et al. 2005, *ApJ*, **629**, 61

Kelner, S. R., Aharonian, F. A., & Bugayov, V. V. 2006, *PhRvD*, **74**, 034018

King, A., & Pounds, K. 2015, *ARA&A*, **53**, 115

King, A. R. 2010, *MNRAS*, **402**, 1516

King, A. R., & Pounds, K. A. 2003, *MNRAS*, **345**, 657

Kleinmann, S. G., Cutri, R. M., Young, E. T., Low, F. J., & Gillett, F. C. 1986, *yCat*, **II**, 126

Kormendy, J., & Ho, L. C. 2013, *ARA&A*, **51**, 511

Koyama, K., Maeda, Y., Sonobe, T., et al. 1996, *PASJ*, **48**, 249

Lamastra, A., Fiore, F., Guetta, D., et al. 2016, *A&A*, **596**, A68

Lamastra, A., Menci, N., Fiore, F., et al. 2017, *A&A*, **607**, A18

Lister, M. L., Aller, M. F., Aller, H. D., et al. 2016, *AJ*, **152**, 12

Liu, R.-Y., Murase, K., Inoue, S., Ge, C., & Wang, X.-Y. 2018, *ApJ*, **858**, 9

Longinotti, A., Vega, O., Krongold, Y., et al. 2018, *ApJL*, **867**, L11

Mannheim, K., & Schlickeiser, R. 1994, *A&A*, **286**, 983

Marinucci, A., Bianchi, S., Nicastro, F., Matt, G., & Goulding, A. D. 2012, *ApJ*, **748**, 130

Mattox, J. R., Bertsch, D. L., Chiang, J., et al. 1996, *ApJ*, **461**, 396

Matzeu, G. A., Reeves, J. N., Braito, V., et al. 2017, *MNRAS*, **472**, L15

McCarthy, I. G., Schaye, J., Ponman, T. J., et al. 2010, *MNRAS*, **406**, 822

McNamara, B. R., & Nulsen, P. E. J. 2007, *ARA&A*, **45**, 117

Mehdipour, M., & Costantini, E. 2019, *A&A*, **625**, A25

Morlino, G., & Caprioli, D. 2012, *A&A*, **538**, A81

Moshir, M., Copan, G., Conrow, T., et al. 1990, *yCat*, **II**, 156A

Nardini, E., Reeves, J. N., Gofford, J., et al. 2015, *Sci*, **347**, 860

Nicastro, F., Martocchia, A., & Matt, G. 2003, *ApJL*, **589**, L13

Onken, C. A., Valluri, M., Brown, J. S., et al. 2014, *ApJ*, **791**, 37

Orienti, M., & Prieto, M. A. 2010, *MNRAS*, **401**, 2599

Padovani, P., Turcati, A., & Resconi, E. 2018, *MNRAS*, **477**, 3469

Paliya, V. S., Dominguez, A., Ajello, M., Franckowiak, A., & Hartmann, D. 2019, *ApJL*, **882**, L3

Panessa, F., Baldi, R. D., Laor, A., et al. 2019, *NatAs*, **3**, 387

Peterson, B. M., Ferrarese, L., Gilbert, K. M., et al. 2004, *ApJ*, **613**, 682

Pounds, K. A., Reeves, J. N., King, A. R., et al. 2003, *MNRAS*, **345**, 705

Predehl, P., Sunyaev, R. A., Becker, W., et al. 2020, *Natur*, **588**, 227

Pushkarev, A. B., Kovalev, Y. Y., Lister, M. L., & Savolainen, T. 2017, *MNRAS*, **468**, 4992

Reeves, J., O'Brien, P. T., & Ward, M. 2003, *ApJL*, **593**, L65

Saez, C., & Chartas, G. 2011, *ApJ*, **737**, 91

Sazonov, S. Y., Ostriker, J. P., & Sunyaev, R. A. 2004, *MNRAS*, **347**, 144

Silk, J., & Rees, M. J. 1998, *A&A*, **331**, L1

Skilling, J. 1975, *MNRAS*, **172**, 557

Smith, K. L., Mushotzky, R. F., Koss, M., et al. 2020, *MNRAS*, **492**, 4216

Somerville, R. S., Hopkins, P. F., Cox, T. J., Robertson, B. E., & Hernquist, L. 2008, *MNRAS*, **391**, 481

Su, M., Slatyer, T. R., & Finkbeiner, D. P. 2010, *ApJ*, **724**, 1044

Sunyaev, R. A., Markevitch, M., & Pavlinsky, M. 1991, *ApJ*, **407**, 606

Tombesi, F., Cappi, M., Reeves, J., & Braito, V. 2012, *MNRAS*, **422**, L1

Tombesi, F., Cappi, M., Reeves, J., et al. 2010a, *A&A*, **521**, A57

Tombesi, F., Sambruna, R., Reeves, J., et al. 2010b, *ApJ*, **719**, 700

Ulvestad, J. S., Wong, D. S., Taylor, G. B., Gallimore, J. F., & Mundell, C. G. 2005, *AJ*, **130**, 936

Vasudevan, R., Fabian, A., Gandhi, P., Winter, L., & Mushotzky, R. 2010, *MNRAS*, **402**, 1081

Vasudevan, R. V., & Fabian, A. C. 2007, *MNRAS*, **381**, 1235

Vasudevan, R. V., & Fabian, A. C. 2009, *MNRAS*, **392**, 1124

Wang, X., & Loeb, A. 2016a, *NatPh*, **12**, 1116

Wang, X., & Loeb, A. 2016b, *JCAP*, **2016**, 012

Williams, D. R. A., McHardy, I. M., Baldi, R. D., et al. 2017, *MNRAS*, **472**, 3842

Wilson, A. S., & Ulvestad, J. S. 1982, *ApJ*, **263**, 576

Woo, J.-H., Treu, T., Barth, A. J., et al. 2010, *ApJ*, **716**, 269

XRISM Science Team 2020, arXiv:2003.04962

Zabalza, V. 2015, Proc. ICRC (The Hague), **34**, 922

Zirakashvili, V. N., & Aharonian, F. 2007, *A&A*, **465**, 695

Zubovas, K., & King, A. R. 2012, *MNRAS*, **426**, 2751