Study of the Very High Energy Emission of M87 through its Broadband Spectral Energy Distribution

Gamma rays constitute the highest-energy electromagnetic radiation tracing the most energetic phenomena in the universe. Active galactic nuclei (AGNs) are important sources of extragalactic gamma rays, and according to the current consensus they are powered by accreting supermassive black holes (SMBHs). Most AGNs that are VHE emitters are classified as blazars, i.e., radio-loud AGNs whose jets are pointing nearly toward the observer. Since particles within the jet travel at nearly the speed of light, relativistic beaming increases the brightness of these objects. According to unification schemes (Urry & Padovani 1995), radio galaxies (RDGs) correspond to the misaligned counterparts of blazars, and so far, VHE emission has been detected from six of them (Rieger & Levinson 2018). Since RDGs are located on average closer than blazars, it is possible obtain more detailed observations to test theoretical models of their emission. The broadband spectral energy distributions (SEDs) of AGN jets, which are prominent in the emission of blazars and RDGs, are globally nonthermal and they are characterized by the presence of two components (generally referred to as "peaks") (Blandford et al. 2019). The low-energy component is usually attributed to synchrotron emission, which is produced when relativistic particles are moving in the presence of a magnetic field (Rybicki & Lightman 1979). On the other hand, the second peak has been explained by many different models. These models can be divided into two types: leptonic, where the high-energy component of the broadband SED is explained as inverse Compton (IC) emission produced by an electron population in the jet (Longair 2011), and hadronic, where the second component is produced by mechanisms involving the collision of accelerated protons with the surrounding environment (Mücke et al. 2003).

M87 (R.A. 12^h30^m472, decl. $+{12}^{^\circ }23^{\prime} 51^{\prime\prime}$), which is classified as a giant Fanaroff–Riley I (FR-I) RDG, is the central dominant galaxy of the Virgo Cluster. It is an elliptical galaxy with a diameter of ∼300 kpc (Doherty et al. 2009), a dynamical mass within 180 kpc estimated as (1.5 ± 0.2) × 10¹³ M_⊙ (Zhu et al. 2014), and a redshift of z = 0.0044. It is located at a distance of 16.4 ± 0.5 Mpc, which is a redshift-independent measurement (Bird et al. 2010). M87 hosts an SMBH, named M87*, whose shadow was the first imaged by the Event Horizon Telescope (The EHT Collaboration et al. 2019). The prominent jet of M87 is one its most noticeable characteristics. This jet has been studied for the last one hundred years (Curtis 1918); it has a length of about 2.5 kpc (Biretta et al. 1999) and has been resolved from radio to X-rays. It presents complex structures such as knots and diffuse emission (Perlman et al. 1999, 2001), apparent superluminal motion (Cheung et al. 2007), and a complex variability (Harris et al. 2006).

M87 was the first RDG detected at VHE (Aharonian et al. 2003). It has been detected by different imaging air Cherenkov telescopes (IACTs) such as HESS (Aharonian et al. 2006), VERITAS (Acciari et al. 2010), and MAGIC (Acciari et al. 2020). Recently, the High Altitude Water Cherenkov (HAWC) Collaboration (Albert et al. 2021) reported weak evidence (3.6σ) of long-term VHE emission from this source. M87 has shown a complex behavior at VHE (Ait Benkhali et al. 2019) with a rapid variability during flaring states (Abramowski et al. 2012). According to available observations, three VHE flares from M87 have been detected in 2005 (Aharonian et al. 2006), 2008 (Albert et al. 2008), and 2010 (Abramowski et al. 2012). Gamma-ray angular resolution is not sufficient to determine the region of the galaxy where this emission is produced. Variability studies have suggested that the innermost jet zone of M87 is most likely the source of its VHE emission. The other candidate location is the jet feature HST-1, but it is disfavored by the timescale of the VHE variability and the lack of correlated activity in other bands during TeV flares (Abramowski et al. 2012; Ait Benkhali et al. 2019).

