© ESO, 2015

1. Introduction

Low-luminosity active galactic nuclei (AGNs) are probably powered by optically thin, hot accretion flows (for recent reviews see, e.g., Poutanen & Veledina 2014; Yuan & Narayan 2014). In the central parts of these flows, protons have energies above the threshold for pion production, and then the flows produce considerable γ-ray fluxes that are due to the decay of neutral pions. This property was noted early in the development of the accretion theory, for instance, by Shapiro et al. (1976), and was reported for the currently applied advection-dominated accretion flow (ADAF) models by Mahadevan et al. (1997). As estimated in Shapiro et al. (1976) and Oka & Manmoto (2003), the γ-ray luminosity, L_γ, of a flow surrounding a rapidly rotating black hole may be similar to its X-ray luminosity, L_X. However, these studies neglected the general relativistic (GR) transfer of γ-ray photons as well as the absorption of these photons in the radiation field of the flow. Improved computations by Niedźwiecki et al. (2013), hereafter N13, showed that L_γ is at least an order of magnitude lower than L_X if these effects are taken into account. The precise γ-ray emission model, however, has not been directly compared with current observational data.

No prominent γ-ray signal was found based on observations in the Fermi/LAT data from radio-quiet AGNs (e.g., Ackermann et al. 2012a, 2015), with a few notable exceptions. These include the three X-ray brightest Seyfert 2 galaxies, NGC 4945, NGC 1068, and Circinus (Lenain et al. 2010; Ackermann et al. 2012b; Hayashida et al. 2013). However, all three are composite AGN plus starburst systems, and it has not yet been definitively established whether their γ-ray emission is related with the starburst or with the AGN activity. The upper limits (UL) on the photon flux in the GeV range, derived in Ackermann et al. (2012a), typically constrain the luminosity ratio to L_γ/L_X< 0.1 and in several Seyferts to L_γ/L_X< 0.01, which shows that the sensitivity of the Fermi/LAT surveys has reached the level at which predictions of the hadronic emission from hot accretion flows can be probed. This motivated the detailed comparison we present here.

We analyzed 6.4 years of Fermi-LAT data of nearby, low-luminosity AGNs, which means that we use a data set that is more than three years longer than the one used by Ackermann et al. (2012a) and is more than two years longer than the recent third catalog of AGNs detected by Fermi/LAT (Ackermann et al. 2015). A detailed comparison with the model requires a precisely determined black hole mass and intrinsic X-ray luminosity, therefore we focus here on several of the best-studied objects. As a particularly interesting case, we thoroughly examine the data from NGC 4151, for which the ratio of L_γ/L_X< 0.0025, found in Ackermann et al. (2012a), is the lowest of the 120 Seyfert galaxies considered in their work.

The physics of particle acceleration and heating in hot flows is not well understood. If non-thermal processes (e.g., magnetic reconnection) take place, a power-law-like distribution of particles may be expected. We here take into account predictions of hot flow models with proton distributions that include thermal and nonthermal components. We compute the γ-ray emission from accretion flows following the approach of N13. We use a GR hydrodynamic description of the flow combined with GR, Monte Carlo (MC) computations of radiative processes. We improve the model of N13 by refining the hydrodynamic solutions for flows whose proton distribution is dominated by a nonthermal component. We obtain solutions for a range of accretion rates that allows comparing our results with objects that have bolometric luminosities of between ~10^-4 and 10^-2 of the Eddington limit. In Sect. 2 we discuss the dependence of the predicted L_γ on various parameters of the hot flow model, this has not been fully explored in precise computations before. In particular, N13 only considered a weak magnetic field and a fixed accretion rate. In Sect. 3 we present the results of our analysis of the Fermi/LAT data together with some results collected from the literature. Our main results are presented in Figs. 1 and 3. In Sect. 4 we briefly consider effects that are related with the nonthermal acceleration of electrons, in particular, their Compton emission in γ-rays. In Sect. 5 we discuss the general implications of our results for the physics of hot flows; we also briefly discuss hard spectral states of black hole binaries, where similar hot flows are most likely present.

2. Model of X- and γ-ray emission from hot flows

We applied the model developed in Niedźwiecki et al. (2012, 2013, 2014, 2015) and Xie et al. (2010). Here we briefly summarize its basic properties and some improvements we added here.

We considered a black hole that is characterized by its mass, M, and angular momentum, J, surrounded by a geometrically thick accretion flow with an accretion rate, Ṁ. We defined the following dimensionless parameters: r = R/R_g, a = J/ (cR_gM), ṁ = Ṁ/Ṁ_Edd, where Ṁ_Edd = L_Edd/c², R_g = GM/c² is the gravitational radius and L_Edd ≡ 4πGMm_pc/σ_T is the Eddington luminosity. We assumed a viscosity parameter, α = 0.3 and the ratio of the gas pressure (electron and ion) to the magnetic pressure, β.

The fraction of the dissipated energy that directly heats electrons is denoted by δ. We considered two limiting cases, one with δ = 10^-3, for which electron heating is dominated by Coulomb interactions at higher ṁ and compression at lower ṁ, and the other with δ = 0.5, for which electron heating is dominated by magnetohydrodynamic (MHD) processes. Our results for δ = 10^-3 are valid for all values of δ ≲ 0.1.

The Eddington ratio for a given energy range is denoted by λ with a relevant subscript, for example, λ_2−10keV = L_2−10keV/L_Edd and λ_1−10GeV = L_1−10GeV/L_Edd, where L_2−10keV and L_1−10GeV are the luminosities in the 2–10 keV and 1–10 GeV range detected far away from the flow, that is, they are affected by GR transfer and γγ absorption effects (only the dashed line in Fig. 1b shows the rest-frame luminosity). Similarly, the photon flux integrated over a given energy range is denoted by F with a relevant subscript.

Fig. 1 a) X-ray photon spectral index as a function of intrinsic λ2−10keV. The observational data correspond to the best fits of the high-quality observations described in Sect. 3. The gray regions indicate the location of solutions for the model with seed photons from thermal synchrotron only lower and including nonthermal synchrotron (of pion-decay or directly accelerated electrons; upper region). Ticks in the top axis show the approximate value of rtr (see text). The thick vertical line indicates the maximum luminosity for the standard ADAF solutions (i.e., unaffected by the Coulomb cooling of protons). b) λ1−10GeV as a function of λ2−10keV, see Tables 1 and 2 for the observational data. The solid lines are for model N with δ = 10-3, a = 0.95, β = 9 (triangles), model N with δ = 0.5, a = 0, β = 9 is plotted with circles), model H0.1 with δ = 10-3, a = 0.95, β = 9 with squares, and all three are the same as in Fig. 2b below. The dashed line with open squares denotes the rest-frame 1–10GeV luminosity in the last model (i.e., H0.1). In panel b), λ2−10keV corresponds to the average X-ray luminosity during the Fermi/LAT observations (see text). This differs in some objects by up to a factor of 2 relative to panel a). The point for Cyg X-1 shows the hard state parameters according to Malyshev et al. (2013).

We found the global hydrodynamical solutions of the GR structure equations following Manmoto (2000) with minor improvements as described in Niedźwiecki et al. (2012). From this, we found the self-consistent electron temperature distribution, T_e(r), for which the electron energy balance (with global Compton cooling) is achieved. This involved several iterations between the solutions of the electron energy equation and the MC Comptonization simulations. The procedure involves the assumption that the flow structure is unaffected by changes in T_e, which constrains the range of exact solutions to λ_2−10keV ≲ (1−2) × 10^-3 (depending on the model, see also Yang et al. 2015); the maximum λ_2−10keV is denoted by λ_ADAF,max and is indicated by the thick vertical line in Fig. 1a. The Coulomb cooling of protons becomes strong at higher luminosities, and the flow is characterized by a dramatic dependence on even small changes of T_e (when global Compton cooling is taken into account, cf. Xie et al. 2010), which should lead to a collapse of the flow and formation of a standard cold disk at r>r_tr, with the transition radius, r_tr, dependent on ṁ.

All model results we present here were obtained for flows extending out to r = 10⁴. We currently investigate the regime of λ_2−10keV>λ_ADAF,max. Although the results are rather preliminary, we find that (1) up to λ_2−10keV ≃ 3 × 10^-3 (indicated by the boundary of gray area in Fig. 1a) the hot flow can occupy the inner region within r_tr ≃ 100 (this accretion flow geometry is supported by observations, as noted below); (2) up to this λ_2−10keV, the X-ray and the rest-frame γ-ray emission is correctly approximated by models with flows extending to large r, for which we used the hydrodynamical solution neglecting the Coulomb cooling of protons; and (3) this model (with a flow extending to large r) may overestimate the γ-ray opacity (see below).

