Neutrino and CR Emission from RIAFs in LLAGN

# Neutrino and Cosmic-Ray Emission and Cumulative Background from Radiatively Inefficient Accretion Flows in Low-Luminosity Active Galactic Nuclei

## Abstract

We study high-energy neutrino and cosmic-ray (CR) emission from the cores of low-luminosity active galactic nuclei (LLAGN). In LLAGN, the thermalization of particles is expected to be incomplete in radiatively inefficient accretion flows (RIAFs), allowing the existence of non-thermal particles. In this work, assuming stochastic particle acceleration due to turbulence in RIAFs, we solve the Fokker-Planck equation and calculate spectra of escaping neutrinos and CRs. The RIAF in LLAGN can emit CR protons with  PeV energies and TeV–PeV neutrinos generated via and/or reactions. We find that, if % of the accretion luminosity is carried away by non-thermal ions, the diffuse neutrino intensity from the cores of LLAGN may be as high as , which can be compatible with the observed IceCube data. This result does not contradict either of the diffuse gamma-ray background observed by Fermi or observed diffuse cosmic-ray flux. Our model suggests that, although very-high-energy gamma rays may not escape, radio-quiet AGN with RIAFs can emit GeV gamma-rays, which could be used for testing the model. We also calculate the neutron luminosity from RIAFs of LLAGN, and discuss a strong constraint on the model of jet mass loading mediated by neutrons from the diffuse neutrino observation.

acceleration of particles — accretion, accretion disks — galaxies: nuclei — neutrinos — diffuse radiation

## 1. Introduction

The IceCube collaboration reported a discovery of extraterrestrial neutrinos with deposited energies ranging from 30 TeV to a few PeV, and the significance now exceeds (Aartsen et al., 2013a, b, 2014). The signals are likely to be astrophysical, and the consistency with isotropic distribution suggests extragalactic components, which is also supported by diffuse gamma-ray data (Murase et al., 2013; Ahlers & Murase, 2014). Significant clustering has not been observed, and the origin of IceCube neutrinos is a new big mystery even though some early models can match the observed data within large model uncertainties (Anchordoqui et al., 2004; Stecker, 2005; Loeb & Waxman, 2006; Murase et al., 2006; Gupta & Zhang, 2007; Murase et al., 2008; Kotera et al., 2009). One of the popular possibilities is neutrino emission from cosmic-ray (CR) reservoirs. Star-forming galaxies, including starbusrt galaxies, may explain the IceCube data, and various possibilities have been speculated to have  PeV CR protons (Murase et al., 2013; Katz et al., 2013; Tamborra et al., 2014). Galaxy groups and clusters may also account for the data without violating gamma-ray limits. While  PeV CR protons can be supplied by active galactic nuclei (AGN), galaxies and galaxy mergers 7 as well as intergalactic shocks, hard spectral indices and contributions from low-mass clusters and groups would be needed  8 (Murase et al., 2013; Kashiyama & Mészáros, 2014). However, due to multi-messenger constraints (Murase et al., 2013), the above scenarios may be challenged by the latest data implying steep indices (Aartsen et al., 2014, 2015) and models that attribute the diffuse gamma-ray background to unresolved blazars.

High-energy neutrino production in AGN has been of interest for many years (e.g., Protheroe & Kazanas, 1983; Kazanas & Ellison, 1986; Stecker et al., 1991; Szabo & Protheroe, 1994; Mannheim, 1995; Atoyan & Dermer, 2001; Alvarez-Muñiz & Mészáros, 2004; Anchordoqui et al., 2004). The most popular possibility is photohadronic () neutrino production in relativistic jets that are established as gamma-ray sites (e.g., Mannheim, 1995; Atoyan & Dermer, 2001; Mücke & Protheroe, 2001). The diffuse neutrino intensity of radio-loud AGN is typically dominated by luminous blazars, in which external radiation fields due to accretion disk, broadline and dust emission are relevant (Murase et al., 2014). However, it has been shown that the simple inner jet model has difficulty in explaining the IceCube data, and additional assumptions are required (Dermer et al., 2014; Tavecchio et al., 2014). Alternatively, high-energy neutrino emission could mainly come from the cores of AGN, and both of (Becker Tjus et al., 2014) and (Stecker, 2013; Winter, 2013; Kalashev et al., 2014) scenarios have been considered. In the latter case, it has been assumed that target photon fields come from the standard, Shakura-Sunyaev disk (Shakura & Sunyaev, 1973), but the disk temperature to account for  TeV neutrinos has to be higher than typical values (Dermer et al., 2014).

