# Discovery of a new extragalactic population of energetic particles

Anthony M. Brown, Céline Bœhm, Jamie Graham, Thomas Lacroix, Paula Chadwick, Joseph Silk

► To cite this version:

Anthony M. Brown, Céline Bœhm, Jamie Graham, Thomas Lacroix, Paula Chadwick, et al. Discovery of a new extragalactic population of energetic particles. Physical Review D, 2017, 95 (6), pp.063018. 10.1103/PhysRevD.95.063018. hal-01554310

HAL Id: hal-01554310

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

Submitted on 5 Oct 2025

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

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

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

# Discovery of a new extragalactic population of energetic particles

Anthony M. Brown,¹ Céline Bœhm,²,³ Jamie Graham,¹ Thomas Lacroix,⁴ Paula Chadwick,¹ and Joseph Silk⁴,⁵

¹Department of Physics and Centre for Advanced Instrumentation, University of Durham,
South Road, Durham DH1 3LE, United Kingdom

²Institute for Particle Physics Phenomenology (IPPP), University of Durham,
Durham DH1 3LE, United Kingdom

³LAPTH, U. de Savoie, CNRS, BP 110, 74941 Annecy-Le-Vieux, France

⁴Institut d'Astrophysique de Paris, UMR 7095, CNRS, UPMC Université Paris 6,
Sorbonne Université, 98 bis boulevard Arago, 75014 Paris, France

⁵Department of Physics & Astronomy, Johns Hopkins University, Baltimore, Maryland 21218, USA
(Received 27 January 2017; published 28 March 2017)

We report the discovery of a statistically significant hardening in the Fermi-LAT $\gamma$-ray spectrum of Centaurus A's (Cen A) core, with the spectral index hardening from $\Gamma_1 = 2.73 \pm 0.02$ to $\Gamma_1 = 2.29 \pm 0.07$ at a break energy of ($2.6 \pm 0.3$) GeV. Using a likelihood analysis, we find no evidence for flux variability in Cen A's core light curve above or below the spectral break when considering the entire 8 year period. Interestingly, however, the first $\sim$3.5 years of the low energy light curve shows evidence of flux variability at the $\sim$3.5$\sigma$ confidence level. To understand the origin of this spectral break, we assume that the low energy component below the break feature originates from leptons in Centaurus A's radio jet, and we investigate the possibility that the high energy component above the spectral break is due to an additional source of very high energy particles near the core of Cen A. We show for the first time that the observed $\gamma$-ray spectrum of an active galactic nucleus is compatible with either a very large localized enhancement (referred to as a spike) in the dark matter halo profile or a population of millisecond pulsars. Our work constitutes the first robust indication that new $\gamma$-ray production mechanisms can explain the emission from active galaxies and could provide tantalizing first evidence for the clustering of heavy dark matter particles around black holes.

DOI: 10.1103/PhysRevD.95.063018

## I. INTRODUCTION

At a distance of 3.7 Mpc, Centaurus A (Cen A) is the closest known $\gamma$-ray emitting active galaxy (e.g. [1]). With a Fanaroff-Riley type I (FR I) radio morphology, Cen A possesses a relativistic jet that is orientated at a large angle relative to Earth's line of sight. As such, Cen A is believed to belong to the parent population of BL Lac objects within the unified model of active galactic nuclei (AGN; [2,3]). The off axis nature of Cen A's jet allows for an easy decoupling of the observed relativistic effects from the intrinsic properties of the jet. Furthermore, since the jet emission is weakly beamed, different emission components which typically are not observed in BL Lacs due to the dominant beamed jet component, may become visible.

Due to its proximity, Cen A affords us an excellent opportunity to study the physics of relativistic outflows. As such, Cen A has been extensively studied at many wavelengths. Radio observations have found Cen A to possess a variety of radio structures, with size scales ranging from arcseconds in the inner jet to giant radio lobes which extend $10^\circ$ across the sky. At x-ray energies, emission has been resolved from the inner jet which appears to be well collimated [4,5]. At $\gamma$-ray wavelengths, the first evidence for emission above 1 MeV was provided by balloon-borne experiments [6], with further evidence from the GRANAT

satellite [7]. Cen A was confirmed as a $\gamma$-ray point source by the OSSE and COMPTEL detectors onboard the Compton Gamma-Ray Observatory in the 0.1-30 MeV energy range, with the EGRET detector also discovering Cen A to be a source of $\gamma$-rays up to GeV energies [8,9]. At very high energies (VHE; >250 GeV), Cen A was discovered by the H.E.S.S. telescope array with a soft spectral index ($\Gamma$) of $-2.73 \pm 0.45_{\text{stat}} \pm 0.2_{\text{sys}}$ and a flux at 0.8% of the Crab nebula [10].

