High-energy Neutrinos from Magnetized Coronae of Active Galactic Nuclei and Prospects for Identification of Seyfert Galaxies and Quasars in Neutrino Telescopes

It is widely believed that turbulent magnetic fields generated by the MRI are responsible for the angular momentum transport in accretion flows, which has been confirmed by MHD simulations (Stone & Pringle 2001; Machida & Matsumoto 2003; Ohsuga & Mineshige 2011; Narayan et al. 2012; Porth et al. 2019). Turbulence can also be generated through magnetic reconnections resulting from magnetic fields amplified by the MRI. CRs in the MHD turbulence randomly change their energies through interactions with MHD waves. If the acceleration in momentum space homogeneously occurs without significant trapping (Lemoine & Malkov 2020), the phenomenon can be described by the following diffusion equation in energy space:

$$\begin{eqnarray}&&\begin{array}{r}\begin{array}{l}\displaystyle \frac{{\rm{\partial }}{{ \mathcal F }}_{{\rm{p}}}}{{\rm{\partial }}t}=\displaystyle \frac{1}{{\varepsilon }_{{\rm{p}}}^{2}}\displaystyle \frac{{\rm{\partial }}}{{\rm{\partial }}{\varepsilon }_{{\rm{p}}}}\left({\varepsilon }_{{\rm{p}}}^{2}{D}_{{\varepsilon }_{{\rm{p}}}}\displaystyle \frac{{\rm{\partial }}{{ \mathcal F }}_{{\rm{p}}}}{{\rm{\partial }}{\varepsilon }_{{\rm{p}}}}+\displaystyle \frac{{\varepsilon }_{{\rm{p}}}^{3}}{{t}_{{\rm{c}}{\rm{o}}{\rm{o}}{\rm{l}}}}{{ \mathcal F }}_{{\rm{p}}}\right)\\ \,-\displaystyle \frac{{{ \mathcal F }}_{{\rm{p}}}}{{t}_{{\rm{e}}{\rm{s}}{\rm{c}}}}+{\dot{{ \mathcal F }}}_{{\rm{p}},{\rm{i}}{\rm{n}}{\rm{j}}},\end{array}\end{array}\end{eqnarray} \tag{ 1 }$$

where ${{ \mathcal F }}_{{\rm{p}}}$ is the CR distribution function ($dN/d{\varepsilon }_{{\rm{p}}}\,=4\pi {{\rm{p}}}^{2}{{ \mathcal F }}_{{\rm{p}}}/c$), ${D}_{{\varepsilon }_{{\rm{p}}}}$ is the diffusion coefficient in energy space, t_cool is the cooling timescale, and t_esc is the escape timescale. The injection function is given by ${\dot{{ \mathcal F }}}_{{\rm{p}},{\rm{i}}{\rm{n}}{\rm{j}}}={f}_{{\rm{i}}{\rm{n}}{\rm{j}}}{L}_{X}\delta ({\varepsilon }_{{\rm{p}}}-{\varepsilon }_{{\rm{i}}{\rm{n}}{\rm{j}}})/[4\pi {({\varepsilon }_{{\rm{i}}{\rm{n}}{\rm{j}}}/c)}^{3}{ \mathcal V }],$ where ε_inj is the injection energy and f_inj is the injection fraction at ε_inj, L_X is the X-ray luminosity, ${ \mathcal V }$ is the volume of the corona, and δ(x) is the delta function. Note that in this scenario, the injection fraction is proportional to the dissipation rate in the corona, and therefore is proportional to L_X . The diffusion coefficient is assumed to scale as ${D}_{{\varepsilon }_{{\rm{p}}}}\propto {\varepsilon }_{{\rm{p}}}^{q}$, where q is the power-law index of the diffusion coefficient in momentum space. If the gyro-resonant scattering is dominant, which is not necessarily the case, it corresponds to the spectral index of the turbulence power spectrum. For demonstrative purposes, we use q = 5/3 motivated by the Kolmogorov turbulence, but q = 3/2 (that can be expected for isotropic MHD turbulence in the magnetically dominated regime) or q = 2 (that is a hard-sphere type; see, e.g., Kimura et al. 2019b) may be possible. Note that in general the CR spectrum can be a power law if the trapping is significant (Lemoine & Malkov 2020). As cooling processes, we take into account Bethe−Heitler, pp inelastic collision, and photomeson (p γ) production processes. We consider advective (infall to a central black hole) and diffusive escapes. The maximum energy of the CRs is set by the balance between the acceleration and these losses. For the latter, the principal component is the Bethe−Heitler process due to abundant disk photons, as emphasized in Murase et al. (2020). This results in a cutoff feature in the CR spectrum below PeV energies. In this scenario, we introduce two parameters, the ratio of the CR pressure to the thermal pressure P_CR/P_th and turbulent strength ${\eta }_{\mathrm{tur}}^{-1}$. The former determines the normalization of CRs, whereas the latter is related to the maximum energy of CRs. Note that we use the thermal pressure at the virial temperature, i.e., P_th = n_p kT_vir, where n_p is the thermal proton density and T_vir is the virial temperature, and the realistic thermal pressure can be reduced if the ion temperature is lower than T_vir. See the Supplemental Material of Murase et al. (2020) for details (see also Section III-B of Kimura et al. (2019a), in the context of radiatively inefficient accretion flows).