Fig. 2 a) Photosphere radius for photons with E = 0.1, 1, and 10 GeV from bottom to top as a function of λ2−10keV in the model with δ = 10-3, a = 0.95, β = 1 (in other models the dependence of rph on λ2−10keV is very similar). b) λ1−10GeV as a function of λ2−10keV in models with ηp> 0 and s = 2.1; the dotted lines plot models N with δ = 10-3, (magenta) circles indicate models N with a = 0.95 and β = 9 with ṁ = 0.05, 0.1, and 0.3, (green) squares stand for a = 0, β = 9, ṁ = 0.1, 0.3 and 0.5; the dashed lines show models H0.1 with δ = 10-3, (black) circles models with a = 0.95, β = 1 with ṁ = 0.033, 0.1, 0.3 and 0.5; (blue) squares indicate a = 0.95, β = 1 with ṁ = 0.1, 0.3, 0.5, 0.8 and 1.2; the solid lines show models N with δ = 0.5, (red) squares depict a = 0, β = 1 with ṁ = 0.004, 0.01, 0.033, 0.1 and 0.3; (orange) circles stand for a = 0.95, β = 9 with ṁ = 0.003, 0.01, 0.033 and 0.1. c) The high-energy part of the observed spectra of leptonic emission. (Black) solid lines indicate model H0.1 with δ = 10-3, a = 0.95, β = 1 and ṁ = 0.033, 0.1, 0.3; (red) dashed lines show model N with δ = 0.5, a = 0, β = 1 and ṁ = 0.004, 0.01, 0.033, 0.1; (blue) dotted line shows model N with δ = 0.5, a = 0.95, β = 9 and ṁ = 0.01.

The minimum luminosity of flows that can be studied with our computational method, determining the lowest λ_2−10keV (denoted by λ_min) in all presented models, is set by the requirement that the energy balance for electrons is determined by radiative cooling (rather than advection). Close to λ_min, the compressive heating of electrons dominates for δ = 10^-3, yielding an approximately constant radiative efficiency. This effect is crucial for assessing the γ-ray emission in objects with λ_2−10keV ≪ 10^-4, and we return to this below.

The magnitude of the hadronic processes involves a significant uncertainty related to the distribution of proton energies. At λ_2−10keV<λ_ADAF,max the thermalization and cooling timescales for protons are much longer than the accretion timescale, and the distribution of protons is determined by poorly understood heating and acceleration processes. We considered three cases:

• a thermal model T in which all protons have a Maxwelliandistribution;

• a nonthermal model N, where the total power is used to accelerate a small fraction of protons, for which we assume a power-law distribution, n_pl(γ) ∝ γ^{− s} up to γ_max = 100, and we assume that the other protons remain cold; the radius-dependent fraction of the protons with the power-law distribution is given by Eq. (7) in N13;

• and a hybrid model H_0.1 that has 10% of the energy content in nonthermal protons, denoted by η_p, and 90% of the energy content in thermal protons; for the nonthermal protons we assume a power-law distribution like that in model N.

The assumptions of these models are identical to models T, N, and H in N13, except for the energy content of nonthermal component in model H, which was radius-dependent in N13, whereas here we assumed a radius-independent η_p = 0.1 (indicated by the subscript in H_0.1). For model N, however, we here take into account the effect neglected in N13, namely, a higher compressibility of the flow supported by the pressure of a nonthermal, relativistic proton gas, corresponding to the decrease of the adiabatic index from ≃5/3 (for a thermal proton gas with typically subrelativistic temperatures) to 4/3. All results for model N in this work use the hydrodynamic solutions including this effect. The higher compressibility of the flow results in its smaller geometrical thickness and, hence, higher density (see also Fig. 4 in Popham & Gammie 1998). For a given set of parameters, the density is higher by a factor of ~1.5 (for β = 1) to ~5 (for β = 9) in model N than in models T or H_0.1.

We note that the dependence of the geometrical structure on β, which was discussed in Niedźwiecki et al. (2014), is almost completely reduced in model N. The results presented for this model, with both components significantly contributing to the total pressure (i.e., protons and magnetic field), described by the adiabatic index 4/3, correspond to the maximum density, and hence the maximum Thomson depth, τ, as well as the maximum L_γ of the flow for a given set of ṁ, a, β, and δ.

To compute the X-ray radiation, we took into account seed photons from nonthermal synchrotron emission, which is crucial to reach agreement for the model spectra with the observational data. As first noted by Yuan & Zdziarski (2004), ADAF models that follow the original formulation by Narayan & Yi (1995), that is, whose seed photons are provided by thermal synchrotron emission, predict too high T_e and disagree with observational data for AGNs. Following the results of Niedźwiecki et al. (2015), we self-consistently included the synchrotron emission of pion-decay electrons for δ = 10^-3, whereas for δ = 0.5 we assumed that a fraction η_e of the electron heating power was used for their nonthermal acceleration. In the latter case, we assumed η_e = 0.1 and a monoenergetic injection of electrons at γ = 100, for which synchrotron self-absorption is relatively unimportant (in Sect. 4.2 we release this assumption and consider the dependence on the electron acceleration distribution). The two gray regions are delineated in Fig. 1a to indicate the location of model solutions that include a significant nonthermal synchrotron radiation (upper region) and solutions with only thermal synchrotron emission (lower region), see Sect. 5 for a more detailed discussion.

Using our MC Comptonization code, we tabulated the angular-, energy-, and location-dependent distribution of photons propagating in the central region. Then, to compute L_γ, we took into account the absorption of escaping γ-ray photons on pair creation in interactions with the (tabulated) target photons, strictly following N13. In Fig. 2a we show the size, r_ph, of the γ-ray photosphere (inside which the flow is opaque to γ-rays), which is determined as the radius of the emission point in the equatorial plane, from which the optical depth along the outward radial trajectory is τ_γγ = 1. For radial directions close to the symmetry axis, values of r_ph determined by a similar condition are lower by a factor of ~1.5. Results presented in this work, except for the dashed line in Fig. 3, correspond to spectra averaged over the observation angle, θ_obs. A strong dependence on θ_obs only occurs for a rapidly rotating black hole and only at low λ_2−10keV, at which the contribution from the innermost few R_g is not attenuated by γγ absorption. In other cases, the difference between γ-ray fluxes received at different θ_obs does not exceed a factor of ~2. The amount of attenuation due to γγ absorption is illustrated by the dashed (rest-frame L_γ) and solid (observed L_γ) lines with squares in Fig. 1b; we note that at high λ_2−10keV the internal L_γ is higher by two orders of magnitude than observed.

The dependence of r_ph on λ_2−10keV is similar in all models (i.e., for different a, β, or δ). However, the scaling presented in Fig. 2a cannot be simply extrapolated to λ_2−10keV ≳ 10^-3, at which this extrapolated r_ph would exceed the flow size of r_tr ≃ 100. In this case, the γ-ray opacity should be comparatively lower, at least in outer parts of the flow (i.e., close to r_tr). In a flow extending to large r, the opacity for a γ-ray photon escaping from r_em is dominated by head-head interactions with X-ray photons produced at r>r_em. For photons escaping from the outer parts of a flow contained within r_tr, the opacity is dominated by much less efficient head-tail interactions with X-ray photons produced at r<r_em.

Below we summarize the dependence of L_γ on various parameters that is relevant for the comparison presented in Sect. 3. For thermal protons, the bulk of the hadronic radiation is emitted in the 0.1–1 GeV range. The presence of nonthermal protons is reflected in the production of photons with E> 1 GeV, with a power-law spectrum for a power-law proton distribution. We consider these two regimes separately.

Fig. 3 λ0.2−1GeV as a function of λ2−10keV, see Tables 1 and 2 for the observational data. All model points are for thermal protons (model T) with δ = 10-3; the solid line with circles shows a = 0.95, β = 1, ṁ = 0.033, 0.1, 0.3, and 0.5; the solid line with full squares depicts a = 0.95, β = 9, ṁ = 0.1, 0.3, 0.5, 0.8, and 1.2; the dashed line with open squares plots the estimate (see text) for an edge-on observer for a = 0.998, β = 9, ṁ = 0.033, and 0.1.

2.1. Thermal protons: λ_0.2−1GeV vs. δ, β, and a

The rate of pion production depends on the number of protons with energies above the pion production threshold. For the thermal distribution of protons, this rate is extremely sensitive to the proton temperature, T_p. The pion production rate is low for T_p ≲ 10¹² K and negligible for T_p ≲ 4 × 10¹¹ K. The highest T_p, and accordingly the highest L_γ, can be expected for flows with large β and a and small δ. The maximum λ_0.2−1GeV that can be observed for thermal protons, with λ_0.2−1GeV ≃ 10^-5 for λ_2−10keV between ~10^-4 and 10^-3, is shown by the line with squares in Fig. 3. In models with a = 0, δ = 0.5 or β = 1, λ_0.2−1GeV is smaller at least by a factor of several, and even by orders of magnitude at λ_2−10keV ≃ 10^-3; in all these cases, the model predictions are below the Fermi ULs discussed in Sect. 3.