One big issue in the AGN core models is how to have non-thermal protons in such inner regions. Electric field acceleration has been discussed (Levinson, 2000), but the formation of a gap above the black hole is not clear in the presence of copious plasma. The shock dissipation may also occur in accretion flows (Begelman et al., 1990; Alvarez-Muñiz & Mészáros, 2004; Becker et al., 2008), although efficient shock acceleration is possible when the shock is not mediated by radiation. When the accretion rate is high enough to form the standard disk (Shakura & Sunyaev, 1973) or the slim disk (Abramowicz et al., 1988), protons and electrons should be thermalized via the Coulomb scattering within the infall time. This indicates that turbulent acceleration is unlikely in the bulk of the accretion flow, although possible non-thermal proton acceleration in the corona of standard disks has also been discussed (e.g., Dermer et al., 1996; Romero et al., 2010).

In this work, we consider the possibility of high-energy neutrino emission from low-luminosity AGN (LLAGN), which has not been discussed in light of the IceCube data. It has been considered that LLAGN do not have standard or slim disks, since their spectra show no blue bump (Ho, 2008). Instead, LLAGN are believed to have the radiation inefficient accretion flows (RIAFs, Narayan & Yi, 1994), in which plasma can be collisionless with low accretion rates, allowing the existence of non-thermal particles (Mahadevan & Quataert, 1997; Toma & Takahara, 2012). In the unified picture of AGN, BL Lac objects correspond to LLAGN viewed from an on-axis observer, whereas quasar-hosted blazars correspond to highly accreted AGN with optically thick disks.

Turbulent magnetic fields play an important role for the angular momentum transport in accretion disks, and the magnetic rotational instability (MRI) has been believed to be responsible for the effective viscosity (e.g., Balbus & Hawley, 1991; Sano et al., 2004). Recent particle-in-cell simulations have shown that the magnetic reconnection occurring in the MRI turbulence generates non-thermal particles (Riquelme et al., 2012; Hoshino, 2013, 2015), although these simulations follow the plasma scale structure that is much smaller than the realistic scale of the accretion flow. It would be also natural to expect stochastic acceleration in the presence of strong turbulence in RIAFs (e.g., Lynn et al., 2014).

Protons accelerated in RIAFs lead to gamma-ray emission via and processes (Mahadevan et al., 1997; Niedźwiecki et al., 2013). In this work, we focus on neutrino and CR emission, motivated by the latest IceCube discovery. The acceleration efficiency is uncertain at present. Kimura et al. (2014) shows that high-energy particles do not affect the dynamical structure, except for the cases that they extract most of the released energy from RIAFs. This implies that the energy loss rate by high-energy protons is limited to several percents of , where is the mass accretion rate, but we show that it is still possible for LLAGN to make significant contributions to the cumulative neutrino background without violating existing observational limits.

In this paper, we consider stochastic acceleration in RIAFs of LLAGN and suggest that they are potential sources of high-energy neutrinos. In Section 2, we set up physical states of RIAFs. In Section 3, we formulate and calculate energy spectrum of the non-thermal protons inside a typical RIAF. The spectra of escaping protons and neutrinos from the RIAF are also presented in Section 3. Then, the diffuse neutrino intensity is estimated and compared to the IceCube data in Section 4. We discuss several related issues such as the detectability in gamma rays in Section 5, and summarize our results in Section 6.

## 2. Physical Setup

We model emission from RIAFs with the one-zone approximation, where it is assumed that particles are accelerated only within some radius . When one considers the structure of accretion disks, the multi-dimensionality is important in general. Nevertheless, we consider that our approach is enough as the first step to consider high-energy neutrino emission from RIAFs.