Recent Particle-In-Cell Simulations revealed that relativistic magnetic reconnections in the ion-electron plasma, where the ion magnetization parameter σ_i = B²/(8π n_p m_p c²) ≳ 1, can accelerate nonthermal relativistic particles very efficiently, given that the plasma β, which is defined by β ≡ 8π n_p kT_p/B², is smaller than ∼0.01 (Guo et al. 2016; Werner et al. 2018; Ball et al. 2018b). For AGN coronae, β can be sufficiently small (e.g., β ≲ 1–3) but we typically expect that the plasma is semi-relativistic, σ_i = 2kT_p/(m_p c² β) < 1 (Werner et al. 2018). It is still uncertain whether protons can be accelerated above ∼ m_p c², and further acceleration via turbulence would be necessary (Comisso & Sironi 2018; Zhdankin et al. 2019; Lemoine & Malkov 2020; Wong et al. 2020). Not only the stochastic acceleration but also the first-order Fermi acceleration may operate due to scatterings with reconnection outflows (e.g., Hoshino 2012; Pisokas et al. 2018). In this work, for the purpose of phenomenological studies, we simply assume that the injection spectrum of protons is given by a single power law with an exponential cutoff, providing the power-law index, s, as a parameter. We assume that the corona can be regarded as an electron-proton plasma. X-ray observational data seem consistent with the electron-proton plasma (Ricci et al. 2018), although better quality data are necessary to confirm it. For electron-positron-proton plasma, the nonthermal proton production efficiency via magnetic reconnections is expected to be lower than that for electron-proton plasma (e.g., Petropoulou et al. 2019).

To calculate neutrino emission, we solve the transport equation for the nonthermal protons with single-zone and steady-state assumptions:

$$\begin{eqnarray}&&-\displaystyle \frac{d}{d{\varepsilon }_{{\rm{p}}}}\left(\displaystyle \frac{{\varepsilon }_{{\rm{p}}}{N}_{{\varepsilon }_{{\rm{p}}}}}{{t}_{{\rm{c}}{\rm{o}}{\rm{o}}{\rm{l}}}}\right)={\dot{N}}_{{\varepsilon }_{{\rm{p}}},{\rm{i}}{\rm{n}}{\rm{j}}}-\displaystyle \frac{{N}_{{\varepsilon }_{{\rm{p}}}}}{{t}_{{\rm{e}}{\rm{s}}{\rm{c}}}},\end{eqnarray} \tag{ 2 }$$

where ${N}_{{\varepsilon }_{{\rm{p}}}}=dN/d{\varepsilon }_{{\rm{p}}}$, and ${\dot{N}}_{{\varepsilon }_{{\rm{p}}},{\rm{i}}{\rm{n}}{\rm{j}}}\propto {\varepsilon }_{{\rm{p}}}^{-s}$ is the injection term.

We consider the same cooling and escape processes as those in the stochastic acceleration scenario. The maximum proton energy of the injection spectrum is determined by either the size of the acceleration region or the balance between acceleration and the total losses via cooling and escape processes, as ${E}_{{\rm{p}}}^{max}=min[{E}_{{\rm{p}}}^{{\rm{r}}{\rm{e}}{\rm{c}}},{E}_{{\rm{p}}}^{{\rm{c}}{\rm{o}}{\rm{o}}{\rm{l}}}]$. In general, particle acceleration consists of multiphase processes including the electric field acceleration and Fermi acceleration (Li et al. 2021). In the absence of losses, the former energy can be written as ${E}_{{\rm{p}}}^{{\rm{r}}{\rm{e}}{\rm{c}}}\sim {eBl}_{{\rm{r}}{\rm{e}}{\rm{c}}}$, where l_rec is the reconnection region thickness that can be smaller than the system size (Hoshino 2012; del Valle et al. 2016). The latter is determined by t_acc = t_cool, where t_acc is the proton acceleration time and t_cool is the proton cooling time. The acceleration time depends on details of the magnetic reconnection process. For example, particle bounces may occur in converging magnetic fluxes with reconnection velocity (de Gouveia dal Pino & Lazarian 2005; del Valle et al. 2016; Medina-Torrejon et al. 2021) and/or in outflows with Alfvén velocity (Drake et al. 2006; Hoshino 2012). Recent simulations have shown that the acceleration is efficient especially for relativistic reconnections (Hoshino & Lyubarsky 2012). For demonstrative purposes, we assume t_acc = η_acc r_L/c, where η_acc is the acceleration efficiency parameter and r_L is the Larmor radius. At sufficient high energies the Fermi-like acceleration and η_acc can be as small as ∼3(c/V_rec); (Giannios 2010; Zhang et al. 2021), although details are uncertain in nonrelativistic reconnections (see Hoshino 2012; del Valle et al. 2016). In this work, we assume η_acc = 300, but we will see that ${E}_{{\rm{p}}}^{max}={E}_{{\rm{p}}}^{{\rm{r}}{\rm{e}}{\rm{c}}}$ to better explain the observation of NGC 1068 (see Section 3.1). The normalization of CRs are provided such that $\int {\dot{N}}_{{\rm{i}}{\rm{n}}{\rm{j}}}{\varepsilon }_{{\rm{p}}}d{\varepsilon }_{{\rm{p}}}={\epsilon }_{{\rm{C}}{\rm{R}}}\dot{M}{c}^{2}$, where _CR is the energy fraction carried by CRs and $\dot{M}$ is the mass accretion rate in the accretion disk. Only the fraction of the accretion energy is spent to accelerate CRs. Therefore, we use the total mass accretion rate onto the black hole to deduce the CR normalization. The total mass accretion rate is proportional to the bolometric luminosity, which is estimated from the intrinsic X-ray luminosity. Note that the corona is connected to the accretion disk through magnetic fields, and a significant fraction of the accretion energy may be dissipated in the disk−corona interface.

The previous works focused on the stochastic acceleration (Murase et al. 2020) or shock acceleration mechanism (see below), and the magnetic reconnection scenario was not considered in detail. In this work, we show that the diffuse neutrino flux can be explained by the magnetic reconnection scenario in the AGN corona model. In Section 6 we discuss the possibility of explaining the diffuse neutrino flux at medium energies with this scenario and comment on the feasibility of identifying this scenario in the near future.