The successful launch of the Fermi $\gamma$-ray Space Telescope, and the unprecedented sensitivity and resolution of the Large Area Telescope (LAT) onboard Fermi, affords us an ideal opportunity to understand the inner working of Cen A. The first ten months of LAT observations confirmed the core of Cen A to be a source of MeV and GeV $\gamma$-rays [11], with $\Gamma = -2.69 \pm 0.10_{\text{stat}} \pm 0.08_{\text{sys}}$ and a flux of ($1.50 \pm 0.25_{\text{stat}} \pm 0.37_{\text{sys}}$) $\times 10^{-7} \text{ ph cm}^{-2} \text{ s}^{-1}$. A deeper 4-year analysis of Cen A's core with PASS7 data revealed a similar spectral index and flux, with the authors noting that the spectrum above $\sim$4 GeV appeared to depart from a power-law [12]. This departure was statistically insignificant, with a broken power-law preferred over a power-law at the $\langle 3\sigma \rangle$ level. Interestingly, early Fermi-LAT observations revealed extended $\gamma$-ray emission spatially coincident with Cen A's giant radio lobes [13]. Both core and lobe
emissions appeared to have similar luminosities, and neither showed variability in their flux. The lack of variability in Cen A's flux is at odds with properties of other prominent γ-ray bright radio galaxies, such as M87 and NGC 1275, whose flux varies on very short timescales (as little as days) [14,15].

Building on earlier work that investigated the contribution to M87's γ-ray emission from dark matter annihilation, [16], in this paper we investigate the spectral and temporal properties of the γ-ray emission from Cen A's core. Using 8 years of Fermi-LAT data, and taking advantage of the improvements in effective area, energy and angular resolution afforded by the PASS8 data analysis, we discover a statistically significant hardening feature in Cen A's MeV-GeV spectrum.

The structure of our paper is as follows. After establishing in Sec. II the statistical significance of the spectral hardening feature, we investigate the temporal characteristics of the flux below and above the onset of the spectral feature. In Sec. III we investigate two possible interpretations; one involves the annihilation products from dark matter with a “spike” density profile and the other one millisecond pulsars.

## II. FERMI-LAT OBSERVATIONS AND DATA ANALYSIS

The LAT detector onboard Fermi, described in detail in [17,18], is a pair-conversion telescope observing a photon energy range from below 20 MeV to above 500 GeV, with unprecedented sensitivity and resolution compared to previous γ-ray space missions. The PASS8 data analysis and instrument response functions (IRFs) of the LAT detector has brought further improvement in detector performance. Since 2008 August 4, the vast majority of data taken by Fermi has been in all-sky-survey mode. This observing mode, coupled with the large effective area of the LAT detector and the long mission lifetime of the Fermi mission, has produced the deepest extragalactic scan ever at MeV-GeV energies. This enables us to study the γ-ray emission from Cen A's core with unprecedented sensitivity and accuracy.

The data used in this study comprise of all Fermi-LAT event and spacecraft data taken during the first 8 years of the Fermi-LAT science mission, from 2008 August 4 to 2016 August 8, equating to a mission elapsed time (MET) period of 239557417 [s] to 492328870 [s]. All $0.1 < E_γ < 300$ GeV SOURCE events, across all point spread function (PSF) classes, within a 15° radius of interest (RoI) centred on the Cen A core position were considered. In accordance with criteria for PASS8 data analysis, a zenith cut of 90° was applied, and good time intervals selected by removing data that did not satisfy the criteria "$DATA Qual>0 & &LAT_CONFIG==1$". A summary of the criteria used in the analysis is given in Table 1.

<table><caption>TABLE I. Summary of the criteria utilised in this analysis. Note that since a summed likelihood analysis was used, with PSF classes 0 through 3 considered, four isotropic diffuse models were used.</caption><tr><td>Science tools version</td><td>v10R0p5</td></tr><tr><td>IRF</td><td>P8_SOURCE_V6</td></tr><tr><td>Event class</td><td>SOURCE, Pass 8</td></tr><tr><td>Photone Energies</td><td>0.1–300 GeV</td></tr><tr><td>RoI</td><td>15°</td></tr><tr><td>Zenith angle cut</td><td>&lt;90°</td></tr><tr><td>Data quality</td><td>&gt;0</td></tr><tr><td>LAT config</td><td>1</td></tr><tr><td>Galactic diffuse model</td><td>gl1_lem_v06.fit</td></tr><tr><td>Isotropic diffuse model</td><td>iso_P8R2_SOURCE_V6_PSF0–3.txt</td></tr></table>