### 2.1. Physical quantities of RIAFs

RIAFs are the hot and rapid infall accretion flows. We set the radial velocity , thermal proton density , thermal pressure , and strength of magnetic fields of our RIAF model as follows,

 vr = αvK, (1) np = ˙M2πR2vrmp, (2) Pth = npGMBH3Rmp, (3) B = √8πPthβ, (4)

where is the alpha parameter (Shakura & Sunyaev, 1973), is the Keplerian velocity, is the mass accretion rate, is the mass of the super massive black hole (SMBH), and is the plasma beta parameter. We assume the scale height of the flow . We normalize the radius and mass accretion rate as and , respectively, where we use the Schwarzschild radius and the Eddington accretion rate . This makes

 R = 2.95×1013 r1MBH,7 cm , (5) vr = 6.7×108 r−1/21α−1 cm s−1 , (6) np = 1.1×109 r−3/21α−1−1M−1BH,7˙m−2 cm−3 , (7) B = 4.9×102 r−5/41α−1/2−1β−1/23M−1/2BH,7˙m1/2−2 Gauss , (8)

where , except and . For reference, the accretion luminosity and Thomson optical depth are computed as

 ˙Mc2 = 1.3×1043MBH,7˙m−2 erg s−1, (9) τT = npσTR=2.2×10−2r−1/21α−1−1˙m−2, (10)

where is the Thomson cross section. If we consider small , and are quite different from above expression because the flow becomes supersonic and particles go into the SMBH quickly (Narayan et al., 1997; Kimura et al., 2014). In this paper, we fix the parameters and for demonstration.

### 2.2. Thermal electrons and target photon fields

The photomeson production is an important process of the neutrino and/or gamma-ray production. The sub-PeV and PeV neutrinos are generated by the interaction between protons and photons of the orders of eV and eV, respectively. Thus, we estimate the target photon spectrum in RIAFs. The estimate of luminosity of LLAGN, , will be also used to calculate the diffuse neutrino flux in Section 4.

It has been suggested that in LLAGN, emission comes from a jet, an outer thin disk, and a RIAF (Nemmen et al., 2006, 2014). In this paper, we consider only radiation from the RIAF because radiation from the jet and thin disk would be sub-dominant. We use the one-zone approximation and calculate the photon spectrum within acceleration radius . Thermal electrons in the RIAF emit radiation through the synchrotron, bremsstrahlung, and inverse Compton scattering. We use fitting formulae of the emissivity of bremsstrahlung and synchrotron (Narayan & Yi, 1995). Assuming the local thermodynamic equilibrium with Eddington approximation (Rybicki & Lightman, 1979), we can get the photon fields from synchrotron and bremsstrahlung. This treatment consistently includes the synchrotron self absorption (Manmoto et al., 1997). Using this photon fields as the seed photons, spectra of inverse Compton scattering are calculated. See the Appendix for details of the calculation of target photon fields.

To obtain the target photon spectra, we need to know the electron temperature. Since the relaxation time between electrons and protons in RIAFs is longer than the infall time (see Section 3.1), electrons would have different temperature from that of the protons (Takahara & Kusunose, 1985). The electron temperature is usually determined by an energy balance, , where and is the heating rate and cooling rate of electrons, respectively. The mechanism of electron heating in RIAFs is determined by details of dissipation in collisionless plasma (Quataert & Gruzinov, 1999; Sharma et al., 2007; Howes, 2010), but the accurate prescription for the turbulent heating is not well understood. Here, we assume a simple heating prescription where the heating rate of electrons are proportional to the accretion luminosity, i.e.,

 Q+=δe˙Mc2, (11)

where is the heating parameter that represents the fraction of energy that directly heats up electrons. We consider the range of following to Nemmen et al. (2014), and use as a fiducial value 9. The cooling rate is estimated from the total luminosity of target photons , where we use the assumption of optically thin limit. We calculate the equilibrium temperature, , from the energy balance using the bisection method (Press et al., 1992). Since the source of thermal energy is the released gravitational energy, the electron temperature cannot exceed the virial temperature of electron-proton gas, . Thus, the electron temperature is written as . Then, we assume the distribution function of thermal electrons is the relativistic Maxwellian,

 Ne(γe)=neγ2eβeexp(−γe/θe)θeK2(1/θe), (12)