The dependence on these three parameters involves the following effects. In flows with smaller β, a larger part of the accretion power is used to build up the magnetic field strength, therefore, the proton heating power is weaker. As a result, for β = 1 the proton temperature is by a factor of ~2 lower than for β = 9. The related decrease of λ_0.2−1GeV is illustrated in Fig. 3; for β = 1 and a = 0.95, T_p only exceeds 4 × 10¹¹ K at r< 6, hence, at λ_2−10keV> 10^-4 the γ-ray emission is strongly attenuated by γγ absorption. The increase of δ yields the increase of electron heating at the expense of the proton heating, and it therefore obviously results in the decrease of the L_γ/L_X ratio. For the decreasing a, a decreasing fraction of the accreted rest-energy is dissipated, therefore, T_p decreases. In the extreme case of a = 0 and β = 1, T_p weakly exceeds 4 × 10¹¹ K only close to the event horizon, and no hadronic emission can be expected.

At λ_2−10keV ≃ 10^-5 photons with E ≲ 1 GeV may escape from the ergosphere (the extremely relativistic region located within r< 2) and some specific effects of the Kerr metric can be observed. These are that gravitational focusing of the photon trajectories toward the equatorial plane enhances the flux received by edge-on observers. Furthermore, the increase of a from 0.95 to 0.998 increases L_γ by a factor of ~2; at r> 2 both values of a give the same γ-ray production rate. Taking these effects into account, we show the estimated maximum λ_0.2−1GeV that can be expected at λ_2−10keV ~ 10^-5 by the dashed line in Fig. 3. For β = 9 this λ_2−10keV is lower than λ_min, meaning that it corresponds to the regime where advection of electron energies dominates. This in turn means that the lowest λ_2−10keV point (for a = 0.998 and ṁ = 0.033) is not based on a precise solution of our model; we used the scaling of radiative efficiency, ∝ṁ^0.7, found in Xie & Yuan (2012) for low ṁ, to roughly estimate the X-ray luminosity.

2.2. Nonthermal protons: λ_1−10GeV vs. δ and η_p

The dependence of λ_1−10GeV on λ_2−10keV in models with η_p> 0 is shown in Fig. 2b. Here the number of protons above the pion production threshold depends only linearly on the proton heating power. This means that L_γ is much less sensitive to both β and a than in model T. The difference between various models results mostly from different radiative efficiencies (of electrons), which establishes the relation between λ_2−10keV and ṁ, and accordingly the rest-frame L_γ ∝ ṁ². The efficiency is determined primarily by the value of δ. Obviously, λ_1−10GeV also depends on η_p. However, we note that the difference between models H_0.1 and N does not simply involve the difference of λ_1−10GeV by a factor of 10, corresponding to an increase of η_p from 0.1 to 1. The difference is more complex because of the difference between the hydrodynamical solutions that underlie these models.

For δ = 10^-3, when the heating of electrons is dominated by Coulomb interactions (at larger ṁ), the rest-frame γ-ray luminosity scales roughly linearly with the X-ray luminosity (both ∝ṁ²). For this scaling, combined with the increase of r_ph with increasing ṁ, the observed λ_1−10GeV is roughly constant while λ_2−10keV changes by more than an order of magnitude (see dotted and dashes lines in Fig. 2b). At low ṁ, when the compressive heating of electrons dominates, leading to λ_2−10keV ∝ ṁ, we note a fast decrease of λ_1−10GeV with decreasing λ_2−10keV; this is most clearly seen for model H_0.1 with β = 9 (dashed line with squares in Fig. 2b), in which compressive heating dominates for ṁ ≲ 0.3, corresponding to λ_2−10keV ≲ 10^-4.

In models with δ = 0.5, the radiative efficiencies are higher, therefore, λ_1−10GeV is lower (because of the lower ṁ) at a given λ_2−10keV than for δ = 10^-3. Hence, for δ = 0.5 we only present results for model N; for H_0.1 the predicted γ-ray flux is below the sensitivity of LAT. For a constant δ, λ_2−10keV ∝ ṁ results in a fast decrease of λ_1−10GeV with decreasing λ_2−10keV. Furthermore, for large δ, the radiative efficiency increases with a, as a result, λ_1−10GeV at a given λ_2−10keV is by a factor of several larger for a = 0 than for a = 0.95 (somewhat counter-intuitively).

The uncertain value of δ introduces the main uncertainty in testing the hadronic emission model by allowing for a range of ṁ at a given λ_2−10keV. However, different ṁ also correspond to different Thomson depths of the flow, which is reflected in the shape of the X-ray spectrum. The above uncertainty can therefore be reduced by using the X-ray spectral information. Figure 2c shows example X-ray spectra for large and small δ. Spectra of large-δ models in general deviate from a power-law shape at λ_2−10keV ≲ 10^-4, that is, for ṁ ≲ 0.01 with a = 0.95 and ṁ ≲ 0.03 with a = 0. For low-δ models, deviations from a power-law shape are weak down to λ_2−10keV ~10^-5 (the bottom solid line in Fig. 2c).

We considered various proton distribution slopes between s = 2.1 and 2.7. In all models, the lowest λ_1−10GeV corresponds to s = 2.1 (for s = 2.6, λ_1−10GeV is larger by up to a factor 2), and conservatively we assumed this s for Fig. 1b.

3. Sample and data analysis

We considered nearby AGNs with spectral properties consistent with the hot flow scenario. This includes several Seyfert 1 galaxies: NGC 4151, NGC 5548, IC 4329a, NGC 6814, NGC 4258, and NGC 7213. We also took into account the FR I galaxy Centaurus A, whose X-ray radiation may be dominated by a Seyfert-like emission. Spectral parameters of the intrinsic X-ray emission, measured in all these objects, agree with our model of thermal Comptonization in hot flows (see, e.g., Fig. 1a). Furthermore, they show direct evidence of the lack of an optically thick disk in the central region. For objects with the lowest luminosities (NGC 4258 and NGC 7213) studies of the Fe Kα line place the inner edge of optically thick disk between r_tr ~ 10³ and ~10⁴ (Reynolds et al. 2009; Lobban et al. 2010). In more luminous objects (NGC 4151, NGC 5548, IC 4329a, and NGC 6814), the measured Fe Kα line widths, ≲100 eV (Lubiński et al. 2010; Walton et al. 2013; Brenneman et al. 2012, 2014), indicate truncation of the disk at r_tr ≳ 100; a similar location of an optically thick disk is indicated by reverberation measurements (Edelson et al. 2015). These accretion geometries are consistent with the model of spectral evolution developed for black hole binaries, for instance, by Esin et al. (1997), with r_tr changing in response to changes in ṁ (for a review see Done et al. 2007).

We also considered three Seyfert 2 galaxies, NGC 1068, Circinus, and NGC 4945. All three are observed to have high bolometric luminosities, above the range that can be precisely studied with the current version of our model. Nevertheless, comparing their γ-ray luminosities with properties of their active nuclei appears to be interesting.

For most AGNs in our sample, direct measurements of the masses of their supermassive black holes are available, only for IC 4329a and NGC 7213 did we use estimates based on the velocity dispersion. The adopted M and distance values are given in Table 1. For NGC 5548 and IC 4329a we assumed the luminosity distances for the following cosmological parameters: Ω_m = 0.27, Ω_Λ = 0.73, and H₀ = 71 km s^-1 Mpc^-1. The error bars on the values of λ in Figs. 1 and 3 are mostly due to uncertainties in the determination of M. Of the measurements of M that are based on H₂0 megamaser kinematics, the uncertainty is reported only for Circinus; for the remaining cases (NGC 1068, NGC 4258, and NGC 4945) we assumed the same relative uncertainties in M for Figs. 1 and 3. For NGC 6814 the difference of measurements of M between Bentz et al. (2009) and Pancoast et al. (2015) amount to an order of magnitude, so we present results for both values; the case with the higher value of M is denoted by a star superscript in Tables 1 and 2 and Fig. 1.

Table 2 shows the results of our analysis of the Fermi/LAT data, described below, except for NGC 1068, Circinus, NGC 4945, and Centaurus A, for which we used the results of previous studies collected from literature. In Figs. 1b and 3 we show the measured λ_0.2−1GeV and λ_1−10GeV, or UL on them, as a function of the intrinsic λ_2−10keV. To find the latter, we used results of detailed X-ray spectral studies (specified below) of high-quality data from Suzaku, NuStar, or simultaneous XMM and INTEGRAL observations, which allow distinguishing the primary X-ray emission and reflection or absorption components. The parameters of these fits are shown in Fig. 1a. Using the Swift/BAT (Krimm et al. 2013) light curves¹, we find that even for the most variable objects (NGC 7213, IC 4329a, and NGC 6814; their light curves are shown in Fig. 4) the adopted X-ray observations correspond to the flux levels close to the average Swift/BAT flux during the analyzed 6.4 years of Fermi/LAT observations; in some objects we used the intrinsic λ_2−10keV rescaled to account for a small difference, see Table 1.

