Hidden Cores of Active Galactic Nuclei as the Origin of Medium-Energy Neutrinos: Critical Tests with the MeV Gamma-Ray Connection

TL;DR

This paper proposes a model linking the cores of active galactic nuclei to the origin of medium-energy neutrinos, predicting correlated gamma-ray emissions and suggesting observational tests with current and future telescopes.

Contribution

The study introduces a phenomenological model connecting AGN coronae to neutrino production, emphasizing the role of plasma turbulence and Bethe-Heitler processes in multimessenger emissions.

Findings

Protons in AGN coronae can be accelerated to PeV energies.

The model explains the diffuse flux of 10-100 TeV neutrinos with a few percent of thermal energy in cosmic rays.

A correlation between MeV gamma-ray background and 10 TeV neutrinos is predicted.

Abstract

Mysteries about the origin of high-energy cosmic neutrinos have deepened by the recent IceCube measurement of a large diffuse flux in the 10-100 TeV range. Based on the standard disk-corona picture of active galactic nuclei (AGN), we present a phenomenological model enabling us to systematically calculate the spectral sequence of multimessenger emission from the AGN coronae. We show that protons in the coronal plasma can be stochastically accelerated up to PeV energies by plasma turbulence, and find that the model explains the large diffuse flux of medium-energy neutrinos if the cosmic rays carry only a few percent of the thermal energy. We find that the Bethe-Heitler process plays a crucial role in connecting these neutrinos and cascaded MeV gamma rays, and point out that the gamma-ray flux can even be enhanced by the reacceleration of secondary pairs. Critical tests of the model are given by its prediction that a significant fraction of the MeV gamma-ray background correlates with about 10 TeV neutrinos, and nearby Seyfert galaxies including NGC 1068 are promising targets for IceCube, KM3Net, IceCube-Gen2, and future MeV gamma-ray telescopes.

Tables

Table 1. Table S1: Physical quantities deduced in our model (for ℛ = 30 ℛ 30 {\mathcal{R}}=30 ) as a function of L X subscript 𝐿 𝑋 L_{X} . Units are [erg s − 1 superscript s 1 \rm s^{-1} ] for L X subscript 𝐿 𝑋 L_{X} and L disk subscript 𝐿 disk L_{\rm disk} , [ M ⊙ subscript 𝑀 direct-product M_{\odot} ] for M 𝑀 M , and [ cm − 3 superscript cm 3 \rm cm^{-3} ] for n p subscript 𝑛 𝑝 n_{p} . Note that θ e ≡ k B ​ T e / ( m e ​ c 2 ) subscript 𝜃 𝑒 subscript 𝑘 𝐵 subscript 𝑇 𝑒 subscript 𝑚 𝑒 superscript 𝑐 2 \theta_{e}\equiv k_{B}T_{e}/(m_{e}c^{2}) is the normalized electron temperature.

$\log L_X$	$\log L_{\mathrm{disk}}$	$\log M$	${\Gamma }_X$	${\theta }_e$	${\tau }_T$	$\log n_p$
42.0	43.0	6.51	1.72	0.27	0.59	10.73
43.0	44.2	7.25	1.80	0.23	0.52	9.93
44.0	45.4	8.00	1.88	0.20	0.46	9.13
45.0	46.6	8.75	1.96	0.16	0.41	8.33
46.0	47.9	9.49	2.06	0.12	0.36	7.53

Equations

$$f_{pp}$$

$$f_{p\gamma }$$

$$f_{\mathrm{BH}}$$

$$E_{\nu }^2{\Phi }_{\nu }$$

$$t_{\mathrm{fall}}\approx \frac{\mathcal{R}{R}_S}{\alpha V_K},$$

$$t_{\mathrm{diss}}\approx \frac{H}{V_A}\approx \frac{\mathcal{R}{R}_S}{V_K}\sqrt{\frac{\beta }{2}},$$

$$t_{C,ee}\approx \frac{4\sqrt{\pi }{\theta }_e^{3/2}}{n_p{\sigma }_{T}c\ln \Lambda },$$

$$t_{C,pp}\approx \frac{4\sqrt{\pi }{\theta }_p^{3/2}}{n_p{\sigma }_{T}c\ln \Lambda }{\left(\frac{m_p}{m_e}\right)}^2,$$

$$t_{C,pe}\approx \sqrt{\frac{\pi }{2}}\frac{{\left({\theta }_p+{\theta }_e\right)}^{3/2}}{n_p{\sigma }_{T}c\ln \Lambda }\left(\frac{m_p}{m_e}\right),$$

$$t_{\mathrm{acc}}\approx \eta {\left(\frac{c}{V_A}\right)}^2\frac{H}{c}{\left(\frac{{\epsilon }_p}{eBH}\right)}^{2-q},$$

$$t_{\mathrm{diff}}\approx \frac{{\eta }^{\prime }}{\eta }\frac{H}{c}{\left(\frac{eBH}{{\epsilon }_p}\right)}^{2-q},$$

$$t_{pp}^{-1}\approx n_p{\kappa }_{pp}{\sigma }_{pp}c,$$

$$t_{p\gamma }^{-1}\approx \frac{c}{2{\gamma }_p^2}{\int }_{{\overline{\epsilon }}_{\mathrm{th}}}^{\infty }d{\overline{\epsilon }}_{\gamma }{\sigma }_{p\gamma }{\kappa }_{p\gamma }{\overline{\epsilon }}_{\gamma }{\int }_{{\overline{\epsilon }}_{\gamma }/(2{\gamma }_p)}^{\infty }d{\epsilon }_{\gamma }{\epsilon }_{\gamma }^{-2}\frac{d{n}_{\gamma }}{d{\epsilon }_{\gamma }},$$

$$t_{p,\mathrm{syn}}\approx \frac{6\pi m_p^4{c}^3}{m_e^2{\sigma }_T{B}^2{\epsilon }_p}.$$

$$\frac{\partial {\mathcal{F}}_p}{\partial t}=\frac{1}{{\epsilon }_p^2}\frac{\partial }{\partial {\epsilon }_p}\left({\epsilon }_p^2{D}_{{\epsilon }_p}\frac{\partial {\mathcal{F}}_p}{\partial {\epsilon }_p}+\frac{{\epsilon }_p^3}{t_{p-\mathrm{cool}}}{\mathcal{F}}_p\right)-\frac{{\mathcal{F}}_p}{t_{\mathrm{esc}}}+{\dot{\mathcal{F}}}_{p,\mathrm{inj}},$$

$${\dot{\mathcal{F}}}_{p,\mathrm{inj}}=\frac{f_{\mathrm{inj}}{L}_X}{4\pi {({\epsilon }_{\mathrm{inj}}/c)}^3\mathcal{V}}\delta ({\epsilon }_p-{\epsilon }_{\mathrm{inj}}),$$

$$t_{e,\mathrm{syn}}\approx \frac{6\pi m_{e}c}{{\sigma }_T{B}^2{\gamma }_e},$$

$$t_{\mathrm{IC}}\approx \frac{3m_{e}c}{4{\sigma }_T{U}_{\gamma }{\gamma }_e},$$

$${\Phi }_{\nu }$$

Full text