where is the electron number density, and are the velocity and the Lorentz factor of the thermal electrons, respectively, and is the second modified Bessel function. We ignore effects of the pair production on the thermal component for simplicity, which gives . We do not consider non-thermal electrons because electrons have much shorter relaxation time than protons. When , they become thermalized within infall time through the synchrotron self absorption process (Mahadevan & Quataert, 1997). For , the electrons seem to be non-thermal. However, we ignore the effect of non-thermal electrons because such low-luminosity objects are less important for the diffuse neutrino flux (see Section 4).

We tabulate the values of for some models in Table 1 (where the other resultant quantities will be introduced later). We fix the parameters , , , and . For the reference model A1, and it does not depend on very much. The high leads to high due to the high heating rate. For the flows with lower density (lower ) or weaker magnetic fields (higher ), the electron temperature is higher because the cooling rate for a given is lower in such flows.

Figure 1 shows target photon spectra in RIAFs for models A1, A2, and A3 (for models A4 and A5, see Section 3 and 4, respectively). The values of and are tabulated in Table 1, where is the X-ray luminosity in the 2-10 keV band. For model A1, the synchrotron component has a peak at eV. The thermal electrons scatter seed synchrotron photons efficiently, and make a few peaks from the infrared to soft X-ray range. Multiple-scattered photons may make an almost flat spectrum for the hard X-ray range. The spectrum has a cutoff corresponding to the electron temperature. The inverse Compton scattering dominates over the bremsstrahlung in all the frequency range for A1.

The efficiency of inverse Compton scattering depends on the parameter, , where is the optical depth for Thomson scattering. Although is higher as is lower, it makes the parameter lower. Thus, the spectrum by inverse Compton scattering is softer for lower . For A2, parameter is less than unity, but it still dominates over the bremsstrahlung in all the range. The parameter is independent of in our formulation. The high makes luminosity higher due to high values of and . It also makes the synchrotron peak frequency low because of weak . The profile of the spectrum for A3 is similar to that for A1 but the luminosity for A3 is about ten times higher than that for A1. When the electron temperature is higher with fixed , the parameter become higher. Thus, the spectrum is harder for higher and higher .

## 3. Spectra of Non-thermal Particles in a Typical RIAF

### 3.1. Plasma in accretion flows

If the infall time is shorter than the relaxation time due to the Coulomb scattering , it allows the existence of non-thermal particles. The infall time for RIAFs is estimated to be

 tfall≃Rvr∼4.4×104r3/21α−1−1MBH,7 s , (13)

whereas the proton-proton relaxation time is estimated as

 trel = Missing or unrecognized delimiter for \left (14) ∼ 2.1×107α−1MBH,7˙m−1−2 s

where is the Coulomb logarithm (e.g., Spitzer, 1962). Thus, RIAFs satisfy , which allows to be non-thermal (cf. Takahara & Kusunose, 1985; Mahadevan & Quataert, 1997). For RIAFs, has the same order as the dissipation time via the viscosity (e.g., Pringle, 1981). Thus, the proton distribution function in RIAFs may not be Maxwellian within the dissipation time.

The protons inside RIAFs are scattered by turbulent magnetic fields. This process changes a momentum of each proton whose distribution function may be different from Maxwellian. In this paper, we consider relativistic protons accelerated through stochastic acceleration in RIAFs. It is expected that the stochastic acceleration leads to a hard spectrum of protons with , where (e.g., Becker et al., 2006; Stawarz & Petrosian, 2008). Thus, most of the accelerated protons accumulate on the high-energy end of proton distribution (see Equation (26)). This implies that it is impossible to accelerate all the protons in RIAFs because the protons are accelerated using the gravitational energy released by accretion, which is typically 0.1 per a proton. We assume only a small fraction of protons are injected to relativistic energy through some plasma processes, such as the magnetic reconnection (Hoshino, 2013, 2015), and those relativistic protons are governed by the Fokker-Plank equation (e.g., Stawarz & Petrosian, 2008),

 ∂∂tF(p) = 1p2∂∂p[p2(Dp∂∂pF(p)+ptcoolF(p))] (15) − F(p)(t−1diff+t−1fall)+˙Finj,