Fig. 4 Small black circles with error bars show Swift/BAT 15–50 keV light curves in seven-day bins for NGC 7213, NGC 4151, IC 4329a, and NGC 6814. The red vertical line in each panel shows the average flux during the 6.4 years of analyzed Fermi/LAT observations. The large orange circles in panels (a) c) d)) show the times of the adopted Suzaku and NuStar observations and the average flux in time bins equal to the duration of these pointing observations. The red vertical line in panel b) indicates the end of the three-year period used for comparison with Ackermann et al. (2012a). The blue line in the bottom of panel b) delineates period S (see Sect. 3.2). The flux levels of dim (D), medium (M), and bright (B) states, defined in Lubiński et al. (2010), are shown in the left part of panel b).

Fig. 5 TS maps for the region around NGC 4151; in all panels the green cross shows the location of NGC 4151 and the orange cross the location of source S. In all panels only the 3FGL sources were subtracted from the maps, except for panel b), where an additional point source was added at the position of NGC 4151. a) The 8° × 8° degrees region with a pixel size of . The map was built for the energy range 0.2–30 GeV, neglecting period S (see text); the red and white crosses indicate the location of 3FGL sources. b) The same as in a), but the model includes a point source at the position of NGC 4151 to compensate for the residual seen in a). c) The degree region with a pixel size of 0.̊075 for 0.2–30 GeV for period S. d) The region with a pixel size of 0.̊075 for 0.2–6 GeV and the full period of 6.4 years; the black cross shows the location of a point source found with gtfindsrc under the assumption that the residual seen in the map is produced by a single source; source S is not seen in this map because photons with E> 10 GeV are excluded.

3.1. LAT data analysis

For NGC 4151, NGC 5548, NGC 4258, NGC 6814, NGC 7213, and IC 4329a we analyzed 6.4 years of the Fermi/LAT data, comprising observations carried out between 2008 August 4 and 2015 January 10. For each object, events were selected from a region with a radius of 15° centered on the position of the analyzed source. We performed the unbinned likelihood analysis using the v9r33p0 Fermi Science Tools with CALDB instrument response functions. We used the standard templates for the Galactic (gll_iem_v05_rev1.fits) and the isotropic (iso_source_v05_rev1.txt) backgrounds. In the likelihood analysis we took into account all sources reported in the Fermi/LAT Third Source Catalog (Acero et al. 2015, hereafter 3FGL) within a radius of 15° around the analyzed object. Each of the catalog sources was modeled with a best-fit spectral function (as specified in the catalog, in most cases a power-law) with parameters left free in the model fitting.

Except for NGC 6814 we did not find statistically significant signals. We derived the 95% confidence level ULs for the photon flux using the Bayesian method.

We first checked how the use of both an extended data set and improved models (based on 3FGL) of the celestial regions of interest affects the UL values. Assuming the γ-ray photon spectral index Γ = 2.5, we found ULs for F_{> 0.1GeV} lower by a factor ~4 (for NGC 7213 and IC 4329a) and ~2 (for NGC 5548) than the corresponding UL in Ackermann et al. (2012a). For NGC 4151 the UL value is somewhat dependent on the approach to modeling (as discussed below), but in general, we were unable to reduce it significantly below the value quoted by Ackermann et al. (2012a).

The pion decay spectra can only be described by a simple power-law in limited energy ranges. For the thermal distribution of protons, the π⁰-decay spectrum can be very roughly approximated by a very soft power-law, with Γ ≃ 4, above ≃0.2 GeV. For the power-law proton distribution, it is a power-law above ~1 GeV. Then, to compare with predictions of model T, we assumed Γ = 4 and found UL for F_0.2−1GeV. To compare with predictions of models with η_p> 0, we found UL for F_1−10GeV assuming different values of Γ between 2.1 and 2.7. For all AGNs, except for NGC 5548, the highest ULs correspond to Γ = 2.1 (the difference for other values of Γ does not exceed a factor of 2), and we used them in Fig. 1b. The comparison in Fig. 1b is presented in the most conservative manner, meaning that we used the lowest λ_1−10GeV predicted by the model (for s = 2.1) and the highest UL value (for Γ = 2.1).

3.2. NGC 4151

At the position of NGC 4151 we found a signal with the test-statistic significance TS ≃ 17 for E = 0.2−6 GeV, see Fig. 5d. It most likely contains a contribution from a new γ-ray source found in our analysis, denoted below by S, shifted by only ~ from NGC 4151. Interestingly, however, the gtfindsrc best-fit location for this signal (shown by the black cross), assuming that it comes from a single source, is shifted from the position of NGC 4151 by only 8′. Source S is strongly variable, which hinders a proper assessment of its contribution to the γ-ray signal. The details of our analysis are as follows.

We first examined the significances of γ-ray signals around NGC 4151 by means of their TS values. This revealed three new sources in the region, which were not reported in the Fermi/LAT Second Source Catalog and hence not included in the model used in Ackermann et al. (2012a). For each of these sources we used the gtlike and gtfindsrc tools to find their significance, best-fit position, and spectral parameters. Two of them (indicated by the red symbols in Fig. 5a) have recently been reported in 3FGL, with parameters approximately consistent with those estimated in our analysis. They are 3FGL J1220.2+3434, shifted from NGC 4151 by ~5°, for which our analysis gives TS ≃ 164, Γ ≃ 2.2 ± 0.1 and F_{> 0.1GeV} ≃ (8.4 ± 1.5) × 10^-9 ph s^-1 cm^-2, and 3FGL J1203.2+3847, shifted from NGC 4151 by ~1.5°, for which we derived TS ≃ 32, Γ ≃ 2.3 ± 0.2 and F_{> 0.1GeV} ≃ (3.5 ± 1.4) × 10^-9 ph s^-1 cm^-2. These sources are included in the model used for all results presented in this work.

The position of the third source, S, ~ from NGC 4151, is determined by several photons with energies between 10 and 20 GeV, which arrived from the same direction (within ~8′) between December 2011 and June 2013; we assume that the source exhibited an outburst activity during this 1.5-year period, which we denote as period S. Source S is clearly seen in the TS map built for period S alone, see Fig. 5c, whereas the map built without this period only shows a weak residual, centered on NGC 4151 and not on S, see Fig. 5a. After subtracting all 3FGL sources, we find for source S (1) TS ≃ 30, Γ = 1.78 ± 0.37, F_{> 0.1GeV} ≃ 1.2 × 10^-9 ph s^-1 cm^-2 using the data from period S, and (2) TS ≃ 22, Γ = 2.17 ± 0.20, F_{> 0.1GeV} ≃ 2.1 × 10^-9 ph s^-1 cm^-2 using the full 6.4-year data set. Source S is not reported in 3FGL, which was built using data up to 2012 July, that is, covering only ~30% of period S. For the data from period S, gtfindsrc gives the location of S at α_J2000 = 12^h11^m27^s, δ_J2000 = 38°56′48″. A possible candidate for this source is a BL Lac object, 2E 1209.0+3917, only ≃4′ from this location.

For a 4.9-year data set, without the period S, gtlike shows no signal (TS ≃ 1.5) at the S position. The residual seen in the TS map for this data set, Fig. 5a, can be fully compensated for by adding the point source at the position of NGC 4151 (compare with panel b), and gtlike gives TS ≃ 8, Γ = 2.7 ± 0.3 and F_{> 0.1GeV} = (1.5 ± 0.6) × 10^-9 for this source.

The presence of a variable source S results in an ambiguity for the UL on the γ-ray flux from NGC 4151, especially below 1 GeV, where the point spread function of LAT is much higher than the angular separation between S and NGC 4151. In our analysis we considered two variants:

• (1)

  assuming that S strongly softened before and after period S, butstill provided some contribution at low energies, we used the full6.4-year data set and included source S in the model in addition tothe 3FGL sources, which yielded the UL values denoted bysuperscript “a” in Table 2;

• (2)

  assuming that S completely faded away before and after period S, we used the data for 4.9 years without period S and included only the 3FGL sources in the model (i.e., without source S), which yielded the UL values denoted by superscript “b” in Table 2.

The thus obtained ULs on F_0.2−1GeV differ by a factor of ~4, both values of the resulting UL on λ_0.2−1GeV are shown in Fig. 3. The UL on F_1−10GeV is much less model-dependent. In Fig. 1b we use the UL obtained for the above case (1) with assumed Γ = 2.1, other ULs for this energy range differ by at most 50% and do not affect our conclusions. We also note that adding source S in the model for the first three years has a negligible effect, see cases “cd” in Table 2. We also checked that our results are unchanged if we introduce additional sources to account for several residuals located at ~1° to ~4° from NGC 4151, seen in Fig. 5a.