The model employed during our likelihood analyses consisted of pointlike γ-ray sources, spatially extended γ-ray sources and diffuse γ-ray emission. In particular, the position and spectral shape of all γ-ray point sources within 25° of Cen A's core were taken from the third Fermi source catalog (3FGL; [19]). The extended sources considered were the two pulsar wind nebulae located within the 15°–25° annulus from Cen A's core, namely HESS J1303-631 and MSH 15-52, and Cen A's radio lobes [13]. The diffuse γ-ray emission detected by the LAT compresses of two components: the galactic diffuse emission and the isotropic diffuse emission. The galactic component of the diffuse emission was modeled with Fermi's gl1_lem_v06.fit model. Since a summed likelihood analysis was used to take advantage of the performance of PASS8's four PSF, the isotropic diffuse emission was defined by the iso_P8R2_SOURCE_V6_PSF0–3.txt files, where 0–3 refers to PSF class 0 through to 3.

To confirm the accuracy of our “diffuse + point + extended” model description, an initial Binned likelihood analysis was performed over the entire data set. Firstly, a likelihood analysis was performed with the normalization of all sources within 15° being left free to vary. Thereafter, a second likelihood fit was performed with all point sources with a test statistic,<fnref n="1" /> TS, greater than 25 able to vary spectrally. Finally, all sources with TS < 1 were removed from our model, and all sources with a 3FGL TSvar > 72.44 were fitted again with the normalization and spectral shape left free to vary. The Fermi science tool GTTSMAP was then used in conjunction with the final best-fit model from the three-step initial analysis, to construct a 17° × 17° TS map centred on Cen A. This TS map was used to identify additional sources of γ-rays that were not accounted for in our best-fit model. As can be seen in Fig. 1, there are three >5σ excesses in the TS map indicating three new point

<footnote n="1">¹The test statistic, TS, is defined as twice the difference between the log-likelihood of two different models, $2[\log L - \log L_0]$, where L and $L_0$ are defined as the likelihoods of individual model fits [20].</footnote>
<figure><img src="image_3.png" /><figcaption>FIG. 1. A 17° × 17° TS map of all 0.1–300 GeV photons that passed the selection criteria, for the entire 8 year data set. The color scale is in units of TS. Three new point sources were found, we note however that two of these apparent point sources are located on the boundary of Cen A’s lobes, shown by the red contour lines.</figcaption></figure>

sources of γ-rays have been found in our data set. These new
point sources were accounted for by a power-law fixed at the
($\alpha_{J2000}$, $\beta_{J2000}$) of the peak excess, and a final likelihood fit
was performed with the normalization and spectral index of
the new point sources free to vary. Interestingly, we note that
two of these point sources are located on the boundary of
Cen A's lobes.

Once the data set was modeled correctly, a summed-
likelihood analysis was used to study the spectrum of Cen A’s
core. The data were binned into ten logarithmically spaced
energy bins, with a likelihood fit being applied to each bin
separately. For each separate likelihood fit, all spectral
parameters were frozen except for the normalization of
Cen A’s core. For an individual spectral bin, if the calculated
flux had a TS ≪ 25, a 2σ upper limit was calculated. The
resulting spectrum can be seen in Fig. 2, with all error bars on
the LAT data points representing a 1σ level of statistical
uncertainty. The derived Cen A core spectrum seen in Fig. 2
shows a clear departure from a simple power-law description
above photon energies of ~1 GeV. To quantify the signifi-
cance of this discrepancy at high energies, a summed-
likelihood analysis of the (0.1–300) GeV spectrum was
undertaken with the Cen A’s core being described by
a broken power-law. A likelihood-ratio test of the broken
power-law and power-law fit to Cen A core’s (0.1–300) GeV
spectrum found that the broken power-law spectral model is
preferred over the power-law spectral model with a test
statistic of TS = 28.6, which equates to a significance slightly
greater than 5σ (even after the increased number of free
parameters of the broken power law). As such, this analysis
represents the first observation of a statistically significant
hardening of an AGN’s γ-ray spectrum above a spectral
feature.

<figure><img src="image_4.png" /><figcaption>FIG. 2. The 0.1 to 300 GeV spectrum of Cen A's core, as seen by the Fermi-LAT detector. The broken power-law model, shown in blue, is preferred over the power law model, shown in red, with a significance >5σ. The grey histogram shows the TS value for each spectral bin. The bins with a TS < 25, are replaced with an upper limit at 95% confidence level.</figcaption></figure>