The broadband emission of M87 has usually been explained by a one-zone synchrotron self-Compton (SSC) scenario (Abdo et al. 2009; De Jong et al. 2015). In these models, the first component is attributed to synchrotron radiation that is produced by electrons moving at relativistic velocity with random orientation with respect to the magnetic field. The second peak is explained by inverse Compton scattering of synchrotron photons to higher energies by the same electron population. However, some authors claim that SSC models are not able to explain VHE emission in M87. Evidence for a possible spectral turnover in the GeV regime E ≳ 10 GeV was found by Ait Benkhali et al. (2019). This was interpreted as due to the presence of an additional physical component in the emission. However, constraining the decrease of this component from only Fermi data is not possible and long-term TeV observations are needed. One-zone SSC models were found to have difficulties in explaining VHE/X-ray correlated variability in M87 (Abramowski et al. 2012), and Fraija & Marinelli (2016) claimed that one-zone SSC models cannot be extended to VHE in FR-I RDGs. This is why alternative ideas have been proposed to explain the SED such as seed photons coming from other regions in the jet (Georganopoulos et al. 2005) and photohadronic interactions (Fraija & Marinelli 2016).

The main goal of this work is to compare the VHE emission of the RDG M87 observed by IACTs during specific epochs (including the 2005 flare) with the long-term quiescent/average emission provided by continuous observation by the HAWC observatory from 2014 to 2019. We used a lepto-hadronic model, which combines SSC and photohadronic scenarios, to explain this emission. We developed a Python code to simulate the broadband emission of M87 and constructed an average SED of M87, collecting multifrequency observations from data archives. Preliminary results of this work were released in Ureña-Mena et al. (2021).

We collected historical archive data to construct the broadband SED of M87. Data sets from radio to X-rays were partially based on those observations of the innermost jet zone used by Abdo et al. (2009), Fraija & Marinelli (2016), and Prieto et al. (2016). MeV–GeV gamma-ray data were obtained from the Fermi Large Area Telescope Fourth Source Catalog (4FGL), which is based on the first eight years of data from the Fermi Gamma-ray Space Telescope (Abdollahi et al. 2020). ²⁷ The 4FGL covers the energy range from 50 MeV to 1 TeV. Four different sets of data were used for the TeV range: (1) H.E.S.S. observations from 2004, which were taken during a TeV quiescent phase (Aharonian et al. 2006); (2) H.E.S.S. observations from 2005, which were taken during a TeV high-activity state but without evidence of inner jet activity in the rest of the broadband inner jet spectrum (Aharonian et al. 2006). That is why we use the same broadband SED as in the nonflaring state case; (3) MAGIC-I observations from 2005 to 2007, which correspond to an observation campaign where no flaring activity was detected (Aleksić et al. 2012); (4) HAWC observations, which cover a quiescent period from 2014 to 2019 corresponding to 1523 days (Albert et al. 2021).

The HAWC array consists of 300 water Cherenkov detectors, each with four photomultiplier tubes. The HAWC data are divided into nine bins according to the fraction of channels hit, which are used to estimate the energy of the events (see Abeysekara et al. 2017 for more details). In Albert et al. (2021) a power-law spectrum, with a spectral index set to 2.5, was fit to a sample of 138 nearby AGNs. For those sources with test statistic > 9, including M87, an optimized spectrum with free normalization and spectral index was obtained. Then, quasi-differential flux limits were obtained in three bands: 0.5–2.0 TeV, 2.0–8.0 TeV, and 8.0–32.0 TeV.

We used a hybrid model in this work. Emission components from radio to GeV gamma rays are explained with an SSC scenario whereas photohadronic interactions are added to explain the VHE emission. Therefore, the broadband SED has been modeled with three components. Due to the low redshift of the source, extragalactic background light (EBL) absorption was not relevant for photon energies ≲10 TeV and we did not consider it in the modeling. HAWC data, which are the only data set affected by some relevant EBL absorption, were already corrected for this effect (Albert et al. 2021) using the EBL model of Domínguez et al. (2011). We used the one-zone SSC code described in Finke et al. (2008), which considers a homogeneous spherical region or blob in the inner jet moving with a Lorentz factor Γ and a randomly oriented magnetic field with mean intensity B. The Doppler factor δ is given by