The origin of cosmic neutrinos observed in IceCube is a major enigma Aartsen et al. (2013a, b), and the latest data of high- and medium-energy starting events and shower events Aartsen et al. (2015a, b, 2020a) are more puzzling. The atmospheric background of high-energy electron neutrinos is lower than that of muon neutrinos, allowing us to analyze the data below 100 TeV Beacom and Candia (2004); Laha et al. (2013). The extragalactic neutrino background (ENB) at these energies has shown a larger flux with a softer spectrum, compared to the $\gtrsim 100$ TeV data Aartsen et al. (2016); Haack et al. (2017). The comparison with the extragalactic gamma-ray background (EGB) measured by Fermi indicates that the 10–100 TeV ENB originates from hidden sources preventing the escape of GeV–TeV gamma rays Murase et al. (2016).

Active galactic nuclei (AGN) are major contributors to the energetics of high-energy cosmic radiations Murase and Fukugita (2019); radio quiet (RQ) AGN are dominant in the extragalactic x-ray sky Fabian and Barcons (1992); Ueda et al. (2003); Hasinger et al. (2005); Ajello et al. (2008); Ueda et al. (2014), and jetted AGN that are typically radio loud (RL) dominantly explain the EGB Costamante (2013); Inoue (2014); Fornasa and Sánchez-Conde (2015). AGN may also explain the MeV gamma-ray background whose origin has been under debate (e.g., Refs. Inoue et al. (2008); Ajello et al. (2009); Lien and Fields (2012)).

High-energy neutrino production in the vicinity of supermassive black holes (SMBHs) were discussed early on Eichler (1979); Berezinskii and Ginzburg (1981); Begelman et al. (1990); Stecker et al. (1991), in particular to explain x-ray emission by cosmic-ray (CR) induced cascades assuming the existence of high Mach number accretion shocks at the inner edge of the disk Kazanas and Ellison (1986); Zdziarski (1986); Sikora et al. (1987); Stecker et al. (1991). However, cutoff features evident in the x-ray spectra of Seyfert galaxies and the absence of electron-positron annihilation lines ruled out the simple cascade scenario for the x-ray origin (e.g., Refs. Madejski et al. (1995); Ricci et al. (2018)). In the standard disk-corona scenario, the observed x rays are attributed to thermal Comptonization of disk photons Shapiro et al. (1976); Sunyaev and Titarchuk (1980); Zdziarski et al. (1996); Poutanen and Svensson (1996); Haardt et al. (1997), and electrons are presumably heated in the coronal region Liang and Price (1977); Galeev et al. (1979). There has been significant progress in our understanding of accretion disks with the identification of the magnetorotational instability (MRI) Balbus and Hawley (1991, 1998), which can result in the formation of a corona above the disk as a direct consequence of the accretion dynamics and magnetic dissipation (e.g., Refs. Miller and Stone (2000); Merloni and Fabian (2001); Liu et al. (2002); Blackman and Pessah (2009); Io and Suzuki (2014); Suzuki and Inutsuka (2014); Jiang et al. (2014)).

Accompanied turbulence and magnetic reconnections are important for particle acceleration Lazarian et al. (2012). The roles of nonthermal particles have been studied in the context of radiatively inefficient accretion flows (RIAFs) Narayan and Yi (1994); Yuan and Narayan (2014), in which the plasma is often collisionless because Coulomb collisions are negligible for protons (e.g., Refs. Takahara and Kusunose (1985); Mahadevan et al. (1997); Mahadevan and Quataert (1997); Kimura et al. (2014); Lynn et al. (2014); Ball et al. (2018)). Recent studies based on numerical simulations of the MRI Kimura et al. (2016, 2019a) support the idea that high-energy ions may be accelerated in the presence of the magnetohydrodynamic (MHD) turbulence.

The vicinity of SMBHs is often optically thick to GeV–TeV gamma rays, so that CR acceleration 111In this work CRs are used for nonthermal ions and nucleons that do not have to be observed on Earth. cannot be directly probed by these photons, but high-energy neutrinos can be used as a unique probe of the physics of AGN cores. In this work, we present a new concrete model for these high-energy emissions (see Fig. 1). Spectral energy distributions (SEDs) are constructed from the data and from empirical relations, and then we compute neutrino and cascade gamma-ray spectra by consistently solving particle transport equations. We demonstrate the importance of future MeV gamma-ray observations for revealing the origin of IceCube neutrinos especially in the medium-energy ($\sim 10-100$ TeV) range and for testing neutrino emission from NGC 1068 and other AGN.

We use a notation with $Q_{x}=Q\times 10^{x}$ in CGS units.

*Phenomenological prescription of AGN disk coronae.—*We begin by providing a phenomenological disk-corona model based on the existing data. Multiwavelength SEDs of Seyfert galaxies have been extensively studied, consisting of several components; radio emission (see Ref. Panessa et al. (2019)), infrared emission from a dust torus Netzer (2015), optical and ultraviolet components from an accretion disk Koratkar and Blaes (1999), and x rays from a corona Sunyaev and Titarchuk (1980). The latter two components are relevant for this work.

The “blue” bump, which has been seen in many AGN, is attributed to multitemperature blackbody emission from a geometrically thin, optically thick disk Shakura and Sunyaev (1973). The averaged SEDs are provided in Ref. Ho (2008) as a function of the Eddington ratio, $\lambda_{\rm Edd}=L_{\rm bol}/L_{\rm Edd}$, where $L_{\rm bol}$ and $L_{\rm Edd}\approx 1.26\times{10}^{45}~{}{\rm erg}~{}{\rm s}^{-1}{(M/10^{7}M_{\odot})}$ are bolometric and Eddington luminosities, respectively, and $M$ is the SMBH mass. The disk component is expected to have a cutoff in the ultraviolet range. Hot thermal electrons in a corona, with an electron temperature of $T_{e}\sim 10^{9}$ K, energize the disk photons by Compton upscattering. The consequent x-ray spectrum can be described by a power law with an exponential cutoff, in which the photon index ($\Gamma_{X}$) and the cutoff energy ($\varepsilon_{X,\rm cut}$) can also be estimated from $\lambda_{\rm Edd}$ Trakhtenbrot et al. (2017); Ricci et al. (2018). Observations have revealed the relationship between the x-ray luminosity $L_{X}$ and $L_{\rm bol}$ Hopkins et al. (2007) [where one typically sees $L_{X}\sim(0.01-0.1)L_{\rm bol}$], by which the disk-corona SEDs can be modeled as a function of $L_{X}$ and $M$. In this work, we consider contributions from AGN with the typical SMBH mass for a given $L_{X}$, using $M\approx 2.0\times 10^{7}~{}M_{\odot}~{}(L_{X}/1.16\times 10^{43}~{}{\rm erg}~{}{\rm s}^{-1})^{0.746}$ Mayers et al. (2018). The resulting disk-corona SED templates in our model are shown in Fig. 2 (see Supplemental Material for details), which enables us to quantitatively evaluate CR, neutrino and cascade gamma-ray emission.