This work focuses on the standard magnetized corona scenario based on recent global MHD simulations, which does not involve accretion shocks.

Alternatively, accretion shocks have been discussed as a CR production site in the accretion flows (Begelman et al. 1990; Stecker et al. 1991; Inoue et al. 2019). Originally, this scenario was proposed to explain X-ray observations of AGNs by CR-induced electromagnetic cascades (Berezinskii & Ginzburg 1981; Zdziarski 1986). However, cutoff and softening features were discovered in bright AGNs in the 1990s (Maisack et al. 1993; Zdziarski et al. 2000), which ruled out the hadron-induced electromagnetic cascade scenario for the observed X-rays. The accretion shock could still be viable as a neutrino production site as long as the CR-induced cascade flux is below the X-ray and gamma-ray data (Stecker 2005, 2013; Inoue et al. 2019). However, there are several problems with this scenario. First, the accretion shocks are not seen in any global MHD simulations dedicated to accretion flows (e.g., Narayan et al. 2012; Kimura et al. 2014; Jiang et al. 2019). Although there are solutions with an accretion shock in the one-dimensional hydrodynamic equation system with a steady-state assumption (Becker et al. 2011; Chattopadhyay & Kumar 2016), such solutions are not realized in these dedicated simulations. The second point is the angular momentum of the accreting matter. For example, Inoue et al. (2019) explicitly assumed the freefall spherical accretion (see their Equation (2)). However, we expect a large specific angular momentum at an outer scale, which prevents the matter from freely falling to the vicinity of the supermassive black hole. An anomalously efficient angular momentum transport mechanism is necessary to form an accretion shock, which contradicts with the freefall assumption made by Inoue et al. (2019, 2020). Lastly, the accretion rate through the accretion shocks is assumed to be extremely high to explain the IceCube data (Inoue et al. 2019). It is comparable to or even higher than the accretion rate expected in the geometrically thin, optically thick disk, which requires two accretion components with totally different angular momentum distributions. Note that P_CR < 0.5P_th is required in this scenario, otherwise hot coronae cannot be maintained.

We also note that the Bethe−Heitler process, which was often ignored, is relevant. The 10–100 TeV neutrino flux from the photomeson production with coronal X-rays is suppressed by Bethe−Heitler interactions with disk photons (Murase et al. 2020). By contrast, given that shock acceleration leads to higher maximum energies, PeV neutrinos are mainly produced via the photomeson production with disk photons (Stecker et al. 1991).

Also, the realistic coronal structure is inhomogeneous. As shown in Gutiérrez et al. (2021), radio emission is more likely to come from outer coronal regions. Such a multizone situation is commonly invoked to explain radio emission from blazars and radio galaxies because it is known that radio emission cannot be explained by single-zone models.

Note that Kalashev et al. (2015) assumed electric acceleration in a gap formed in the magnetosphere. The gap formation is possible only in low-luminosity AGNs such as M87, but it is likely to be screened in the case of Seyferts and quasars because the accretion rate is so high that the plasma satisfies the quasi-neutrality condition (Levinson & Rieger 2011).

NGC 1068 is a Seyfert 2 galaxy at z = 0.00303 (Tully 1988) with a heavily obscured nucleus. It is considered to be one of the best studied AGNs, thanks to its role in the AGN unification scenarios, which was initially proposed for explaining the broad optical lines in its polarized light. The column density of NGC 1068 is ≳ 10²⁵ cm⁻², implying a Compton thick environment. The high-density environment at the core of NGC 1068 and the high level of mass accretion rate provide favorable environments for the efficient production of high-energy neutrinos.

In order to estimate the neutrino spectrum in the AGN corona model for NGC 1068, we use the intrinsic X-ray luminosity, L_X , suggested by Swift-BASS (Ricci et al. 2017). The BASS catalog provides X-ray fluxes and luminosities in a wide range, covering both soft (below 2 keV) and hard (above 10 keV) X-rays. Here, we use the measured spectrum and luminosity in the 2–10 keV band, which provides a medium range. This choice suits the disk-corona modeling, as the abundance of data in this energy band enables us to use an empirical correlation.

The intrinsic X-ray luminosity, L_X, measured from the direction of NGC 1068 by BAT is ∼ 10⁴³ erg s⁻¹ for the 2–10 keV band, although it could be larger depending on estimates of the column density. We use this luminosity to model the neutrino spectrum from NGC 1068 in three different scenarios.

Assuming a single power-law spectrum, the IceCube Collaboration has reported a best-fit flux of 3 × 10⁻¹¹ TeV⁻¹ cm⁻² s⁻¹ at 1 TeV with an index of ≃3.2 (Aartsen et al. 2020b). As a criterion, we adopt parameters that maintain the flux within the best-fit reported flux and its uncertainties reported by the IceCube Collaboration.

First, we model neutrino production assuming the stochastic acceleration scenario. As mentioned earlier, in this scenario, the neutrino spectrum has a more complicated shape than a single power law. Accommodating the IceCube flux at TeV energies requires a relatively high normalization, while the spectrum has to cut off fast enough that the spectrum drops around 100 TeV. Such conditions would result in a high level of CR pressure in the corona model.

In order to maintain realistic scenarios, we restrict ourselves to the range of parameters for which the ratio of the CR pressure (P_CR) to the thermal pressure (P_th) is bound to less than 0.5. In this limit, the nonthermal energy is equal to half of the gravitational binding energy at the coronal radius without leaving room for thermal particles. Although the coronal plasma may be heated more through magnetic fields connected to the inner disk, we assume 0.5 as the maximal case in this work, and the neutrino spectrum peaks at ∼5 TeV and falls sharply around 20 TeV. We refer to this scenario as "High CR pressure."

