Soft gamma rays from low accreting supermassive black holes and connection to energetic neutrinos

First, we describe the properties of the thermal plasma in RIAFs (see subsection Emission from thermal electrons in RIAFs in Methods and Ref. [29] for technical details). Thermal electrons in RIAFs emit broadband photons through synchrotron radiation, bremsstrahlung, and inverse Compton scattering. For the parameter set shown in Table 1, values of the magnetic field, $B$, the Thomson optical depth, $\tau_{T}=n_{p}\sigma_{T}R$ ($n_{p}$ is number density, $\sigma_{T}$ is Thomson cross section, and $R$ is the size of accreting plasma), and the normalized electron temperature, $\Theta_{e}=k_{B}T_{e}/(m_{e}c^{2})$ ($m_{e}$ is electron mass, $k_{B}$ is Boltzmann constant, and $c$ is speed of light), are tabulated in Table 2. RIAFs are optically thick to synchrotron self-absorption at the synchrotron characteristic energy, $\varepsilon_{\rm syn,ch}\approx{3h_{p}eB\Theta_{e}^{2}}/(4\pi{m_{e}}c)$, where $h_{p}$ is Planck constant and $e$ is the elementary charge. The absorption energy above which RIAFs become optically thin is determined by equating the synchrotron emissivity to the blackbody radiation: $\varepsilon_{\rm syn,ab}\approx 3x_{M}h_{p}eB\Theta_{e}^{2}/(4\pi m_{e}c)\simeq 7.0\times 10^{-3}(B/{0.18\rm~{}kG})(\Theta_{e}/1.5)^{2}(x_{M}/10^{3})\rm~{}eV$, where $x_{M}=\varepsilon_{\rm syn,ab}/\varepsilon_{\rm syn,ch}\sim 10^{3}$ represents the correction from $\varepsilon_{\rm syn,ch}$ to $\varepsilon_{\rm syn,ab}$. See Ref. [30] for the values and the method to estimate $x_{M}$. The synchrotron luminosity from an optically thick source is given by

$$L_{\rm syn}\approx\frac{4\varepsilon_{\rm syn,ab}^{3}k_{B}T_{e}}{h_{p}^{3}c^{2}}\pi^{2}R^{2}\simeq 2.3\times 10^{40}{\rm~{}erg~{}s}^{-1}\left(\frac{M}{10^{8}\rm~{}M_{\odot}}\right)^{2}\left(\frac{\mathcal{R}}{10}\right)^{2}\left(\frac{B}{0.18\rm~{}kG}\right)^{3}\left(\frac{\Theta_{e}}{1.5}\right)^{7}\left(\frac{x_{M}}{10^{3}}\right)^{3},$$

(1)

where $M$ is the mass of the supermassive black hole (SMBH), $\mathcal{R}=R/R_{S}$, and $R_{S}=2GM/c^{2}$ is the Schwarzschild radius. The thermal electrons up-scatter the synchrotron photons, and the resulting spectrum can be approximated by a power-law form, $\varepsilon_{\gamma}L_{\varepsilon_{\gamma}}\propto\varepsilon_{\gamma}^{1-\alpha_{\rm IC}}$, where $\alpha_{\rm IC}=-\ln\tau_{T}/\ln{A_{\rm IC}}$ and $A_{\rm IC}=4\Theta_{e}+16\Theta_{e}^{2}$ [31]. If $\alpha_{\rm IC}<1$, the Comptonization dominates over other cooling processes, which is satisfied for $\dot{m}\gtrsim 2\times 10^{-3}$ as seen in Table 2. In this case, the luminosity of the Comptonized photons is

$$L_{\rm IC}\approx{L}_{\rm syn}\left(\frac{3k_{B}T_{e}}{\varepsilon_{\rm syn,ab}}\right)^{1-\alpha_{\rm IC}}\propto\Theta_{e}^{6+\alpha_{\rm IC}}.$$

(2)

Owing to the strong $\Theta_{e}$ dependence, our model predicts RIAF electron temperatures which lie in a narrow range $\Theta_{e}\sim 1-3$, unavoidably leading to a peak in the MeV range for $\dot{m}\gtrsim 2\times 10^{-3}$. The normalization is determined by the balance with the Coulomb heating:

$$L_{\rm IC}\approx(1-\alpha_{\rm IC})L_{\rm bol}\simeq 1.4\times 10^{42}{\rm~{}erg~{}s^{-1}}\left(\frac{M}{\rm 10^{8}~{}M_{\odot}}\right)\left(\frac{\dot{m}}{0.01}\right)^{2}\left(\frac{\eta_{\rm rad,sd}}{0.1}\right)\left(\frac{\alpha}{0.1}\right)^{2}\left(\frac{1-\alpha_{\rm IC}}{0.33}\right),$$

(3)

where $\dot{m}$ is the normalized accretion rate, $\eta_{\rm rad,sd}$ is radiation efficiency for a standard disk, and $\alpha$ is the viscous parameter (See subsection Emission from thermal electrons in RIAFs in Methods for the Coulomb heating rate and definition of some parameters).

The broadband photon spectra from thermal electrons are shown in Figure 1 for various values of $\dot{m}$. The synchrotron emission produces a peak at $0.001-0.01$ eV depending on $\dot{m}$, and the Comptonization of the synchrotron photons creates higher-energy photons up to 1–10 MeV. Cases with higher $\dot{m}$ have harder spectra because of their higher Thomson optical depths (see also Refs. [32, 33]), making their spectral peaks in the MeV range. These features are quantitatively consistent with the analytic estimates in Equations (1) and (3).

Our model is consistent with observations of nearby LLAGN. Ref. [34] reported a softening feature in the hard X-ray band in NGC 3998, from which they claimed that the electron temperature is $\simeq 30-40$ keV. Our RIAF model can reproduce the softening feature in the NuSTAR band, as well as the Swift BAT data shown in Fig. 2, despite a higher electron temperature. Ref. [34] also provided the X-ray spectrum for NGC 4579, which has a higher $\dot{m}$ and does not show any softening feature. Our model also produces a hard power-law spectrum consistent with the NuSTAR data (see Fig. 2). In our RIAF model, the resulting spectra for NGC 3998 and NGC 4579 are relatively hard, and well below the longer wavelength data (radio, infrared/optical/ultraviolet, and soft X-rays). These should be attributed to other emission components, such as compact jets or outer accretion disks [35]. Indeed, radio jets are observed in both objects [36, 37].