Next we estimate the nucleon density $n_{p}$ and coronal magnetic field strength $B$. Let us consider a corona with the radius $R\equiv{\mathcal{R}}R_{S}$ and the scale height $H$, where ${\mathcal{R}}$ is the normalized coronal radius and $R_{S}=2GM/c^{2}$ is the Schwarzschild radius. Then the nucleon density is expressed by $n_{p}\approx\tau_{T}/(\sigma_{T}H)$, where $\tau_{T}$ is the Thomson optical depth that is typically $\sim 0.1-1$. The standard accretion theory Pringle (1981); Kato et al. (2008) gives the coronal scale height $H\approx(C_{s}/V_{K}){\mathcal{R}}R_{S}={\mathcal{R}}R_{S}/\sqrt{3}$, where $C_{s}=\sqrt{k_{B}T_{p}/m_{p}}=c/\sqrt{6\mathcal{R}}$ is the sound velocity, and $V_{K}=\sqrt{GM/R}=c/\sqrt{2\mathcal{R}}$ is the Keplerian velocity. For an optically thin corona, the electron temperature is estimated by $T_{e}\approx\varepsilon_{X,\rm cut}/(2k_{B})$, and $\tau_{T}$ is empirically determined from $\Gamma_{X}$ and $k_{B}T_{e}$ Ricci et al. (2018). We expect that thermal protons are at the virial temperature $T_{p}={GMm_{p}/(3{\mathcal{R}}R_{S}k_{B})}=m_{p}c^{2}/(6{\mathcal{R}}k_{B})$, implying that the corona may be characterized by two temperatures, i.e., $T_{p}>T_{e}$ Di Matteo et al. (1997); Cao (2009). Finally, the magnetic field is given by $B=\sqrt{8\pi n_{p}k_{B}T_{p}/\beta}$ with plasma beta ($\beta$).

Many physical quantities (including the SEDs) can be estimated observationally and empirically. Thus, for a given $L_{X}$, parameters characterizing the corona (${\mathcal{R}}$, $\beta$, $\alpha$) are remaining. They are also constrained in a certain range by observations Jin et al. (2012); Morgan et al. (2012) and numerical simulations Io and Suzuki (2014); Jiang et al. (2014). For example, recent MHD simulations show that $\beta$ in the coronae can be as low as 0.1–10 (e.g., Refs. Miller and Stone (2000); Suzuki and Inutsuka (2014)). We assume $\beta\lesssim 1-3$ and $\alpha=0.1$ for the viscosity parameter Shakura and Sunyaev (1973), and adopt ${\mathcal{R}}=30$.

*Stochastic proton acceleration in coronae.—*Standard AGN coronae are magnetized and turbulent, in which it is natural that protons are stochastically accelerated via plasma turbulence or magnetic reconnections. In this work, we solve the known Fokker-Planck equation that can describe the second order Fermi acceleration process (e.g., Refs. Becker et al. (2006); Stawarz and Petrosian (2008); Chang and Cooper (1970); Park and Petrosian (1996)). Here we describe key points in the calculations of CR spectra (see Supplemental Material or an accompanying paper Kimura et al. (2019b) for technical details). The stochastic acceleration time is given by $t_{\rm acc}\approx\eta{(c/V_{A})}^{2}(H/c){(\varepsilon_{p}/eBH)}^{2-q}$, where $V_{A}$ is the Alfvén velocity and $\eta$ is the inverse of the turbulence strength Dermer et al. (1996, 2014). We consider $q\sim 3/2-5/3$, which is not inconsistent with the recent simulations Kimura et al. (2019a), together with $\eta\sim 10$. The stochastic acceleration process is typically slower than the first order Fermi acceleration, which competes with cooling and escape processes. We find that for luminous AGN the Bethe-Heitler pair production ($p\gamma\rightarrow pe^{+}e^{-}$) is the most important cooling process because of copious disk photons, which determines the proton maximum energy. For our model parameters, the CR spectrum has a cutoff at $\varepsilon_{p}\sim 0.1-1$ PeV, leading to a cutoff at $\varepsilon_{\nu}\sim 5-50$ TeV in the neutrino spectrum. Note that all the loss timescales can uniquely be evaluated within our disk-corona model, and this result is only sensitive to $\eta$ and $q$ for a given set of coronal parameters. Although the resulting CR spectra (that are known to be hard) are numerically obtained in this work, we stress that spectra of $p\gamma$ neutrinos are independently predicted to be hard, because the photomeson production occurs only for protons whose energies exceed the pion production threshold Murase et al. (2016); Kimura et al. (2019b). The CR pressure to explain the neutrino data turns out to be $\sim(1-10)$% of the thermal pressure, by which the normalization of CRs is set.

For coronae considered here, the infall and dissipation times are $t_{\rm fall}\simeq 2.5\times{10}^{6}~{}{\rm s}~{}\alpha_{-1}^{-1}{({\mathcal{R}}/30)}^{3/2}R_{S,13.5}$ and $t_{\rm diss}\simeq 1.7\times{10}^{5}~{}{\rm s}~{}{(\mathcal{R}/30)}^{3/2}R_{S,13.5}\beta^{1/2}$, respectively. The Coulomb relaxation timescales for protons [e.g., $t_{C,pe}\sim 4\times{10}^{5}~{}{\rm s}~{}({\mathcal{R}}/30)R_{S,13.5}{(\tau_{T}/0.5)}^{-1}{(k_{B}T_{e}/0.1~{}{\rm MeV})}^{3/2}$] are longer than $t_{\rm diss}$ (especially for $\beta\lesssim 1$), so turbulent acceleration may operate for protons rather than electrons (and acceleration by small-scale magnetic reconnections may occur Hoshino (2015); Li et al. (2015)). This justifies our assumption on CR acceleration (cf. Refs. Ozel et al. (2000); Kimura et al. (2015); Ball et al. (2016); Kimura et al. (2019b) for RIAFs).

Connection between 10–100 TeV neutrinos and MeV gamma rays.— Accelerated CR protons interact with matter and radiation modeled in the previous section, producing secondary particles. We compute neutrino and gamma-ray spectra as a function of $L_{X}$, by utilizing the code to solve kinetic equations with electromagnetic cascades taken into account Murase (2018); Murase et al. (2019). Secondary injections by the Bethe-Heitler and $p\gamma$ processes are approximately treated as $\varepsilon_{e}^{2}(d\dot{N}_{e}^{\rm BH}/d\varepsilon_{e})|_{\varepsilon_{e}=(m_{e}/m_{p})\varepsilon_{p}}\approx t_{\rm BH}^{-1}\varepsilon_{p}^{2}(dN_{\rm CR}/d\varepsilon_{p})$ Chodorowski et al. (1992); Stepney and Guilbert (1983); Murase and Beacom (2010), $\varepsilon_{e}^{2}(d\dot{N}_{e}^{p\gamma}/d\varepsilon_{e})|_{\varepsilon_{e}=0.05\varepsilon_{p}}\approx(1/3)\varepsilon_{\nu}^{2}(d\dot{N}_{\nu}^{p\gamma}/d\varepsilon_{\nu})|_{\varepsilon_{\nu}=0.05\varepsilon_{p}}\approx(1/8)t_{p\gamma}^{-1}\varepsilon_{p}^{2}(dN_{\rm CR}/d\varepsilon_{p})$, and $\varepsilon_{\gamma}^{2}(d\dot{N}_{\gamma}^{p\gamma}/d\varepsilon_{\gamma})|_{\varepsilon_{\gamma}=0.1\varepsilon_{p}}\approx(1/2)t_{p\gamma}^{-1}\varepsilon_{p}^{2}(dN_{\rm CR}/d\varepsilon_{p})$. The cascade photon spectra are broad, being determined by the energy reprocessing via two-photon annihilation, synchrotron radiation, and inverse Compton emission.