From the best-fit broken power-law, the spectral break
occurs at a photon energy of $E_{\text{break}} = 2.6 \pm 0.3$ GeV,
with a spectral index of $\Gamma_1 = -2.73 \pm 0.02$ below the
spectral break, hardening to a spectral index of
$\Gamma_2 = -2.29 \pm 0.07$ above the break energy. The total
energy flux from Cen A's core in the 0.1–300 GeV energy
range is $\sim 1.1 \times 10^{-12}$ ergs cm$^{-2}$ s$^{-1}$. Assuming isotropic
emission for the $\gamma$-ray flux and a luminosity distance of
3.7 million parsecs, the total luminosity of Cen A's core in
the 0.1–300 GeV is $1.8 \times 10^{39}$ ergs s$^{-1}$. Interestingly, if we
deconvolve the two spectral components, we find that the
luminosity of the low and high energy components are com-
parable, with $9.9 \times 10^{38}$ ergs s$^{-1}$ and $8.4 \times 10^{38}$ ergs s$^{-1}$
respectively.

To investigate the temporal characteristics of the two
spectral components, we constructed individual light
curves for the γ-ray flux above and below the break energy.
To maintain significant statistics for each temporal bin in
the high-energy light curve, the 8-year data set was binned
into 6-month temporal bins. Using the global best-fit model
with the normalization of Cen A's core left free, a like-
lihood analysis was applied to each bin separately. The
resulting light curves can be seen in Fig. 3.

To determine if there was evidence for flux variability in
the resultant light curves, both a $\chi^2$ fit of the light curves to
a constant flux value and a full likelihood evaluation, $\text{TS}_{\text{var}}$
[21], were used. When considering the 8-year data set,
neither test found significant evidence of variability in
either the low or high energy light curve. It is interesting to
note that when considering the first 3.5 years of the low-
energy light curve, there is evidence of variability at the
$3.3\sigma$ level, with the remaining 4.5 years of the light curve
<figure><img src="image_5.png" /><figcaption>FIG. 3. Above: light curve of 0.1–2.6 GeV flux Below: light curve of 2.6–300 GeV flux, binned in six month temporal bins. The red horizontal band shows the 8-year averaged flux for each energy range.</figcaption></figure>

being consistent with a constant flux. As such, while the flux above the spectral break is statistically consistent with being constant, there is evidence of variability in the low-energy spectral component during some of the time period we studied Cen A.

### III. INTERPRETATION

Our analysis shows a hardening of the $\gamma$-ray spectrum above 2.6 GeV at a significance level of 5$\sigma$. Importantly, while an extrapolation of the power-law description of the *Fermi*-LAT spectrum to very high energies would underpredict the flux observed by the H.E.S.S. telescope array by an order of magnitude, the spectral hardening above 2.6 GeV allows us to reconcile our *Fermi*-LAT spectrum with previous H.E.S.S. observations. Furthermore, the nonvariable flux we observe with *Fermi*-LAT above the spectral break is also in agreement with the H.E.S.S. observations of Cen A finding no evidence of variability in the $E_{\gamma} > 250$ GeV flux. The combined *Fermi*-LAT and H.E.S.S. spectrum is shown in Fig. 4.

With the existence of the hardening feature confirmed, we now turn our attention to understanding the origin of said feature. Traditionally, the $\gamma$-ray emission from radio galaxies has been attributed to a single zone population of relativistic electrons within the radio galaxy’s jet that up-scatters the synchrotron photon field via the inverse Compton process [11], usually cited as the synchrotron self-Compton (SSC) model. In the SSC model, the high-energy tail of the electron population is responsible for both the x-ray and $\gamma$-ray emission and as such, if the SSC model was an accurate description of Cen A’s core, the spectral break we observe in the $\gamma$-ray spectrum, should also be present in the x-ray spectrum, which is not seen in the latest x-ray observations [22]. Furthermore, the lack of variability above the spectral break, combined with the evidence for variability below the spectral break, is difficult to reconcile with a single-zone SSC model as high energy electrons are expected to cool

<figure><img src="image_6.png" /><figcaption>FIG. 4. A combined spectrum of our LAT analysis, shown with black circles, and the H.E.S.S. spectrum above $E_{\gamma} = 250$ GeV, taken from [10], shown in red. The broken power-law model, shown in blue, is preferred over the power law model, shown in red, with a significance >5$\sigma$. The grey histogram, with the right-hand y-axis, shows the TS value for each spectral bin for the LAT data points. The bin with a TS < 25 is replaced with an upper limit at 95% confidence level. Note that the last upper limit is not shown since it overlaps in energy with the H.E.S.S. data points. The broken power law fit allows for a smooth transition with the $\gamma$-ray spectrum reported here and the H.E.S.S. spectrum.</figcaption></figure>

quicker than the low energy electrons, thus generating flux variability above the break energy rather than below.