We consider the second scenario for neutrino emission from NGC 1068 assuming coronal emission from stochastically accelerated particles, where instead of matching the flux at TeVs, we match the diffuse neutrino flux at tens of TeV, motivated by the medium-energy excess in the neutrino spectrum. In this case, as shown previously (Murase et al. 2020), we adopt parameters that can explain the high-energy neutrino flux excess observed at medium energies (Aartsen et al. 2020a). In this case, P_CR/P_th is set to ≃0.01. Here, the neutrino spectrum peaks at ∼40 TeV, which corresponds to a lower level of neutrino flux compared to the previous scenario. We refer to this case as "Modest CR pressure" hereafter.

These results are compatible with the spectra presented previously by Murase et al. (2020) where the CR pressure considered to explain the medium-energy neutrino flux and NGC 1068 are found at the level of ∼1% and ∼30% of the thermal pressure, respectively. Here, we allow the pressure ratio to be as high as 50% to explain the soft spectrum reported for NGC 1068 by the IceCube Collaboration (Aartsen et al. 2020b). Note that, in principle, both the High CR pressure and Modest CR pressure cases can be viable within the same stochastic acceleration scenario. For example, Modest CR pressure may be realized in an average AGN, whereas some sources such as NGC 1068 may have a large CR pressure.

Finally, we consider the magnetic reconnection scenario for particle acceleration. In this case, the neutrino flux approximately follows mainly the initial CR spectrum until the p γ process becomes the dominant channel for the production of pions. Therefore, this scenario leads to the spectrum having a shape close to that of a power-law spectrum with a cutoff at high energies. For the injected CR spectrum, we assume a spectral index of 2. The normalization and CR maximum energy are set such that the modeled flux is constrained to the IceCube steep spectrum reported for NGC 1068 while P_CR/P_th is bound to be smaller than 0.5. We find ${E}_{{\rm{p}}}^{{\rm{r}}{\rm{e}}{\rm{c}}}\approx 5\,{\rm{P}}{\rm{e}}{\rm{V}}$ for this purpose. Smaller values of ${E}_{{\rm{p}}}^{{\rm{r}}{\rm{e}}{\rm{c}}}$ cannot accommodate the IceCube flux without violating the CR to thermal pressure maximum band. Larger values, however, would create an excess at high energies that is disfavored by the steep spectra reported for NGC 1068. As described in Section 2, we set η_acc = 300 for magnetic reconnection acceleration. For NGC 1068, ${E}_{{\rm{p}}}^{{\rm{c}}{\rm{o}}{\rm{o}}{\rm{l}}}$ is too high to match the IceCube data.

Figure 1 shows the three modeled neutrino fluxes from NGC 1068. We also projected the best-fit spectrum reported by the IceCube Collaboration. The best-fit power-law spectrum corresponds to the ∼51 excess neutrinos found from the direction of NGC 1068. The shaded area shows the uncertainty on the fitted spectrum as reported by IceCube. As shown, all modeled neutrino spectra are within the 68% uncertainty of the measured spectrum. The parameters that we adapt in each scenario for particle acceleration and interaction efficiency are presented in Table 1. The common parameters among different scenarios are the same as in Murase et al. (2020). The injected CR, i.e., proton, differential luminosity for the three scenarios shown in Figure 1 is presented in the Appendix (see Figure 14).

Zoom In Zoom Out Reset image size

Table 1. Model-dependent Parameters and Quantities Adapted for Estimating the Neutrino Flux from NGC 1068 in each Scenario of Particle Acceleration shown in Figure 1

Model	PCR/Pth (CR)	q	ηtur	s	ηacc	${E}_{{\rm{p}}}^{{\rm{r}}{\rm{e}}{\rm{c}}}$
Stochastic acceleration with High CR pressure	0.5 (0.03)	5/3	50	⋯	⋯	⋯
Stochastic acceleration with Modest CR pressure	0.008 (0.0009)	5/3	10	⋯	⋯	⋯
Magnetic Reconnection	0.5 (0.01)	⋯	⋯	2.0	300	5 PeV

Download table as:  ASCIITypeset image

We should note that a single power-law spectrum is not a realistic spectral energy distribution for neutrino emission from individual astrophysical objects. While neutrino and γ-ray spectra may, in general, reflect the initial CR spectrum, the shape of neutrino and γ-ray fluxes depends on the nature of the interaction, thresholds, and the opacity of the source. The neutrino spectra provided in this study take all of these into account. That said, the diffuse flux of high-energy neutrinos (or γ-rays) over a specific range of energies may be explained by a power law since the superposition of the individual sources would wash out the features.

We use the modeled neutrino spectra for NGC 1068 to compare with the findings of the IceCube 10 yr point-source study. In addition, we investigate the prospects for identification of each neutrino emission scenario in the next decade of IceCube operation.

In order to find the p-value for the observation of neutrinos from NGC 1068 over the background of atmospheric neutrinos, we calculate the number of signal neutrinos using the publicly available effective area for the IceCube point-source selection (Aartsen et al. 2017). We also estimate the expected number of background atmospheric neutrinos using the zenith-dependent atmospheric neutrino flux reported by Honda et al. (2007), assuming a resolution of 07. This is indeed larger than IceCube's nominal angular uncertainty of 05 at high energies (≳100 TeV). We chose this value because most neutrinos in our predicted spectra are found at the 1–10 TeV range. The energy distribution of the events per year corresponding to the three neutrino scenarios we consider here is shown in Figure 2.

Zoom In Zoom Out Reset image size

The energy distribution of the events is a convolution of the spectrum and the IceCube effective area in the direction of NGC 1068. The stochastic acceleration scenario with High CR pressure case results in an event distribution that peaks around 10 TeV. The Modest CR pressure case for the stochastic acceleration yields a lower rate, and the majority of the events are found with energies around 30 TeV. That said, the harder spectrum at higher energies in the magnetic reconnection scenario leads to a relatively flat distribution of events and creates an excess beyond 100 TeV where atmospheric neutrinos are scarce.