The Event Horizon Telescope Collaboration (EHT) reported a horizon-scale image of the SMBH in M87 [38]. Its brightness temperature is $\sim 5\times 10^{9}$ K (equivalent to $\Theta_{e}\sim 0.8$), while the real temperature should be $\Theta_{e}\simeq 3-10$, because the image is beam-smeared and the RIAF is likely optically thin at the observed frequency. Our model predicts $\Theta_{e}\simeq 3.5$ with the parameters appropriate for M87 ($M=6.3\times 10^{9}{\rm~{}M_{\odot}},~{}\dot{m}=6.1\times 10^{-4})$, which matches the expected temperature. The electron temperature in RIAFs also affects the interpretation of the photon ring observed by EHT [38, 39]. The emission region of the photon ring is determined by the electron temperature and magnetization, which should be clarified through the future multi-wavelength modeling of nearby LLAGN.

Protons in RIAFs are accelerated by magnetohydrodynamic (MHD) turbulence and/or magnetic reconnection generated by the magnetorotational instability (MRI). Here, we focus on the stochastic proton acceleration mechanism, in which non-thermal particles randomly gain or lose their energy via interactions with turbulent MHD waves. The accelerated protons, or cosmic-rays (CRs), produce neutrinos and gamma-rays via hadronuclear and photohadronic interactions with thermal protons and photons inside the RIAFs. Neutrinos freely escape from the system, whereas the gamma-rays create electron/positron pairs via $\gamma\gamma\rightarrow e^{+}e^{-}$, which initiates proton-induced electromagnetic cascades.

Figure 2 shows the resulting proton, neutrino, and proton-induced cascade gamma-ray spectra for NGC 3998 and NGC4579. The proton spectrum is hard because of the stochastic acceleration and has a cutoff around $10-100$ PeV due to photohadronic interactions. The neutrinos are mainly produced by $pp$ interactions for $\varepsilon_{\nu}\lesssim 10^{5}-10^{6}$ GeV, where $\varepsilon_{\nu}$ is the neutrino energy, while $p\gamma$ interactions are more efficient around the cutoff energy (See subsection Non-thermal particles in RIAFs in Methods for details). The resulting cascade gamma-ray spectrum is flat for $\varepsilon_{\gamma}<\varepsilon_{\gamma\gamma}$, where $\varepsilon_{\gamma}$ is the gamma-ray energy and $\varepsilon_{\gamma\gamma}$ is the cutoff energy above which gamma-rays are attenuated. The gamma-ray flux decreases rapidly above the energy. This feature is commonly seen in photon spectra from well-developed electromagnetic cascades [23]. The cascade gamma-ray spectra are well below the upper limit from the Fermi Large Area Telescope (LAT) [40] and the design sensitivity of future MeV satellites [41, 42, 43].

We calculate the gamma-ray and neutrino background intensities from the RIAFs in LLAGN (see subsection Cumulative background intensities in Methods for details). Fig. 3 shows the resulting gamma-ray and neutrino intensities from LLAGN. In our RIAF model, the Comptonized emission from the thermal electrons naturally accounts for the soft gamma-ray background in the 0.3–3 MeV range, below which canonical AGN coronae explain the X-ray background. Furthermore, non-thermal protons in RIAFs can simultaneously reproduce the IceCube data above $\sim 0.3$ PeV, below which hot coronae in luminous AGN can account for the neutrino data in the 10–100 TeV range. Hence, our synthesized AGN core scenario provides an attractive unified explanation for a wide energy range of the cosmic keV-MeV photon and TeV-PeV neutrino backgrounds, as demonstrated in Figure 3. Notably, AGN accretion flows can do this using a reasonable set of plasma parameters of the non-thermal proton component: $\beta\sim 1-10$, $\eta_{\rm tur}\sim 10-20$, and $P_{\rm CR}/P_{\rm th}\sim 0.01$ [23] (see Table 2). This value of $P_{\rm CR}/P_{\rm th}$ is reasonable in the sense that the CR energy density is lower than the magnetic field energy density. Also, the parameters related to non-thermal proton production, such as $\eta_{\rm tur}$ and $P_{\rm CR}/P_{\rm th}$, are expected to be similar in both AGN coronae and RIAFs, because they share the same CR acceleration and turbulence generation mechanisms. In this sense, the parameters of the non-thermal particles in our RIAF model are effectively calibrated by the 10–100 TeV neutrino data and the AGN corona model.

As shown in panel (a) of Figure 4, it is the relatively more luminous LLAGN that mainly contribute to the MeV intensity. In particular, LLAGN with $L_{\rm H\alpha}\sim L_{\rm crit}$ provide the dominant contribution, where $L_{\rm crit}$ is the ${\rm H\alpha}$ luminosity for $\dot{m}=\dot{m}_{\rm crit}$, where $\dot{m}_{\rm crit}$ is the critical mass accretion rate above which RIAFs no longer exist (See subsection Emission from thermal electrons in RIAFs in Methods for the value of $\dot{m}_{\rm crit}$). Thus, we can analytically estimate the gamma-ray background intensity to be [23]

		$\displaystyle E_{\gamma}^{2}\Phi_{\gamma}$	$\displaystyle\sim\frac{c}{4\pi H_{0}}\xi_{z}\rho_{*}\left(\frac{L_{*}}{L_{\rm crit}}\right)^{s_{1}-1}\varepsilon_{\gamma}L_{\varepsilon_{\gamma}}$
		$\displaystyle\sim$	$\displaystyle 3\times 10^{-6}{\rm~{}GeV~{}s^{-1}~{}cm^{-2}~{}sr^{-1}}\left(\frac{\xi_{z}}{0.6}\right)\left(\frac{\varepsilon_{\gamma}L_{\varepsilon_{\gamma}}}{45L_{\rm crit}}\right)\left(\frac{\dot{m}_{\rm crit}}{0.03}\right)^{-0.05},$