A number of possible explanations for a spectral hardening, combined with nonvariable emission above 2.6 GeV, are mentioned in [23]. One such model is a pulsarlike magnetospheric acceleration mechanism of electrons, with the entire MeV-GeV $\gamma$-ray flux attributed to inverse-Compton emission from the same lepton population [24]. However, this is inconsistent with both the spectral break and the lack of flux variability above 2.6 GeV. Moreover these magnetospheric models predict that inverse-Compton luminosity, $L_{\text{IC}}$, is proportional to the supermassive black hole mass, $M_{\text{BH}}$, which is not observed across the mass range of prominent $\gamma$-ray bright radio galaxies [15].

Another possible explanation is that the hardening feature is associated with the decay of the pion by-products caused by a hadronic population within the jet interacting with an ambient photon field (eg. [25,26]), as has been done by [27]. However, given the lack of variability above the spectral break, this explanation would require the pion population to be in a “steady-state”.

Building on the apparently smooth transition between the data points above 2.6 GeV derived in this work and the H.E.S.S. observations, [10], we discuss the validity of models that can account for the spectral hardening, by jointly fitting these data points. More specifically, we use the first three H.E.S.S. points, but exclude the last three which are less statistically significant, and more likely to be modified with new observations and an updated analysis. In any event, the last three points do not affect the best-fit spectrum due to their large error bars, but simply lead to a slight increase in the $\chi^2$.
## A. Dark matter

Dark matter (DM), as yet undetected, is a central explanation for structure formation, the stability of galaxies and the acoustic peaks of the cosmic microwave background [28]. Most DM models assume annihilation into lighter (Standard Model) particles but the nature of DM is still an open question. While it is established that DM agglomerates in the cores of galaxies, its density profile near the central black hole (BH) has yet to be characterized. Regardless, DM annihilations should produce high energy cosmic rays in the central parts of galaxies. The decay or hadronization of the particles injected by the DM is expected to produce $\gamma$-rays at an energy smaller than or equal to the DM mass (assuming they are charged and unstable). Here we consider only prompt emission. We disregard any acceleration process [29] or secondary emission [16,30] that could distort the prompt spectrum. Here we examine whether this prompt emission is sufficient to explain the anomalous spectral component between 2.6 GeV and ~5 TeV.

Throughout this work, we have assumed the existence of a spike in the distribution of the DM—induced by the adiabatic growth of a supermassive black hole at the center of the galaxy—with a power-law index of $\gamma_{sp} = 7/3$. Such a spiky profile has been suggested in Ref. [31], although it is very debated. Indeed, a plethora of astrophysical phenomena could lead to a much shallower profile, such as nonadiabatic BH growth (as expected if the BH seed was brought in by a merger) [32,33] or an off-center position of the central BH, while others may have the opposite effect of making the formation or survival of spikes more likely, such as for instance enhanced accretion of DM to counteract the depopulation of chaotic orbits in a triaxial DM halo [34].

It is also important to remember that even if a spike could form, the process of dynamical relaxation by DM scattering off stars could smooth down the spike and lead to a DM halo profile of the form $\rho \propto r^{-3/2}$ instead [32,34]. However, Cen A is dynamically young: its relaxation time is estimated to be $t_r \sim 10^2$ Gyr (compared to ~2.5 Gyr for the Milky Way) due to the dependence on the mass of the central black hole. Hence we will assume that a spike formed in the core of Cen A at early times has survived dynamical processes.

In our best-fit analysis, we consider DM self-annihilations into leptons or quarks and leave the DM mass and value of the annihilation cross section (as well as the normalization and slope of the spectral power law emission below 2.6 GeV) free to vary. As shown in Fig. 5, our best fit favors a DM candidate with a mass of 3 TeV, annihilating into a pair of top and antitop quarks ($t\bar{t}$) with cross section $\langle\sigma v\rangle \approx 1.6 \times 10^{-32} \text{ cm}^2 \text{ s}^{-1}$, and a spike in the density profile. While the impact of a spike on the $\gamma$-ray emission from DM prompt emission in an AGN was first studied in [16], our work constitutes the first evidence

<figure><img src="image_7.png" /><figcaption>FIG. 5. Best fit to the $\gamma$-ray spectrum of Cen A, obtained by assuming a single power law plus the prompt emission from 3 TeV DM particles annihilating into $t\bar{t}$ with a cross section of $1.6 \times 10^{-32} \text{ cm}^2 \text{ s}^{-1}$, and a spike in the DM density profile, with slope $\gamma_{sp} = 7/3$.</figcaption></figure>

that $\gamma$-ray observations have the means to probe an anomalously high concentration of DM in the very inner core of AGNs.