We present the p-value ⁵ for the observation of each scenario in 5–15 yr of operation of IceCube in Figure 3. We use the method described in ATLAS (2011) to estimate the statistical significance for identification of the neutrino emission; for details, see the Appendix. The stochastic acceleration scenario with High CR pressure, compatible with the IceCube 10 yr flux measurements, provides the most likely scenario for observation of NGC 1068 in IceCube, reaching to 5σ level in less than 12 yr of data. This is in accordance with the reported local p-value in the IceCube 10 yr point-source study. The magnetic reconnection scenario is expected to be identified at better than 3σ within the 10 yr, while nearly a decade more of observations might be required to establish a significance at the discovery level. Finally, the stochastic acceleration scenario with Modest CR pressure seems difficult to identify with enough significance. However, the prospects for this scenario may be conservative given uncertainties on the X-ray luminosity. As demonstrated in Figure 4, the neutrino flux would be higher for larger intrinsic X-ray luminosities.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The key parameter in modeling the neutrino emission from NGC 1068 is the intrinsic X-ray luminosity. The column density (N_H) for NGC 1068 is very high. Therefore, it is difficult to estimate the intrinsic X-ray luminosity, and measurements often carry large uncertainties. The dominant source of uncertainty in measuring the X-ray luminosity from NGC 1068 is the fraction of X-rays that are scattered into our line of sight (Janssen et al. 2015). X-ray luminosities as high as 10⁴⁴ erg s⁻¹ have been considered for NGC 1068. NuSTAR and XMM-Newton monitoring campaigns (Marinucci et al. 2016) found ${L}_{X}\simeq {7}_{-4}^{+7}\,\times {10}^{43}$ erg s⁻¹ with N_H ≃ 10²⁵ cm⁻². Here, we used a central value of 10⁴³ erg s⁻¹, which is a conservative choice. We further investigate the variation of the measured intrinsic X-ray luminosity for the first scenario. Figure 4 shows how the spectrum changes when the X-ray luminosity is varied between 10^42.8 and 10^43.8 erg s⁻¹. Here, we set the CR to thermal pressure ratio to the maximal value of 0.5. The larger values of X-ray luminosity will increase the flux by almost an order of magnitude, while the peak of the spectrum is slightly shifted toward higher energies. With a larger value of L_X , smaller values of η_tur and P_CR/P_th are demanded to match the IceCube data, which may help in relaxing the extreme CR production efficiency for the High CR pressure case.

As we discussed earlier, in the magnetic reconnection case, the neutrino spectrum follows the CRs' flat spectrum with a cutoff at very high energies. IceCube's 10 yr observation disfavors a hard spectrum and suggests that the flux is suppressed at energies beyond 100 TeV. We evaluate the neutrino detection rate using a neutrino spectrum of a simple power law with an exponential cutoff to constrain the normalization and cutoff energy for the neutrino flux. We present in Figure 5 the 2σ and 3σ exclusion regions for 10 years of IceCube observations. As such, magnetic reconnection scenarios with large normalization and/or large values of maximum energies are not preferred. The constraints on E_ν,cut will also be useful for testing the accretion shock scenario. The diffusive shock acceleration theory predicts t_acc ≈ η_acc r_L/c with ${\eta }_{{\rm{a}}{\rm{c}}{\rm{c}}}\approx (20/3){(c/{V}_{{\rm{s}}})}^{2}$, where V_s is the shock velocity. This leads to ${E}_{{\rm{p}}}^{max}\sim {\rm{a}}\,{\rm{f}}{\rm{e}}{\rm{w}}\,$ PeV in NGC 1068 (Inoue et al. 2020). Thus, our results imply that fast acceleration models including both magnetic reconnection and accretion shock scenarios can be critically tested with near-future neutrino observations.

Zoom In Zoom Out Reset image size

We now use the parameters of the modeled neutrino flux from NGC 1068 to estimate neutrino spectra and prospects for observations of other bright Seyfert galaxies. We first focus on IceCube as it is currently the major operating neutrino telescope. We will extend our study to KM3NeT and the next generation of IceCube, IceCube-Gen2, in the next section.

For this purpose, we rank Seyfert galaxies by their intrinsic X-ray fluxes measured by the BAT AGN Spectroscopic Survey (Ricci et al. 2017). We select the 10 brightest sources for our study and extract their intrinsic X-ray luminosities from the BASS catalog. Table 2 presents the list of sources and their decl., distance, redshift, intrinsic X-ray flux, and luminosity in the 2–10 keV band. In order to estimate the luminosity distance for a source, we incorporate the redshift reported in BASS together with the cosmological parameters introduced in the Introduction, with the exception of Cen A, Circinus Galaxy, and NGC 4945, for which we employ astronomical measurements for the evaluation of the distance as reported by Harris et al. (2010), Karachentsev et al. (2013), and Tully et al. (2013). We should note that the neutrino flux is insensitive to the distance, because the neutrino flux is roughly proportional to the X-ray flux.