The EGB and ENB are numerically calculated via the line-of-sight integral with the convolution of the x-ray luminosity function given by Ref. Ueda et al. (2014) (see also Supplemental Material, which includes Refs. Murase et al. (2014); Ajello et al. (2015); Murase et al. (2008a); Kotera et al. (2009); Fang and Murase (2018); Inoue (2011); Chen et al. (2015); Palladino and Winter (2018)). Note that the luminosity density of AGN evolves as redshift $z$, with a peak around $z\sim 1-2$, and our prescription enables us to simultaneously predict the x-ray background, EGB and ENB. The results are shown in Fig. S5, and our AGN corona model can explain the ENB at $\sim 30$ TeV energies with a steep spectrum at higher energies (due to different proton maximum energies), possibly simultaneously with the MeV EGB. We find that the required CR pressure ($P_{\rm CR}$) is only $\sim 1$% of the thermal pressure ($P_{\rm th}$), so the energetics requirement is not demanding in our AGN corona model (see Supplemental Material).

Remarkably, we find that high-energy neutrinos are produced by both $pp$ and $p\gamma$ interactions. The disk-corona model indicates $\tau_{T}\approx n_{p}\sigma_{T}{\mathcal{R}}R_{S}/\sqrt{3}\sim 0.1-1$, leading to the effective $pp$ optical depth

[TABLE]

where $\sigma_{pp}\sim 4\times{10}^{-26}~{}{\rm cm}^{2}$ is the $pp$ cross section, $\kappa_{pp}\sim 0.5$ is the proton inelasticity, and $V_{\rm fall}=\alpha V_{K}$ is the infall velocity. Coronal x rays provide target photons for the photomeson production, whose effective optical depth Murase et al. (2008b, 2016) for $\tau_{T}\lesssim 1$ is

[TABLE]

where $\eta_{p\gamma}\approx 2/(1+\Gamma_{X})$, $\hat{\sigma}_{p\gamma}\sim 0.7\times{10}^{-28}~{}{\rm cm}^{2}$ is the attenuation $p\gamma$ cross section, $\bar{\varepsilon}_{\Delta}\sim 0.3$ GeV, $\tilde{\varepsilon}_{p\gamma-X}=0.5m_{p}c^{2}\bar{\varepsilon}_{\Delta}/\varepsilon_{X}\simeq 0.14~{}{\rm PeV}~{}{(\varepsilon_{X}/1~{}{\rm keV})}^{-1}$, and $n_{X}\sim L_{X}/(2\pi R^{2}c\varepsilon_{X})$ is used. The total meson production optical depth is given by $f_{\rm mes}=f_{p\gamma}+f_{pp}$, which always exceeds unity in our model. Note that the spectrum of $p\gamma$ neutrinos should be hard at low energies, because only sufficiently high-energy protons can produce pions via $p\gamma$ interactions with x-ray photons.

Note that $\sim 10-100$ TeV neutrinos originate from $\sim 0.2-2$ PeV CRs. Unlike in previous studies explaining the IceCube data Stecker (2013); Kalashev et al. (2015), here in fact the disk photons are not much relevant for the photomeson production because its threshold energy is $\tilde{\varepsilon}_{p\gamma-\rm th}\simeq 3.4~{}{\rm PeV}~{}{(\varepsilon_{\rm disk}/10~{}{\rm eV})}^{-1}$. Rather, CR protons responsible for the medium-energy neutrinos should efficiently interact via the Bethe-Heitler process because the characteristic energy is $\tilde{\varepsilon}_{\rm BH-disk}\approx 0.5m_{p}c^{2}\bar{\varepsilon}_{\rm BH}/\varepsilon_{\rm disk}\simeq 0.47~{}{\rm PeV}~{}{(\varepsilon_{\rm disk}/10~{}{\rm eV})}^{-1}$, where $\bar{\varepsilon}_{\rm BH}\sim 10(2m_{e}c^{2})\sim 10$ MeV (Chodorowski et al., 1992; Stepney and Guilbert, 1983; Murase and Beacom, 2010). With the disk photon density $n_{\rm disk}\sim L_{\rm disk}/(2\pi R^{2}c\varepsilon_{\rm disk})$ for $\tau_{T}\lesssim 1$, the effective Bethe-Heitler optical depth (with $\hat{\sigma}_{\rm BH}\sim 0.8\times{10}^{-30}~{}{\rm cm}^{2}$) is

[TABLE]

which is much larger than $f_{p\gamma}$. The dominance of the Bethe-Heitler cooling is a direct consequence of the observed disk-corona SEDs. The 10–100 TeV neutrino flux is suppressed by $\sim f_{\rm mes}/f_{\rm BH}$, predicting the tight relationship with the MeV gamma-ray flux.

Analytically, the medium-energy ENB flux is given by

[TABLE]

which is indeed consistent with the numerical results shown in Fig. S5. Here $K=1$ and $K=2$ for $p\gamma$ and $pp$ interactions, respectively, $\xi_{z}\sim 3$ due to the redshift evolution of the AGN luminosity density Waxman and Bahcall (1998); Murase and Waxman (2016), ${\mathcal{R}}_{p}$ is the conversion factor from bolometric to differential luminosities, and $\xi_{\rm CR}$ is the CR loading parameter defined against the x-ray luminosity, where $P_{\rm CR}/P_{\rm th}\sim 0.01$ corresponds to $\xi_{\rm CR}\sim 0.1$ in our model. The ENB and EGB are dominated by AGN with $L_{X}\sim{10}^{44}~{}{\rm erg}~{}{\rm s}^{-1}$ Ueda et al. (2014), for which the effective local number density is $\rho_{X}\sim 5\times{10}^{-6}~{}{\rm Mpc}^{-3}$ Murase and Waxman (2016).

The $pp$, $p\gamma$ and Bethe-Heitler processes all initiate cascades, whose emission appears in the MeV range. Thanks to the dominance of the Bethe-Heitler process, AGN responsible for the medium-energy ENB should contribute a large fraction $\gtrsim 10-30$% of the MeV EGB.

When turbulent acceleration operates, the reacceleration of secondary pairs populated by cascades Murase et al. (2012) can naturally enhance the gamma-ray flux. The critical energy of the pairs, $\varepsilon_{e,\rm cl}$, is determined by the balance between the acceleration time $t_{\rm acc}$ and the electron cooling time $t_{\rm e-cool}$ (see Supplemental Material and Refs. Murase et al. (2012); Howes et al. (2011)). We find that the condition for the reacceleration is rather sensitive to $B$ and $t_{\rm acc}$. For example, with $\beta=3$ and $q=1.5$, the reaccelerated pairs can upscatter x-ray photons up to $\sim{(\varepsilon_{e,\rm cl}/m_{e}c^{2})}^{2}\varepsilon_{X}\simeq 3.4~{}{\rm MeV}~{}{(\varepsilon_{e,\rm cl}/30~{}{\rm MeV})}^{2}(\varepsilon_{X}/1~{}{\rm keV})$, which may lead to the MeV gamma-ray tail. This possibility is demonstrated in Fig. S5, and the effective number fraction of reaccelerated pairs is constrained as $\lesssim 0.1$%.