For a DM spike with slope $\gamma_{sp} = 7/3$, the best fit shown in Fig. 5, and corresponding to annihilations into $t\bar{t}$ gives $\chi^2 = 1.7$, for 11 spectral data points. The points are taken from our Fermi-LAT analysis and the first three data points from H.E.S.S. observations [10] and four free parameters ($m_{DM}$, $\langle\sigma v\rangle$ and the normalization and the slope of the power law spectrum below 2.6 GeV). We find a $\chi^2$/d.o.f = 0.24 which illustrates a remarkably good fit and suggests that it is dominated by statistical errors. Annihilations into $b\bar{b}$ give $\chi^2 = 3$, i.e. $\chi^2$/d.o.f = 0.43, which also corresponds to a very good fit (see Fig. 6). If we include the last

<figure><img src="image_8.png" /><figcaption>FIG. 6. Best fit to the $\gamma$-ray spectrum of Cen A, obtained by assuming a single power law plus the prompt emission from 3 TeV DM particles annihilating into $b\bar{b}$ with a cross section of $1.4 \times 10^{-32} \text{ cm}^2 \text{ s}^{-1}$, and a spike in the DM density profile, with slope $\gamma_{sp} = 7/3$.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>FIG. 7. Best fit to the γ-ray spectrum of Cen A, obtained by assuming a single power law plus the prompt emission from 400 GeV DM particles annihilating into $\tau^+\tau^-$ with a cross section of $4 \times 10^{-33} \text{ cm}^2 \text{ s}^{-1}$, and a spike in the DM density profile, with slope $\gamma_{sp} = 7/3$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIG. 8. Contribution to the γ-ray spectrum of Cen A from a single power law plus the prompt emission from 3 TeV DM particles annihilating into $b\bar{b}$ with the canonical cross section of $3 \times 10^{-26} \text{ cm}^2 \text{ s}^{-1}$, for a DM density following the NFW profile. The DM-induced emission is several orders of magnitude below the data.</figcaption></figure>

three H.E.S.S. points, we obtain $\chi^2$/d.o.f. = 0.62 for $\bar{t}\bar{t}$ and $\chi^2$/d.o.f. = 0.87 for $b\bar{b}$.

In Fig. 7, we show the best fit when considering DM particles annihilating into $\tau^+\tau^-$. Now the best-fit mass is ~400 GeV, smaller than for hadronic channels. In this case, we find $\chi^2$/d.o.f. = 2.14, so this channel gives a less good fit than the previous cases, but it still reproduces the spectral hardening reported in this work. For simplicity we have focused only on single annihilation channels; depending on the underlying model, DM particles may annihilate into several different final states.

In the absence of a spike, DM annihilation for a maximum cross section of $\langle\sigma v\rangle \approx 3 \times 10^{-26} \text{ cm}^2 \text{ s}^{-1}$ cannot account for the spectral hardening, since in that case the DM flux is several orders of magnitude smaller than the observed flux. This is illustrated in Fig. 8, for the NFW profile. In practice a spike with slope $\gtrsim 2$ is needed to account for the observed flux.

Our best-fit annihilation cross section of the order of $\langle\sigma v\rangle \approx 1.6 \times 10^{-32} \text{ cm}^2 \text{ s}^{-1}$ is far too small to explain the observed fraction of DM in the Universe. However, this might simply be revealing the existence of a rich dark sector with several (nonthermal) DM particles [35,36],<fnref n="2" /> or a combination of velocity-dependent and independent terms in the annihilation cross section. We note that there is a degeneracy between the annihilation cross section and the spike characteristics (normalization, size, power law index) which could affect our estimates. Either way, our findings would suggest the existence of heavy DM particles exchanging at least one heavy particle mediator.

## B. Millisecond pulsars

An alternative explanation to the DM model is a population of millisecond pulsars (MSPs). These objects are rapidly rotating neutron stars which have been observed mostly in globular clusters. They are believed to spin up to millisecond periods due to frequent interactions with neighboring stars (e.g. by angular momentum accretion from a binary companion). For this very reason, they are expected to form in high stellar density environments, including the central parsec around the Milky Way Galactic center (GC), where the density is considerably higher than in globular clusters.

MSPs, along with DM, are currently the leading explanations for the low energy γ-ray excess which has been observed in the central region of the Milky Way (see e.g. Ref. [38]). It is estimated that around $10^3-10^4$ MSPs are needed to explain the GeV excess in the GC [39]. The production mechanism of high energy γ-rays from MSPs involves electron-positron pairs that may eventually be accelerated, radiate and produce more electron-positron pairs [40]. This should lead to a spectral feature at GeV energies. However the resulting γ-ray signature can be broadened up to TeV energies by inverse Compton processes [41], when the electrons accelerated by MSP winds up-scatter the ambient soft photon field (from e.g. UV and IR bands). Two critical assumptions for these estimates are an electron injection spectrum extending to a few tens of TeV and a large enough interstellar radiation field for the inverse Compton losses to dominate over synchrotron losses. Here we use the same propagation technique as for the Milky Way (see e.g. Ref. [42]). For the magnetic field, we assume a constant value of $10\mu\text{G}$ constrained by Chandra limits on synchrotron radiation from high energy electron-positron pairs and use a conservative power law,