where $H_{0}\sim 70\rm~{}km~{}s^{-1}~{}Mpc^{-3}$ is the Hubble constant, $L_{*}$, $\rho_{*}$, and $s_{1}$ are the break luminosity, break density, and power-law index for the luminosity function [44] (see subsection Cumulative background intensities in Methods), respectively, and $\xi_{z}$ is the redshift evolution factor of the luminosity function. This is consistent with our numerical results and the observed MeV background from COMPTEL. We stress that the resulting MeV gamma-ray intensity does not depend strongly on the parameters related to the RIAF, such as $\alpha$, $B$, $M$, $\dot{m}_{\rm crit}$, nor on the electron heating prescription, as long as we can use the luminosity function and bolometric correlation. This is because LLAGN of $\dot{m}\sim\dot{m}_{\rm crit}$ provide the dominant contribution, and the energy budget of such LLAGN does not strongly depend on the value of $\dot{m}_{\rm crit}$.

On the other hand, the relatively faint LLAGN mainly contribute to the neutrino background (see panel (b) of Fig. 4). The high-energy gamma-rays accompanying the neutrinos are considerably attenuated inside RIAFs (see Table 2 for values of the cutoff energy due to $\gamma\gamma\rightarrow e^{+}e^{-}$ , $\varepsilon_{\gamma\gamma}$), making the GeV gamma-ray intensity well below the Fermi data. Thus, LLAGN can be regarded as gamma-ray hidden neutrino sources [45].

Our RIAF model provides a $\sim 5-10$% contribution to the cosmic soft X-ray (0.5–8 keV) background, which is dominated by canonical AGN and star-forming galaxies [46, 6]. In observations of distant LLAGN, fainter ones are likely to be classified as star-forming galaxies, while relatively luminous ones will be indistinguishable from the faint end of canonical AGN. Thus, it is difficult to determine the contribution by LLAGN accurately. Here, we make a rough estimate of the LLAGN contribution. The luminous objects of $L_{X}>3.2\times 10^{42}\rm~{}erg~{}s^{-1}$ provides 63% of the CXB flux. Our LLAGN should have X-ray luminosities below this range, and thus, 37% contribution can be regarded as an upper limit for the LLAGN contribution. Ref. [26] reported that 43% of nearby galaxies have signatures of AGN activity. Then, the LLAGN contribution can be as high as $0.43\times 0.37\simeq 0.16$. Our model predicts a lower LLAGN contribution than our rough estimate. Future X-ray and optical spectroscopic surveys with better sensitivities are necessary to unravel the X-ray contribution from LLAGN.

High-energy multimessenger observations are essential for identifying the origin of the gamma-ray and neutrino backgrounds. We have estimated the detectability of neutrinos from RIAFs (see subsection Neutrino detectability from nearby LLAGN in Methods for details). The expected number of through-going muon track events from nearby LLAGN, $N_{\mu}$, is shown in Figure 5. Current facilities are not sufficient to detect the signal, even using stacking techniques. However, stacking $\sim 10$ LLAGN will enable the planned IceCube-Gen2 facility to detect a few neutrino events by means of its larger effective area. Its better angular resolution will reduce the atmospheric background, making individual source detections possible.

Future MeV satellites, such as e-ASTROGAM [41], the All-sky Medium Energy Gamma-ray Observatory (AMEGO) [47], and Gamma-Ray and AntiMatter Survey (GRAMS) [43], will be able to measure the electron temperature in RIAFs by detecting a clear cutoff feature (see Fig. 2), which will shed light on the electron heating mechanisms in the collisionless plasma. In addition, anisotropy tests for the MeV gamma-ray background are promising [48]. Our RIAF model predicts a smaller anisotropy than those from Seyfert and blazar models, owing to the higher source number density of LLAGN.

If protons are accelerated by another mechanism, such as magnetic reconnection, the resulting proton spectrum is usually described by a power-law form. In the upper panel of Fig. 6, we plot the diffuse neutrino and gamma-ray intensities for power-law injection models with the parameter sets tabulated in Table 3 (see subsection Power-law injection models for CRs in Methods). The power-law injection models can also account for the PeV neutrino data with $P_{\rm CR}/P_{\rm th}\sim 0.01$ for model B and $P_{\rm CR}/P_{\rm th}\sim 0.1$ for model C. A higher $P_{\rm CR}/P_{\rm th}$ is demanded for a higher $s_{\rm inj}$, which may lead some feedback to the MHD turbulence (see discussion in the next subsection). The cascade gamma-ray emissions are well below the Fermi data owing to the low $\gamma\gamma$ cutoff energy in the RIAFs. Note that the resulting diffuse MeV gamma-ray intensities for models B and C are the same as that for model A, because the thermal electrons in the RIAFs are identical.

For the sake of completeness, we also investigate the possibility of whether our model can reproduce the 10–100 TeV neutrino background. The observed intensity of the 10–100 TeV neutrinos is higher than that of the 100 GeV gamma-rays. Hence, the emission region of origin of the 10–100 TeV neutrinos is expected to be opaque to GeV–TeV gamma-rays [45].

The lower panel of Fig. 6 depicts the diffuse gamma-ray and neutrino intensities for models D (stochastic acceleration), E (power-law injection with $s_{\rm inj}=1$), and F (power-law injection with $s_{\rm inj}=2$), whose parameters are tabulated in Table 3. Since the thermal component is identical to that for model A, we only show the non-thermal components.

The neutrino intensity is normalized so that it matches the observed 10–100 TeV neutrinos, which requires higher values of $\eta_{\rm CR}$ and $\eta_{\rm tur}$ or $\eta_{\rm acc}$. The gamma-ray intensity from proton-induced cascades is also higher than that of PeV neutrinos, but still considerably lower than the Fermi data. Therefore, our LLAGN model could in principle be the source of the mysterious 10-100 TeV neutrinos, although the contribution from Seyferts and quasars is typically more important in view of the energetics.

The resulting pressure ratio of the CRs to the thermal component is about 10% for models D and E and 50% for model F. These are much higher than those for models for PeV neutrinos. Nevertheless, even with such a high value of $P_{\rm CR}/P_{\rm th}$, the global dynamical structure of RIAFs are very similar as long as CRs are confined inside the flow as shown in Ref. [49]. However, such a high value of $P_{\rm CR}/P_{\rm th}$ is at odds with the picture of turbulent acceleration, and the feedback from CRs to MHD turbulence may be significant. In Ref. [29], we estimated the point-source detectability of the neutrinos from nearby LLAGN using the method in Ref. [50] with the parameter sets similar to those for models D, E, and F. The diffuse neutrino spectra for the models in Ref. [29] are almost identical to those for models D, E, and F. Ref. [29] showed that the TeV–PeV neutrinos and MeV gamma-rays are detectable with the planned neutrino and gamma-ray experiments, IceCube-Gen2 and AMEGO/eASTROGAM/GRAMS, respectively.