*Multimessenger tests.—*Our corona model robustly predicts $\sim 0.1-10$ MeV gamma-ray emission in either a synchrotron or an inverse Compton cascade scenario, without any primary electron acceleration (see Fig. 4). A large flux of 10–100 TeV neutrinos should be accompanied by the injection of Bethe-Heitler pairs in the 100–300 GeV range (see Supplemental Material for details) and form a fast cooling $\varepsilon_{e}^{-2}$ spectrum down to MeV energies in the steady state. In the simple inverse Compton cascade scenario, the cascade spectrum is extended up to a break energy at $\sim 1-10$ MeV, above which gamma rays are suppressed by $\gamma\gamma\rightarrow e^{+}e^{-}$. In reality, both synchrotron and inverse Compton processes can be important. The characteristic energy of synchrotron emission from Bethe-Heitler pairs is $\varepsilon_{\rm syn}^{\rm BH}\sim 1~{}{\rm MeV}~{}{B}_{2.5}{(\varepsilon_{p}/0.5~{}{\rm PeV})}^{2}$ Murase and Beacom (2010). Because disk photons lie in the $\sim 1-10$ eV range, the Klein-Nishina effect is important for the Bethe-Heitler pairs. Synchrotron cascades occur if the photon energy density is smaller than $\sim 10B^{2}/(8\pi)$, i.e., $B\gtrsim 170~{}{\rm G}~{}L_{\rm disk,45.3}^{1/2}{(\mathcal{R}/30)}^{-1}R_{S,13.5}^{-1}$.

The detectability of nearby Seyferts such as NGC 1068 and ESO 138-G001 is crucial for testing the model. MeV gamma-ray detection is promising with future telescopes like eASTROGAM De Angelis et al. (2017), GRAMS Aramaki et al. (2020), and AMEGO Moiseev et al. (2017), e.g., AMEGO’s differential sensitivity suggests that point sources with $L_{X}\sim 10^{44}~{}{\rm erg}~{}{\rm s}^{-1}$ are detectable up to $d\sim 70-150$ Mpc. At least a few of the brightest sources will be detected, and detections or nondetections of the MeV gamma-ray counterparts will support or falsify our corona model as the origin of $\sim 30$ TeV neutrinos. Interestingly, as demonstrated in Fig. 4, our corona model can explain the $\sim 3\sigma$ excess neutrino flux from NGC 1068 Aartsen et al. (2020b). It also predicts that the x-ray brightest Seyferts (that are more in the southern sky) can be detected as neutrino point sources by IceCube-Gen2 and KM3Net (see also Supplemental Material, which includes Refs. Aartsen et al. (2019, 2014); Adrian-Martinez et al. (2016)).

*Summary and discussion.—*We presented a new AGN corona model that can explain the medium-energy neutrino data. The observed disk-corona SEDs and known empirical relations enabled us to estimate model parameters, with which we solved particle transport equations and consistently computed subsequent electromagnetic cascades. Our coronal emission model provides clear, testable predictions relying on the neutrino–gamma-ray relationship that are largely independent of CR spectra. In particular, the dominance of the Bethe-Heitler pair production process is a direct consequence of the observed SEDs, leading to a natural connection with MeV gamma rays. Nearby Seyferts such as NGC 1068 and ESO 138-G001 will be promising targets for future MeV gamma-ray telescopes such as eASTROGAM and AMEGO. A good fraction of the MeV EGB may come from RQ AGN especially with secondary reacceleration, in which gamma-ray anisotropy searches should also be powerful Inoue et al. (2013). Neutrino multiplet Murase and Waxman (2016) and stacking searches with x-ray bright AGN are also promising, as encouraged by the latest neutrino source searches Aartsen et al. (2020b).

The proposed tests are crucial for unveiling nonthermal phenomena in the vicinity of SMBHs. In Seyferts and quasars, the plasma density is so high that a vacuum polar gap is diminished and GeV–TeV gamma rays are blocked. MeV gamma rays and neutrinos can escape and serve as a smoking gun of hidden CR acceleration that cannot be probed by x rays and lower-energy photons. Our results also strengthen the importance of further theoretical studies of disk-corona systems. Simulations on turbulent acceleration in coronae and particle-in-cell computations of magnetic reconnections are encouraged to understand the CR acceleration in such systems. Global MHD simulations will also be relevant to examine other speculative postulates such as accretion shocks Kazanas and Ellison (1986); Stecker et al. (1991); Szabo and Protheroe (1994); Stecker and Salamon (1996) or colliding blobs Alvarez-Muniz and Mészáros (2004), and to reveal the origin of low-frequency emission Inoue and Doi (2014, 2018).

I Supplemental Material

I.1 Disk-corona spectra and properties of coronae

Here we describe details on the modeling of disk-corona spectral energy distributions (SEDs) that are constructed from observational data and empirical relationships.

For a disk component, we use the averaged SEDs in Ref. Ho (2008), which are given as a function of the Eddington ratio, $\lambda_{\rm Edd}=L_{\rm bol}/L_{\rm Edd}$, where $L_{\rm bol}$ and $L_{\rm Edd}$ are bolometric and Eddington luminosities, respectively. The disk spectrum is expected to have a cutoff at $\varepsilon_{\rm disk,cut}\approx 2.8k_{B}T_{\rm disk}$, where $T_{\rm disk}\approx 0.49(GM\dot{M}/72\pi\sigma_{\rm SB}R_{S}^{3})^{1/4}$ is the maximum effective temperature of the disk near the central supermassive black hole (SMBH) (e.g., Ref. Pringle (1981)). Here, $M$ is the SMBH mass, $\dot{M}$ is the mass accretion rate, $R_{S}=2GM/c^{2}$ is the Schwarzschild radius, and $\sigma_{\rm SB}$ is the Stefan-Boltzmann constant. For a standard disk, one may use $\dot{M}\approx L_{\rm bol}/(\eta_{\rm rad}c^{2})=\lambda_{\rm Edd}L_{\rm Edd}/(\eta_{\rm rad}c^{2})$ with a radiative efficiency of $\eta_{\rm rad}=0.1$ Kato et al. (2008). Although the SEDs in Ref. Ho (2008) extend to low energies, we only consider disk photons with $\varepsilon_{\rm disk}>2$ eV because infrared photons would come from a dust torus. The ultraviolet and X-ray data are interpolated with an exponential cutoff taken into account, which gives conservative estimates on target photon densities for cosmic-ray interactions.

The coronal spectrum can be modeled by a power law with an exponential cutoff. The photon index, $\Gamma_{X}$, is correlated with $\lambda_{\rm Edd}$ as $\Gamma_{X}\approx 0.167\times\log(\lambda_{\rm Edd})+2.0$ (see Eq. 2 of Ref. Trakhtenbrot et al. (2017)), and the cutoff energy is given by $\varepsilon_{X,\rm cut}\sim[-74\log(\lambda_{\rm Edd})+1.5\times 10^{2}]$ keV Ricci et al. (2018). The x-ray luminosity $L_{X}$ can be converted into $L_{\rm bol}$ following Ref. Hopkins et al. (2007), by which the disk-corona SEDs are constructed as a function of $L_{X}$ and $M$.