where is the distribution function of the non-thermal protons () , is the momentum of the protons, is the diffusion coefficient for the momentum space, is the injection term, is the cooling time, is the diffusion time, and is the infall time.

When we consider the relativistic particles, we should compare the Coulomb loss time for relativistic particles to . The Coulomb loss time is estimated to be (e.g., Dermer et al., 1996)

 tCoul≃1225(γp−1)(3.8θ3/2e+1.0)τTlnΛRc ∼7×107r3/21α−1MBH,7˙m−1−2θ3/2eγp,1 s (16)

where is the Lorentz factor of the proton. Since is satisfied for RIAFs, we can neglect the Coulomb loss in RIAFs.

It is considered that quasars have standard disks, in which the physical quantities are much different from those in RIAFs. For the Shakura-Sunyaev disks in the gas pressure dominant regime (gas-SSD, Shakura & Sunyaev, 1973), we have longer ( sec), and shorter () than those of RIAFs. The dissipation time is the same as that of RIAFs (see Equation [13]). Thus, is satisfied in gas-SSDs. The distribution function is expected to be Maxwellian due to the efficient Coulomb scattering. Even for the relativistic particles, the Coulomb loss time is much shorter than the dissipation time for because they have large optical depth (for the value of , see Equation (2.16) of Shakura & Sunyaev, 1973). Therefore, it seems difficult to accelerate the particles in gas-SSDs. For other solutions, such as standard disks in the radiation pressure dominant regime (Shakura & Sunyaev, 1973) and magnetically arrested disks (Bisnovatyi-Kogan & Ruzmaikin, 1974), the Thomson optical depth may not be as large as gas-SSDs, and it might be possible to satisfy .

### 3.2. Timescales

Equation (15) involves four important timescales, the acceleration time , the diffusion time , the infall time , and the cooling time .

In this paper, we assume a power spectrum , and fix the index of the power spectrum for simplicity. This value is motivated by the Alfvénic turbulence (Goldreich & Sridhar, 1995), although other modes may also play an important role on particle acceleration. According to the quasi-linear theorem, the diffusion coefficient is (e.g., Dermer et al., 1996)

 Dp≃(mpc)2(ckmin)(vAc)2ζ(rLkmin)q−2γqp, (17)

where is the minimum wave number of the turbulence, is the Alfven speed, , and is the ratio of the strength of turbulent fields to that of the non-turbulent fields. Then, the acceleration time is

 taccel ≃ p2Dp≃1ζ(vAc)−2Rc(rLR)2−qγ2−qp (18) ∼ 1.1×103r25/121α1/6−1β7/63M5/6BH,7˙m−1/6−2ζ−1−1 × γ1/3p,1 s.

We consider the diffusive escape and infall as the sink term of Equation (15). The particles fall to the SMBH in the infall time, given by Equation (13). For isotropically turbulent magnetic fields, the diffusion time is (e.g., Stawarz & Petrosian, 2008)

 tdiff ≃ 9Rcζ(rLR)q−2γq−2p (19) ∼ 6.7×106r11/121α−1/6−1β−1/63M7/6BH,7˙m1/6−2ζ1−1 × γ−1/3p,1 s.

For the cooling time, we consider inelastic and reactions, and the proton synchrotron emission process. The total cooling rate is given as

 t−1cool=t−1pp+t−1pγ+t−1sync, (20)

where , , and are cooling time scales for each process. We neglect the inverse Compton scattering by protons and the Bethe-Heitler process because they are typically sub-dominant. The synchrotron cooling rate is

 t−1sync=43(memp)3cσTUBmec2γp, (21)