$$\begin{eqnarray}&&\delta ={\left[{\rm{\Gamma }}(1-\beta \mu )\right]}^{-1},\end{eqnarray} \tag{ 1 }$$

where β is the ratio of the speed of the jet to the speed of light and $\mu =\cos \theta$ where θ is the angle of the jet with respect to the observer's line of sight.

The minimum variability timescale is

$$\begin{eqnarray}&&{t}_{v,\min }=\displaystyle \frac{(1+z){R}_{b}^{{\prime} }}{c\delta },\end{eqnarray} \tag{ 2 }$$

where ${R}_{b}^{{\prime} }$ corresponds to the comoving radius of the region, c to the speed of light, δ to the Doppler factor, and z to the redshift of the source. Comoving quantities are primed following the convention used by Finke et al. (2008).

The electron population of the region, which follows an energy distribution ${N}^{{\prime} }({\gamma }^{{\prime} })$, is moving in a randomly oriented magnetic field and producing synchrotron radiation. The electron energy distribution for this model was assumed to be a broken power law given by

$$\begin{eqnarray}&&{N}_{e}({\gamma }^{{\prime} })={K}_{e}\left\{\begin{array}{l}{\left(\displaystyle \frac{{\gamma }^{{\prime} }}{{\gamma }_{\mathrm{break}}^{{\prime} }}\right)}^{-{p}_{1}}\,\mathrm{for}\,{\gamma }_{1}^{{\prime} }\lt {\gamma }^{{\prime} }\lt {\gamma }_{\mathrm{break}}^{{\prime} }\\ {\left(\displaystyle \frac{{\gamma }^{{\prime} }}{{\gamma }_{\mathrm{break}}^{{\prime} }}\right)}^{-{p}_{2}}\,\mathrm{for}\,{\gamma }_{\mathrm{break}}^{{\prime} }\lt {\gamma }^{{\prime} }\lt {\gamma }_{2}^{{\prime} }\end{array}\right.,\end{eqnarray} \tag{ 3 }$$

where ${\gamma }^{{\prime} }$ is the electron Lorentz factor, p₁ and p₂ are the power-law indices, ${\gamma }_{\mathrm{break}}^{{\prime} }$ is the break electron Lorentz factor, K_e is a normalization constant, and ${\gamma }_{1}^{{\prime} }$ and ${\gamma }_{2}^{{\prime} }$ are the minimum and maximum electron Lorentz factors.

The photohadronic component is based on the model presented by Sahu (2019). Recent evidence of neutrino emission from AGNs seems to support the relevance of photohadronic interactions (Aartsen et al. 2018). An accelerated proton population is assumed to be contained in a spherical volume of radius ${R}_{f}^{{\rm{{\prime} }}}$ inside the blob of radius ${R}^{{\rm{{\prime} }}}$ (SSC blob) with ${R}_{f}^{{\rm{{\prime} }}}\lt {R}^{{\rm{{\prime} }}}$. The inner region is also assumed to have a higher seed photon density because the low density in the SSC emission zone makes the photohadronic process inefficient. The proton population has a power-law energy distribution (Fraija & Marinelli 2016; Sahu et al. 2019):

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}{N}_{p}}{{\rm{d}}{E}_{p}}\propto {E}_{p}^{-\alpha },\end{eqnarray} \tag{ 4 }$$

where the spectral index is α > 2.

Due to the higher photon density in this inner volume, protons interact with the background photons through the following mechanism (Dermer & Menon 2009):