We implement the parameters that we adopted in the previous section for modeling the neutrino spectrum from NGC 1068, listed in Table 1. Again, we consider three scenarios: the stochastic acceleration scenario compatible with the 10 yr IceCube observation of NGC 1068 (High CR pressure), the stochastic acceleration scenario compatible with the medium-energy excess in the diffuse neutrino flux (Modest CR pressure), and the magnetic reconnection scenario. Similar to NGC 1068, we impose the physical constraint on the ratio of CR pressure to the thermal pressure of 0.5. Therefore, the normalization of the CR spectrum is adjusted to maintain this bound whenever this ratio exceeds 0.5 with the NGC 1068 parameters. Figure 6 demonstrates the spectrum for the source list considered in this study. For the stochastic acceleration scenario with High CR pressure, we have adjusted the flux to match the ratio of CR pressure to the thermal pressure to 0.5 by rescaling f_inj, which is used in our analyses hereafter. The measured intrinsic X-ray luminosity and the distance of the source define the shape and magnitude of the spectra. In the stochastic acceleration scenario, sources with L_X ≳ 10^43.5 erg s⁻¹ peak at relatively higher energies with a narrow width while sources with a smaller intrinsic X-ray luminosity demonstrate a broader spectrum peaking at lower energies. In this scenario, Cen A has the highest level of flux in the list of bright Seyfert galaxies. For the stochastic acceleration with Modest CR pressure, we employ f_inj value that yields a CR to thermal pressure ratio of ∼0.008 for NGC 1068. That corresponds to f_inj ≃ 10⁻⁷ (5 × 10⁻⁵) at 2m_p c² (1000m_p c²). Finally, in the magnetic reconnection scenario, the larger X-ray luminosity leads to a larger contribution of photohadronic processes. For each source, as the maximum energy in this scenario depends on the size of the system as well as the cooling effects, we evaluate ${E}_{{\rm{p}}}^{{\rm{r}}{\rm{e}}{\rm{c}}}$ by rescaling with M_BH as ${E}_{{\rm{p}}}^{{\rm{r}}{\rm{e}}{\rm{c}}}\approx {E}_{{\rm{p}}}^{{\rm{r}}{\rm{e}}{\rm{c}}-{\rm{N}}{\rm{G}}{\rm{C}}\,1068}{\left(M/{M}_{{\rm{N}}{\rm{G}}{\rm{C}}1068}\right)}^{1/2}$. This scaling is motivated by the model assumptions with l_rec/(2GMc⁻²) and β constant.

Zoom In Zoom Out Reset image size

The majority of bright nearby Seyfert galaxies in Table 2 are located in the Southern Hemisphere, as pointed out by Murase et al. (2020). IceCube's event selection is optimal for the Northern sky, where the Earth acts as a shield for the atmospheric muons. In the Southern Hemisphere, the event selection imposes a higher-energy threshold on the energy of the neutrinos to suppress the atmospheric muon background. This feature suppresses the event rate for the majority of the luminous Seyfert galaxies. We show the expected events from the sources in this list in Figure 15 in Appendix. Except for the sources NGC 1275, UGC 11910, and CGCG 164-019 that are in the Northern Hemisphere, the event rates for the rest of the sources are low, weakening the likelihood of identifying individual sources in IceCube.

Using the expected signal and background rates, we estimate the likelihood for observations of these sources in IceCube. Table 3 summarizes the expected p-values under each emission scenario for 10 years of IceCube operation. The listed p-values show that for all three acceleration scenarios, NGC 1068 is the brightest source in IceCube. While the prospects for the identification of most sources are not very promising due to the suppression of events in the Southern Hemisphere, with continued data collection CGCG 164-019 and NGC 1275 are likely to be observed at the 3σ level in 20 years of IceCube operations. However, we should note that the likelihood of observations depends on the neutrino emission scenario: the stochastic acceleration scenario with High CR pressure, compatible with NGC 1068 parameters, would yield ∼3σ.

Another source in the list worth discussing is NGC 4945. The IceCube 10 yr analysis found an excess of ∼1 event in its direction, corresponding to a p-value of 0.48, which is consistent with the expectations found in our study. NGC 4945 is a starburst galaxy. We should note that, similar to NGC 1068, the neutrino flux from this source cannot be explained by the starburst scenarios such as proposed by Eichmann & Becker Tjus (2016).

Other than NGC 1068, our predictions indicate that even optimistic scenarios are not strong enough to yield a statistically significant measurement of the neutrino emission from the rest of the bright nearby sources. As such, a stacking search for neutrino emission from the bright Seyfert galaxies is going to offer the best chance for identifying these sources in IceCube. Stacking analyses are widely used to study the correlation of the arrival direction of high-energy neutrinos and a catalog of sources. The cumulative signal from all of the sources improves the sensitivity and makes it possible to examine a class of sources with lower fluxes.

We project the prospects for identifying each scenario of neutrino emission in a stacking search for 5–20 yr of IceCube in Figure 7. The solid lines show the p-value for the duration of IceCube observations for each scenario. The most promising scenario is the stochastic acceleration scenario with High CR pressure case. The high-level neutrino flux in this scenario would result in an excess of signal that can reach a discovery level of 5σ in less than 10 years of IceCube operation. The hard spectrum in the magnetic reconnection scenario would result in an excess distinguishable from the background at high energies that can yield a 3σ significance in about 7 years of IceCube. The significance can reach a discovery level of 5σ within the 20 years of operation for this scenario.

Zoom In Zoom Out Reset image size

The most conservative scenario here is the stochastic acceleration scenario with Modest CR pressure. This scenario would not yield a signal at a significant level with the current generation of IceCube, mainly due to the fact that most sources are in the Southern Hemisphere. We investigate the likelihood of identifying this scenario in the next generation of neutrino telescopes in the next section.

We should note that the significance for the stacking analysis for the stochastic acceleration with High CR pressure scenario is dominated by the neutrino flux in the direction of NGC 1068. However, it is worth mentioning that, for this scenario, a stacking analysis of nearby bright Seyfert galaxies when NGC 1068 is excluded is expected to identify neutrino emission at the level of 3σ with 15 years of IceCube data. Such an analysis can test this scenario independently, without relying on the neutrino emission reported from the direction NGC 1068, which this scenario is based upon.