We used the recent dust-parallax distance measurement and the implied stellar-velocity-based mass (see Table 1), both higher by ~50% than assumed in previous works on high-energy emission from NGC 4151. NGC 4151 shows a moderate X-ray variability with changes of the X-ray flux by up to a factor of ~4 between the dim and bright states. Using the adopted distance and M values and the power-law fits of Lubiński et al. (2010), we found the intrinsic λ_2−10keV ≃ (0.5,1,2) × 10^-3 for the (dim, medium, and bright) state; the parameters of these three states are shown in Fig. 1a). The average Swift/BAT flux during the considered Fermi/LAT observations is higher by a factor 1.2 than the flux of the medium state, see Fig. 4b.

As we see in Figs 1b and 3, the ULs on F_0.2−1GeV and F_1−10GeV derived for NGC 4151 yield the tightest constraints on both λ_0.2−1GeV and λ_1−10GeV of the bright Seyfert 1 galaxies (although IC 4329a has the lowest L_γ/L_X ratio in our sample). For δ = 10^-3, model N overpredicts λ_1−10GeV by a factor of several, regardless of a or β values. Then, if heating of electrons is dominated by Coulomb interactions, the energy content of nonthermal protons is limited to η_p ≲ 0.2.

For δ = 0.5, the model predictions, except for model N with a = 0, are below the UL value. However, as noted before, the much lower L_γ predicted by large-δ models is accompanied by much lower τ. The X-ray spectrum of NGC 4151 is precisely measured, at least in the bright state, and allows testing this property. Fits of the bright-state spectrum with the slab compPS model in Lubiński et al. (2010) give τ^PS ≃ 1.3. To compare it with our model predictions, we simulated the X-ray spectra for our models with λ_2−10keV ≃ 0.002; we used the INTEGRAL/ISGRI response function and assumed the normalization corresponding to the bright state of NGC 4151. Then, the simulated spectra were fitted with the slab compPS model. For δ = 10^-3, the best fit has τ^PS ≃ 1.1 (model H_0.1 with a = 0.95, β = 1, ṁ = 0.65). For δ = 0.5, all simulated spectra were fitted with τ^PS ≲ 0.5. We conclude that although large-δ models are essentially unconstrained by the Fermi UL, they also appear to be inconsistent with the X-ray data.

λ_0.2−1GeV predicted by model T, compared with with the Fermi UL for the above case (1), disfavors models with both a large a and a large β, in a manner consistent with NGC 7213. For case (2), however, the UL is above the model prediction. Stronger constraints can be obtained if NGC 4151 enters the flux level of the dim state for a longer period, as at λ_2−10keV ≃ 10^-3, r_ph is close to the extent of the γ-ray emitting region for thermal protons.

3.3. NGC 6814

The γ-ray signal from the direction of NGC 6814 with TS = 25.6 was found by Ackermann et al. (2012a). However, the source is not reported in 3FGL. It may not have passed the TS > 25 criterion, assumed in 3FGL, as the above value is close to the limit. Our analysis, using data for 6.4 years, confirms the signal with TS ≃ 32.3, and our best-fit parameters for a power-law spectrum (see Table 2) are approximately consistent with those of Ackermann et al. (2012a).

In Fig. 1a we show the parameters of the intrinsic X-ray emission from spectral fits to Suzaku observations in 2011 (Walton et al. 2013). Interestingly, the source is located in the region for an inefficient (thermal synchrotron) source of seed photons. During this observation the source was in a relatively dim state, with the average Swift/BAT flux lower by a factor of ~1.5 than the average during the Fermi/LAT observations, see Fig. 4d.

As noted above, we considered two values of M for NGC 6814. For both, the measured F_1−10GeV gives λ_1−10GeV exceeding the maximum value predicted for the hot flow (in models N with δ = 10^-3), which rules out an origin of the observed γ-ray signal in an inner hot flow. For the higher value of M, NGC 6814 has a similar λ_1−10GeV to the γ-ray-loud Seyfert 2 galaxies, a possible explanation is discussed in Sect. 5.

3.4. NGC 5548 and IC 4329a

The likelihood analysis of the Fermi/LAT data does not reveal any significant signal around these two X-ray-bright Seyfert galaxies (in both TS < 1) and therefore we computed the ULs given in Table 2.

For NGC 5548, we show in Fig. 1a the parameters of the power-law fits from Brenneman et al. (2012) for seven Suzaku observations in 2007. The average Swift/BAT flux during the 6.4 years of Fermi/LAT observations is equal to the average flux during the Suzaku observations, and we used the intrinsic λ_2−10keV from the fit of Brenneman et al. (2012) for the average Suzaku spectrum.

For IC 4329a, we show in Fig. 1a the parameters of the power-law fits from Brenneman et al. (2014) for lower and higher flux states during Suzaku and NuSTAR observations in 2012. We used λ_2−10keV larger by a factor 1.3 to account for the difference between the average Swift/BAT flux for 6.4 years and the flux from measurement simultaneous with Suzaku/NuSTAR observations, see Fig. 4c.

Both sources are observed above λ_ADAF,max but within the range of approximate applicability of our model. During one observation, NGC 5548 is seen in the area of ADAF solutions with inefficient cooling. In the remaining observations both sources show interesting hints of modest hardening below or softening above λ_2−10keV ≃ 0.002 with increasing λ_2−10keV.

For both objects the ULs on the γ-ray photon flux is lower by a factor of 2–3 times than for NGC 4151. In addition, both sources have a similar L_γ/L_X ratio to NGC 4151 (for IC 4329a it is even lower). However, we see in Figs. 1b and 3 that they insignificantly constrain the hot flow models.

3.5. NGC 7213

Despite the large uncertainty in M, NGC 7213 gives interesting constraints on the hot flow model; crucially, it is observed at λ_2−10keV ≃ 10^-4, between the X-ray bright Seyfert 1 galaxies discussed above and NGC 4258. The likelihood analysis of the Fermi/LAT data does not reveal any significant signal around NGC 7213 (TS < 1), and we computed the ULs given in Table 2. In Fig. 1a we show the parameters from a Suzaku observation in 2006 (Lobban et al. 2010) during which the average flux level was similar to the average during 6.4 years, see Fig. 4a.

Figure 1b shows that the derived UL for F_1−10GeV constrains the nonthermal energy content to η_p ≲ 0.1 for δ = 10^-3. Models with δ = 0.5 predicts λ_1−10GeV below the Fermi limit, except for a = 0 and η_p = 1. However, the X-ray spectra predicted for δ = 0.5 appear to be qualitatively inconsistent with the spectrum measured in NGC 7213. The spectra at λ_2−10keV ≃ 10^-4 (the dotted line for a = 0.95 and the overlying dashed line for a = 0 in Fig. 2c) deviate significantly from a power-law shape, whereas the Suzaku spectrum between 0.6 and 50 keV is described by a simple power-law (Lobban et al. 2010). The UL for F_0.2−1GeV rules out a high β flow around a high-a black hole.

3.6. NGC 4258

The precisely determined parameters of this AGN made it an essential object for the development of the hot flow models (e.g., Lasota et al. 1996). However, it has not been included in previous analyses of the Fermi/LAT data. At the position of NGC 4258 we found a weak γ-ray excess above the background, with TS = 9.4. A likelihood analysis of the signal with gtlike gives Γ = 2.5 ± 0.2, F ≃ (2.8 ± 1.2) × 10^-9 ph s^-1 cm^-2. Despite the weak residual, the derived ULs on the photon flux yield the tightest constraint on λ_1−10GeV and λ_0.2−1GeV of all considered objects, see Table 2.

In Fig. 1a we show parameters of the analysis reported by Reynolds et al. (2009) of the Suzaku observation in 2006, who also noted that during this observation NGC 4258 increased its intrinsic X-ray luminosity by a factor of ~2 relative to the average from Swift/BAT measurements. We used λ_2−10keV from the Suzaku fit reduced by a factor of 2.

At λ_2−10keV ≃ 10^-5, the 0.2–1 GeV emission is weakly affected by γγ absorption, and we can observe radiation from the ergospheric region (see Sect. 2). However, as Fig. 3 shows, the strongest expected signal is only at the level of the Fermi UL.

For δ = 10^-3, model N predicts a flux strongly exceeding the derived UL for F_1−10GeV, ruling out η_p ~ 1. For δ = 0.5, the predicted γ-ray flux is an order of magnitude below the UL value. However, the corresponding X-ray spectrum (bottom dashed line in Fig. 2c) seems to rule out this (large-δ) case if the X-ray radiation observed from NGC 4258 is produced by thermal Comptonization.

3.7. Circinus, NGC 1068, and NGC 4945

For the Seyfert 2 galaxies we used the parameters reported in 3FGL, which are approximately consistent (F_1−10GeV within 30%) with previous analyses (Lenain et al. 2010; Ackermann et al. 2012b; Hayashida et al. 2013).