Here, we describe the properties of thermal plasma in RIAFs [24, 58] in detail. Hereafter, we use the notation of $Q_{X}=Q/10^{X}$ in cgs unit unless otherwise noted. We consider an accreting plasma of size $R$ around a supermassive black hole (SMBH) of mass $M$ with an accretion rate $\dot{M}$. We use the normalized radius and mass accretion rate, $\mathcal{R}=R/R_{S}$ and $\dot{m}=\dot{M}c^{2}/L_{\rm Edd}$, where $R_{S}$ is the Schwarzschild radius and $L_{\rm Edd}$ is the Eddington luminosity. Plasma quantities in RIAFs are described by two parameters, the viscosity parameter, $\alpha$, and the pressure ratio of gas to magnetic field, $\beta$. Based on recent numerical simulations (e.g., Refs. [59, 60, 61, 62, 63, 64]), the radial velocity, number density, proton thermal temperature, magnetic field, Thomson optical depth, and Alfvén velocity in the RIAF are analytically approximated to be

	$\displaystyle V_{R}\approx\alpha V_{K}/2\simeq 3.4\times 10^{8}~{}\mathcal{R}_{1}^{-1/2}\alpha_{-1}\rm~{}cm~{}s^{-1},$		(5)
	$\displaystyle n_{p}\approx\frac{\dot{M}}{4\pi m_{p}RHV_{R}}\simeq 4.6\times 10^{8}~{}\mathcal{R}_{1}^{-3/2}\alpha_{-1}^{-1}M_{8}^{-1}\dot{m}_{-2}\rm~{}cm^{-3},$		(6)
	$\displaystyle k_{B}T_{p}\approx\frac{GMm_{p}}{4R}\simeq 12~{}\mathcal{R}_{1}^{-1}\rm~{}MeV,$		(7)
	$\displaystyle B\approx\sqrt{\frac{8\pi n_{p}k_{B}T_{p}}{\beta}}\simeq 1.5\times 10^{2}~{}\mathcal{R}_{1}^{-5/4}\alpha_{-1}^{-1/2}M_{8}^{-1/2}\dot{m}_{-2}^{1/2}\beta_{1}^{-1/2}\rm~{}G,$		(8)
	$\displaystyle\tau_{T}\approx n_{p}\sigma_{T}R\simeq 0.090~{}\mathcal{R}_{1}^{-1/2}\dot{m}_{-2}\alpha_{-1}^{-1},$		(9)
	$\displaystyle\beta_{A}\approx\frac{B}{\sqrt{4\pi n_{p}m_{p}c^{2}}}\simeq 0.050~{}\mathcal{R}_{1}^{-1/2}\beta_{1}^{-1/2},$		(10)

where $V_{K}=\sqrt{GM/R}$ is the Keplerian velocity and $H\approx R/2$ is the scale height. Observations of X-ray binaries and AGN demand $\alpha\sim 0.1-1$ [65], while the global MHD simulations result in $\alpha\sim 0.01-0.1$ and $\beta\sim 3-30$ [61, 66]. Hence, we set $\alpha=0.1$ and $\beta=7.0$ as their reference values.

Thermal electrons in RIAFs emit broadband photons through synchrotron radiation, bremsstrahlung, and inverse Compton scattering. The electron temperature is determined so that the resulting photon luminosity is equal to the bolometric luminosity estimated by $\dot{m}$. Assuming that thermal electrons are heated by Coulomb collisions with protons, the bolometric luminosity is estimated to be $L_{\rm bol}\approx\eta_{\rm rad,sd}\dot{m}_{\rm crit}L_{\rm Edd}(\dot{m}/\dot{m}_{\rm crit})^{2}$, where $\eta_{\rm rad,sd}\sim 0.1$ is the radiation efficiency for the standard disk, and $\dot{m}_{\rm crit}\approx 0.03(\alpha/0.1)^{2}$ is the critical mass accretion rate above which RIAFs no longer exist [30]. We calculate the photon spectra by synchrotron radiation, bremsstrahlung, and inverse Compton scattering by thermal electrons with the steady-state and one-zone approximations. The synchrotron and bremsstrahlung spectra are calculated by the method given in Appendix of Ref. [27], where we use the fitting formulae for the emissivity of these processes and Eddington approximation to take into account the effects of the radiative transfer. For the inverse Compton scattering spectrum, we utilize the corrected delta-function method given in Ref. [67]. In this method, the distribution of the scattered photon energy is approximated to be a delta function, but the mean energy of the scattered photon is calculated using the exact kernel. This method approximately takes into account the electron recoil effect for $\varepsilon_{\gamma}\gtrsim k_{B}T_{e}$. The error of the method is about 50%. This is sufficiently accurate for our purpose, considering significant uncertainty in the MeV gamma-ray data. Given the uncertainty, our calculation results are consistent with those by the Monte Carlo simulations [68]. The spectral decline due to the cutoff in this method is somewhat stronger than that in the exact method, implying that our results on the MeV fluxes are regarded as conservative. In this work, we assume that the Coulomb heating is the dominant electron heating mechanism, and $L_{\rm bol}\propto\dot{m}^{2}$ is used. This treatment is qualitatively different from the previous work [27], where $L_{\rm bol}\propto\dot{m}$ is assumed. Such a treatment may be more appropriate if the electrons are directly heated by the plasma dissipation process  [69, 70, 71, 72]. Nevertheless, these details will not change our conclusions that LLAGN are bright in soft MeV gamma-rays. Also, the electron heating prescription do not strongly affect the critical mass accretion rate above which RIAFs no longer exist [73, 28, 30, 74].

Here, we describe the details of the stochastic acceleration model, where protons are accelerated by MRI turbulence [75, 64]. To obtain the CR spectrum, we solve the transport equation for CR protons, which is a diffusion equation in the momentum space [76]:

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