where is the energy density of the magnetic fields. The cooling rate is

 t−1pp=npσppcKpp, (22)

where is the proton inelasticity of the process. The total cross section of this process is represented as a function of the proton energy ,

 σpp≃(34.3+1.88L+0.25L2)[1−(Epp,thrEp)4]2 mb (23)

for , where and 1.22 GeV (Kelner et al., 2006). The cooling rate is

 t−1pγ=c2γ2p∫∞¯εthrd¯εσpγ(¯ε)Kpγ(¯ε)¯ε ×∫∞¯ε/(2γp)dEγNγ(Eγ)E2γ, (24)

where and are the photon energy in the proton rest frame and the black hole frame, respectively, is the photon occupation number, and 145 MeV. We use the rectangular approximation for this process (Stecker, 1968). Assuming , we write as

 t−1pγ=c2γ2p¯ϵpkΔ¯ϵpkσpkKpk ×∫∞¯ϵpk/(2γp)dEγNγ(Eγ)E2γ, (25)

where GeV, , , GeV.

Figure 2 shows the timescales for models A1, A2, A3, and A4, whose parameters are tabulated in Table 1. For low values of , is the shortest for all the models. At some energy , is satisfied. Above the energy, the acceleration is limited by the diffusive escape. Equating and , we can estimate to be

 γp,eq ∼ (3ζvAc)3(RrL) (26) ∼ 1.4×105˙m1/2−2M1/2BH,7α1/2−1ζ3−1β−23r−7/41.

We tabulate the value of in Table 1. This characteristic energy strongly depends on . The higher or lower makes the magnetic fields stronger, so that the is higher. The larger weakens , which leads to lower . The estimation in Equation (26) is correct as long as we choose . As becomes higher, becomes longer but does not change. Thus, the infall limits the acceleration for . We note that does not correspond to the peak energy of the spectrum (see the next subsection). When diffusive escape limits acceleration, the distribution function declines gradually above , whose asymptote is for (see Equation (56) of Becker et al., 2006). This allows the protons to have about 10 times higher energy than the estimate in Equation (26). Thus, LLAGN can have the protons up to eV when .

For all the models, at low energies, inelastic collisions dominate over the synchrotron and photomeson production processes. At high energies around GeV, the photomeson production becomes relevant although the synchrotron cooling is comparable to it. For large or large cases, the reaction is more efficient than the synchrotron owing to the high target photon density.

### 3.3. Spectra of Non-thermal Particles

When we solve Equation (15), we treat the injection term as a delta-function , where is the injection proton momentum and is the normalization factor of injection. We fix because little affects the profile of distribution function as long as we choose . We assume that the total luminosity expended to inject and accelerate relativistic protons is proportional to the accretion luminosity, . As seen in the previous subsection, the proton acceleration is limited by escape. We determine the normalization of relativistic protons such that the luminosity of injection and acceleration balances with the escape luminosity, i.e.,

 ηcr˙Mc2=∫dV∫dp4πp2F(p)Ep(t−1fall+t−1diff), (27)

where is a parameter of injection efficiency. This parameter determines the normalization of the non-thermal protons, not affecting the shapes of the spectra. Kimura et al. (2014) shows that the non-thermal particles do not substantially affect the dynamical structure if . We use as a fiducial value.

We solve Equation (15) until steady solutions are realized by using the Chang-Cooper method (Chang & Cooper, 1970). We set the computational region from to GeV and divide the grids so that they are uniform in the logarithmic space. The number of the grid points is . We calculate some models with and find that the results are unchanged by the number of grids.

From the calculation results, we estimate the cosmic-ray pressure defined as

 Pcr=4π∫dpp2fcp3. (28)

We tabulate the ratio of cosmic-ray pressure to thermal pressure in Table 1. We find that , and is almost independent of the other parameters. For all the models, is satisfied, and thus, it is justified that non-thermal protons do not affect the dynamical structure of RIAFs (Kimura et al., 2014).