Of these three AGNs, only NGC 4945 shows rapid variability in the hard X-rays, which indicates that it is a transmission-dominated source. We can therefore directly probe the nuclear emission. In Fig. 1a we show fits from Puccetti et al. (2014) for NuSTAR observations in 2013. The average Swift/BAT flux during the four years of Fermi observations used in 3FGL agrees with the average flux during the NuSTAR observations.

In NGC 1068, the direct X-ray emission of is completely obscured along our line of sight and we only see the reflected component. Then, the assessed intrinsic emission is strongly model dependent. We used λ_2−10keV corresponding to the best-fit of Bauer et al. (2015). This assessment, based on the observed reflected component, gives the level of emission from the active nucleus averaged over a long time, possibly over hundreds of years, because a significant fraction of reflection arises at a ~100 pc scale (see Bauer et al. 2015). Similarly, the X-ray emission from Circinus is reflection dominated; we used results from Arévalo et al. (2014). NGC 1068 and Circinus are not shown in Fig. 1a because the estimates for their intrinsic Γ_X are model dependent.

Figure 1b shows that all three γ-ray-loud Seyfert 2 galaxies agree remarkably well after scaling by their central black hole masses, with λ_2−10keV ≃ 0.02 and λ_1−10GeV ≃ 3 × 10^-5 in all three. We also note that their L_1−10GeV/L_2−10keV ratio does not exceed the limit on the luminosity ratio in other AGNs, except for that corresponding to the UL in IC 4329a for Γ = 2.7.

3.8. Centaurus A

The interpretation of the high-energy spectrum of Cen A has some open questions. Its X-ray emission may contain contribution from both the jet and the accretion flow. The latter origin is supported by the thermal-like cutoff, claimed in Beckmann et al. (2011), and by the location in the area of the λ_2−10keV–Γ_X plane occupied by Seyfert galaxies. In Fig. 1a we show the parameters from three Suzaku observations in 2009 (Fukazawa et al. 2011). Fukazawa et al. (2011) also noted hints for a harder power-law component, with Γ_X< 1.6, weakly contributing below 100 keV, which may represent emission from the jet. The main contribution to the γ-ray emission of Cen A most likely comes from the jet. However, Sahakyan et al. (2013) found evidence for a second component, an order of magnitude less luminous than the main γ-ray component. This weaker component has a hard spectrum and dominates above ~4 GeV. Assuming that this component extends down to 1 GeV with the same slope, we find λ_1−10GeV ≃ 1.7 × 10^-6. Based on a comparison with predictions of model H_0.1, this component can be emitted from a flow with δ = 10^-3 and η_p ≃ 0.1−0.2 (slightly dependent on plasma magnetization).

The Swift/BAT flux level during the 2009 Suzaku observations is equal to the average flux during the Fermi/LAT data-taking period used for the analysis in Sahakyan et al. (2013).

4. Nonthermal electrons

Here we consider effects related with nonthermal electrons that are produced by various mechanisms, specifically, their contribution to the observed γ-ray luminosity and to the optical depth.

4.1. Pion decay

If nonthermal protons are present, electrons with γ ≳ 10³ may be produced by charged pion decay. The nonthermal Compton radiation of these electrons may extend to the LAT energy range. The power injected in π^± decay electrons is similar to the luminosity in π⁰ decay γ-rays. However, electrons with γ ≳ 10³ lose most of their energy in synchrotron radiation, therefore their nonthermal Compton emission in the γ-ray range is always much weaker than the π⁰ decay component.

The density of π^±-decay e^±, determined by the pair equilibrium condition (cf. Niedźwiecki et al. 2015), is at least two orders of magnitude lower than the density of the ionization electrons. The production rate of π⁰ is similar to that of π^±, therefore e^± pairs produced by absorption of π⁰-decay γ-rays give at most a similar contribution to the total density as those from π^±-decay. Production of further generations of e^± pairs is unlikely for the luminosities described by our model, see below. Overall, the products of hadronic processes contribute negligibly to the optical depth.

If a strong γγ absorption affects the π⁰-decay γ-rays, the synchrotron radiation of thus produced e^± pairs increases the seed photon input, again at most by a factor of ~2. We neglected this effect in computations presented in this work; we checked that this twice higher seed photon input yields an electron temperature that is reduced by up to 10% and only insignificantly affects our results.

4.2. Direct acceleration

We briefly discuss constraints on the acceleration of electrons resulting from the Fermi/LAT ULs. For the results presented in this section we do not use the simplified description (i.e., with a monoenergetic injection) of the accelerated electrons, and we assume that they are accelerated to the power-law distribution with the acceleration index s_e. As above, the amount of nonthermal electrons is parametrized by η_e. Specifically, the power used for the nonthermal acceleration is Q_acc = η_eδQ_diss, where Q_diss is the total power dissipated in the flow.

Conditions in the inner parts of hot flows surrounding supermassive black holes allow for acceleration of electrons to GeV energies; for instance, using Eq. (7) from Zdziarski et al. (2009), we found the maximum Lorentz factors (limited by the synchrotron losses) of γ_max ~ 10⁵−10⁶ for the typical magnetic field strength of B ~ (10²−10⁴) G. The synchrotron emission of these electrons extends up to ~10 MeV, so they can produce photons in the LAT energy range only by Compton scattering. We implemented the description of nonthermal Compton emission in our MC code and performed simulations with γ_max = 2 × 10⁴, for which the nonthermal Compton spectrum extends up to ~10 GeV.

The γ-ray luminosity produced by the nonthermal Compton scattering primarily depends on η_e and s_e. We found that for s_e ≳ 2.4 and η_e = 1, the values of λ_0.2−1GeV and λ_1−10GeV for the nonthermal Compton are below any UL derived from the LAT data for objects with λ_2−10keV ≲ 10^-3. For s_e ≲ 2.4, the γ-ray Compton component for η_e = 1 exceeds the ULs for NGC 4258, NGC 7213 and NGC 4151, then, for such s_e the fraction of accretion power used for nonthermal acceleration is constrained. We found for instance that for s_e ≃ 2 it is limited to ≲0.05Q_diss (corresponding to δ = 0.5 and η_e ≲ 0.1).

In our computations we assumed that synchrotron radiation emitted below the self-absorption frequency is thermalized, as described in several works on the synchrotron boiler mechanism (e.g., Malzac & Belmont 2009), meaning that the effective heating of thermal electrons is higher than (1−η_e)δQ_diss. The power in the self-absorbed radiation is ≳ 0.2Q_acc for s_e ≳ 2. For s_e = 2.6, it is ~0.8Q_acc. Then, if δ is large, the heating of thermal electrons is always strong, with the synchrotron boiler dominating for η_e ~ 1. This also implies that our conclusions regarding the allowed values of η_e are weakly dependent on other parameters of the flow. Although Q_diss shows a significant dependence on a, for example, both the nonthermal Compton L_γ and the thermal Compton L_X scale roughly linearly with Q_diss (the latter either through the synchrotron boiler or direct MHD heating), making the L_γ/L_X ratio weakly dependent on the other parameters.

4.3. γγ absorption

Secondary e^± pairs, produced by γγ absorption, in general affect radiative properties of flows in the considered range of λ_2−10keV only insignificantly. The first generation of pairs contributes at most a few per cent to the Thomson optical depth. Copious creation of secondary pairs through electromagnetic cascades requires (1) the Compton process to dominate over synchrotron, and (2) Compton photons produced above the threshold for pair production. Typically, at least one of these conditions is not fulfilled. In the innermost parts of the flow, Compton energy losses are comparable to synchrotron losses only for electrons with Lorentz factors not exceeding a few hundred. At higher energies, synchrotron losses strongly dominate as a result of the Klein-Nishina decline in the Compton scattering rate, hence, the development of a cascade is suppressed. In turn, the lower energy electrons, for which the more efficient Compton losses are mostly due to interactions with X-ray photons from thermal Comptonization, mainly produce photons with E ≲ 10 MeV. However, at λ_2−10keV ≲ 10^-3, the flows are mostly transparent to such photons. Therefore, we do not expect a significant contribution of secondary pairs to the Thomson depth of the flow. We note, however, that at λ_2−10keV ≳ 10^-3 the innermost region should be opaque to photons with E ≳ 1 MeV, and in this regime the secondary pairs may then contribute considerably to the optical depth.

We also note that spatially extended pair cascades, such as discussed by Zdziarski et al. (2009), for instance, are unlikely to strengthen the observed γ-ray flux. Such extended cascades can only develop in the region with τ_γγ ~ 1. In a hot flow, τ_γγ increases quickly toward the center, see Fig. 7 in N13, therefore the free path of γ-rays produced within r_ph is typically very short. Furthermore, as noted above, the absorbed energy is lost to synchrotron radiation rather than used for cascade development.

5. Discussion