(11)

where $\mathcal{F}_{p}$ is the momentum distribution function for protons ($dN/d\varepsilon_{p}=4\pi p^{2}\mathcal{F}_{p}/c$), $D_{\varepsilon_{p}}$ is the diffusion coefficient that mimics the stochastic particle acceleration, $t_{\rm cool}$ is the cooling time for protons, $t_{\rm esc}$ is the escape time, $\dot{\mathcal{F}}_{p,\rm inj}=\dot{\mathcal{F}}_{0}\delta(\varepsilon_{p}-\varepsilon_{p,\rm inj})$ is the injection function, and $\varepsilon_{p,\rm inj}$ is the initial energy of the particles injected to the stochastic acceleration process. We assume that the particles are injected to the stochastic acceleration process by fast acceleration processes such as magnetic reconnections, which are induced by MRI [77, 78, 79]. We consider a delta-function injection term with $\varepsilon_{p,\rm inj}$ much higher than the thermal proton energy, which may mimic the injection by the reconnection. The value of $\varepsilon_{p,\rm inj}$ has no influence on the resulting spectrum as long as $\varepsilon_{p,\rm inj}$ is much lower than the cutoff energy. We consider resonant scatterings between MHD waves and CR particles, where CR particles interact with the turbulent eddy of their gyration radius, $r_{L}=\varepsilon_{p}/(eB)$. Then, diffusion coefficient in energy space can be written as

$$D_{\varepsilon_{p}}\approx\frac{c\beta_{A}^{2}}{\eta_{\rm tur}H}\left(\frac{r_{L}}{H}\right)^{q-2}\varepsilon_{p}^{2},$$

(12)

where $\eta_{\rm tur}=B^{2}/(8\pi\int P_{k}dk)$ is the turbulence parameter ($P_{k}$ is the turbulence power spectrum), $q$ is the power-law index of $P_{k}$, and we set the scale height, $H$, to the injection scale of the MHD turbulence. We assume a Kolmogorov turbulence, $q=5/3$. The acceleration time is given by $t_{\rm acc}\approx\varepsilon_{p}^{2}/D_{\varepsilon_{p}}\approx\eta_{\rm tur}\beta_{A}^{-2}(H/c)[\varepsilon_{p}/(eBH)]^{1/3}$. $\eta_{\rm tur}$ is lower for a stronger turbulence, and a low value of $\eta_{\rm tur}$ results in a higher maximum energy of the CR protons. Since RIAFs are expected to be turbulent due to MRI [80, 81], the value of $\eta_{\rm tur}$ should be small. The turbulent strength is also related to the value of $\alpha$, as stronger turbulence leads to a higher $\alpha$. Our fiducial value of $\eta_{\rm tur}\sim 15$ is reasonable in the sense that $\eta_{\rm tur}$ is close to $\alpha^{-1}$. The amount of CRs is determined so that $\int L_{\varepsilon_{p}}d\varepsilon_{p}=\eta_{\rm CR}\dot{m}L_{\rm Edd}$ is satisfied, where $L_{\varepsilon_{p}}=t_{\rm loss}^{-1}\varepsilon_{p}dN_{p}/d\varepsilon_{p}$ is the differential proton luminosity [23], $\eta_{\rm CR}$ is the CR production efficiency, and $t_{\rm loss}^{-1}=t_{\rm cool}^{-1}+t_{\rm esc}^{-1}$ is the total energy loss rate including cooling and escape processes.

We solve the transport equation until the steady state is reached using the Chang-Cooper method [82, 83]. As the proton cooling mechanism, we consider the proton synchrotron, Bethe-Heitler ($p+\gamma\rightarrow p+e^{+}+e^{-}$), photomeson ($p+\gamma\rightarrow p+\pi$), and $pp$ inelastic collision ($p+p\rightarrow p+p+\pi$) processes, The timescale of the $pp$ inelastic collisions is given as $t_{pp}^{-1}\approx n_{p}\sigma_{pp}\kappa_{pp}c$, where $\sigma_{pp}$ and $\kappa_{pp}$ are the cross section and inelasticity of $pp$ interactions [84]. The photomeson production and Bethe-Heitler cooling timescales are estimated to be

$$t_{i}^{-1}=\frac{c}{2\gamma_{p}^{2}}\int_{\varepsilon_{\rm th}}^{\infty}d\overline{\varepsilon}_{\gamma}\sigma_{i}\kappa_{i}\overline{\varepsilon}_{\gamma}\int_{\overline{\varepsilon}_{\gamma}/(2\gamma_{p})}^{\infty}\frac{d\varepsilon_{\gamma}}{\varepsilon_{\gamma}^{2}}n_{\varepsilon_{\gamma}},,$$

(13)

where $\gamma_{p}=\varepsilon_{p}/(m_{p}c^{2})$, $\overline{\varepsilon}_{\gamma}$ is the photon energy in the proton rest frame, and $\sigma_{i}$ and $\kappa_{i}$ are the cross section and inelasticity for the process ($i=p\gamma$ for photomeson production [85] and $i={\rm BH}$ for Bethe-Heitler process [86, 87]). The cooling time by the proton synchrotron is $t_{p,\rm syn}\approx 6\pi m_{p}^{4}c^{3}/(\sigma_{T}B^{2}\varepsilon_{p})$. The total cooling rate is then given by $t_{\rm cool}^{-1}=t_{pp}^{-1}+t_{p\gamma}^{-1}+t_{\rm BH}^{-1}+t_{p,\rm syn}^{-1}$. Regarding the escape process, we only consider the advective escape, i.e., infall to the SMBH. We write the escape timescale as $t_{\rm esc}=t_{\rm adv}\approx R/V_{R}$. We consider that the CR component has the same bulk velocity as that for the thermal component owing to efficient interactions through the turbulent magnetic field. The diffusive escape would be inefficient because the high-energy protons tend to move in the azimuthal direction, which is the direction of the background magnetic field in differentially rotating accretion flows [75, 64]. See also Refs. [23, 29] for technical details of the calculation methods for the cooling and escape timescales.