<footnote n="2">²Our estimate assumes a unique DM candidate. Subcomponent DM particles might require efficient coannihilation processes, see e.g. [37].</footnote>
<figure><img src="image_11.png" /><figcaption>FIG. 9. Best fit for a population of MSPs plus a power law.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>FIG. 10. Best fit for a power law, a population of MSPs and a DM candidate of 30 TeV plus a spike.</figcaption></figure>

$\rho_{\text{MSP}} \propto r^{-2.4}$, for the MSP density profile (consistent with the GC $\gamma$-ray excess). As for the DM analysis, we keep the normalization and slope of the jet contribution as free parameters.

Assuming the existence of such a population of MSPs
leads to the best fit shown in Fig. 9, corresponding to
$\chi^2 \approx 10$, that is $\chi^2/\text{d.o.f} \approx 1.4$. This fit improves with the
addition of a DM component at very high energy, as shown
in Fig. 10 where we assumed a 30 TeV DM candidate
annihilating into tau leptons. While the reduced $\chi^2$ is
good, the poor knowledge of the MSPs density profile,
the soft-photon field and our rather crude model of the
magnetic field severely limits our interpretation of the
goodness of the fit.

IV. CONCLUSION

In this paper, we report a 5σ evidence for a hardening of
the Fermi-LAT γ-ray spectrum and show that either heavy
DM particles or a population of MSPs could explain this
high-energy spectral feature. While we cannot rule out that

the jet itself is at the origin of the hardening of the
spectrum, the lack of variability of the emission above
2.6 GeV, both within the Fermi-LAT and H.E.S.S. spectra,
argues against jet-induced leptonic models (such as SSC).
Hence, at the very least, our results are a strong indication
that the modeling of γ-ray production mechanisms in active
galaxies needs to be modified.

Our findings hint at new physics inside objects like Cen A or astrophysical objects that are rarely detected outside our Galaxy. Therefore the precise modeling of cosmic ray propagation, together with observations of Cen A's γ-ray spectrum at TeV energies by H.E.S.S. and the forthcoming Cherenkov Telescope Array [43], will be critical for determining the origin of the spectral break.

A DM explanation would reveal the existence of
particles beyond the standard model of particle physics
as well as spikes in the DM profile. While the existence of
spikes is highly debated, Cen A is dynamically young, so
such an enhancement may survive nuclear star cluster
dynamics. The fit to the combination of Fermi and H.E.S.S.
data in Fig. 10 could provide tantalizing first evidence for
the clustering of heavy dark matter particles around black
holes. If confirmed by other observations, our case for a
boost of the DM annihilation signal constitutes the first
evidence that DM clusters around black holes, and have
important implications for our understanding of the feed-
back mechanisms between ordinary matter and DM in
galactic cores.

Likewise, if the spectral hardening of Cen A's core
reported here is due to a population of MSPs, then the γ-ray
spectrum of Cen A constitutes the first insights into the
pulsar population in another galaxy. Either way, these
findings open up a new window on the physics of the
cores of active galaxies and provide the prospect of more
exciting discoveries.

ACKNOWLEDGMENTS

The authors would like to thank the referee for their
comments which improved the quality and clarity of the
paper. This work has made use of public Fermi-LAT data
obtained from the High Energy Astrophysics Science
Archive Research Centre (HEASARC), provided by
NASA's Goddard Space Flight Center. We thank the
Fermi-LAT collaboration for the quality of the data and
analysis tools that were used in this study. A. M. B.
acknowledges the financial support of Durham
University. J. G. acknowledges the financial support of
his STFC Ph.D. Award No. ST/N50404X/1. The research
of C. B., T. L., and J. S. has received support at IAP from
ERC Project No. 267117 (DARK) hosted by UPMC. T. L.
research has been carried out in the ILP LABEX
(ANR-10-LABX-63) and supported by French state funds
managed by the ANR, within the Investissements d'Avenir
programme (ANR-11-IDEX-0004-02).
[1] L. Ferrarese, J. R. Mould, P. B. Stetson, J. L. Tonry, J. P. Blakeslee, and E. A. Ajhar, Astrophys. J. **654**, 186 (2007).