We first briefly comment on the X-ray spectral properties of the considered AGNs with λ₂₋₁₀ ≲ 10^-3. At these luminosities, the seed photons for Comptonization must be produced internally to the flow, and nonthermal synchrotron radiation is an almost obvious process providing the amount of seed photons needed to explain the AGN data (Poutanen & Veledina 2014; Niedźwiecki et al. 2015). If δ is small, the nonthermal electrons are provided by the decay of charged pions, unless both a and β are small and protons are thermal, see Sect. 2.1. We emphasize that the weakness of γ-ray signals related with π⁰ decay does not contradict a copious injection of nonthermal electrons by π^± decay because the most efficient pion production occurs in the region opaque to γ-rays. If δ is large, the π^± decay does not provide a sufficient amount of nonthermal electrons because the flow density is too low (here inefficient e^± production is strictly related with very low level of γ-ray activity). Then, direct acceleration must take place if δ is large. These models (i.e., with either a significant π^± production or η_e ≳ 0.1) fill the upper gray region in Fig. 1a. Typical Seyfert spectra clearly agree with these solutions.

Solutions with seed photons from thermal synchrotron alone (including models with δ = 0.5 and η_e = 0, as well as model T with δ = 10^-3, a = 0 and β = 1) fill the lower gray region. Obviously, intermediate cases are also possible, for instance, with δ = 0.5 and η_e ≲ 0.01. Seyferts are occasionally observed with parameters corresponding to inefficient cooling (NGC 6814 and one of NGC 5548 data points). These harder-than-typical spectra may result from a decline of nonthermal electron injection. For large δ it may correspond to the decrease of η_e, and for small δ to the decrease of β (again, for small a and thermal protons).

Observational assessment of nonthermal acceleration processes in astrophysical plasmas is an essential question. Hot flows in which ultra-relativistic acceleration of protons takes place can be considered a source of neutrinos detected by IceCube (e.g., Kimura et al. 2015). Accretion flows should also be strong acceleration sites according to hybrid thermal and nonthermal Comptonization models (e.g., Vurm & Poutanen 2009; Malzac & Belmont 2009; Veledina et al. 2011) that are widely applied to explain X-ray spectra of accreting black holes; here at least ~50% of the energy must be provided to electrons by nonthermal processes if thermalization of electrons via Coulomb collision and synchrotron self-absorption is taken into account (Del Santo et al. 2013). Interestingly, these models are typically fitted with s_e ≳ 2.4. As discussed in Sect. 4.2, the γ-ray luminosities from nonthermal Compton emission of flows with such s_e are below the LAT ULs derived in Sect. 3 for Seyfert galaxies, even if η_e = 1.

Observational evidence of the presence of non-thermal particles in accretion flows includes detection of MeV tails (e.g., Droulans et al. 2010; Zdziarski et al. 2012) and patterns of optical/infrared evolution (Poutanen et al. 2014), both observed in stellar black-hole systems. However, these features are produced by nonthermal electrons, which may come from pion decay, and in principle, they do not require direct acceleration.

Observations in the GeV range provide the means to directly measure the efficiency of proton acceleration in hot flows. The obtained ULs for F_1−10GeV in NGC 4258, NGC 7213 and NGC 4151 constrain the fraction of accretion power used for relativistic acceleration of protons to at most η_p ≃ 0.1 (if δ ≲ 0.1). Remarkably, these three AGNs allow probing various parts of the flow because at λ_2−10keV ~ 10^-5, 10^-4 and 10^-3, radiation produced at a few, ten, and several tens of R_g, respectively, dominates the observed 1–10 GeV flux. The component dominating the core emission of Centaurus A above ~4 GeV may come from a hot flow in which η_p ≃ 0.1−0.2. In this scenario, its η_p should be higher at least by a factor of ~2 than in NGC 7213 (see Fig. 1b), which has similar λ_2−10keV. This might reflect differences between flows powering radio-loud and radio-quiet AGNs. We also note that in this hot flow scenario for Centaurus A, the observed γ-rays (above 4 GeV) provide information about the region close to the jet formation site².

We now briefly discuss hard spectral states of black hole binaries in which hot accretion flows, similar to those powering AGNs, are most likely present. In the innermost parts of flows powering both classes of objects the same physical processes are expected. This view is supported by the similarity between the X-ray spectra of Seyferts and black hole binaries in the hard spectral states (e.g., Zdziarski et al. 1996), with slight differences simply explained by scaling of seed photon energies with the black hole mass (cf. Niedźwiecki et al. 2015). These differences do not affect the γ-ray emission, and the relation between λ_{2−10 keV} and λ_{0.2−1 GeV} or λ_{1−10 GeV} shown in Figs. 1–3 is probably also relevant for flows around stellar-mass black holes.

An estimate of the γ-ray activity level, which can be compared with the hot flow model, has only been obtained for the hard state of Cyg X-1 by Malyshev et al. (2013). Figure 1b shows that it gives a slightly tighter constraint than Seyfert galaxies, limiting η_p to at most a few per cent for δ = 10^-3, and at most ~10% for δ = 0.5 and a = 0. A significant nonthermal proton component, with η_p ~ 0.5, is only allowed if δ ~ 0.5 and a is high. This finding places interesting constraints on hybrid models (see above) if protons are accelerated with at least a similar efficiency to electrons (as expected in many acceleration scenarios). It may be also interesting to note that similar values of λ_{1−10 GeV} correspond to the Cyg X-1 hard state (≃6 × 10^-7) and our gtlike fits to signals at the positions of NGC 4151 (≃5 × 10^-7) and NGC 4258 (≃3 × 10^-7), although the significance of the γ-ray signal is rather weak (TS = 15.6) in Cyg X-1 and marginal in NGC 4258 and NGC 4151.

Finally, we comment on the γ-ray emission of NGC 4945, NGC 1068, and Circinus. All three exhibit both an AGN and starburst activities. NGC 4945 and NGC 1068 approximately follow some multiwavelength correlations established in other star-forming galaxies (Ackermann et al. 2012b), which suggests that their γ-ray emission may be related to the latter activity. However, Circinus does not agree with these correlations, furthermore, its L_γ exceeds the calorimetric limit (where all cosmic rays produced by supernovae interact with interstellar matter), see Hayashida et al. (2013). Furthermore, Lenain et al. (2010) noted that L_γ measured in NGC 1068 exceeds by a factor of ~10 the level expected for the total gas mass and supernova rate estimated in this galaxy, which led them to suggest that the γ-rays may be produced in a jet. Similarly, Yoast-Hull et al. (2014) found that their detailed starburst model underestimates the observed γ-ray flux and overestimates the radio flux for NGC 1068, and they argued that these problems would be resolved if the AGN were the primary source of γ-rays. This shows that the origin of γ-rays still seems rather unclear.

Seyfert galaxies show a strong correlation between nuclear star-formation and the AGN luminosity (e.g., Diamond-Stanic & Rieke 2012). Masses of supermassive black holes in Circinus, NGC 1068, and NGC 4945 are precisely determined by H₂0 megamaser measurements. Their values imply an intriguing similarity of the values of both λ_2−10keV and λ_1−10GeV in all three objects. The question then arises whether their γ-ray emission might be related with accretion processes in their active nuclei. For the intrinsic λ_2−10keV ≃ 0.02, the geometry of accretion flow is very uncertain. It cannot be studied through Fe Kα line distortions because the X-ray spectra of all three AGNs are completely absorbed below 10 keV. Some analogies with black hole binaries, in which (1) luminous hard states have been observed at such λ_2−10keV (e.g., Zdziarski et al. 2004); (2) softening of intrinsic X-ray slopes above λ_2−10keV ≃ 10^-3 is observed that closely resembles the behavior seen in NGC 4945 (see, e.g., Fig. 1 in Yang et al. 2015); and (3) strong thermal components are seen at λ_2−10keV ~ 0.01, indicating that the cold disk extends down to r_tr ≲ 10 (e.g., Cabanac et al. 2009) suggest that inner hot flows may be present in these AGNs, but only within the innermost several R_g. This value of r_tr ~ 10 is tentatively shown in Fig. 1a. For such a small, strongly energized (the main release of gravitational energy occurs in this region) and relatively dense flow, we cannot exclude the possibility that a sufficient γ-ray flux is produced in an outer layer and escapes avoiding γγ absorption (taking into account a strong local anisotropy of the X-ray radiation field, see Sect. 2). Any stronger conclusions for this scenario would be premature; a self-consistent model, compatible with the X-ray data, needs to be developed before a quantitative assessment of the γ-ray flux can be made.

We do not expect direct γ-ray signals from central parts of accretion flows in black hole binaries even if they have the same geometry at this λ_2−10keV. Any emission in the GeV range would be absorbed in interaction with soft X-ray emission of a cold disk extending to small r, see Bednarek (1993) for quantitative estimates. In AGNs, a UV-emitting cold disk only attenuates γ-rays at ~100 GeV.