In Fig. 7, we plot the acceleration and loss rates of CRs for model A (reference model) for NGC 4579, whose parameters are shown in Table 1 and caption. The dominant loss process is the advective escape for $\varepsilon_{p}\lesssim 1\times 10^{8}$ GeV. Although the acceleration time is longer than the advective escape time for $\varepsilon_{p}\gtrsim 2\times 10^{5}$ GeV, the CR proton spectrum continues to a higher-energy because the weak $\varepsilon_{p}$ dependence of the advective escape leads to a very gradual cutoff in the proton spectrum [88, 27]. At $\varepsilon_{p}\sim 1\times 10^{7}$ GeV, the photomeson production becomes more efficient than the acceleration, which makes a sharp cutoff owing to a strong $\varepsilon_{p}$ dependence (see Fig. 2).

The CRs produce pions via both $pp$ and $p\gamma$ interactions, and pions decay to gamma-rays, electron/positron pairs, and neutrinos ($\pi^{0}\rightarrow 2\gamma$; $\pi^{\pm}\rightarrow e^{\pm}+3\nu$). These neutrinos are believed to explain IceCube neutrinos [89, 90, 91]. We calculate neutrino spectra by $pp$ collisions using the formalism given by Ref. [92]. For the neutrinos by $p\gamma$ interactions, we use a semi-analytic prescription given in Refs. [93, 94] (see, e.g., Ref. [85] for numerical results). Since the effect of meson cooling is negligible in the RIAFs, neutrino flavor ratio is $\nu_{e}:\nu_{\mu}:\nu_{\tau}=1:2:0$ at the source, which becomes $\sim 1:1:1:$ on Earth after the flavor mixing. As shown in Figure 2, neutrinos are mainly produced by the $pp$ collisions for $\varepsilon_{\nu}\lesssim 4\times 10^{5}$ GeV, and the photomeson production is effective above the energy. The Bethe-Heitler and proton synchrotron processes are subdominant in the range we investigated. For higher $\dot{m}$ cases, the photomeson production is more efficient, and hence, the peak energy of the proton spectrum is lower. In our model, the neutrino background is dominated by relatively faint LLAGN, and their target photon spectra are not hard. Then, the multi-pion production channel is subdominant, and our approximation provides reasonably accurate results.

The hadronic interactions also produce gamma-rays and electron/positron pairs. The gamma-rays are absorbed by two-photon annihilation, and create electron/positron pairs. These pairs also emit gamma-rays, and electromagnetic cascades are initiated. We calculate the cascade emission by solving the kinetic equations of electron/positron pairs and photons [95, 96]:

$\displaystyle\frac{\partial n_{\varepsilon_{e}}}{\partial t}+\frac{\partial}{\partial\varepsilon_{e}}\left[\left(P_{\rm IC}+P_{\rm syn}+P_{\rm ff}+P_{\rm Cou}\right)n_{\varepsilon_{e}}\right]=\dot{n}_{\varepsilon_{e}}^{(\gamma\gamma)}-\frac{n_{\varepsilon_{e}}}{t_{\rm esc}}+\dot{n}_{\varepsilon_{e}}^{\rm inj},$

(14)

$$\frac{\partial n_{\varepsilon_{\gamma}}}{\partial t}=-\frac{n_{\varepsilon_{\gamma}}}{t_{\gamma\gamma}}-\frac{n_{\varepsilon_{\gamma}}}{t_{\gamma,\rm esc}}+\dot{n}_{\varepsilon_{\gamma}}^{(\rm IC)}+\dot{n}_{\varepsilon_{\gamma}}^{(\rm ff)}+\dot{n}_{\varepsilon_{\gamma}}^{(\rm syn)}+\dot{n}_{\varepsilon_{\gamma}}^{\rm inj},$$

(15)

where $n_{\varepsilon_{i}}$ is the differential number density ($i=e$ or $\gamma$), $\dot{n}_{\varepsilon_{i}}^{(XX)}$ is the particle source term from the process $XX$ ($XX=\rm IC$ (inverse Compton scattering), $\gamma\gamma$ ($e^{+}e^{-}$ pair production by $\gamma\gamma$ interactions), syn (synchrotron), or ff (bremsstrahlung)), $\dot{n}_{\varepsilon_{i}}^{\rm inj}$ is the injection term from the hadronic interaction, and $P_{YY}$ is the energy loss rate for the electrons from the process $YY$ ($YY=\rm IC$ (inverse Compton scattering), syn (synchrotron), ff (bremsstrahlung), or Cou (Coulomb collision)). See also Refs. [95, 96] for technical details. We approximately treat the pair injection processes by the Bethe-Heitler process and photomeson production as in Refs. [23, 29].

In our RIAF models, secondary pairs do not contribute to the the high-energy gamma-ray background, even for the case that the secondary pairs are re-energized by MHD turbulence. The secondary pairs suffer from a strong cooling by the synchrotron emission, whose timescale is estimated to be $t_{e,\rm syn}=6\pi m_{e}^{2}c^{3}/(\sigma_{T}B^{2}\varepsilon_{e})$. When their energy becomes sufficiently low, they may be re-energized by MHD turbulence, as suggested by Ref. [23]. The energization timescale is the same with the proton acceleration timescale: $t_{e,\rm acc}\approx\eta_{\rm tur}\beta_{A}^{-2}(H/c)[\varepsilon_{e}/(eBH)]^{2-q}$. Equating these two timescales, we obtain the critical energy at which the secondary pairs piles up:

$\displaystyle\varepsilon_{e,\rm crit}\approx\left(\frac{6\pi m_{e}^{2}c^{4}\beta_{A}^{2}}{\sigma_{T}H\eta_{\rm tur}B^{2}}\right)^{\frac{1}{3-q}}\left(eBH\right)^{\frac{2-q}{3-q}}\simeq 4.2~{}\mathcal{R}_{1}^{5/16}\alpha_{-1}^{5/8}M_{8}^{1/8}\dot{m}_{-2}^{-5/8}\beta_{0.5}^{-1/8}\eta_{\rm tur,1.5}^{-3/4}\rm~{}MeV,$

(16)