[2] B. L. Fanaroff and J. M. Riley, Mon. Not. R. Astron. Soc. **167**, 31P (1974).

[3] C. M. Urry and P. Padovani, Publ. Astron. Soc. Pac. **107**, 803 (1995).

[4] I. N. Evans and A. P. Koratkar, Astrophys. J. **617**, 209 (2004).

[5] M. J. Hardcastle, D. M. Worrall, R. P. Kraft, W. R. Forman, C. Jones, and S. S. Murray, Astrophys. J. **593**, 169 (2003).

[6] P. von Ballmoos, R. Diehl, and V. Schoenfelder, Astrophys. J. **312**, 134 (1987).

[7] E. Jourdain et al., Astrophys. J. **412**, 586 (1993).

[8] P. Sreekumar, D. L. Bertsch, R. C. Hartman, P. L. Nolan, and D. J. Thompson, Astropart. Phys. **11**, 221 (1999).

[9] H. Steinle et al., Astron. Astrophys. **330**, 97 (1998).

[10] F. A. Aharonian et al., Astrophys. J. Lett. **695**, L40 (2009).

[11] A. A. Abdo et al., Astrophys. J. **719**, 1433 (2010).

[12] N. Sahakyan, R. Yang, F. A. Aharonian, and F. Reiger, Astrophys. J. Lett. **770**, L6 (2013).

[13] A. A. Abdo et al., Science **328**, 725 (2010).

[14] F. A. Aharonian et al., Science **314**, 1424 (2006).

[15] A. M. Brown and J. Adams, Mon. Not. R. Astron. Soc. **413**, 2785 (2011).

[16] T. Lacroix, C. Bœhm, and J. Silk, Phys. Rev. D **92**, 043510 (2015).

[17] M. Ackermann et al., Astrophys. J. Suppl. Ser. **203**, 4 (2012).

[18] W. B. Atwood et al., Astrophys. J. **697**, 1071 (2009).

[19] F. Acero et al., Astrophys. J. Suppl. Ser. **218**, 23 (2015).

[20] J. R. Mattox et al., Astrophys. J. **461**, 396 (1996).

[21] P. L. Nolan et al., Astrophys. J. Suppl. Ser. **199**, 31 (2012).

[22] F. Fuerst et al., Astrophys. J. **819**, 150 (2016).

[23] F. Rieger, Mem. Soc. Astron. Ital. **83**, 127 (2012).

[24] A. Neronov and F. A. Aharonian, Astrophys. J. **671**, 85 (2007).

[25] A. Mücke, R. Engel, J. P. Rachen, R. J. Protheroe, and T. Stanev, Comb. Probab. Comput. **124**, 290 (2000).

[26] A. Mücke et al., Nucl. Phys. B, Proc. Suppl. **80**, 810 (2000).

[27] M. Cerruti, A. Zech, G. Emergy, and D. Guarin, AIP Conf. Proc. **1792**, 050027 (2017).

[28] Planck Collaboration, Astron. Astrophys. **571**, A15 (2014).

[29] J. H. Davis, J. Silk, C. Bœhm, and K. Kotera, Phys. Rev. D **93**, 103523 (2016).

[30] T. Lacroix, C. Bœhm, and J. Silk, Phys. Rev. D **89**, 063534 (2014).

[31] P. Gondolo and J. Silk, Phys. Rev. Lett. **83**, 1719 (1999).

[32] O. Y. Gnedin and J. R. Primack, Phys. Rev. Lett. **93**, 061302 (2004).

[33] P. Ullio, H. Zhao, and M. Kamionkowski, Phys. Rev. D **64**, 043504 (2001).

[34] D. Merritt and M. Y. Poon, Astrophys. J. **606**, 788 (2004).

[35] C. Bœhm, P. Fayet, and J. Silk, Phys. Rev. D **69**, 101302 (2004).

[36] K. M. Zurek, Phys. Rev. D **79**, 115002 (2009).

[37] C. Bœhm, Phys. Rev. D **62**, 035012 (2000).

[38] L. Goodenough and D. Hooper, arXiv:0910.2998.

[39] J. Petrović, P. D. Serpico, and G. Zaharijas, J. Cosmol. Astropart. Phys. **2** (2015) 023.

[40] A. K. Harding, I. A. Grenier, and P. L. Gonthier, Astrophys. Space Sci. **309**, 221 (2007).

[41] W. Bednarek and T. Sobczak, Mon. Not. R. Astron. Soc. **435**, L14 (2013).

[42] T. Lacroix, C. Bœhm, and J. Silk, Phys. Rev. D **90**, 043508 (2014).

[43] B. S. Acharya et al., Astropart. Phys. **43**, 3 (2013).