To evaluate the cumulative neutrino and gamma-ray intensities from all active galactic nuclei (AGN) along the line-of-sight, we convolute the x-ray luminosity function with neutrino and gamma-ray spectra calculated as a function of $L_{X}$. In this work, instead of using the SMBH mass function, we assume the typical SMBH mass estimated by $M\approx 2.0\times 10^{7}~{}M_{\odot}~{}(L_{X}/1.16\times 10^{43}~{}{\rm erg}~{}{\rm s}^{-1})^{0.746}$ Mayers et al. (2018). The resulting SMBH mass $M$ and disk luminosity $L_{\rm disk}(\sim 0.5L_{\rm bol})$ are given in Table S1.

The x-ray cutoff energy is also used to estimate the coronal electron temperature $T_{e}\approx\varepsilon_{X,\rm cut}/(2k_{B})$. Ref. Ricci et al. (2018) gave $\tau_{T}\approx 10^{(2.16-\Gamma_{X})/1.06}{(k_{B}T_{e}/{\rm keV})}^{-0.3}$ in a slab geometry, by which the Thomson optical depth $\tau_{T}$ is determined. For a given SMBH mass, the coronal radius $R$ is parameterized by the normalized coronal radius $\mathcal{R}$ (that is set to 30 in this work). Then the nucleon density $n_{p}$ is determined by $\tau_{T}$ and ${\mathcal{R}}$. In Table S1, we summarize $\Gamma_{X}$, $kT_{e}$, $\tau_{T}$ and $n_{p}$, which are obtained in the model.

I.2 Timescales for thermal particles

To accelerate particles through stochastic acceleration in turbulence, the Coulomb relaxation timescales at their injection energy should be longer than the dissipation timescale. We discuss these plasma timescales in this subsection. The infall timescale is expected to be similar to that of the advection dominated accretion flow Narayan and Yi (1994); Yuan and Narayan (2014):

[TABLE]

where $\alpha$ is the viscous parameter in the accretion flow Shakura and Sunyaev (1973), ${\mathcal{R}}$ is the normalized coronal radius, and $R_{S}$ is the Schwarzschild radius. Assuming that the coronal dissipation is related to some magnetic process like reconnections, the dissipation timescale is expressed as

[TABLE]

which can be short especially for low $\beta$. The relaxation timescales for electrons and protons are Takahara and Kusunose (1985); Kimura et al. (2014)

[TABLE]

[TABLE]

[TABLE]

where $\theta_{i}=k_{B}T_{i}/(m_{i}c^{2})$ ($i=p$ or $e$), $\ln\Lambda\sim 20$ is the Coulomb logarithm, and we may consider a proton-electron plasma Ricci et al. (2018). We plot these timescales as a function of $L_{X}$ in Fig. S1 for $\beta=1$ and $\alpha=0.1$. We see that among the five timescales $t_{C,ee}$ and $t_{C,pp}$ are the shortest and longest, respectively. This means that electrons are easily thermalized while nonthermal protons are naturally expected and could be accelerated through stochastic acceleration. In order to form a two-temperature corona that is often discussed in the literature (e.g., Refs. Di Matteo et al. (1997); Cao (2009)), the dissipation timescale should be shorter than the proton-electron relaxation time. Because $t_{\rm diss}\lesssim t_{C,pe}$ is satisfied for the range of our interest, one may expect the two-temperature corona to be formed. Note that electrons in coronae would be heated by magnetic fields from the disk to explain $L_{X}$ and $kT_{e}$.

I.3 Timescales for high-energy protons

Nonthermal proton spectra are determined by the balance among particle acceleration, cooling, and escape processes. We consider stochastic acceleration by turbulence, and take account of infall and diffusion as escape processes. We also treat inelastic $pp$ collisions, photomeson production, Bethe-Heitler pair production, and synchrotron radiation as relevant cooling processes.

It is well known that stochastic acceleration is modeled as a diffusion phenomenon in momentum or energy space (e.g., Refs. Becker et al. (2006); Stawarz and Petrosian (2008); Lynn et al. (2014); Kimura et al. (2016, 2019a)). Assuming gyro-resonant scattering through turbulence with a power spectrum of $P_{k}\propto k^{-q}$, the acceleration time is written as Dermer et al. (1996); Murase et al. (2012); Dermer et al. (2014); Kimura et al. (2015)

[TABLE]

where $\eta^{-1}=8\pi\int dk\,P_{k}/B^{2}$ is the turbulence strength parameter, $V_{A}=B/\sqrt{4\pi m_{p}n_{p}}$ is the Alfvén velocity, and $R$ is the effective size of the coronal region. The infall time is given by Eq. (S1). Using the same scattering process for the stochastic acceleration, the diffusive escape time is estimated to be Dermer et al. (1996); Stawarz and Petrosian (2008); Kimura et al. (2015)

[TABLE]

where $\eta^{\prime}$ is a prefactor, which we set to $\eta^{\prime}=9$ Stawarz and Petrosian (2008) within uncertainties of our model with $H\lesssim R$. The cooling rate by $pp$ inelastic collisions is given by

[TABLE]

where $\sigma_{pp}$ and $\kappa_{pp}\approx 0.5$ are the cross section and inelasticity for $pp$ interactions, as implemented in Refs. Murase (2018); Murase et al. (2019). The photomeson production energy loss rate is calculated by

[TABLE]

where $\gamma_{p}=\varepsilon_{p}/(m_{p}c^{2})$ is the proton Lorentz factor, $\overline{\varepsilon}_{\rm{th}}\approx 145$ MeV is the threshold energy for the photomeson production, $\overline{\varepsilon}_{\gamma}$ is the photon energy in the proton rest frame, $\sigma_{p\gamma}$ and $\kappa_{p\gamma}$ are the cross section and inelasticity, respectively, and the normalization is given by $\int d\varepsilon_{\gamma}\varepsilon_{\gamma}(dn_{\gamma}/d\varepsilon_{\gamma})\approx(1+\tau_{T})(L_{\rm disk}+L_{X})/(2\pi R^{2}c)$, where $\mathcal{V}=(4\pi/3)r^{3}=2\pi{(\mathcal{R}R_{S})}^{2}H$ is the volume of the coronal region and $r/c$ is used for the photon and neutrino escape time. We utilize the fitting formula based on GEANT4 for $\sigma_{p\gamma}$ and $\kappa_{p\gamma}$, which are used in Ref. Murase et al. (2008b). The Bethe-Heitler energy loss rate ($t_{\rm BH}^{-1}$) is written in the same form of Eq. (S9) by replacing the cross section and inelasticity with $\sigma_{\rm BH}$ and $\kappa_{\rm BH}$, respectively, where we use the fitting formula given in Refs. Chodorowski et al. (1992) and Stepney and Guilbert (1983). Finally, the synchrotron timescale for protons is given by

[TABLE]