Our predicted neutrino emission was built upon the disk-corona model of AGNs. Two of the sources in the list of bright Seyfert galaxies, Cen A and NGC 1275, are seen with a high jet activity and it is likely that the X-ray emission arises from the jet rather than the disk. We will further investigate the prospects for observations of bright Seyfert galaxies if the disk-corona emission would not be dominant. The dashed lines in Figure 7 show the p-value when Cen A and NGC 1275 are not considered in the source list. While the likelihood of observation is decreased, the cumulative neutrino emission is strong and the stacking analysis of the rest of the bright Seyfert galaxies can still reach a significant level in the lifetime of the IceCube detector.

For Cen A, we should note that although one-zone models typically attribute the high energy to the jet (Abdo et al. 2010), it may be difficult to explain some of the X-ray properties with such a scenario (Tachibana et al. 2016). The observed soft lags in X-rays from Cen A may indicate their coronal origin. We should note that despite the large magnitude of the flux predicted under both stochastic acceleration and magnetic reconnection scenarios for Cen A in our study, our prediction is compatible with the current upper limits imposed by the ANTARES and IceCube Collaborations. In Figure 8 we compare the predicted neutrino flux with the current limits from the IceCube cascade source search in the Southern sky (Albert et al. 2020). The upper limit is obtained for a power law with an exponential cutoff, which is the closest shape to the predicted spectra in our study.

Zoom In Zoom Out Reset image size

In summary, observations of stacked neutrino emission from nearby, bright Seyfert galaxies is promising with IceCube, and could reveal the dominant sources responsible for the medium-energy excess in the spectrum of high-energy cosmic neutrinos. The likelihood for observations of these sources will be enhanced by the improvements in the event selection in the Southern Hemisphere. We will address this in Section 6.

Highly magnetized and turbulent coronae can be possible sites of particle acceleration. The system is calorimetric in the sense that sufficiently high-energy CRs are depleted via hadronuclear and photohadronic interactions. The large magnitude of the neutrino flux at 10–100 TeV makes this scenario a primary candidate for the medium-energy neutrino flux observed in IceCube at the level of ${E}_{\nu }^{2}{{\rm{\Phi }}}_{\nu }\sim {10}^{-7}\,\mathrm{GeV}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1}\,{\mathrm{sr}}^{-1}$ (Murase et al. 2020). The diffuse flux mainly originates from AGNs at high redshifts (with z ∼ 1−2), which are too far to detect as individual sources. The contribution from local sources is small, but it is still of interest to evaluate their aggregated flux.

Figure 11 shows the individual (thin lines) and sum (thick line) of the neutrino fluxes from nearby, bright Seyfert galaxies for different acceleration scenarios considered in this study. We have divided the fluxes by 4π in order to compare with the total neutrino flux from the 6 yr cascade analysis of IceCube (Aartsen et al. 2020a). Overall, each scenario predicts the contribution of the cataloged nearby sources to the total neutrino flux at 10 TeV to be within 2%–10%.

Zoom In Zoom Out Reset image size

The stochastic acceleration scenario with Modest CR pressure would mainly contribute to the 10–100 TeV region. However, the High CR pressure case would generate a significant excess of the flux below 10 TeV. This region is hard to investigate with the overwhelming flux of atmospheric neutrinos, and detailed veto techniques are required to distinguish the flux at TeV energies with good accuracy. The magnetic reconnection scenario has the highest contribution to the flux at ≳100 TeV. Distinguishing this scenario from the one responsible for the flux above 100 TeV would be difficult because of the scarcity of the data at high energies. While the extension of the flux to higher energies could make it easier to distinguish this flux from the background, such a possibility is less motivated by the compelling signal from NGC 1068. We discuss this further later in the next subsection.

It has been already shown that Seyfert galaxies with a CR pressure at the level of 1%–10% of the thermal pressure at the virial temperature can explain the magnitude of the diffuse neutrino flux at medium energies (Murase et al. 2020).

One of the novel points in this work is that we suggest the magnetic reconnection scenario as an alternative source of CR acceleration. In this scenario, we expect a power-law spectrum for CRs. The diffuse neutrino flux in this scenario is estimated to be

$$\begin{eqnarray}&&\begin{array}{l}{E}_{\nu }^{2}{{\rm{\Phi }}}_{\nu }\sim {10}^{-7}\,\displaystyle \frac{\mathrm{GeV}}{{\mathrm{cm}}^{2}\,{\rm{s}}\,\mathrm{sr}}\,\left(\displaystyle \frac{2K}{1+K}\right)\\ \quad \times \ \left(\displaystyle \frac{10}{{{ \mathcal R }}_{p}}\right)\left(\displaystyle \frac{15{f}_{\mathrm{mes}}}{1+{f}_{\mathrm{BH}}+{f}_{\mathrm{mes}}}\right)\\ \quad \times \ \left(\displaystyle \frac{{\xi }_{z}}{3}\right){\xi }_{\mathrm{CR}}\left(\displaystyle \frac{{L}_{X}{\rho }_{X}}{2\times {10}^{47}\,\mathrm{erg}\,{\mathrm{Mpc}}^{-3}\,{\mathrm{yr}}^{-1}}\right),\end{array}\end{eqnarray} \tag{ 3 }$$

where K = 1(2) for p γ(pp) interactions, ξ_z ∼ 3 summarizes the redshift evolution of the AGN luminosity density (Waxman & Bahcall 1999; Murase & Waxman 2016), ${{ \mathcal R }}_{p}$ is the conversion factor from bolometric to differential luminosities, ξ_CR = L_CR/L_X = (_CR/_rad)(L_disk/L_X ) is the CR loading factor defined against the X-ray luminosity, L_X is the X-ray luminosity, and ρ_X is the local density of X-ray selected AGN. Therefore, the total flux is found to be 10⁻⁷ GeV cm⁻² s⁻¹ sr⁻¹.