Since the peak energy is determined by CR escape for all the models, the profiles of the distribution functions are quite similar to each other. They show a power law for low (see Equation 56 of Becker et al., 2006). For , they deviate from the power-law and their asymptote is for (Becker et al., 2006). After obtaining , we estimate the differential luminosity spectra of the escaping protons to be

 EpLEp=∫dV4πp3F(p)Eptdiff=4π2cR3p4F(p)tdiff. (29)

We plot in Figure 3. We tabulate parameter sets in Table 1, fixing the parameters , , , , , and . For A1, the total luminosity of protons is , and the differential luminosity has a peak at PeV. Since the total luminosity of escaping protons is proportional to the released energy , the peak luminosity of escaping protons is almost proportional to and (see A2 and A3 in Figure 3), while it is almost independent of other parameters. All the models have a power law, , for . For , the spectra deviate from the power law and their asymptote is for (see Equation (70) of Becker et al., 2006). This makes a peak at the energy . The parameter dependence of is consistent with the estimation by Equation (26).

The neutrino spectrum is estimated to be

 EνLEν=(12tpp+38tpγ)4π2R3p3EpF(p), (30)

where is the neutrino energy. As long as the reaction is the dominant process of neutrino production, this treatment becomes invalid for spectra that are harder than (e.g., Kelner et al., 2006). Since we expect hard proton spectra with , our analytical method to calculate neutrino spectra will not be accurate at low energies. Thus, we show neutrino spectra only at TeV energies.

Figure 4 depicts spectra of neutrinos, The neutrinos are mainly made via the collisions for A1, A2, A3, because for . The cooling rate is almost independent of the proton energy. Thus, neutrino spectra are similar to those of protons unless proton spectra are too hard. The neutrino luminosity at the peak is estimated to be , where is the peak neutrino energy. We can see this feature in Figure 4 by comparing the dashed lines. On the other hand, both the and processes are important for A4. The photomeson production is dominant for GeV. This makes another peak in the spectra because the neutrino spectrum by reactions reflects the target photon spectrum. For example, in A4, the target photon field has a bump made by the inverse Compton scattering at eV, which leads to a peak in the neutrino spectrum at GeV. The total neutrino luminosity, , is tabulated in Table 1.

The gamma rays are also emitted through pion decay. The total luminosity of the gamma rays would be the same order as the neutrino luminosity, although the spectrum is different because of the pair production process (see Section 5.2).

Since proton acceleration is limited by escape in our models, the total injection luminosity is almost the same as the proton escape luminosity. This allows us to write the efficiency of pion production as

 fπ≈t−1pp+t−1pγt−1fall+t−1diff. (31)

If and at , we can write the parameter dependence of pion production efficiency at as . Thus, LLAGN with high can emit neutrinos more efficiently than those with low . On the other hand, we cannot simply write down the parameter dependence of neutrino luminosity for dominant cases because it depends on the target photon spectrum. The efficiency of pion production is tabulated in Table 1. For all models, the pion production efficiency is . Thus, most of the high-energy protons escape from RIAFs without losing their energies.

If we consider models that have high , , and , compared to A1, CR acceleration is limited by the photomeson production because high increases , and high and/or decrease . In this case, the scaling of the peak energy is different from the one obtained with Equation (26), and proton spectra could change (see, e.g., Stawarz & Petrosian, 2008). However, this parameter range looks extreme in our model. High and lead to high photon luminosities, which are inconsistent with the concept of RIAFs. In addition, should be less than unity for the validity of the quasi-linear theory. Thus, we can focus on the models where the diffusive escape limits the acceleration.

## 4. Diffuse Intensities of Neutrinos and Cosmic-Ray Protons

The diffuse neutrino intensity from extragalactic sources is given by (e.g., Alvarez-Muñiz & Mészáros, 2004; Murase et al., 2014)

 Φν=c4πH0∫zmax0dz√ΩM(1+z)3+ΩΛ ×∫LmaxLmindLbolϕ(Lbol, z)LE′ν(Lbol)E′ν, (32)

where is the luminosity function, is the neutrino energy at the rest flame of LLAGN. We assume that LLAGN exist from to and from to .