We plot the times scales in Fig. S2 with a parameter set of ${\mathcal{R}}=30$, $\beta=1$, $\alpha=0.1$, $\eta=10$, and $q=5/3$. We can see that particle acceleration is limited by interactions with photons except for $L_{X}=10^{42}\rm~{}erg~{}s^{-1}$. In the lowest-luminosity case, the photomeson production, the Bethe-Heitler process, the $pp$ reaction, and the diffusive escape rates are comparable to the acceleration rate around $10^{6}$ GeV. In the other cases, the Bethe-Heitler process hinders the acceleration and the maximum energy is reduced to $\sim 10^{5}-10^{6}$ GeV due to larger photon number densities of more luminous Seyferts.

I.4 Spectra of nonthermal protons

It is well known that the second order Fermi acceleration can be described by the Fokker-Planck equation (e.g., Refs. Becker et al. (2006); Stawarz and Petrosian (2008); Kimura et al. (2015)). To obtain spectra of nonthermal cosmic-ray (CR) protons, using the standard Chang-Cooper method Chang and Cooper (1970); Park and Petrosian (1996), we solve the following equation,

[TABLE]

where ${\mathcal{F}}_{p}$ is the CR distribution function, $D_{\varepsilon_{p}}\approx\varepsilon_{p}^{2}/t_{\rm acc}$ is the diffusion coefficient in energy space, $t_{p-\rm cool}^{-1}=t_{pp}^{-1}+t_{p\gamma}^{-1}+t_{\rm BH}^{-1}+t_{p,\rm syn}^{-1}$ is the total cooling rate, $t_{\rm esc}^{-1}=t_{\rm fall}^{-1}+t_{\rm diff}^{-1}$ is the CR escape rate (for the coronal region), and $\dot{\mathcal{F}}_{p,\rm inj}$ is the injection function that is given by

[TABLE]

where $\varepsilon_{\rm inj}$ is the injection energy and $f_{\rm inj}$ is the injection fraction at $\varepsilon_{\rm inj}$. When the normalization is set by the diffuse neutrino flux, we only need $P_{\rm CR}/P_{\rm th}\sim 0.1-1$%. Importantly, this CR energy fraction is only $\sim 1-10$% of the turbulent energy, which is energetically reasonable. The values of $\varepsilon_{\rm inj}$ and $f_{\rm inj}$ do not affect the resulting spectral shape as long as $\varepsilon_{\rm inj}\ll\varepsilon_{p,\rm cl}$. Note that we do not have to specify both of $\varepsilon_{\rm inj}$ and $f_{\rm inj}$, because $f_{\rm inj}$ is larger as $\varepsilon_{\rm inj}$ is higher. For $\varepsilon_{\rm inj}<\varepsilon_{p}\ll\varepsilon_{p,\rm cl}$, it is known that the stochastic acceleration mechanism predicts a very hard spectrum, $dN_{\rm CR}/d\varepsilon_{p}\propto\varepsilon_{p}^{1-q}$ (e.g., Refs. Becker et al. (2006); Stawarz and Petrosian (2008); Kimura et al. (2015)). All details are presented in our accompanying paper Kimura et al. (2019b).

We plot $\varepsilon_{p}L_{\varepsilon_{p}}\equiv 4\pi(\varepsilon_{p}^{4}/c^{3})F_{p}{\mathcal{V}}(t_{\rm esc}^{-1}+t_{p-\rm cool}^{-1})$ in Fig. S3. In Table S2, we tabulate the critical energy for CR protons, $\varepsilon_{p,\rm cl}$, at which the stochastic acceleration balances with energy and escape loss processes. For lower values of $L_{X}$, the critical energy is higher owing to their lower loss rates. We consider two cases, $q=5/3$ and $q=3/2$, which give similar CR proton spectra. For $L_{X}\gtrsim 10^{43}\rm~{}erg~{}s^{-1}$, owing to inefficient escape, the accelerated protons pile up around $\varepsilon_{p}\sim\varepsilon_{p,\rm cl}$, which creates a hardening feature around $\varepsilon_{p}\sim 0.1\varepsilon_{p,\rm cl}$ and a strong cutoff above $\varepsilon_{p}\sim 2\varepsilon_{p,\rm cl}$. Note that the resulting spectra of $p\gamma$ neutrinos are rather insensitive to such hard CR spectra, because they cannot be harder than the neutrino spectrum from pion decay Murase et al. (2016). The energy density of CRs is typically lower than that of the thermal protons: $P_{\rm CR}/P_{\rm th}$ is of the order of $10^{-2}$ for all the cases that allow us to explain the diffuse neutrino flux (see Table S2), where $P_{\rm CR}=\int\,d\varepsilon_{p}(dn_{\rm CR}/d\varepsilon_{p})\varepsilon_{p}/3$ and $P_{\rm th}=n_{p}k_{B}T_{p}$. These CRs do not affect the dynamical structure Kimura et al. (2014).

I.5 Timescales for high-energy pairs

Even if primary electrons are not accelerated through the turbulence due to their efficient Coulomb losses, sufficiently high-energy electron-positron pairs that are injected via hadronic processes and populated through electromagnetic cascades can be accelerated to higher energies without suffering from Coulomb losses, as studied in the literature of gamma-ray bursts Murase et al. (2012). In the meantime, such high-energy pairs also cool down by synchrotron radiation and inverse Compton emission. In Seyferts, the energy of the dominant target photons for inverse Compton scatterings is around 1–10 eV, implying that the Klein-Nishina effect is unimportant at $\varepsilon_{e}\lesssim 100-300$ GeV. If the secondary pairs can be accelerated by the turbulence, the acceleration and cooling processes are balanced at the critical energy, $\varepsilon_{e,\rm cl}$, and this effect is relevant if $\varepsilon_{e,\rm cl}$ is higher than the energy of the thermal protons. The timescale of the turbulent reacceleration is given by Eq. (S6) above the thermal energy of protons, $k_{B}T_{p}\simeq 5.2~{}{({\mathcal{R}}/30)}^{-1}$ MeV. Below this energy, the turbulent power spectrum should become steeper due to kinetic effects Howes et al. (2011), and the reacceleration time is considerably longer than that by Eq. (S6).

The electron synchrotron cooling timescale is

[TABLE]

where $\gamma_{e}=\varepsilon_{e}/(m_{e}c^{2})$ is the electron Lorentz factor. Then, the electron inverse Compton cooling time in the Thomson limit is estimated to be

[TABLE]

where $U_{\gamma}$ is the target photon energy density. Note that the Klein-Nishina effect is accurately taken into account in the calculations of cascaded photon spectra Murase (2018); Murase et al. (2019); Kimura et al. (2019b).