$$\begin{eqnarray}&&\begin{array}{c}{\rm{p}}+\gamma \,\to \,{{\rm{\Delta }}}^{+}\\ \to \,\left\{\begin{array}{l}{\rm{p}}+{\pi }^{0}\,\to \,{\rm{p}}+2\gamma \\ {\rm{n}}+{\pi }^{+}\,\to \,{\rm{n}}+{\rm{e}}+3\nu \,\to \,{\rm{p}}+2{\rm{e}}+4\nu .\end{array}\right.\end{array}\end{eqnarray} \tag{ 5 }$$

This process requires the center-of-mass energy of the interaction to exceed the Δ-mass (Sahu 2019; Fraija & Marinelli 2016),

$$\begin{eqnarray}&&{E}_{p}^{{\prime} }{\epsilon }_{\gamma }^{{\prime} }=\displaystyle \frac{({m}_{{\rm{\Delta }}}^{2}-{m}_{p}^{2})}{2(1-{\beta }_{p}\cos \theta )}\cong 0.32\,{{\rm{G}}{\rm{e}}{\rm{V}}}^{2},\end{eqnarray} \tag{ 6 }$$

where E_p is the energy of the proton and _γ is the energy of the target photon. Considering collisions with SSC photons from all directions, β_p ≈ 1, and viewing from the observer frame,

$$\begin{eqnarray}&&{\epsilon }_{{\rm{\Gamma }}}{\epsilon }_{\gamma }\cong 0.32\displaystyle \frac{{\delta }^{2}}{{(1+z)}^{2}}\,{{\rm{G}}{\rm{e}}{\rm{V}}}^{2},\end{eqnarray} \tag{ 7 }$$

where _Γ is the energy of the emitted photon. According to Sahu (2019) the π⁰ decay photon flux is given by

$$\begin{eqnarray}&&{f}^{p\gamma }({\epsilon }_{{\rm{\Gamma }}})={A}_{\gamma }{f}^{{\rm{ICS}}}({\nu }_{\gamma }){\left(\displaystyle \frac{{\epsilon }_{{\rm{\Gamma }}}}{{\rm{TeV}}}\right)}^{-\alpha +3},\end{eqnarray} \tag{ 8 }$$

where ν_γ is the frequency of a photon with energy _γ and f^ICS(ν_γ ) is the flux at ν_γ . Thus, the total emitted flux at VHE energies is given by the sum of the inverse Compton flux and photohadronic flux:

$$\begin{eqnarray}&&{f}^{{\rm{VHE}}}({\nu }_{{\rm{\Gamma }}})={f}^{{\rm{ICS}}}({\nu }_{{\rm{\Gamma }}})+{f}^{p\gamma }({\epsilon }_{{\rm{\Gamma }}}({\nu }_{{\rm{\Gamma }}})),\end{eqnarray} \tag{ 9 }$$

where ν_Γ is the frequency of the emitted photon and _Γ its energy (_Γ = h ν_Γ where h is the Planck constant).

The contribution of proton–proton (pp) interactions would require a very proton-loaded jet to be significant (Reynoso et al. 2011). According to Atoyan & Dermer (2003), for an emission region with Doppler factor δ ≈ 10 and variability timescale t_var ≈ 1 day, which is the case for this source, the addition of the accelerated protons necessary for the pp interactions would imply a super-Eddington jet power of ≈10⁵⁰ erg s⁻¹. This power is six orders of magnitude larger than the estimates for M87, which is why we did not consider pp interactions in this work. Synchrotron emission from protons and muons would need a much higher magnetic field intensity to be relevant. We also neglected synchrotron from pions because of their short lifetime.

We developed a Python code to reproduce the lepto-hadronic model described above. The broadband SED was fit with this emission model. We obtained the best fit values for the physical parameters and estimated their errors using Monte Carlo simulations.

4.1. Fitting Technique: SSC Model

4.2. Fitting Technique: Photohadronic Component

We carried out the SSC model fitting following the procedure described in Section 4.1. The best fit values of the physical parameters are presented in Table 1. The best fit model is plotted together with the observational data and residuals in Figure 1.

Zoom In Zoom Out Reset image size

Table 1. Best Fit Values for the SSC Model Parameters with Estimated Errors

Parameter	Value
Magnetic field B (G)	0.046 ± 0.003
Doppler factor δ	4.3 ± 0.2
Electron energy distribution parameters
Power-law index (lower energies) p1	1.52 ± 0.02
Power-law index (higher energies) p2	3.53 ± 0.02
Break Lorentz factor ${\gamma }_{c}^{{\prime} }$	${3.80}_{-0.05}^{+0.06}\times {10}^{3}$
${\chi }_{\nu }^{2}$(d.o.f)	26.1 (20)

Download table as:  ASCIITypeset image

As mentioned before, we considered four different VHE data sets for modeling the VHE emission of M87. The best fit values for the photohadronic model parameters are presented in Table 2. The best fit models, data, and residuals are plotted in Figure 2.

Zoom In Zoom Out Reset image size

Table 2. Best Fit Values for the Photohadronic Component Fitting Parameters with Their Error Estimates

	α	Aγ	${\chi }_{\nu }^{2}$(d.o.f)
H.E.S.S.: 2004 observations a	${3.2}_{-0.4}^{+0.2}$	${0.1}_{-0.1}^{+0.2}$	25.5 (22)
H.E.S.S.: 2005 observations a	2.8 ± 0.2	${0.6}_{-0.2}^{+0.4}$	22.5 (22)
MAGIC: 2005–2007 observations b	3.0 ± 0.2	${0.2}_{-0.1}^{+0.2}$	23.7 (26)
HAWC observations (1523 days) c	3.1 ± 0.2	0.2 ± 0.1	25.8 (22)

Notes. α, proton energy distribution index; A_γ , normalization constant.

^a Data taken from Aharonian et al. (2006). ^b Data taken from Aleksić et al. (2012). ^c Data taken from Albert et al. (2021).

Download table as:  ASCIITypeset image

We constructed a broadband SED using data from steady-state epochs. It was fit with a lepto-hadronic model in order to explain its VHE emission. The lepto-hadronic model is a combination of the one-zone SSC model proposed by Finke et al. (2008) and the photohadronic model used by Sahu (2019). The model proposes the existence of three components: a synchrotron-dominated component from radio to X-rays, an inverse Compton-dominated component from X-rays to GeV gamma rays, and a photohadronic-dominated component in the range of TeV gamma rays.

As shown in Figure 1, the SSC model provides a good fit to the SED up to GeV gamma rays, but it underestimates the VHE emission, consistently with what is found in the literature (e.g., Ait Benkhali et al. 2019). The X-ray emission was of particular interest for corresponding to the transition between the two components of the SSC scenario. X-rays are also important for producing the VHE emission in the photohadronic scenario, because the typical energies of the target photons for the photohadronic interactions fall in the X-ray range. The only X-ray data point that could not be well fit was the NuSTAR 20–40 KeV observation. In Wong et al. (2017), which reported this observation, they also discussed its inconsistency with inverse Compton models to explain gamma-ray emission. According to them, the observation uncertainties are limited by the statistical power of their data, and deeper observations are necessary to resolve the tension with the NuSTAR data.

Table 3 shows a comparison of the SSC parameter values obtained by different studies. Two of them (the magnetic field intensity and the break Lorentz factor) present a dispersion of several orders of magnitude. This dispersion can be caused by degeneracy in the SED models, different assumptions in the emission zone geometry, or different levels of completeness in the data sets. In the case of magnetic field intensity (B), most of the degeneracy comes from the relation $B{\delta }^{2}{R}_{b}^{{\prime} }\approx$ constant (Finke et al. 2008), which was first derived by Tavecchio et al. (1998). As this relation holds in every case of Table 3 ($B{\delta }^{2}{R}_{b}^{{\prime} }\approx (2-4)\times {10}^{16}$), variations in the magnetic field intensity (B) can be caused by different assumptions regarding the radius of the emission zone ${R}_{b}^{{\prime} }$. The other parameters have less dispersion. It is important to remark that the lack of error estimates for most of the parameters reported in the literature prevents a more precise comparison of these results.

Table 3. Comparison of Results for the SSC Parameters from Different Studies

Parameter	This Study	Abdo et al. (2009)	De Jong et al. (2015)	Fraija & Marinelli (2016)	Acciari et al. (2020)
B(G)	0.046 ± 0.003	0.055	0.002	1.61	0.0031
δ	4.3 ± 0.2	3.9	5	2.8	5.3
p1	1.52 ± 0.02	1.6	−1.8	3.21 ± 0.02	1.9
p2	3.53 ± 0.02	3.6	3.4	4.21	3.2
${\gamma }_{c}^{{\prime} }$	${3.80}_{-0.05}^{+0.06}\times {10}^{3}$	4 × 103	4.0 × 102	1.7 × 103	1.4 × 104
${R}_{b}^{{\prime} }$ (cm)	(1.54 ± 0.07) × 1016	1.4 × 1016	5.6 × 1017	2.1 × 1015	4.0 × 1017

Note. All the SEDs were constructed to model the nonflaring state of M87.

Download table as:  ASCIITypeset image

The viewing angle (θ) plays an important role in AGN properties. According to unification schemes, the transition between blazars and RDGs is produced around θ ∼ 10°. The Doppler factor can be constrained from θ by Abdo et al. (2010):

$$\begin{eqnarray}&&\delta \leqslant \displaystyle \frac{1}{\sin (\theta )}.\end{eqnarray} \tag{ 10 }$$

Higher Doppler factors enhance the HE and VHE emission, which explains why many more blazars have been detected at gamma-ray energies than RDGs.

The Doppler factor obtained in this work (δ = 4.3 ± 0.2) is consistent with the lowest estimates for the viewing angle (θ ≈ 10°–20°) (Biretta et al. 1999), which have been based mainly on optical observations of the jet feature HST-1. However, it is in disagreement with other estimates (θ ≈ 30°–45°) (Ly et al. 2007) based on very long baseline interferometry (VLBI) observations of the jet base of M87. In fact, all the values for δ presented in Table 3 are inconsistent with the VLBI measurements. One way to solve this tension is to take into account the width of the jet base. According to Hada et al. (2016) the jet base of M87 has an apparent opening angle of ∼100°, which would correspond to an intrinsic opening angle of ∼50° (if θ ∼ 30°). Therefore, a VHE emission zone located in the outer zones of the jet base could have a viewing angle as low as ∼5°, which would be consistent with a Doppler factor δ ≤ 11.8. It is important to mention that the Doppler factor can have spatial and temporal variations along the jet (Hada et al. 2016), which are not considered in this modeling. However, the high degeneracy of the multizone models makes it difficult to constrain these variations.

As shown in Figure 2, the photohadronic component is able to explain the VHE emission in both flaring and quiescent states. This component is produced by the interaction between inverse Compton photons and accelerated protons from the jet. We studied the quiescent state using three different sets of TeV data: H.E.S.S. observations from 2004, MAGIC observations from 2005 to 2007, and HAWC observations from the 2014 to 2019 period (all three of them without VHE flares). The high-activity state was studied using H.E.S.S. data corresponding to the 2005 VHE flare. There were no reports of high activity in the inner jet in other bands during that flare. That is why we used the same broadband SED as in the nonflaring state. HAWC results were in agreement with the 2004 H.E.S.S. observations and the MAGIC observations, but HAWC fluxes are lower than the fluxes of the 2005 flaring state. This indicates that HAWC observations constrain the average VHE emission from M87 during quiescent periods.

The proton energy distribution index α was estimated with the four VHE data sets. Those measurements agree with the result obtained by Fraija & Marinelli (2016), α = 2.80 ± 0.02, where a similar lepto-hadronic model was used to fit the 2004 H.E.S.S. data. It is also important to remark the change in α observed in the 2005 H.E.S.S. observation with respect to other TeV results and the lack of high activity in the other bands during this flare. Those results indicate that the flare could have been caused by an increase in energy of the accelerated proton population.

It is important to mention that HAWC data provide the most continuous and uniform TeV measurements of M87. The IACT campaign observations may have a good cadence, but they are made with exposure times corresponding to a few hours that could be affected by rapid VHE flux variations. Therefore, HAWC results represent a steadier constraint to the mean VHE emission of M87 during the HAWC's period of observation.

With these results, we calculated an electron luminosity of L_e ∼ 7 × 10⁴² erg s⁻¹. Accelerated proton flux (F_p ) and luminosity (L_p ) can also be estimated using the following equation (Sahu 2019):

$$\begin{eqnarray}&&{F}_{p}({E}_{p})\approx 10\times \displaystyle \frac{{f}^{p\gamma }({\epsilon }_{{\rm{\Gamma }}})}{{\tau }_{p\gamma }({E}_{p})},\end{eqnarray} \tag{ 11 }$$

where f^{p γ} (_Γ) is given by Equation (8) and τ_{p γ} is the optical depth of the Δ resonance process in the inner jet region. As mentioned before, this model assumes the photohadronic interactions to occur in an inner compact region of the blob with a smaller size and a higher photon density. Unfortunately, these quantities are not directly observable and the value of τ_{p γ} cannot be calculated. However, Sahu (2019) gives two prescriptions, τ_{p γ} < 2 and τ_{p γ} > f^{p γ} (_Γ)/f_Edd, where f_Edd is the Eddington flux. As in this case f^{p γ} (_Γ)/f_Edd ≈ 10⁻⁷, we assumed the intermediate value τ_{p γ} ≈ 10⁻² and we obtained a proton luminosity of L_p ∼ 6 × 10⁴³ erg s⁻¹. These results are in agreement with the estimates of total jet power L_j ∼ 10⁴⁴ erg s⁻¹ (Owen et al. 2000).

The decay of charged pions produces neutrinos. The neutrino flux (f^ν ) can be estimated assuming (Aartsen et al. 2020)

$$\begin{eqnarray}&&{E}_{\nu }\approx {\epsilon }_{{\rm{\Gamma }}}/2,\end{eqnarray} \tag{ 12 }$$

where E_ν is the emitted neutrino energy, and also (Murase et al. 2016)

$$\begin{eqnarray}&&{f}^{\nu }({E}_{\nu })\approx \displaystyle \frac{3}{4}{f}^{p\gamma }({\epsilon }_{{\rm{\Gamma }}}/2).\end{eqnarray} \tag{ 13 }$$

The estimated neutrino flux for E_ν = 5 TeV corresponds to _Γ = 10 TeV and is f^ν ∼ 1 × 10⁻¹³ TeV cm⁻² s⁻¹, which is below IceCube upper limits (Aartsen et al. 2020).

With regard to the results of some other alternative models, in Fraija & Marinelli (2016) the SED of M87 was fit with a very similar lepto-hadronic model. In this case, VHE emission was represented only by the 2004 H.E.S.S. results, which were obtained with just ∼50 hr of observation. However, the best fitting values of the photohadronic parameters were in agreement, within their uncertainties, with those obtained in this work. In Sahu & Palacios (2015) a lepto-hadronic scenario was used to explain the 2010 TeV flare, and their results were also in agreement with those obtained in this work. However, this flare had an X-ray counterpart (Abramowski et al. 2012), which may be interpreted as purely leptonic. A leptonic and a hybrid model were used to model MAGIC results from a 2012–2015 campaign (Acciari et al. 2020), where the hybrid model was found to be more consistent with gamma-ray data. Finally, according to Ait Benkhali et al. (2019), extended gamma-ray production scenarios such as Compton scattering in the kiloparsec-scale jet (e.g., Hardcastle & Croston 2011) are disfavored by gamma-ray variability.

We cannot rule a purely leptonic scenario out. The reason is that the single-zone model used here does not take into account how emissions from different source regions are related. Actually, multizone structured leptonic models are necessary to explain specific features in the whole multiwavelength SED (Algaba et al. 2021). However, these models have a large number of parameters, which introduces a high degeneracy, making it difficult to derive firm conclusions from their fit. As the aim of this paper is to explain the VHE emission, and we are therefore limited by the quality of these observations, we decided to use a one-zone SSC scenario to model the leptonic contribution. However, we are confident that other physical models can be applied to future releases of HAWC data that will provide higher signal-to-noise ratio.

M87 is considered a laboratory for understanding AGN properties since it is the only source that has been mapped from the SMBH shadow (∼0.005 pc) to the outer jet (∼25 kpc). Understanding its gamma-ray emission, which is practically unaffected by EBL absorption up to ∼10 TeV, could be the key to understanding gamma-ray emission in the rest of the AGNs (not only RDGs). We fit a broadband SED of M87 with a lepto-hadronic model with the aim of explaining its VHE emission. Emission from radio to GeV gamma rays has been modeled with an SSC scenario. The best fit values for SSC model parameters were δ = 4.3 ± 0.2 for the Doppler factor, B = 0.046 ± 0.003 G for the mean magnetic field intensity, and p₁ = 1.52 ± 0.02, p₂ = 3.53 ± 0.02, and ${\gamma }_{c}^{{\prime} }=$ ${3.80}_{-0.5}^{+0.6}\times {10}^{3}$ for the electron energy distribution parameters. The value of the Doppler factor is in agreement with a low viewing angle of the jet base (θ ∼ 13°). However, a high viewing angle is also possible if the opening angle of the jet base is wide enough to place the emission zone closer to the observer's line of sight.

A photohadronic model was fit to the VHE emission. Results show that this model is able to explain the quiescent VHE emission represented by H.E.S.S., MAGIC, and HAWC observations. H.E.S.S. data corresponding to the 2005 VHE flare were also fit using this model. The results show that the model can explain the so-called orphan flares, which are only detected at VHE bands, such as the flare observed in 2005. If the value of the proton spectral index decreases, a harder VHE spectrum would be obtained, therefore the resulting VHE spectrum would take the shape described by the equation ${E}^{2}({dN}/{dE})\approx {E}^{-\alpha +2.6}$, where α is the proton spectral index. Hence, changes in the proton energy distribution would produce those flares.

HAWC observations constrained the VHE emission from M87 for the 2014 to 2019 period in which no evidence of VHE flares was reported. We obtained a power-law index for the proton energy distribution of α = 3.1 ± 0.2 and a TeV gamma-ray flux normalization constant of A_γ = 0.2 ± 0.1. HAWC will be taking data for a few more years. Therefore, the significance of the M87 detection will probably be improved, allowing a better estimation of the photohadronic model parameters.

We acknowledge the support from: the US National Science Foundation (NSF); the US Department of Energy Office of High-Energy Physics; the Laboratory Directed Research and Development (LDRD) program of Los Alamos National Laboratory; Consejo Nacional de Ciencia y Tecnología (CONACyT), Mexico, grants 271051, 232656, 260378, 179588, 254964, 258865, 243290, 132197, A1-S-46288, A1-S-22784, cátedras 873, 1563, 341, 323, Red HAWC, Mexico; DGAPA-UNAM grants IG101320, IN111716-3, IN111419, IA102019, IN110621, IN110521; VIEP-BUAP; PIFI 2012, 2013, PROFOCIE 2014, 2015; the University of Wisconsin Alumni Research Foundation; the Institute of Geophysics, Planetary Physics, and Signatures at Los Alamos National Laboratory; Polish Science Centre grant, DEC-2017/27/B/ST9/02272; Coordinación de la Investigación Científica de la Universidad Michoacana; Royal Society–Newton Advanced Fellowship 180385; Generalitat Valenciana, grant CIDEGENT/2018/034; Chulalongkorn University's CUniverse (CUAASC) grant; Coordinación General Académica e Innovación (CGAI-UdeG), PRODEP-SEP UDG-CA-499; Institute of Cosmic Ray Research (ICRR), University of Tokyo, H.F. acknowledges support by NASA under award number 80GSFC21M0002. We also acknowledge the significant contributions over many years of Stefan Westerhoff, Gaurang Yodh, and Arnulfo Zepeda Dominguez, all deceased members of the HAWC collaboration. Thanks to Scott Delay, Luciano Díaz, and Eduardo Murrieta for technical support.

Software: Python (Van Rossum & Drake 1995), NumPy (Harris et al. 2020), Astropy (Astropy Collaboration et al. 2018), Matplotlib (Hunter 2007).