In the magnetic reconnection scenario, we assumed an injected spectral index of 2 for the injected CRs. This assumption yields a relatively flat spectrum of high-energy cosmic neutrinos with a cutoff around 100 TeV for η_acc that could describe the flux reported for NGC 1068. In addition, depending on the intrinsic X-ray luminosity of the source, a significant contribution from p γ interactions would emerge before the cutoff. The normalization in this scenario is set by P_CR/P_th ≤ 0.5 in our analysis for NGC 1068. While this scenario yields a smaller flux in the 1–10 TeV energy range, the flux becomes dominant over the background at higher energies, providing an excess that could be distinguished over the background in a sufficient number of years of detector operations. A higher cutoff energy in the spectrum would be therefore constrained by the lack of observation of such a hard neutrino flux from these sources. Therefore, both magnetic reconnection and shock acceleration scenarios (Inoue et al. 2020) would be disfavored unless η_acc is fine-tuned.

We further investigate the expected flux in this scenario by considering a steeper CR spectrum. In Figure 12, we show the neutrino flux corresponding to injected CR indices of 2.3 and 2.5 in addition to the spectral index of 2 that we used before. As can be seen, when restricting the pressure of CRs to less than 50% of the thermal pressure, it is difficult to explain the measured flux NGC 1068 by the IceCube Collaboration. We should note that although a spectral index of ∼2.5 or softer is motivated for CR acceleration in magnetic reconnection (Ball et al. 2018a), such an index is obtained for electrons, not protons, and it has not been supported by simulations thus far. For magnetic reconnections, harder spectra with an index ∼1 may be motivated; see, e.g., Hoshino (2015) and Guo et al. (2016). However, this case would induce a high σ regime and a super-relativistic plasma, which is not compatible with our modeling that considers the low σ regime. Nevertheless, we have shown the case for s = 1 in Figure 12 for parameters that we adopted for NGC 1068 in Section 3.1. As shown, in this scenario, we can explain the IceCube data if we change the maximum energy.

Zoom In Zoom Out Reset image size

The differential CR luminosity for the three scenarios for particle acceleration in NGC 1068 are shown in Figure 14. The CR luminosity is deduced from the measured intrinsic X-ray luminosity by incorporating the escape and cooling timescales, the injection function for the CRs, and the volume of the coronal region. For more details, see the Supplemental Material of Murase et al. (2020).

Zoom In Zoom Out Reset image size

To calculate the number of background atmospheric neutrino events, we have integrated the atmospheric flux (Honda et al. 2007 and Sanuki et al. 2007) over an opening angle of Ω = π σ² around the direction of the source, where the angle σ is the angular resolution of the detector. Given that the majority of the events in our models appear below ∼10 TeV, we set the angular resolution for the IceCube and KM3NeT to 07 as specified in Figure 2 of Aartsen et al. (2017) and Figure 22 of Adrian-Martinez et al. (2016), respectively.

To calculate the number of cosmic neutrino events for the flux predicted in each scenario, we integrate the product of the neutrino spectrum and effective area over energy and time. For IceCube, we use the zenith-dependent effective area for the point-source analysis (Aartsen et al. 2017). Similarly, we use this effective area, together with the atmospheric neutrino flux of Honda et al. (2007) and Sanuki et al. (2007) to estimate the number of background events in the direction of each source within the opening angle of Ω = π σ². The signal and background event distributions for NGC 1068 are presented in Figure 2 in the main text and in Figure 15 for the rest of the bright Seyfert galaxies in Table 2. Given that the majority of the sources in Figure 15 are located in the Southern Hemisphere, the corresponding event distribution is limited by the strong event selection cuts imposed to remove the atmospheric muon background penetrating the detector.

Zoom In Zoom Out Reset image size

The event distribution for KM3NeT is presented in Figure 16. We incorporate the average effective area for upgoing events in KM3NeT according to Figure 19 of Adrian-Martinez et al. (2016) and apply the trigger efficiency for the point-source searches. To take into account the fraction of the duration for which the source is below the horizon for KM3NeT, we rescale the number of events by the visibility of the sources, given by their decl., as is defined in Figure 37 of Adrian-Martinez et al. (2016).

Zoom In Zoom Out Reset image size

We estimate the statistical significance for observing the sources using the analytic expression introduced by ATLAS (2011). This method has been used extensively to evaluate the prospects for observation of Galactic sources of high-energy neutrinos; see, e.g., Gonzalez-Garcia et al. (2014) and Halzen et al. (2017). In this method, the p-value of identifying signal events from a source is given by

$$\begin{eqnarray}&&{p}_{{\rm{v}}{\rm{a}}{\rm{l}}{\rm{u}}{\rm{e}}}=\displaystyle \frac{1}{2}\left[1-{\rm{e}}{\rm{r}}{\rm{f}}\left(\sqrt{{q}_{0}^{{\rm{o}}{\rm{b}}{\rm{s}}}/2}\right)\right],\end{eqnarray} \tag{ A1 }$$

where q₀ ^obs is defined as

$$\begin{eqnarray}&&\begin{array}{r}\begin{array}{l}{q}_{0}^{{\rm{o}}{\rm{b}}{\rm{s}}}\equiv -2{\rm{l}}{\rm{n}}{{ \mathcal L }}_{b,D}\\ \,=\,2\displaystyle \sum _{i}\left({Y}_{b,i}-{N}_{D,i}+{N}_{D,i}{\rm{l}}{\rm{n}}\left(\displaystyle \frac{{N}_{D,i}}{{Y}_{b,i}}\right)\right),\end{array}\end{array}\end{eqnarray} \tag{ A2 }$$

with i running over the different energy bins. Here, Y_b,i is the theoretical expectation for the background hypothesis, while N_D,i is the estimated signal generated as the median of events Poisson distributed around the signal plus background. The background expectation here is the number of atmospheric neutrinos in the solid angle set by the angular resolution. This solid angle corresponds to roughly 72% of the signal events from the source; see Alexandreas et al. (1993) for a discussion.