We plot timescales for high-energy pairs in Fig. S4, where thin (thick) lines are for the case with reacceleration allowed (forbidden), respectively. In the case of no secondary reacceleration, we have $\varepsilon_{e,\rm cl}\ll 5$ MeV, which is the critical energy is lower than the thermal energy of protons. The waves are expected to be dissipated at the relevant scales, and the stochastic acceleration of pairs is unlikely in all the range of $L_{X}$. In the other case allowing reacceleration, electrons and positrons injected via hadronic cascades can be maintained with energies around $\sim 10-30$ MeV, depending on $L_{X}$ (see Table S2), and the reaccelerated pairs can naturally enhance the MeV gamma-ray emission. To see the effect of the reacceleration, we take the approximate approach, in which we calculate electromagnetic cascades assuming that the pairs do not cool down below the critical energy that is determined by the balance between the acceleration time and cooling time. For $L_{X}=10^{42}\rm~{}erg~{}s^{-1}$, the synchrotron cooling is more likely to be important, while for $L_{X}\gtrsim 10^{43}\rm~{}erg~{}s^{-1}$, the inverse Compton cooling is dominant. The critical electron energy, $\varepsilon_{e,\rm cl}$, at which the cooling and reacceleration balance with each other is lower for a higher value of $L_{X}$ due to the more efficient inverse Compton cooling. Only when the critical energy is higher than the thermal proton temperature, the secondary reacceleration can occur. As shown in the main text, the effective number fraction of the reaccelerated pairs is constrained to be less than $\sim 0.1$% not to overshoot the MeV extragalactic gamma-ray background (EGB). This can be consistent with a prediction of stochastic acceleration (see, e.g., Ref. Stawarz and Petrosian (2008)). If pairs above $\varepsilon_{e,\rm reacc}\sim{\rm a~{}few}$ MeV can be reaccelerated, the number of pairs in the steady state is reduced by $\sim t_{e-\rm cool}(\varepsilon_{e,\rm reacc})/t_{\rm fall}$ (e.g., $\sim 0.05$% for $L_{X}=10^{43}\rm~{}erg~{}s^{-1}$) compared to the total number of injected pairs during $t_{\rm fall}$. The reacceleration is more difficult in luminous AGN due to small values of $t_{\rm IC}$.

I.6 Extragalactic neutrino background

The extragalactic neutrino background (ENB) flux (mean neutrino intensity) is calculated by (e.g., Ref. Murase et al. (2014))

[TABLE]

where $d\rho_{X}/dL_{X}$ is the x-ray luminosity function of AGN per comoving volume per luminosity, $L_{E^{\prime}_{\nu}}(L_{X})$ is the neutrino luminosity per energy, and $z_{\rm max}$ is the maximum value of redshift $z$. Also, $H_{0}$ is the local Hubble constant, whereas $\Omega_{m}=0.3$ and $\Omega_{\Lambda}=0.7$ are cosmological parameters. The differential neutrino luminosity, including spectral information, is calculated as a function of the x-ray luminosity, following our model described in this work. The x-ray luminosity function of AGN can be described by a smoothly connected double power-law model, and we use Eq. (14) and Table 4 of Ref. Ueda et al. (2014). The luminosity function evolves as redshift, for which Eqs. (16)–(19) and Table 4 of Ref. Ueda et al. (2014) are adopted. It is known that the x-ray luminosity density peaks around a redshift of $z\sim 1-2$ with a fast evolution. With the typical redshift evolution of x-ray AGN, the dimensionless parameter $\xi_{z}$, which represents the redshift evolution of the sources, gives $\xi_{z}\sim 3$ Murase and Waxman (2016); Waxman and Bahcall (1998).

The results for the ENB flux from AGN coronae are shown in Fig. 3 in the main text. We consider the sum of all three neutrino flavors. Note that these AGN are dominated by radio quiet (RQ) AGN that are mostly Seyfert galaxies with a standard accretion disk. On the other hand, some AGN are radio loud (RL), which are characterized by the radio loudness and attributed to jet activities. AGN jets are known to be powerful CR accelerators, and blazars and their misaligned counterparts – radio galaxies – account for the dominant fraction of the EGB in the 0.1–100 GeV range (e.g., Refs. Ajello et al. (2009, 2015)). Although the blazar contribution is likely to be subdominant in the neutrino sky Murase and Waxman (2016), RL AGN are among the most promising candidates as the sources of ultrahigh-energy cosmic rays (UHECRs). The CRs escaping from the jets can be confined in magnetized environments such as galaxy clusters and groups, generating high-energy neutrinos and gamma rays Murase et al. (2008a); Kotera et al. (2009). It is suggested that the jetted AGN embedded in clusters and groups give a “grand unified” explanation for UHECRs, sub-PeV neutrinos, and the nonblazar component of the EGB Murase and Waxman (2016); Fang and Murase (2018). Fig. S5 demonstrates that total contributions from AGN coronae, jetted AGN, and cosmogenic emission can explain the EGB from MeV to TeV energies as well as the ENB from TeV to PeV energies. The EGB from blazars Ajello et al. (2015) and radio galaxies Inoue (2011), and the ENB from galaxy clusters and groups with jetted AGN Fang and Murase (2018) are shown in Fig. S5. The latest IceCube data suggest a soft spectrum below 100 TeV with a large ENB flux Aartsen et al. (2020a), whereas the higher-energy data above 100 TeV can be described by a harder spectrum with $s\sim 2.2-2.3$ Haack et al. (2017). Although the current IceCube data do not exclude the single component scenario, they intriguingly hint at a structure in the neutrino spectrum Aartsen et al. (2020a), which is compatible with two- or multi-component models Chen et al. (2015); Palladino and Winter (2018). Indeed, the x-ray luminosity density of AGN is comparable to the CR luminosity of jetted AGN, which are $\sim 3\times{10}^{46}~{}{\rm erg}~{}{\rm Mpc}^{-3}~{}{\rm yr}^{-1}$ Murase and Fukugita (2019), and it is reasonable that the AGN coronae and the jetted AGN embedded in gaseous environments give comparable contributions to the ENB flux with different spectral indices. In addition, it was shown that the 10–100 TeV neutrino data can be explained only by hidden CR accelerators Murase et al. (2016), for which our AGN corona model provides a viable explanation and has reasonable energetics consistent with both theory and observations.

I.7 Detectability of the most promising sources

Our AGN corona model has a prediction for the relationship among the neutrino luminosity, soft-gamma-ray luminosity, and x-ray luminosity. AGN with higher “intrinsic” luminosities are better as the sources of high-energy neutrinos, and our model can be uniquely tested by future MeV gamma-ray satellites as well as neutrino experiments such as IceCube-Gen2 and KM3Net. In Fig. S6, we show the case for $L_{X}=10^{44}~{}{\rm erg}~{}{\rm s}^{-1}$ at $d=40$ Mpc (cf. ESO 138-G001). As one can see, the Bethe-Heitler process is crucial, and cascade gamma rays appear in the MeV range via electromagnetic cascades.

Stacking and cross correlation studies are also powerful, for which the x-ray catalogues, e.g., one from the Swift BAT AGN Spectroscopic Survey (BASS), should be useful BAS (BASS Project). Among Seyfert galaxies listed in BASS, based on the brightest x-ray objects in the 2–10 keV band, we find that the seven brightest targets are Circinus Galaxy, ESO 138-G001, NGC 7582, Cen A, NGC 1068, NGC 424, and CGCG 164-019. In particular, NGC 1068 and CGCG 164-019 are the brightest for upgoing muon neutrino observations by IceCube and IceCube-Gen2. On the other hand, the other Seyferts should be more interesting targets for KM3Net because they are located in the southern sky. Note that the merging galaxy NGC 6240 is also an interesting object, but the x-ray emission may not come from the AGN corona.