where we use $q=5/3$ for the last equation. The turbulence power below the mean thermal proton energy, $3k_{B}T_{p}\sim 35\mathcal{R}^{-1}$ MeV, is significantly reduced due to dissipation by plasma kinetic effects, and the turbulent re-energization is not expected when $\varepsilon_{e,\rm crit}<3k_{B}T_{p}$. In our model A (fiducial parameters), we cannot expect re-energization of secondary pairs even for the case with $L_{\rm H\alpha}=L_{\rm min}$. For a further lower $\dot{m}$ case, re-energization may occur. In such a case, we can ignore the inverse Compton component by the secondary pairs owing to its lower photon energy density ($B^{2}\propto\dot{m}$, $L_{\rm bol}\propto\dot{m}^{2}$). Then, the synchrotron peak energy for the re-energized pairs are estimated to be

$$\varepsilon_{\gamma,\rm syn}\approx\frac{3h_{p}\varepsilon_{e,\rm crit}^{2}eB}{4\pi m_{e}^{3}c^{5}}\simeq 0.07~{}\left(\frac{\varepsilon_{e,\rm crit}}{100\rm~{}MeV}\right)^{2}B_{2}\rm~{}eV.$$

(17)

Hence, we conclude that the re-energized pairs cannot contribute to the MeV gamma-ray background for all the $\dot{m}$ range in our model. Also, primary electrons are not expected to be produced efficiently in RIAFs, because of their rapid thermalization in the range of our interest [28]. Thus, they do not contribute to the MeV gamma-ray background.

Here we describe the method to obtain the background intensities. Since the $\rm H\alpha$ luminosity functions include much fainter sources than the X-ray luminosity functions, we use the luminosity function for type-1 Seyfert galaxies provided by Ref. [44]: $d\rho/dL_{\rm H\alpha}\approx(\rho_{*}/L_{*})/[(L_{\rm H\alpha}/L_{*})^{s_{1}}+(L_{\rm H\alpha}/L_{*})^{s_{2}}]$, where $\rho_{*}\approx 1.2\times 10^{-6}\rm~{}Mpc^{-3}$, $L_{*}=3.7\times 10^{42}\rm~{}erg~{}s^{-1}$, $s_{1}=2.05$, and $s_{2}=5.12$ (the values of $L_{*}$ and $\rho_{*}$ are corrected using Hubble constant of $H_{0}=70\rm~{}km~{}s^{-1}~{}Mpc^{-1}$). We extrapolate this luminosity function to $L_{\rm min}=10^{38}\rm~{}erg~{}s^{-1}$, below which the Palomar survey finds a hint of a break [26]. The survey also indicates a correlation between X-ray luminosity, $L_{X}$, and $L_{\rm H\alpha}$ for LLAGN. The ratio, $\kappa_{X/\rm H\alpha}=L_{X}/L_{\rm H\alpha}$ ranges $5\lesssim\kappa_{X/\rm H\alpha}\lesssim 7$ in the luminosity range of our interest for type-1 AGN. We use $\kappa_{X/\rm H\alpha}=6.0$ for simplicity [26], but the difference from the cases with $\kappa_{X/\rm H\alpha}=5$ or $\kappa_{X/\rm H\alpha}=7$ is less than a factor of 1.2.

Observationally, the X-ray luminosity at the 2–10 keV band can be converted to the bolometric luminosity using the bolometric correction factor, $\kappa_{{\rm bol}/X}=L_{\rm bol}/L_{X}\simeq 15$ for LLAGN [97, 98, 99, 100]. Using the two correction factor, $\kappa_{{\rm bol}/X}$ and $\kappa_{X/\rm~{}H\alpha}$, we can convert $\dot{m}$ to $L_{\rm H\alpha}$ if we fix a SMBH mass, $M$. Ref. [101] provided a sample of LLAGN, and the mean and median values of $\log(M/\rm M_{\odot})$ are 8.0 and 8.1, respectively. Also, the X-ray luminosity density is dominated by AGN with $M\sim 10^{8}-3\times 10^{8}\rm~{}M_{\odot}$ if the Eddington ratio function is independent of the SMBH mass [97, 102, 103]. Thus, we set $M=10^{8}\rm~{}M_{\odot}$ as a reference value. With this prescription, the critical $\rm H\alpha$ luminosity above which RIAFs no longer exist is estimated to be $L_{\rm crit}=\eta_{\rm rad,sd}\dot{m}_{\rm crit}L_{\rm Edd}/(\kappa_{X/\rm H\alpha}\kappa_{{\rm bol}/X})\simeq 4.2\times 10^{41}\rm~{}erg~{}s^{-1}$. Finally, we integrate the gamma-ray and neutrino fluxes over the range of $L_{\rm min}\leq L_{\rm H\alpha}\leq L_{\rm crit}$, which corresponds to $4.6\times 10^{-4}<\dot{m}<0.03$. The background intensity is written as [104]

$$E_{i}^{2}\Phi_{i}=\frac{c}{4\pi H_{0}}\int dz\frac{1}{\sqrt{(1+z)^{3}\Omega_{m}+\Omega_{\Lambda}}}\int dL_{\rm H\alpha}\frac{d\rho}{dL_{\rm H\alpha}}E_{i}L_{E_{i}},$$

(18)

where we use the cosmological parameters of $\Omega_{m}\approx 0.3$ and $\Omega_{\Lambda}\approx 0.7$. Since dimmer AGN tend to have weaker redshift evolution based on the X-ray, gamma-ray, and radio luminosity functions [105, 103, 106], no redshift evolution of the $\rm H\alpha$ luminosity function is used. With this treatment, LLAGN at $z<0.5$ contribute about a half of the diffuse intensity, objects at $0.5<z<2$ provide the other half, and the contribution from $z>2$ is negligible. We consider the gamma-ray attenuation by the extragalactic background light, and include the exponential suppression with the optical depth given in Ref. [107]. The two-photon pair annihilation initiates intergalactic electromagnetic cascades [108], but it has little influence on the GeV gamma-ray intensity because TeV gamma-rays cannot escape from the RIAFs due to two-photon annihilation inside the RIAFs (see Table 2 for the value of $\varepsilon_{\gamma\gamma}$, the $\gamma\gamma$ cutoff energy above which two-photon annihilation attenuates the gamma-rays inside the RIAFs).