The luminosity function of from nearby LLAGN is plotted in Figure 8 of Ho (2008). Here we assume a broken power law shape,

 ϕ0(LHα)=n∗/L∗(LHα/L∗)s1+(LHα/L∗)s2. (33)

From Figure 8 of Ho (2008), we find , , , between . For LLAGN, is related to the bolometric luminosity as (Ho, 2008). Since the redshift evolution is poorly known, we assume no evolution of the luminosity function . This is because LLAGN are similar to the BL Lac objects in the sense that they have a faint disk component. The luminosity function of BL Lac objects is nearly consistent with no evolution (Ajello et al., 2014).

In our RIAF model, we can calculate for given , , , and as described in Section 2.2. We assume that all the LLAGN have the same values of , , . Then, we can integrate Equation (32), using the relationship of to , where is a monotonically increasing function of as shown in the upper panel of Figure 5. The parameters are tabulated in Table 2. The break of each line corresponds to the change of from to . In reality, other components such as the radiation from jets may contribute to from LLAGN, but we ignore the other components for simplicity. The bottom panel shows the X-ray luminosity in the 2-10 keV band, , as a function of . We can see that the deviation between our model and observation is not large in , although the faint part is quite different. We set but the detailed value of little affects the results. A RIAF is expected to change to a standard thin disk above the maximum accretion rate . We use following to Xie & Yuan (2012). We tabulated the physical quantities for LLAGN with as model A5 in Table 1. Then, the maximum luminosity is written as

 Lmax=δe˙MEdd˙mmaxc2≃7.6×1041MBH,7δe,−2 erg s−1. (34)

The maximum luminosity for each model is tabulated in Table 2. We calculate the diffuse spectra with some values of and confirm that does not affect the results if we use a sufficiently high value . Recent studies show that the number density of is high at and monotonically decreases with (e.g., Li et al., 2011). On the other hand, it seems that the average of of nearby LLAGN whose multi-band spectra are observed is (Eracleous et al., 2010). This suggests that the average mass of SMBHs in LLAGN is not so clear. We calculate two cases for and . Fixing , , , and , we search suitable , , and to see whether the calculated intensity amounts to the observed one.

### 4.1. Diffuse intensity of neutrinos

In this subsection, we show that our models can fit spectra of neutrinos observed by IceCube (Aartsen et al., 2014). Recently, IceCube reported neutrino spectra around 10 TeV (Aartsen et al., 2015). The flux around 10 TeV is higher than that at PeV energies. Although it may be premature to discuss the origin of this low-energy excess, we show that it is possible for LLAGN to explain the excess.

The results are plotted in upper panel of Figure 6, whose parameter sets are tabulated in Table 2. LLAGN can emit enough neutrinos to amount to the observed flux with reasonable parameter sets. In view of the PeV neutrino observation, protons must be accelerated up to around several tens of PeV energies, and their spectra cannot be extended to higher energies. This is feasible unless is somehow very low. The spectral shape is affected by and because these parameters determine whether the photomeson production is important or not. The injection efficiency just determines the normalization of the diffuse neutrino flux, and we find that to account for the PeV neutrinos. The spectra for B1 and B3 can fit the data at 0.1–1 PeV energies, although we have difficulty in explaining the 10–100 TeV neutrino flux at the same time. They are almost flat for 0.1 PeV 1 PeV and have a cutoff at a few PeV energies. We can also fit the data of 10 TeV neutrinos with lower values of (B2 and B4). They have a peak at TeV and gradually decrease for TeV. The photomeson production is ineffective in these models because of the lack of target photons. Although we need a higher for 10 TeV neutrinos, the required injection efficiency is lower than that for other AGN models (cf., Alvarez-Muñiz & Mészáros, 2004; Murase et al., 2014). Even if each LLAGN is much fainter than quasars, the number density of LLAGN is so high that they can significantly contribute to the diffuse neutrino flux in principle.

We show four models in the upper panel of Figure 6 for demonstration, but it is possible to fit the data with other sets of parameters. For the models with higher , higher is needed to achieve the observed flux. This is because higher makes higher, and thereby the corresponding number density of LLAGN decreases. The higher