Another possibility involves escape of protons accelerated in the flow and their interaction with circumnuclear matter. About 0.1% of the accretion power would need to be carried away by escaping protons to explain the γ-ray measurements. This would also allow explaining the γ-ray signal in NGC 6814 by the same emission mechanism, assuming that (1) its supermassive black hole has M ≃ 2 × 10⁷M_⊙ (the higher of the two values considered above), for which its λ_1−10GeV is similar to those of NGC 4549, Circinus, and NGC 1068; and (2) protons were injected into the circumnuclear region during an earlier phase of higher (by almost 2 orders of magnitude relative to current) activity of NGC 6814, in the scenario similar to that proposed for the Galactic center (e.g., Chernyakova et al. 2011). In this scenario, an energy-dependent diffusive loss of protons may explain the slightly softer γ-ray spectrum of NGC 6814.

6. Summary

Hot flows are supported by the proton pressure, and the production of pions in collisions of the energetic protons should be their generic property. If this is true, a lack of detection of a related γ-ray signal could be regarded as an argument against this class of models. However, the expected flux depends on several parameters, and in the extreme case of a strongly magnetized flow around a nonrotating black hole, the conditions for pion production are not achieved if protons are thermal.

We compared the prediction of the model with several well-studied AGNs, for which the available data allowed a robust determination of the nuclear luminosity scaled by the Eddington value. At the positions of NGC 4151 and NGC 4258, we found weak residuals, whose nature (background fluctuation or an actual emission) could not be assessed because of their low statistical significance.

The highest γ-ray luminosities are predicted if electrons are heated mostly by Coulomb interactions. If additionally most of the accretion power is used for the relativistic acceleration of a small fraction of protons, the highest level of λ_1−10GeV ≃ 10^-5 is expected for λ_2−10keV between ~10^-5 and 10^-3. This exceeds the Fermi limit on λ_1−10GeV for NGC 4258, NGC 7213, and NGC 4151 at least by a factor of several. Thus, the Fermi ULs provide an interesting constraint on the MHD processes that convert the accretion power into the kinetic energy of protons: a rather uniform heating of all protons would need to take place. Given this constraint, we used λ_0.2−1GeV predicted for thermal distribution of protons to rule out a flow with a high (equipartition) value of β around a high-a black hole for NGC 7213 and and NGC 4151; for the latter this result is model dependent.

If a large part of the accretion power is used for a direct heating or acceleration of electrons, the expected γ-ray fluxes are below the Fermi ULs, except for NGC 4151 and NGC 7213, for which η_p ≳ 0.5 was ruled out with a small a (yielding some constraints on the acceleration efficiency, which may be relevant for hybrid models). However, this (large δ) class of models appears to be disfavored by the predicted properties of X-ray spectra, which significantly deviate from a power law at λ_2−10keV ≲ 10^-4 and are described by (too) small τ at larger λ_2−10keV.

The γ-ray-loud Seyfert 2 galaxies radiate at a much higher Eddington ratio than other nearby AGNs, possibly as a result of enhanced fueling that is related with starburst activity in their nuclear regions. At this λ_2−10keV, the transition between the hot flow and the cold disk possibly takes place within the central ~10R_g. Whether it creates conditions suitable for the escape of γ-rays or relativistic protons remains an open question.

Acknowledgments

We thank F. Longo for help with the Fermi/LAT data analysis, P. Lubiński for help with the INTEGRAL software and the referee for motivating us to include the results for nonthermal Compton emission. We made use of data and software provided by the Fermi Science Support Center, managed by the HEASARC at the Goddard Space Flight Center, and Swift/BAT transient monitor results provided by the Swift/BAT team. R.W. was supported by the Polish MNiSW grant for young scientists B1511500001045.02. A.N. was supported by the Polish NCN grant DEC-2011/03/B/ST9/03459. F.G.X. was in part supported by the Strategic Priority Research Program “The Emergence of Cosmological Structures” of CAS (grant XDB09000000), and the Natural Science Foundation of China (grants 11203057, 11133005).

References

All Tables

All Figures

Fig. 1 a) X-ray photon spectral index as a function of intrinsic λ2−10keV. The observational data correspond to the best fits of the high-quality observations described in Sect. 3. The gray regions indicate the location of solutions for the model with seed photons from thermal synchrotron only lower and including nonthermal synchrotron (of pion-decay or directly accelerated electrons; upper region). Ticks in the top axis show the approximate value of rtr (see text). The thick vertical line indicates the maximum luminosity for the standard ADAF solutions (i.e., unaffected by the Coulomb cooling of protons). b) λ1−10GeV as a function of λ2−10keV, see Tables 1 and 2 for the observational data. The solid lines are for model N with δ = 10-3, a = 0.95, β = 9 (triangles), model N with δ = 0.5, a = 0, β = 9 is plotted with circles), model H0.1 with δ = 10-3, a = 0.95, β = 9 with squares, and all three are the same as in Fig. 2b below. The dashed line with open squares denotes the rest-frame 1–10GeV luminosity in the last model (i.e., H0.1). In panel b), λ2−10keV corresponds to the average X-ray luminosity during the Fermi/LAT observations (see text). This differs in some objects by up to a factor of 2 relative to panel a). The point for Cyg X-1 shows the hard state parameters according to Malyshev et al. (2013).

In the text

Fig. 2 a) Photosphere radius for photons with E = 0.1, 1, and 10 GeV from bottom to top as a function of λ2−10keV in the model with δ = 10-3, a = 0.95, β = 1 (in other models the dependence of rph on λ2−10keV is very similar). b) λ1−10GeV as a function of λ2−10keV in models with ηp> 0 and s = 2.1; the dotted lines plot models N with δ = 10-3, (magenta) circles indicate models N with a = 0.95 and β = 9 with ṁ = 0.05, 0.1, and 0.3, (green) squares stand for a = 0, β = 9, ṁ = 0.1, 0.3 and 0.5; the dashed lines show models H0.1 with δ = 10-3, (black) circles models with a = 0.95, β = 1 with ṁ = 0.033, 0.1, 0.3 and 0.5; (blue) squares indicate a = 0.95, β = 1 with ṁ = 0.1, 0.3, 0.5, 0.8 and 1.2; the solid lines show models N with δ = 0.5, (red) squares depict a = 0, β = 1 with ṁ = 0.004, 0.01, 0.033, 0.1 and 0.3; (orange) circles stand for a = 0.95, β = 9 with ṁ = 0.003, 0.01, 0.033 and 0.1. c) The high-energy part of the observed spectra of leptonic emission. (Black) solid lines indicate model H0.1 with δ = 10-3, a = 0.95, β = 1 and ṁ = 0.033, 0.1, 0.3; (red) dashed lines show model N with δ = 0.5, a = 0, β = 1 and ṁ = 0.004, 0.01, 0.033, 0.1; (blue) dotted line shows model N with δ = 0.5, a = 0.95, β = 9 and ṁ = 0.01.

In the text

Fig. 3 λ0.2−1GeV as a function of λ2−10keV, see Tables 1 and 2 for the observational data. All model points are for thermal protons (model T) with δ = 10-3; the solid line with circles shows a = 0.95, β = 1, ṁ = 0.033, 0.1, 0.3, and 0.5; the solid line with full squares depicts a = 0.95, β = 9, ṁ = 0.1, 0.3, 0.5, 0.8, and 1.2; the dashed line with open squares plots the estimate (see text) for an edge-on observer for a = 0.998, β = 9, ṁ = 0.033, and 0.1.

In the text

Fig. 4 Small black circles with error bars show Swift/BAT 15–50 keV light curves in seven-day bins for NGC 7213, NGC 4151, IC 4329a, and NGC 6814. The red vertical line in each panel shows the average flux during the 6.4 years of analyzed Fermi/LAT observations. The large orange circles in panels (a) c) d)) show the times of the adopted Suzaku and NuStar observations and the average flux in time bins equal to the duration of these pointing observations. The red vertical line in panel b) indicates the end of the three-year period used for comparison with Ackermann et al. (2012a). The blue line in the bottom of panel b) delineates period S (see Sect. 3.2). The flux levels of dim (D), medium (M), and bright (B) states, defined in Lubiński et al. (2010), are shown in the left part of panel b).

In the text

Fig. 5 TS maps for the region around NGC 4151; in all panels the green cross shows the location of NGC 4151 and the orange cross the location of source S. In all panels only the 3FGL sources were subtracted from the maps, except for panel b), where an additional point source was added at the position of NGC 4151. a) The 8° × 8° degrees region with a pixel size of . The map was built for the energy range 0.2–30 GeV, neglecting period S (see text); the red and white crosses indicate the location of 3FGL sources. b) The same as in a), but the model includes a point source at the position of NGC 4151 to compensate for the residual seen in a). c) The degree region with a pixel size of 0.̊075 for 0.2–30 GeV for period S. d) The region with a pixel size of 0.̊075 for 0.2–6 GeV and the full period of 6.4 years; the black cross shows the location of a point source found with gtfindsrc under the assumption that the residual seen in the map is produced by a single source; source S is not seen in this map because photons with E> 10 GeV are excluded.

In the text