Ref. [44] also provided the luminosity function for type-2 AGN, but we find that their contributions to the MeV gamma-ray and PeV neutrino backgrounds are negligible because of two reasons. One is the value of $\kappa_{X/\rm H\alpha}$. Type-2 AGN have $\kappa_{X/\rm H\alpha}\sim 1$ [26], and neutrino and gamma-ray luminosity scales with $L_{\nu}\propto\dot{m}^{2}\propto L_{\rm bol}\propto\kappa_{X/\rm H\alpha}$. The other is the shape of the luminosity function. The break luminosity for type-2 AGN is lower than the critical luminosity, $L_{*}\approx 2.8\times 10^{40}{\rm~{}erg~{}s^{-1}}<L_{\rm crit}$, resulting in a lower MeV gamma-ray intensity. The slope in the low-luminosity range is lower than 2, $s_{1}=1.77$, which makes the contribution by relatively faint LLAGN smaller, leading to a lower PeV neutrino intensity. In addition, type-2 objects may be contaminated by non-AGN objects, the fraction of which is still uncertain. Therefore, we focus on contribution by type-1 AGN in this study.

Note that the $\rm H\alpha$ luminosity function does not match the X-ray luminosity function with a constant $\kappa_{X/\rm H\alpha}$. The conversion factor depends on the luminosity, and there is some degree of uncertainty in the conversion factor (see Ref. [52]). The $\rm H\alpha$ luminosity function by Ref. [52] matches the observed X-ray luminosity function at a high luminosity range, while it underestimates the number density for $L_{X}\lesssim 10^{42}\rm~{}erg~{}s^{-1}$. This range is the most relevant to the MeV gamma-ray and PeV neutrino backgrounds, and the luminosity function by Ref. [44] gives a higher number density at the range. From this reason, we use one by Ref. [44] as a reference case. Future optical spectroscopic and hard X-ray surveys are necessary to reduce this uncertainty.

The number of through-going muon track events is calculated by [109, 50]

$$N_{\mu}(>E_{\mu})=\int_{E_{\mu}}^{\infty}dE^{\prime}_{\mu}\frac{\mathcal{N}_{A}\mathcal{A}_{\rm det}}{\alpha_{\mu}+\beta_{\mu}E^{\prime}_{\mu}}\int_{E^{\prime}_{\mu}}^{\infty}dE_{\nu}\phi_{\nu}\sigma_{\rm CC}\exp(-\tau_{\nu N}),$$

(19)

where $E_{\nu}$ is the incoming neutrino energy, $E_{\mu}$ is the muon energy, $\phi_{\nu}=\Delta tL_{E_{\nu}}/(4\pi d_{L}^{2}E_{\nu})$ is the neutrino fluence from a LLAGN for a time interval of $\Delta t$, $\mathcal{N}_{A}$ is the Avogadro Number, $\sigma_{\rm CC}$ is the charged-current cross section, $\tau_{\nu N}$ is the optical depth for the neutrino scattering in the Earth, and $\alpha_{\mu}+\beta_{\mu}E_{\mu}$ represents the energy loss rate of muons. We use $\Delta t=10$ yr and the list of LLAGN given by Ref. [101] (see Ref. [29] for the brightest 10 LLAGN in the list). This method can reproduce the effective area by Ref. [110]. For IceCube-Gen2, we use a $10^{2/3}$ times bigger $A_{\rm det}$ than that for IceCube [111]. The main background of the astrophysical neutrino is atmospheric background, whose intensity depends on the declination. We appropriately take into account both conventional and prompt atmospheric muon neutrinos as well as the attenuation in the Earth [50].

The CR spectrum inside the RIAFs can be obtained by solving the transport equation with a power-law injection term:

$$\frac{d}{d\varepsilon_{p}}\left(-\frac{\varepsilon_{p}}{t_{\rm cool}}N_{\varepsilon_{p}}\right)=\dot{N}_{\varepsilon_{p},\rm inj}-\frac{N_{\varepsilon_{p}}}{t_{\rm esc}},$$

(20)

$$\dot{N}_{\varepsilon_{p},\rm inj}=\dot{N}_{0}\left(\frac{\varepsilon_{p}}{\varepsilon_{p,\rm cut}}\right)^{-s_{\rm inj}}\exp\left(-\frac{\varepsilon_{p}}{\varepsilon_{p,\rm cut}}\right),$$

(21)

where $\varepsilon_{p,\rm cut}$ is the cutoff energy for the injected protons and $\dot{N}_{0}$ is the normalization factor. The injection is normalized such that $\int\varepsilon_{p}\dot{N}_{\varepsilon_{p},\rm inj}d\varepsilon_{p}=\eta_{\rm CR}\dot{m}L_{\rm Edd}$ is satisfied. We estimate $\varepsilon_{p,\rm max}$ by equating the infall timescale to the acceleration timescale, $t_{\rm acc}$. We set $t_{\rm acc}=\eta_{\rm acc}r_{L}/c$, where $\eta_{\rm acc}$ is the acceleration efficiency parameter. The value of $\eta_{\rm acc}$ is quite uncertain, but $\eta_{\rm acc}\gg 1$ is expected in the RIAFs because $\eta_{\rm acc}\sim 1$ is achieved only for relativistic shocks or relativistic reconnections (i.e., magnetic energy density is higher than the rest mass energy). Here, we tune it so that our RIAF models can explain IceCube neutrino data. Equation (20) has an analytic solution:

$$N_{\varepsilon_{p}}=\frac{t_{\rm cool}}{\varepsilon_{p}}\int_{\varepsilon_{p}}^{\infty}d\varepsilon_{p}^{\prime}\dot{N}_{\varepsilon_{p}^{\prime},\rm inj}\exp\left(-\mathcal{G}(\varepsilon_{p},~{}\varepsilon_{p}^{\prime})\right),$$

(22)

$$\mathcal{G}(\varepsilon_{1},~{}\varepsilon_{2})=\int_{\varepsilon_{1}}^{\varepsilon_{2}}\frac{t_{\rm cool}}{t_{\rm esc}}\frac{d\varepsilon_{p}^{\prime}}{\varepsilon_{p}^{\prime}}.$$

(23)

The cooling and loss processes are the same with the stochastic acceleration model. See Ref. [29] for details.

Data Availability.    The data generated in this study have been deposited in
https://doi.org/10.6084/m9.figshare.15170790. The observational data are obtained from the tables and fitting results provided in the references shown in Figure captions. The extracted data files are available upon reasonable requests.

Code Availability.    The codes used for this study are available from S.S.K. and K.M. upon reasonable request.