Neutrino and Cosmic-Ray Emission and Cumulative Background from Radiatively Inefficient Accretion Flows in Low-Luminosity Active Galactic Nuclei
aa r X i v : . [ a s t r o - ph . H E ] J un Draft version June 12, 2018
Preprint typeset using L A TEX style emulateapj v. 5/2/11
NEUTRINO AND COSMIC-RAY EMISSION AND CUMULATIVE BACKGROUND FROMRADIATIVELY INEFFICIENT ACCRETION FLOWS IN LOW-LUMINOSITY ACTIVE GALACTIC NUCLEI
Shigeo S. Kimura , Kohta Murase , and Kenji Toma Astronomical Institute, Tohoku University, Sendai 980-8578, Japan Frontier Research Institute for Interdisciplinary Sciences, Tohoku University, Sendai 980-8578, Japan Hubble Fellow—Institute for Advanced Study, Princeton, New Jersey 08540, USA and Center for Particle and Gravitational Astrophysics; Department of Physics; Department of Astronomy & Astrophysics; ThePennsylvania State University, University Park, Pennsylvania, 16802, USA
Draft version June 12, 2018
ABSTRACTWe study high-energy neutrino and cosmic-ray (CR) emission from the cores of low-luminosity activegalactic nuclei (LLAGN). In LLAGN, the thermalization of particles is expected to be incomplete inradiatively inefficient accretion flows (RIAFs), allowing the existence of non-thermal particles. Inthis 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 emitCR protons with &
10 PeV energies and TeV–PeV neutrinos generated via pp and/or pγ reactions.We find that, if ∼
1% of the accretion luminosity is carried away by non-thermal ions, the diffuseneutrino intensity from the cores of LLAGN may be as high as E ν Φ ν ∼ × − GeV cm − s − sr − ,which can be compatible with the observed IceCube data. This result does not contradict either ofthe diffuse gamma-ray background observed by Fermi or observed diffuse cosmic-ray flux. Our modelsuggests that, although very-high-energy gamma rays may not escape, radio-quiet AGN with RIAFscan emit GeV gamma-rays, which could be used for testing the model. We also calculate the neutronluminosity from RIAFs of LLAGN, and discuss a strong constraint on the model of jet mass loadingmediated by neutrons from the diffuse neutrino observation.
Subject headings: acceleration of particles — accretion, accretion disks — galaxies: nuclei — neutrinos— diffuse radiation INTRODUCTION
The IceCube collaboration reported a discovery of ex-traterrestrial neutrinos with deposited energies rangingfrom 30 TeV to a few PeV, and the significance nowexceeds 5 σ (Aartsen et al. 2013a,b, 2014). The sig-nals are likely to be astrophysical, and the consistencywith isotropic distribution suggests extragalactic com-ponents, which is also supported by diffuse gamma-raydata (Murase et al. 2013; Ahlers & Murase 2014). Sig-nificant clustering has not been observed, and the originof IceCube neutrinos is a new big mystery even thoughsome early models can match the observed data withinlarge 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 neutrinoemission from cosmic-ray (CR) reservoirs. Star-forminggalaxies, including starbusrt galaxies, may explain theIceCube data, and various possibilities have been specu-lated to have & −
100 PeV CR protons (Murase et al.2013; Katz et al. 2013; Tamborra et al. 2014). Galaxygroups and clusters may also account for the data with-out violating gamma-ray limits. While & −
100 PeVCR protons can be supplied by active galactic nuclei(AGN), galaxies and galaxy mergers as well as inter-galactic shocks, hard spectral indices s ν ∼ [email protected]@astr.tohoku.ac.jp CR sources like AGN and galaxies strongly evolve as redshifts. needed (Murase et al. 2013; Kashiyama & M´esz´aros2014). However, due to multi-messenger constraints(Murase et al. 2013), the above scenarios may be chal-lenged by the latest data implying steep indices s ν ∼ . − . pγ ) neutrino productionin relativistic jets that are established as gamma-raysites (e.g., Mannheim 1995; Atoyan & Dermer 2001;M¨ucke & Protheroe 2001). The diffuse neutrino in-tensity of radio-loud AGN is typically dominated byluminous blazars, in which external radiation fieldsdue to accretion disk, broadline and dust emission arerelevant (Murase et al. 2014). However, it has beenshown that the simple inner jet model has difficulty inexplaining the IceCube data, and additional assump-tions are required (Dermer et al. 2014; Tavecchio et al. While neutrino and gamma-ray flux calculations by the re-cent work by Zandanel et al. (2014) is actually consistent withMurase et al. (2008, 2009) for the massive cluster shock case, thesetup noted in Murase et al. (2013) has not been tested. Con-necting individual massive cluster emission to diffuse backgroundsdepends on models and underlying assumptions, and astrophysicaluncertainty is still too large to cover all the relevant parameterspace.
Kimura, Murase, & Toma2014). Alternatively, high-energy neutrino emissioncould mainly come from the cores of AGN, and bothof pp (Becker Tjus et al. 2014) and pγ (Stecker 2013;Winter 2013; Kalashev et al. 2014) scenarios have beenconsidered. In the latter case, it has been assumed thattarget photon fields come from the standard, Shakura-Sunyaev disk (Shakura & Sunyaev 1973), but the disktemperature to account for .
300 TeV neutrinos has tobe higher than typical values (Dermer et al. 2014).One big issue in the AGN core models is how tohave non-thermal protons in such inner regions. Elec-tric field acceleration has been discussed (Levinson 2000),but the formation of a gap above the black hole isnot clear in the presence of copious plasma. Theshock dissipation may also occur in accretion flows(Begelman et al. 1990; Alvarez-Mu˜niz & M´esz´aros 2004;Becker et al. 2008), although efficient shock accelerationis possible when the shock is not mediated by radia-tion. When the accretion rate is high enough to formthe standard disk (Shakura & Sunyaev 1973) or the slimdisk (Abramowicz et al. 1988), protons and electronsshould be thermalized via the Coulomb scattering withinthe infall time. This indicates that turbulent acceler-ation is unlikely in the bulk of the accretion flow, al-though possible non-thermal proton acceleration in thecorona 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-energyneutrino 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 stan-dard or slim disks, since their spectra show no blue bump(Ho 2008). Instead, LLAGN are believed to have the ra-diation inefficient accretion flows (RIAFs, Narayan & Yi1994), in which plasma can be collisionless with low ac-cretion rates, allowing the existence of non-thermal par-ticles (Mahadevan & Quataert 1997; Toma & Takahara2012). In the unified picture of AGN, BL Lac objectscorrespond to LLAGN viewed from an on-axis observer,whereas quasar-hosted blazars correspond to highly ac-creted AGN with optically thick disks.Turbulent magnetic fields play an important role forthe angular momentum transport in accretion disks,and the magnetic rotational instability (MRI) has beenbelieved to be responsible for the effective viscosity(e.g., Balbus & Hawley 1991; Sano et al. 2004). Re-cent particle-in-cell simulations have shown that themagnetic reconnection occurring in the MRI turbulencegenerates non-thermal particles (Riquelme et al. 2012;Hoshino 2013, 2015), although these simulations followthe plasma scale structure that is much smaller than therealistic scale of the accretion flow. It would be also nat-ural to expect stochastic acceleration in the presence ofstrong turbulence in RIAFs (e.g., Lynn et al. 2014).Protons accelerated in RIAFs lead to gamma-ray emis-sion via pp and pγ processes (Mahadevan et al. 1997;Nied´zwiecki et al. 2013). In this work, we focus on neu-trino and CR emission, motivated by the latest IceCubediscovery. The acceleration efficiency is uncertain atpresent. Kimura et al. (2014) shows that high-energyparticles do not affect the dynamical structure, except forthe cases that they extract most of the released energyfrom RIAFs. This implies that the energy loss rate by high-energy protons is limited to several percents of ˙ M c ,where ˙ M is the mass accretion rate, but we show thatit is still possible for LLAGN to make significant contri-butions to the cumulative neutrino background withoutviolating existing observational limits.In this paper, we consider stochastic acceleration inRIAFs of LLAGN and suggest that they are potentialsources of high-energy neutrinos. In Section 2, we setup physical states of RIAFs. In Section 3, we formu-late and calculate energy spectrum of the non-thermalprotons inside a typical RIAF. The spectra of escapingprotons and neutrinos from the RIAF are also presentedin Section 3. Then, the diffuse neutrino intensity is esti-mated and compared to the IceCube data in Section 4.We discuss several related issues such as the detectabilityin gamma rays in Section 5, and summarize our resultsin Section 6. PHYSICAL SETUP
We model emission from RIAFs with the one-zone ap-proximation, where it is assumed that particles are accel-erated only within some radius R . When one considersthe structure of accretion disks, the multi-dimensionalityis important in general. Nevertheless, we consider thatour approach is enough as the first step to consider high-energy neutrino emission from RIAFs. Physical quantities of RIAFs
RIAFs are the hot and rapid infall accretion flows. Weset the radial velocity v r , thermal proton density n p ,thermal pressure P th , and strength of magnetic fields B of our RIAF model as follows, v r = αv K , (1) n p = ˙ M πR v r m p , (2) P th = n p GM BH R m p , (3) B = s πP th β , (4)where α is the alpha parameter (Shakura & Sunyaev1973), v K = p GM BH /R is the Keplerian velocity, ˙ M is the mass accretion rate, M BH is the mass of the su-per massive black hole (SMBH), and β is the plasmabeta parameter. We assume the scale height of the flow H ∼ R . We normalize the radius and mass accretion rateas r = R/R S and ˙ m = ˙ M / ˙ M Edd , respectively, where weuse the Schwarzschild radius R S = 2 GM BH /c and theEddington accretion rate ˙ M Edd = L Edd /c . This makes R = 2 . × r M BH , cm , (5) v r = 6 . × r − / α − cm s − , (6) n p = 1 . × r − / α − − M − , ˙ m − cm − , (7) B = 4 . × r − / α − / − β − / M − / , ˙ m / − Gauss , (8)where A n = A/ n , except M BH ,n = M BH / (10 n M ⊙ )and β = β/
3. For reference, the accretion luminosityand Thomson optical depth are computed as˙
M c = 1 . × M BH , ˙ m − erg s − , (9)eutrino and CR Emission from RIAFs in LLAGN 3 τ T = n p σ T R = 2 . × − r − / α − − ˙ m − , (10)where σ T is the Thomson cross section. If we considersmall r . v r and n p are quite different from aboveexpression because the flow becomes supersonic and par-ticles go into the SMBH quickly (Narayan et al. 1997;Kimura et al. 2014). In this paper, we fix the parame-ters α = 0 . r = 10 for demonstration. Thermal electrons and target photon fields
The photomeson production is an important processof the neutrino and/or gamma-ray production. Thesub-PeV and PeV neutrinos are generated by the in-teraction between protons and photons of the orders of E p ∼ − eV and E γ ∼ −
10 eV, respectively.Thus, we estimate the target photon spectrum in RIAFs.The estimate of luminosity of LLAGN, L E γ , will be alsoused to calculate the diffuse neutrino flux in Section 4.It has been suggested that in LLAGN, emissioncomes from a jet, an outer thin disk, and a RIAF(Nemmen et al. 2006, 2014). In this paper, we consideronly radiation from the RIAF because radiation from thejet and thin disk would be sub-dominant. We use theone-zone approximation and calculate the photon spec-trum within acceleration radius R . Thermal electronsin the RIAF emit radiation through the synchrotron,bremsstrahlung, and inverse Compton scattering. Weuse fitting formulae of the emissivity of bremsstrahlungand synchrotron (Narayan & Yi 1995). Assuming the lo-cal thermodynamic equilibrium with Eddington approxi-mation (Rybicki & Lightman 1979), we can get the pho-ton fields from synchrotron and bremsstrahlung. Thistreatment consistently includes the synchrotron self ab-sorption (Manmoto et al. 1997). Using this photon fieldsas the seed photons, spectra of inverse Compton scatter-ing are calculated. See the Appendix for details of thecalculation of target photon fields.To obtain the target photon spectra, we need to knowthe electron temperature. Since the relaxation time be-tween electrons and protons in RIAFs is longer thanthe infall time t fall (see Section 3.1), electrons wouldhave different temperature from that of the protons(Takahara & Kusunose 1985). The electron temperatureis usually determined by an energy balance, Q + = Q − ,where Q + and Q − is the heating rate and cooling rateof electrons, respectively. The mechanism of electronheating in RIAFs is determined by details of dissipa-tion in collisionless plasma (Quataert & Gruzinov 1999;Sharma et al. 2007; Howes 2010), but the accurate pre-scription for the turbulent heating is not well understood.Here, we assume a simple heating prescription where theheating rate of electrons are proportional to the accretionluminosity, i.e., Q + = δ e ˙ M c , (11)where δ e is the heating parameter that represents thefraction of energy that directly heats up electrons. Weconsider the range of 10 − ≤ δ e ≤ × − following toNemmen et al. (2014), and use δ e ∼ × − as a fidu-cial value . The cooling rate is estimated from the totalluminosity of target photons Q − = R L E γ dE γ , where This definition of δ e is different from that in the previousworks including Nemmen et al. (2014). They define the heatingparameter δ as Q + = δQ vis , where Q vis is the viscous dissipation we use the assumption of optically thin limit. We cal-culate the equilibrium temperature, θ e, eq , from the en-ergy balance Q + = Q − ( θ e, eq ) using the bisection method(Press et al. 1992). Since the source of thermal energyis the released gravitational energy, the electron temper-ature cannot exceed the virial temperature of electron-proton gas, θ e, vir = GM m p / (9 R ). Thus, the electrontemperature is written as θ e = min( θ e, eq , θ e, vir ). Then,we assume the distribution function of thermal electronsis the relativistic Maxwellian, N e ( γ e ) = n e γ e β e exp( − γ e /θ e ) θ e K (1 /θ e ) , (12)where n e is the electron number density, β e and γ e are thevelocity and the Lorentz factor of the thermal electrons,respectively, and K ( x ) is the second modified Besselfunction. We ignore effects of the pair production on thethermal component for simplicity, which gives n p = n e .We do not consider non-thermal electrons because elec-trons have much shorter relaxation time than protons.When ˙ m & − , they become thermalized within in-fall time through the synchrotron self absorption pro-cess (Mahadevan & Quataert 1997). For ˙ m . − , theelectrons seem to be non-thermal. However, we ignorethe effect of non-thermal electrons because such low-luminosity objects are less important for the diffuse neu-trino flux (see Section 4).We tabulate the values of θ e for some models in Table1 (where the other resultant quantities will be introducedlater). We fix the parameters α = 0 . β = 3, r = 10, and δ e = 3 × − . For the reference model A1, θ e ∼ . M BH very much. The high δ e leadsto high θ e due to the high heating rate. For the flowswith lower density (lower ˙ m ) or weaker magnetic fields(higher β ), the electron temperature is higher becausethe cooling rate for a given θ e is lower in such flows.Figure 1 shows target photon spectra in RIAFs formodels A1, A2, and A3 (for models A4 and A5, see Sec-tion 3 and 4, respectively). The values of τ T and L X are tabulated in Table 1, where L X is the X-ray lumi-nosity in the 2-10 keV band. For model A1, the syn-chrotron component has a peak at E γ ∼ .
03 eV. Thethermal electrons scatter seed synchrotron photons effi-ciently, and make a few peaks from the infrared to softX-ray range. Multiple-scattered photons may make analmost flat spectrum for the hard X-ray range. The spec-trum has a cutoff corresponding to the electron temper-ature. The inverse Compton scattering dominates overthe bremsstrahlung in all the frequency range for A1.The efficiency of inverse Compton scattering dependson the y parameter, y ∼ τ T θ e ∝ ˙ mr − / α − θ e , where τ T = n e σ T R is the optical depth for Thomson scatter-ing. Although θ e is higher as ˙ m is lower, it makes the y parameter lower. Thus, the spectrum by inverse Comp-ton scattering is softer for lower ˙ m . For A2, y param-eter is less than unity, but it still dominates over thebremsstrahlung in all the range. The y parameter is in-dependent of M BH in our formulation. The high M BH makes luminosity higher due to high values of R and ˙ M . rate. Since the viscous dissipation rate is typically the order of Q vis ∼ . Mc , the value of δ e = 0 .
05 in our model corresponds to δ = 0 . Kimura, Murase, & Toma
Fig. 1.—
Target photon spectra emitted by thermal electronsin RIAFs. The red-solid, the green-dashed, and blue-dotted linesshow models A1 (reference), A2 (low ˙ m ), A3 (high M BH ), respec-tively. The target photon spectrum for model A4 is the same withthat for A1. It also makes the synchrotron peak frequency low be-cause of weak B . The profile of the spectrum for A3 issimilar to that for A1 but the luminosity for A3 is aboutten times higher than that for A1. When the electrontemperature is higher with fixed ˙ m , the y parameter be-come higher. Thus, the spectrum is harder for higher β and higher δ e . SPECTRA OF NON-THERMAL PARTICLES IN ATYPICAL RIAF
Plasma in accretion flows
If the infall time t fall is shorter than the relaxationtime due to the Coulomb scattering t rel , it allows theexistence of non-thermal particles. The infall time forRIAFs is estimated to be t fall ≃ Rv r ∼ . × r / α − − M BH , s , (13)whereas the proton-proton relaxation time is estimatedas t rel = 4 √ π ln Λ 1 n p σ T c (cid:18) m p m e (cid:19) (cid:18) k B T p m p c (cid:19) / ∼ . × α − M BH , ˙ m − − s (14)where ln Λ is the Coulomb logarithm (e.g., Spitzer1962). Thus, RIAFs satisfy t rel ≫ t fall , which allows F ( p ) to be non-thermal (cf. Takahara & Kusunose 1985;Mahadevan & Quataert 1997). For RIAFs, t fall has thesame order as the dissipation time via the α viscosity t dis (e.g., Pringle 1981). Thus, the proton distributionfunction in RIAFs may not be Maxwellian within thedissipation time.The protons inside RIAFs are scattered by turbulentmagnetic fields. This process changes a momentum ofeach proton whose distribution function may be differ-ent from Maxwellian. In this paper, we consider rel-ativistic protons accelerated through stochastic accel-eration in RIAFs. It is expected that the stochasticacceleration leads to a hard spectrum of protons with s p <
1, where dN p /dE p ∝ E − s p p (e.g., Becker et al. 2006;Stawarz & Petrosian 2008). Thus, most of the acceler- ated protons accumulate on the high-energy end of pro-ton distribution (see Equation (26)). This implies that itis impossible to accelerate all the protons in RIAFs be-cause the protons are accelerated using the gravitationalenergy released by accretion, which is typically 0.1 m p c per a proton. We assume only a small fraction of protonsare injected to relativistic energy through some plasmaprocesses, such as the magnetic reconnection (Hoshino2013, 2015), and those relativistic protons are governedby the Fokker-Plank equation (e.g., Stawarz & Petrosian2008), ∂∂t F ( p ) = 1 p ∂∂p (cid:20) p (cid:18) D p ∂∂p F ( p ) + pt cool F ( p ) (cid:19)(cid:21) − F ( p ) (cid:0) t − + t − (cid:1) + ˙ F inj , (15)where F ( p ) is the distribution function of the non-thermal protons ( dN p /dE p = 4 πE p F ( p ) c ) , p is the mo-mentum of the protons, D p is the diffusion coefficient forthe momentum space, ˙ F inj is the injection term, t cool isthe cooling time, t diff is the diffusion time, and t fall is theinfall time.When we consider the relativistic particles, we shouldcompare the Coulomb loss time for relativistic particles t Coul to t fall . The Coulomb loss time is estimated to be(e.g., Dermer et al. 1996) t Coul ≃ γ p − . θ / e + 1 . τ T ln Λ Rc ∼ × r / α − M BH , ˙ m − − θ / e γ p, s (16)where γ p is the Lorentz factor of the proton. Since t Coul > t fall is satisfied for RIAFs, we can neglect theCoulomb loss in RIAFs.It is considered that quasars have standard disks,in which the physical quantities are much differentfrom those in RIAFs. For the Shakura-Sunyaevdisks in the gas pressure dominant regime (gas-SSD,Shakura & Sunyaev 1973), we have longer t fall ( t fall = R/v r ≃ R/ ( αv K )( R/H ) ∼ × sec), and shorter t rel ( ∼ × − sec ≪ t dis ) than those of RIAFs. Thedissipation time t dis is the same as that of RIAFs (seeEquation [13]). Thus, t rel ≪ t dis ≪ t fall is satisfied ingas-SSDs. The distribution function F ( p ) is expected tobe Maxwellian due to the efficient Coulomb scattering.Even for the relativistic particles, the Coulomb loss timeis much shorter than the dissipation time for γ p . because they have large optical depth τ T ∼ (for thevalue of τ T , see Equation (2.16) of Shakura & Sunyaev1973). Therefore, it seems difficult to accelerate the par-ticles in gas-SSDs. For other solutions, such as stan-dard disks in the radiation pressure dominant regime(Shakura & Sunyaev 1973) and magnetically arresteddisks (Bisnovatyi-Kogan & Ruzmaikin 1974), the Thom-son optical depth may not be as large as gas-SSDs, andit might be possible to satisfy t dis < t Coul . Timescales
Equation (15) involves four important timescales, theacceleration time t accel ≡ p /D p , the diffusion time t diff ,the infall time t fall , and the cooling time t cool .In this paper, we assume a power spectrum P ( k ) ∝ k − q , and fix the index of the power spectrum q = 5 / TABLE 1Model parameters and resultant physical quantities for the spectrum from a LLAGN model ˙ m M
BHa ζ θ e τ T L X b γ p, eq P cr /P th L ν, totc f π E γ, cutd A1 (reference) 1 × − . × − . × . × . × . × − × − . × − . × . × . × . × − × − . × − . × . × . × . × − × − . × − . × . × . × . × − × − . × − . × . × . × . × − a in unit of M ⊙ b X-ray luminosity in the 2-10 keV band in unit of erg s − Total neutrino luminosity in unit of erg s − in unit of GeV for simplicity. This value is motivated by the Alfv´enicturbulence (Goldreich & Sridhar 1995), although othermodes may also play an important role on particle ac-celeration. According to the quasi-linear theorem, thediffusion coefficient is (e.g., Dermer et al. 1996) D p ≃ ( m p c ) ( ck min ) (cid:16) v A c (cid:17) ζ ( r L k min ) q − γ qp , (17)where k min ∼ R − is the minimum wave number of theturbulence, v A = B/ p πm p n p is the Alfven speed, r L = m p c / ( eB ), and ζ = 8 π R P ( k ) dk/B is the ratio of thestrength of turbulent fields to that of the non-turbulentfields. Then, the acceleration time is t accel ≃ p D p ≃ ζ (cid:16) v A c (cid:17) − Rc (cid:16) r L R (cid:17) − q γ − qp ∼ . × r / α / − β / M / , ˙ m − / − ζ − − × γ / p, s . (18)We consider the diffusive escape and infall as the sinkterm of Equation (15). The particles fall to the SMBHin the infall time, given by Equation (13). For isotropi-cally turbulent magnetic fields, the diffusion time is (e.g.,Stawarz & Petrosian 2008) t diff ≃ Rc ζ (cid:16) r L R (cid:17) q − γ q − p ∼ . × r / α − / − β − / M / , ˙ m / − ζ − × γ − / p, s . (19)For the cooling time, we consider inelastic pp and pγ reactions, and the proton synchrotron emission process.The total cooling rate is given as t − = t − pp + t − pγ + t − , (20)where t pp , t pγ , and t sync are cooling time scales for eachprocess. We neglect the inverse Compton scattering byprotons and the Bethe-Heitler process because they aretypically sub-dominant. The synchrotron cooling rate is t − = 43 (cid:18) m e m p (cid:19) cσ T U B m e c γ p , (21)where U B = B / (8 π ) is the energy density of the mag-netic fields. The pp cooling rate is t − pp = n p σ pp cK pp , (22) where K pp ∼ . σ pp is representedas a function of the proton energy E p , σ pp ≃ (34 . . L + 0 . L ) " − (cid:18) E pp, thr E p (cid:19) mb(23)for E p ≥ E pp, thr , where L = log( E p / E pp, thr =1.22 GeV (Kelner et al. 2006). The pγ coolingrate is t − pγ = c γ p Z ∞ ¯ ε thr d ¯ εσ pγ (¯ ε ) K pγ (¯ ε )¯ ε × Z ∞ ¯ ε/ (2 γ p ) dE γ N γ ( E γ ) E γ , (24)where ¯ ε and E γ are the photon energy in the pro-ton rest frame and the black hole frame, respectively, N γ ( E γ ) is the photon occupation number, and ¯ ε thr =145MeV. We use the rectangular approximation for this pro-cess (Stecker 1968). Assuming σ pγ (¯ ε ) K pγ (¯ ε ) = δ (¯ ε − ¯ ǫ pk ) σ pk K pk ∆¯ ǫ pk , we write t pγ as t − pγ = c γ p ¯ ǫ pk ∆¯ ǫ pk σ pk K pk × Z ∞ ¯ ǫ pk / (2 γ p ) dE γ N γ ( E γ ) E γ , (25)where ¯ ǫ pk ∼ . σ pk ∼ × − cm , K pk ∼ . ǫ pk ∼ . E p , t accel is the shortest for all the models.At some energy E p, eq , t accel = t diff is satisfied. Abovethe energy, the acceleration is limited by the diffusiveescape. Equating t diff and t accel , we can estimate γ p, eq ≡ E p, eq / ( m p c ) to be γ p, eq ∼ (cid:18) ζv A c (cid:19) (cid:18) Rr L (cid:19) ∼ . × ˙ m / − M / , α / − ζ − β − r − / . (26)We tabulate the value of γ p, eq in Table 1. This char-acteristic energy strongly depends on ζ . The higher ˙ m or lower β makes the magnetic fields stronger, so thatthe γ p, eq is higher. The larger r weakens B , which leadsto lower γ p, eq . The estimation in Equation (26) is cor-rect as long as we choose β .
10. As β becomes higher, Kimura, Murase, & Toma Fig. 2.—
Energy dependence of the timescales. We plot the cooling time (thick-solid), the diffusion time (thick-dashed), the infall time(thick-dotted), and the acceleration time (dotted-dashed). The thin-solid, thin-dashed, and thin-dotted lines show the t pγ , t pp , and t sync ,respectively. Panels (a)–(d) show the cases for models A1 (reference), A2 (low ˙ m ), A3 (high M BH ), and A4 (high ζ ), respectively. t accel becomes longer but t fall does not change. Thus,the infall limits the acceleration for β &
10. We notethat E p, eq does not correspond to the peak energy of the E p L E p spectrum (see the next subsection). When dif-fusive escape limits acceleration, the distribution func-tion declines gradually above E p, eq , whose asymptote is F ( p ) ∝ E − / p exp( − (27 E p /E p, eq ) / ) for q = 5 / E p & eV when ζ & . pp inelastic col-lisions dominate over the synchrotron and photomesonproduction processes. At high energies around E p & − GeV, the photomeson production becomes rel-evant although the synchrotron cooling is comparable toit. For large β or large δ e cases, the pγ reaction is moreefficient than the synchrotron owing to the high targetphoton density. Spectra of Non-thermal Particles
When we solve Equation (15), we treat the injectionterm as a delta-function ˙ F inj = F δ ( p − p inj ), where p inj isthe injection proton momentum and F is the normaliza-tion factor of injection. We fix p inj = 2 m p c because p inj little affects the profile of distribution function as long as we choose p inj c ≪ E p, eq . We assume that the total lumi-nosity expended to inject and accelerate relativistic pro-tons is proportional to the accretion luminosity, ˙ M c . Asseen in the previous subsection, the proton accelerationis limited by escape. We determine the normalization ofrelativistic protons such that the luminosity of injectionand acceleration balances with the escape luminosity, i.e., η cr ˙ M c = Z dV Z dp πp F ( p ) E p (cid:0) t − + t − (cid:1) , (27)where η cr is a parameter of injection efficiency. Thisparameter determines the normalization of the non-thermal protons, not affecting the shapes of the spectra.Kimura et al. (2014) shows that the non-thermal parti-cles do not substantially affect the dynamical structureif η cr . .
1. We use η cr = 0 .
01 as a fiducial value.We solve Equation (15) until steady solutionsare realized by using the Chang-Cooper method(Chang & Cooper 1970). We set the computational re-gion from E p = 1 . GeV and divide the gridsso that they are uniform in the logarithmic space. Thenumber of the grid points is N = 500. We calculatesome models with N = 1000 and find that the resultsare unchanged by the number of grids.From the calculation results, we estimate the cosmic-eutrino and CR Emission from RIAFs in LLAGN 7ray pressure defined as P cr = 4 π Z dpp f cp . (28)We tabulate the ratio of cosmic-ray pressure to thermalpressure in Table 1. We find that P cr ∝ η cr , and P cr /P th is almost independent of the other parameters. For allthe models, P cr < P th is satisfied, and thus, it is justifiedthat non-thermal protons do not affect the dynamicalstructure of RIAFs (Kimura et al. 2014).Since the peak energy is determined by CR escapefor all the models, the profiles of the distribution func-tions are quite similar to each other. They show apower law F ( p ) ∝ E − (1+ q ) p for low E p (see Equation56 of Becker et al. 2006). For E p & E p, eq , they devi-ate from the power-law and their asymptote is F ( p ) ∝ E − / p exp( − (27 E p /E p, eq ) / ) for q = 5 / F ( p ), we estimate the differentialluminosity spectra of the escaping protons to be E p L E p = Z dV πp F ( p ) E p t diff = 4 π cR p F ( p ) t diff . (29)We plot E p L E p in Figure 3. We tabulate parameter setsin Table 1, fixing the parameters α = 0 . β = 3, r = 10, δ e = 3 × − , q = 5 /
3, and η cr = 0 .
01. For A1, the totalluminosity of protons is ∼ × erg s − , and the dif-ferential luminosity has a peak E p L E p ∼ × erg s − at E p ∼ . M c ,the peak luminosity of escaping protons is almost pro-portional to ˙ m and M BH (see A2 and A3 in Figure3), while it is almost independent of other parameters.All the models have a power law, E p L E p ∝ E − qp ,for E p < E p, eq . For E p > E p, eq , the spectra deviatefrom the power law and their asymptote is E p L E p ∝ E / p exp( − (27 E p /E p, eq ) / ) for q = 5 / E p, pk ∼ E p, eq . The parameter dependence of E p, pk is consistent with the estimation by Equation (26).The neutrino spectrum is estimated to be E ν L E ν = (cid:18) t pp + 38 t pγ (cid:19) π R p E p F ( p ) , (30)where E ν = 0 . E p is the neutrino energy. As longas the pp reaction is the dominant process of neutrinoproduction, this treatment becomes invalid for spec-tra that are harder than F ( p ) ∝ p − . − p − . (e.g.,Kelner et al. 2006). Since we expect hard proton spec-tra F ( p ) ∝ p − (1+ q ) with q = 5 /
3, our analytical methodto calculate neutrino spectra will not be accurate at lowenergies. Thus, we show neutrino spectra only at E ν > pp collisions for A1, A2, A3, because t pp < t pγ for E p . E p, pk . The pp cooling rate is almostindependent of the proton energy. Thus, neutrino spec-tra are similar to those of protons unless proton spectraare too hard. The neutrino luminosity at the peak is es-timated to be E ν L E ν | E ν, pk ∝ η cr ˙ m M BH α − β / , where E ν, pk = 0 . E p, pk is the peak neutrino energy. We cansee this feature in Figure 4 by comparing the dashed lines. On the other hand, both the pp and pγ processesare important for A4. The photomeson production isdominant for E ν & GeV. This makes another peakin the spectra because the neutrino spectrum by pγ reac-tions reflects the target photon spectrum. For example,in A4, the target photon field has a bump made by theinverse Compton scattering at E γ ∼ E ν ∼ × GeV. The total neutrino luminosity, L ν, tot , is tabulatedin Table 1.The gamma rays are also emitted through pion decay.The total luminosity of the gamma rays would be thesame order as the neutrino luminosity, although the spec-trum is different because of the pair production process(see Section 5.2).Since proton acceleration is limited by escape in ourmodels, the total injection luminosity is almost the sameas the proton escape luminosity. This allows us to writethe efficiency of pion production as f π ≈ t − pp + t − pγ t − + t − . (31)If t pp < t pγ and t diff < t fall at E p, eq , we can write theparameter dependence of pion production efficiency at E p, eq as f π ∝ ˙ mα − β / . Thus, LLAGN with high ˙ m can emit neutrinos more efficiently than those with low˙ m . On the other hand, we cannot simply write downthe parameter dependence of neutrino luminosity for pγ dominant cases because it depends on the target photonspectrum. The efficiency of pion production is tabulatedin Table 1. For all models, the pion production efficiencyis f π . .
1. Thus, most of the high-energy protons es-cape from RIAFs without losing their energies.If we consider models that have high δ e , ζ , and ˙ m ,compared to A1, CR acceleration is limited by the pho-tomeson production because high ζ increases t diff , andhigh ˙ m and/or δ e decrease t pγ . In this case, the scalingof the peak energy is different from the one obtainedwith Equation (26), and proton spectra could change(see, e.g., Stawarz & Petrosian 2008). However, this pa-rameter range looks extreme in our model. High ˙ m and δ e lead to high photon luminosities, which are inconsis-tent with the concept of RIAFs. In addition, ζ should beless than unity for the validity of the quasi-linear theory.Thus, we can focus on the models where the diffusiveescape limits the acceleration. DIFFUSE INTENSITIES OF NEUTRINOS ANDCOSMIC-RAY PROTONS
The diffuse neutrino intensity from extragalacticsources is given by (e.g., Alvarez-Mu˜niz & M´esz´aros2004; Murase et al. 2014)Φ ν = c πH Z z max dz p Ω M (1 + z ) + Ω Λ × Z L max L min dL bol φ ( L bol , z ) L E ′ ν ( L bol ) E ′ ν , (32)where φ ( L bol , z ) is the luminosity function, E ′ ν = (1 + z ) E ν is the neutrino energy at the rest flame of LLAGN.We assume that LLAGN exist from z = 0 to z = z max and from L min to L max . Kimura, Murase, & Toma Fig. 3.—
Differential luminosity spectra of escaping protons formodels A1 (red-solid), A2 (green-dashed), A3 (blue-dotted), andA4 (magenta dotted-dashed), respectively.
Fig. 4.—
Differential luminosity spectra of neutrinos. The solidlines represent the total neutrino spectra. The dashed and dot-ted lines show the neutrino spectra from pp and pγ interactions,respectively. Panel (a) shows the models A1 (red-thick lines) andA2 (blue-thin lines). Panel (b) shows the models A3 (red-thicklines) and A4 (blue-thin lines). The thin dotted-dashed line (totalspectra for A1) is also plotted for comparison in Panel (b). The luminosity function of H α from nearby LLAGNis plotted in Figure 8 of Ho (2008). Here we assume abroken power law shape, φ ( L H α ) = n ∗ /L ∗ ( L H α /L ∗ ) s + ( L H α /L ∗ ) s . (33)From Figure 8 of Ho (2008), we find L ∗ = 10 erg s − , n ∗ ∼ . × − Mpc − , s ∼ . s ∼ × erg s − < L H α < × erg s − . ForLLAGN, L H α is related to the bolometric luminosity L bol as L bol ∼ L H α (Ho 2008). Since the redshift evolutionis poorly known, we assume no evolution of the lumi-nosity function ϕ ( L bol , z ) = ϕ ( L bol ). This is becauseLLAGN are similar to the BL Lac objects in the sensethat they have a faint disk component. The luminosityfunction of BL Lac objects is nearly consistent with noevolution (Ajello et al. 2014).In our RIAF model, we can calculate L bol for given M BH , ˙ m , δ e , and β as described in Section 2.2. We as-sume that all the LLAGN have the same values of M BH , δ e , β . Then, we can integrate Equation (32), using therelationship of ˙ m to L bol , where L bol is a monotonicallyincreasing function of ˙ m as shown in the upper panel ofFigure 5. The parameters are tabulated in Table 2. Thebreak of each line corresponds to the change of θ e from θ e, vir to θ e, eq . In reality, other components such as theradiation from jets may contribute to L bol from LLAGN,but we ignore the other components for simplicity. Thebottom panel shows the X-ray luminosity in the 2-10 keVband, L X , as a function of L bol . We can see that the de-viation between our model and observation is not largein L bol & erg s − , although the faint part is quitedifferent. We set ˙ m min = 10 − but the detailed valueof ˙ m min little affects the results. A RIAF is expectedto change to a standard thin disk above the maximumaccretion rate ˙ m max . We use ˙ m max = 0 .
06 following toXie & Yuan (2012). We tabulated the physical quantitiesfor LLAGN with ˙ m max as model A5 in Table 1. Then,the maximum luminosity is written as L max = δ e ˙ M Edd ˙ m max c ≃ . × M BH , δ e, − erg s − . (34)The maximum luminosity for each model is tabulated inTable 2. We calculate the diffuse spectra with some val-ues of z max and confirm that z max does not affect theresults if we use a sufficiently high value z max &
4. Re-cent studies show that the number density of M BH ishigh at M BH ∼ M ⊙ and monotonically decreases with M BH (e.g., Li et al. 2011). On the other hand, it seemsthat the average of M BH of nearby LLAGN whose multi-band spectra are observed is ∼ M ⊙ (Eracleous et al.2010). This suggests that the average mass of SMBHsin LLAGN is not so clear. We calculate two cases for M BH = 10 M ⊙ and M BH = 10 M ⊙ . Fixing q = 5 / α = 0 . β = 3, and r = 10, we search suitable ζ , δ e ,and η cr to see whether the calculated intensity amountsto the observed one. Diffuse intensity of neutrinos
In this subsection, we show that our models can fitspectra of neutrinos observed by IceCube (Aartsen et al.2014). Recently, IceCube reported neutrino spectraaround 10 TeV (Aartsen et al. 2015). The flux aroundeutrino and CR Emission from RIAFs in LLAGN 9
TABLE 2Models for diffuse neutrino flux: parameters and resultant quantities model M BHa δ e ζ η cr f π, max L max b L mid b L max , X b L mid , Xb B1 10 . × − . × . × . × . × B2 10 . × − . × . × . × . × B3 10 . × − . × . × . × . × B4 10 . × − . × . × . × . × in unit of M ⊙ b in unit of erg s − Fig. 5.—
The bolometric and X-ray luminosities in our models.The top panel shows L bol as a function of ˙ m for each model tab-ulated in Table 2. The bottom panel shows L X as a function of L bol . The thick lines are the resultant values for our models. Thethin-solid line in bottom panel is the averaged value based on ob-servations (Ho 2008). The lines for model B4 are not shown, sincethey are completely the same as those for model B3.
10 TeV is higher than that at PeV energies. Although itmay be premature to discuss the origin of this low-energyexcess, we show that it is possible for LLAGN to explainthe excess.The results are plotted in upper panel of Figure 6,whose parameter sets are tabulated in Table 2. LLAGNcan emit enough neutrinos to amount to the observedflux with reasonable parameter sets. In view of the PeVneutrino observation, protons must be accelerated up toaround several tens of PeV energies, and their spectracannot be extended to higher energies. This is feasi-ble unless ζ is somehow very low. The spectral shapeis affected by δ e and ζ because these parameters deter-mine whether the photomeson production is important or not. The injection efficiency η cr just determines thenormalization of the diffuse neutrino flux, and we findthat η cr ∼ − − − to account for the PeV neutri-nos. The spectra for B1 and B3 can fit the data at 0.1–1PeV energies, although we have difficulty in explainingthe 10–100 TeV neutrino flux at the same time. They arealmost flat for 0.1 PeV . E ν . ζ (B2 and B4). Theyhave a peak at E ν ∼
10 TeV and gradually decrease for E ν >
10 TeV. The photomeson production is ineffectivein these models because of the lack of target photons.Although we need a higher η cr for 10 TeV neutrinos, therequired injection efficiency η cr ∼ .
01 is lower than thatfor other AGN models (cf., Alvarez-Mu˜niz & M´esz´aros2004; Murase et al. 2014). Even if each LLAGN is muchfainter than quasars, the number density of LLAGN is sohigh that they can significantly contribute to the diffuseneutrino flux in principle.We show four models in the upper panel of Figure 6for demonstration, but it is possible to fit the data withother sets of parameters. For the models with higher δ e ,higher η cr is needed to achieve the observed flux. Thisis because higher δ e makes L max higher, and thereby thecorresponding number density of LLAGN decreases. Thehigher δ e also causes more efficient pγ reaction. Thus,the maximum f π for model B1 has 0.19, which is largerthan other models (see Table 2). When we considerhigher β , the fitting requires higher ζ . For example, themodels with β = 10 and ζ = 0 . ζ ∼
1, we can use thequasi-linear theory for describing the stochastic acceler-ation if v A /c & . v A /c ∼ . r − / β − / for our model, it is acceptable touse ζ . . η cr = 2 . × − .The others have a parameter set of model B3 except for η cr = 5 . × − . The bottom panel of Figure 6 showsthe result of this two-component model. The resultingspectrum could fit all the four data points at the sametime. Although this two-component model is a possibil-ity to explain the IceCube events, this parameter choiceis ad hoc. The combination of LLAGN model with othersources would be more natural. For example, the low-0 Kimura, Murase, & Toma Fig. 6.—
The diffuse neutrino intensity (per flavor) from RIAFsin the LLAGN model. The top panel shows the diffuse neutrinointensity for each model tabulated in Table 2. The dashed line (B2)almost overlaps the dotted-dashed line (B4). The bottom panelshows the diffuse intensity from two-component model (see text fordetail). The red-solid, green-dashed, and blue-dotted lines showthe total intensity, intensity from low-energy part, and intensityfrom high-energy part, respectively. The green triangles representthe atmospheric muon neutrino background produced by CRs. Theblack squares show the observed data of neutrino signals. energy part could come from LLAGN, while the high-energy part may come from CR reservoirs such as star-burst galaxies and galaxy clusters.The diffuse neutrino flux is dominated by LLAGNwith high ˙ m in our model. The neutrino luminosityis higher as ˙ m is higher, while the number density ofLLAGN is lower for higher ˙ m . The former is more ef-ficient than the latter for the neutrino luminosity. Weshow the contribution to the total intensity from dif-ferent luminosity bins in Figure 7. We set the lumi-nosity bins as a faint part L min < L bol < L ∗ , amiddle part 80 L ∗ < L bol < L mid , and a bright part L mid < L bol < L max , where L mid = √ L ∗ L max whosevalues are tabulated in Table 2. We also tabulate the cor-responding values of L X to L max and L mid . The brightpart emits almost all the neutrinos for all models. Thefaint and middle parts little contribute to the diffuse neu-trino flux due to the low pion production efficiency. Diffuse intensity of cosmic-ray protons
In our model, most of the injected protons escape fromthe accretion flow without depletion due to the low effi-
Fig. 7.—
The contribution to the total intensity (red-thick lines)from different luminosity bins (thin lines). The blue-dashed, green-dotted, and magenta dotted-dashed lines show the fluxes frombright, middle, and faint parts, respectively. See text for defini-tion of the each part. The black squares show the observed dataof neutrino signals. The top and bottom panels show the intensityfor B2 and B3, respectively. ciency of pion production f π . .
2. Here, we discuss theeffects of escaping protons.Assuming that the Universe is filled with CR protons,we can estimate the CR flux as in the neutrino flux.Figure 8 shows the estimated flux of CR protons formodels B1, B2, B3, and B4. This flux of the escap-ing protons is much lower than observed CR flux for10 . eV < E p < eV for all the models. Althoughthe escaping proton luminosity has weaker dependenceon ˙ m than that of neutrino luminosity, the bright part isdominant for the CR proton flux.We note that it is unclear whether CRs of E p ∼ eV are able to arrive at the Earth from LLAGN. Infact, the magnetic fields in the intergalactic medium(IGM) prevent the protons from traveling straightly, sothat the distant sources cannot contribute to the CRflux. The diffusion length of CR protons during the cos-mic time is estimated to be ∼ B − / − E / p, l / , Mpc( E p . eV), where we use B − = B/ (10 − Gauss), E p, = E p / (10 PeV), and the coherence length l coh , = l coh / (100 kpc) (e.g., Ryu et al. 2008). We consider thatthe CRs are in cosmic filaments and/or the galaxy groupswith Kolmogorov turbulence, and ignore the cosmic ex-pansion. In addition, our Galaxy is located in the localeutrino and CR Emission from RIAFs in LLAGN 11group, where the magnetic fields are probably strongerthan the usual IGM. These magnetic fields can poten-tially reduce the UHECR flux of E p ∼ eV arrivingat the Earth (Takami et al. 2014). We should take theeffects of these magnetic fields into account to discuss thearrival CR flux in detail.The escaping protons would diffuse in host galaxiesof LLAGN, and interact with gas in the interstellarmedium (ISM) inside the galaxies. The pion produc-tion efficiency of pp inelastic collisions in the ISM is esti-mated to be f π, gal ≃ K pp n p, gal σ pp ct trap ∼ × − E − . p, ,where n p, gal ∼ − is the mean nucleon densityin the host galaxy, t trap = h / κ is the trapping timein the galaxy. We use the scale height h ∼ κ ∼ × ( E p / . cm s − . Note that the centralregions of galaxies have much denser gases than the meanvalues of ISM that we use for the estimation. The inter-action of escaping protons in such dense cores of galaxiesmight be important. The escaping protons are expectedto be confined in IGM. These protons are likely to inter-act with the protons or photons. The efficiency of pionproduction in IGM is not low, typically ∼ − below100 PeV (Murase et al. 2013), which is likely to be moreimportant than the reactions in ISM. These processesmight affect the diffuse neutrino flux. Constraints on neutron loading in the jet
Toma & Takahara (2012) proposed a mass loadingmodel to relativistic jets by relativistic neutrons madein the accretion flows. They consider that the rela-tivistic neutrons whose Lorentz factor γ n ∼ L jet . × − ˙ M c and the massof ˙ M jet . × − ˙ M . This estimate results from theassumption of L n ∼ .
03 ˙
M c , where L n is the total lu-minosity of neutrons from the accretion flow. The totalluminosity of neutrons in our model is estimated as L n ∼ f n η cr ˙ M c , (35)where f n is the neutron generation efficiency. The neu-tron generation efficiency is the same order of the pionproduction efficiency, f n ∼ f π . .
2. From the fittingof the diffuse neutrino flux, we obtain η cr ∼ .
01. Theseresults restrict L n . × − ˙ M c , which is much lowerthan their assumption. In addition, resultant spectraof relativistic protons that are accelerated via stochasticacceleration are quite hard. This causes the differen-tial luminosity and mass of the neutrons with γ n ∼ ζ ≤ .
03 and q = 5 /
3, LLAGN can-not emit the neutrinos of E ν &
30 TeV. However, it is stillnot easy to achieve the required value of f n η cr ∼ . η cr should be less than 0.1 in or- Fig. 8.—
The maximum flux of the diffuse CR protons. The thicklines show the CR flux for B1 (red-solid), B2 (green-dashed), B3(blue-dotted), and B4 (magenta dotted-dashed). The thin-black-solid line shows the observed CR flux (e.g., Becker 2008). Thedashed line (B2) almost overlaps the dot-dashed line (B4). der to keep the structure of the RIAFs (Kimura et al.2014). Another reason is that the nature of collisionlessplasma requires that the density should be low and limit f n . .
3. Thus, we need an optimized situation for neu-tron generation in order that the neutron injection modelworks as a jet mass loading mechanism. DISCUSSION
Comparison to other AGN models
In this work, we considered one of the AGN coremodels, in which CR acceleration and neutrino pro-duction occur in the vicinity of SMBH. Contrary tothis work, in the previous literature, CR accelerationaround the standard thin disk is assumed. However,since the disk plasma is typically collisional, faster dis-sipation is needed. The shock in accretion flows (e.g.,Protheroe & Kazanas 1983; Stecker et al. 1991) or be-tween blobs (Alvarez-Mu˜niz & M´esz´aros 2004), and elec-tric field acceleration (Kalashev et al. 2014) have beenspeculated as underlying acceleration mechanisms. Theacceleration mechanism at such inner regions is very un-certain. For the efficient shock acceleration mechanismto work, τ T . T max , thetypical neutrino energy is estimated to be E ν ≈ .
05 0 . m p c ¯ ǫ pk k B T max ∼
400 TeV (cid:18) k B T max
20 eV (cid:19) − . (36)Hence, for an averaged accretion disk spectrum observedin quasars, a suppression around sub-PeV energies isexpected (Dermer et al. 2014). In principle, it is pos-sible to have lower-energy neutrinos by assuming high-temperature disks ad hoc. However, in such models, all2 Kimura, Murase, & Tomathe relevant parameters (the CR normalization, spectralindex, maximum energy and disk temperature) are essen-tially free parameters. Also, since gamma rays shouldnot escape because of the large optical depths for pairproduction, these models should be regraded as hid-den neutrino source models . Note that the neutrino lu-minosity will be higher for AGN with higher disk lu-minosities. Then, the well-observed X-ray luminosityfunction (Ueda et al. 2003) suggests that neutrino emis-sion is dominated by AGN with L X & erg s − (Murase et al. 2014).In the vicinity of the standard disk, pγ interac-tions are usually the most important process (e.g.,Szabo & Protheroe 1994; Alvarez-Mu˜niz & M´esz´aros2004). Also, the heavy jet has a problem in its ener-getics (Atoyan & Dermer 2003; Murase et al. 2014). Onthe other hand, there are some discussions on pp scenar-ios in radio galaxies (Becker Tjus et al. 2014), assumingthe existence of dense knots with N H ∼ cm − . IfCRs are supplied by jets, the efficient neutrino produc-tion would significantly be diluted by their volume fillingfactor. Also, note that steep spectra s ν & . ≫ pp scenarios, which has been constrained bymulti-messenger data (Murase et al. 2013). Gamma rays from RIAFs
If neutrinos are produced by pion decay, gamma raysare also inevitably produced. The generated spectrumand luminosity of these gamma rays are similar to those of the neutrinos. However, high-energy gamma rays areabsorbed by soft photons through γ + γ → e + + e − ,so that the observed spectra of the gamma rays can bedifferent from those of the neutrinos. In our model, in-ternal absorption inside sources is relevant, and electro-magnetic cascades are initiated. The emergent spectraare expected to have a break at the energy where theoptical depth of pair production τ γγ = 1. We estimatethe optical depth by τ γγ ( E γ ) ∼ . σ T RN γ ( ǫ t ) ǫ t , (37)where ǫ t ≃ ( m e c ) /E γ is the energy of the soft pho-tons (e.g., Coppi & Blandford 1990). In our model,bright LLAGN, such as A3 model, have the cutoff energy E γ, cut ∼ . E γ, cut ∼
24 GeV, as tabulated in Table 1. Thismeans that bright LLAGN emit only multi-GeV photonsand cannot emit TeV photons. In this subsection, al-though we defer detailed studies of cascade emission, wehere give the order-of-magnitude estimate on expectedgamma-ray signatures, which may serve as tests of theLLAGN model.Recent observations by
Fermi show that the dif-fuse isotropic γ -ray background (IGB) intensity is ∼ × − GeV cm − s − sr − at E γ ∼ ∼ × − GeV cm − s − sr − at E γ ∼
100 GeV(Ackermann et al. 2014). This sets model-independentstrong bounds on pp scenarios (Murase et al. 2013), andspectral indices should be harder than s ν ∼ .
2. TheLLAGN model can avoid these constraints due to tworeasons. First, the neutrino spectrum is harder than s ν ∼ .
0, since stochastic acceleration or magnetic re-connection mechanisms predict hard spectra comparedto the diffusive shock acceleration mechanism. Thus, di-rect gamma rays do not contribute to the IGB. Second,GeV–TeV gamma rays may not escape, and LLAGN canbe regarded as hidden neutrino sources. In reality, thesituation depends on ˙ m , and GeV–TeV gamma rays canbe produced via cascades. The maximum GeV–TeV fluxcan be estimated by assuming that all gamma rays es-cape and get cascaded in intergalactic space. Even inthis case, noting that the gamma-ray flux is comparableto the neutrino flux, the estimated IGB flux is expectedto be . − GeV cm − s − sr − for B2 and B4, and . × − GeV cm − s − sr − for B1 and B3. How-ever, in RIAFs, pairs produced via γ + γ → e + + e − would lose mainly via synchrotron emission rather thaninverse-Compton emission, so the IGB at sub-TeV en-ergies is expected to be sufficiently lower than the ob-served intensity. More accurate calculation and examinethe contribution to the IGB remains as a future work.How could we test the LLAGN model presented here?Unfortunately, the averaged neutrino luminosity persource is too dim to detect individual sources. Gamma-ray detections are also difficult although we expect thatfaint LLAGN can emit TeV photons. One of the ex-amples of LLAGN with RIAFs is Sgr A* at the Galac-tic Center, where thermal electrons at the inner regionemit the radio band photons of L γ ∼ erg s − at E γ ∼ × − eV (Yuan et al. 2003). We calculatethe inner part of the accretion flow of Sgr A* suchthat our model accounts for the observed radio lumi-nosity. We use the same parameters as those for modeleutrino and CR Emission from RIAFs in LLAGN 13A1 except for M BH = 4 × M ⊙ and ˙ m = 1 × − .Then, we find that protons are accelerated up to 10TeV. The differential proton luminosity in this modelis E p L E p ∼ × erg s − at E p ∼
10 TeV and E p L E p ∼ × erg s − at E p ∼
10 GeV. These escap-ing protons are expected to emit GeV–TeV gamma raysvia pp reaction with surrounding materials, but this lu-minosity seems too low to explain the observed GeV–TeVgamma-ray fluxes (Liu et al. 2006; Chernyakova et al.2011). These protons do not contribute to the observedflux of CR protons either. The neutrino flux from SgrA* in this model is E ν F E ν ∼ × − erg cm − s − with peak energy E ν ∼ . E γ ∼ d ≃ . L bol ∼ × erg s − with M BH ∼ × M ⊙ (Eracleous et al. 2010). This lumi-nosity can be obtained by our model with reference pa-rameters except M BH = 6 . × M ⊙ . The peak neutrinoflux in this model is E ν F E ν ∼ × − GeV cm − s − at E ν ∼ . P eV , which is too dim to be detectedby IceCube. The pair-production cutoff energy for thismodel is E γ, cut ∼ E γ F E γ ∼ × − erg cm − s − . Thus, M81 couldbe detectable by Fermi or Cherenkov telescopes withlow thresholds such as Major Atmospheric Gamma-rayImaging Cherenkov Telescope (MAGIC) and CherenkovTelescope Array (CTA), although the estimate above istoo simple to predict the correct flux. Multimessengerstudies are relevant to test the model and more precisecalculations of photon spectra will be presented as a fu-ture work. SUMMARY
We studied particle acceleration and associated neu-trino emission from RIAFs of LLAGN. Various accelera-tion mechanisms have been suggested. In this work, fordemonstration, we consider stochastic acceleration, forwhich we can calculate spectra of escaping particles bysolving the Fokker-Planck equation. We modeled target photon fields in RIAFs by calculating inverse-Comptonemission, based on the one-zone approximation. Then wecompared acceleration, escape, and cooling time scales,and found that in LLAGN, proton acceleration is typi-cally suppressed by diffusive escape rather than coolingprocesses. We also found that LLAGN can have the pro-tons up to more than 10 PeV for reasonable ranges of˙ m , M BH , ζ and q . Then, the pp or pγ production maylead to PeV neutrinos. Note that production of ultra-high-energy CRs is not expected in this model. The CRacceleration efficiency is highly uncertain, so we treatedit as a free parameter assuming that the total luminosityof escaping CRs is equal to η cr ˙ M c , including both diffu-sive and advective escape. Then, the luminosity of CRsescape via diffusion is estimated to be around 1 × erg s − in our reference model. We calculate associatedneutrino emission, and found that high-energy neutrinoproduction occurs mainly via pp interactions, and themeson production efficiency is typically the order of 1%.The neutrino spectrum is hard since CRs are assumed tobe accelerated via the stochastic acceleration mechanism.We also calculated the diffuse neutrino intensity by us-ing the H α luminosity function of LLAGN and assumingno redshift evolution. Interestingly, we found that theobserved IceCube data can be fitted for reasonable pa-rameters if ∼
1% of the accretion luminosity is carriedby CRs. This fraction gurantees our assumption thatthe CRs do not affect the dynamical structure of RIAFs(Kimura et al. 2014). The number density of LLAGN isthe order of ∼ − − − Mpc − , which is much higherthan those of radio-loud AGN including blazars. Sincethe spectrum is hard, this result does not contradict thediffuse gamma-ray bound (Murase et al. 2013) and ob-served CR flux.Whereas RIAFs of LLAGN can provide interestingtargets of high-energy neutrino and gamma-ray obser-vations, unfortunately, there are many uncertainties inthe model. First, parameters related to acceleration areuncertain, although values we adopt are often used inthe different literature such as gamma-ray bursts. How-ever, although we considered stochastic acceleration, oneshould keep in mind that our neutrino flux calculationscan be applied to different possibilities such as accelera-tion via magnetic reconnections. For more reliable pre-dictions, we need better knowledge on the distributionof non-thermal particles, which could be achieved by fu-ture particle-in-cell simulations. Second, the luminosityfunction of LLAGN is quite uncertain due to their faint-ness. Obviously, to estimate the diffuse neutrino inten-sity, more observational data on the shape of the luminos-ity function in the faint end and their redshift evolutionare needed, as well as the mass function of SMBHs hostedby LLAGN. In addition, contributions from RIAFs withthe critical mass accretion rate, at which RIAFs changeto the standard disk, may also be relevant. This impliesthe importance of understanding the physical relation-ship between LLAGN and Seyferts.One of the potentially interesting points of the LLAGNmodel is that one could explain the latest IceCube dataaround 10 TeV. The latest data suggest steeper in-dices of s ν ∼ . − .
5, which seems challenging formany models. Galactic sources may be responsible for .
100 TeV neutrinos, but it is premature to discuss such4 Kimura, Murase, & Tomaa two-component scenario due to the lack of compellinganisotropy. It could be explained by an exponential cut-off or spectral break of starburst galaxies, but hard in-dices of s ν ∼ ζ is very low.We thank Katsuaki Asano, Jun Kakuwa, Matt Kistler,Fumio Takahara, and Shuta Tanaka for fruitful discus-sion. S.S.K. thanks Yutaka Fujita, Kentaro Nagamine,and Hideyuki Tagoshi for continuous encouragement. K.M. thanks Rashid Sunyaev for discussion about the diskphoton field. This work is partly supported by Grant-in-Aid for JSPS Fellowships No. 251784 (S.S.K.). Thiswork is supported by NASA through Hubble FellowshipGrant No. 51310.01 awarded by the STScI, which is op-erated by the Association of Universities for Research inAstronomy, Inc., for NASA, under Contract No. NAS5-26555 (K.M.) The result of this work was presented byK. M. at the JSI workshop in November 2014. APPENDIX
CALCULATION METHOD FOR THE SPECTRUM FROM THERMAL ELECTRONS
Synchrotron and Bremsstrahlung
We use a fitting formula for the synchrotron and bremsstrahlung (Narayan & Yi 1995). Hear, we use the cgs units.The cooling rates by bremsstrahlung of electron-electron q ei and ion-electron q ee are represented as follows, q ei = 1 . × − n e F ei ( θ e ) , (A1) q ee = ( . × − n e θ / e (cid:16) . θ e + θ e − . θ / e (cid:17) ( θ e < . × − n e θ e [ln(1 . θ e + 1 . θ e > , (A2)where we use F ei = ( (cid:0) θ e π (cid:1) . (1 + 1 . θ . e ) ( θ e < θ e π [ln(1 . θ e + 0 .
48) + 1 .
5] ( θ e > . (A3)The emissivity of bremsstrahlung is related to the cooling rates as q ei + q ee = Z ∞ dνj ν, br . (A4)Approximating j ν, br = j exp( − hν/ ( k B T e )), where j does not depend on ν , we can write j ν, br = q br hk B T e exp (cid:18) − hνk B T e (cid:19) (A5)where q br = q ei + q ee is the cooling rate per unit volume. We assume the gaunt factor is unity.The synchrotron emissivity is j ν,sy = 4 πe n e ν √ cK (1 /θ e ) I ′ (cid:18) πm e cν eBθ e (cid:19) , (A6)where I ′ ( x ) = 4 . x / (cid:18) . x / + 0 . x / (cid:19) exp( − . x / ) . (A7) Radiative transfer
We define the optical depth for absorption as τ ν ≡ √ π κ ν H ∼ √ π κ ν R, (A8)where κ ν = j ν,br + j ν,sy πB ν (A9)eutrino and CR Emission from RIAFs in LLAGN 15is the absorption coefficient and B ν is the Plank function. The photon energy flux from a RIAF is represented as(Manmoto et al. 1997) F ν = 2 π √ B ν h − exp( − √ τ ν ) i , (A10)where we use Eddington approximation (Rybicki & Lightman 1979) when estimating the vertical energy flux. Theluminosity by the synchrotron and bremsstrahlung is estimated as L ν, = 2 πR F ν . (A11) Inverse Compton scattering
We calculate the spectrum of the inverse Compton scattering. Seed photons are the photon field by bremsstrahlungand synchrotron. Assuming homogeneous and isotropic distribution, the photon occupation number evolves by (cf.Coppi & Blandford 1990), dN γ ( ǫ ) dt ∆ ǫ = − N γ ( ǫ )∆ ǫ Z dγN e ( γ ) R c ( ǫ, γ )+ Z dγ Z dǫ ′ N e ( γ ) N γ ( ǫ ′ ) R c ( ǫ ′ , γ ) P c ( ǫ ; ǫ ′ , γ ) + ˙ N γ, ∆ ǫ − ˙ N γ,esc ∆ ǫ, (A12)where ǫ = hν/ ( m e c ), N γ ( ǫ ) is the differential number density of photon, R c ( ǫ, γ )[cm s − ] is the reaction rate of theelectrons of Lorentz factor γ and the photons of energy ǫ , and P c ( ǫ ; ǫ ′ , γ ) is the probability that the reactions by thephotons of energy ǫ ′ and electrons of energy γ create the photons of energy ǫ . We add the injection rate of seed photons˙ N γ, and escape rate of photons ˙ N γ,esc = N γ ( ǫ ) / ( R/c ). We calculate the steady state solution of this equation. Thedifferential number density can be expanded by the number of Compton scattering as N γ ( ǫ ) = N γ, + N γ, + ... , where N γ, is determine by the balance of the escape and the injection N γ, / ( R/c ) = (1 − τ T ) ˙ N γ, = L ǫ, / ( πR hν ) , (A13)where L ǫ, = ( m e c /h ) L ν, . Since optical depth of the flow is less than unity, we approximate the effect of first termof right-hand side of Equation (A12) by multiplying the factor 1 − τ T . We obtain the N γ,i by solving the followingequation, N γ,i ( ǫ ) R/c ∆ ǫ = (1 − τ T ) Z dγ Z dǫ ′ N e ( γ ) N γ,i − ( ǫ ′ ) R c ( ǫ ′ , γ ) P c ( ǫ ; ǫ ′ , γ ) , (A14)where we use the approximation of (1 − τ T ) again. We use the fitting formula for R c (Coppi & Blandford 1990) R ( ǫ, γ ) = 3 cσ T γǫ (cid:20)(cid:18) − γǫ − γ ǫ (cid:19) ln (1 + 2 γǫ ) + 12 + 4 γǫ − γǫ ) (cid:21) . (A15)From the energy conservation, the reaction with 4 γ ǫ/ > γ + ǫ does not occur. In this case, we take R ( ǫ, γ ) = 0. Weuse delta function approximation for P c , N γ ( ǫ ′ ) P c ( ǫ ; ǫ ′ , γ ) = N γ ( ǫ ′ )∆ ǫ ′ δ (cid:18) ǫ − γ ǫ ′ (cid:19) . (A16)From the treatment described above, we obtain N γ,i ( ǫ ) = Rc Z dγ γ N e ( γ ) N γ,i − (cid:18) ǫ γ (cid:19) R c (cid:18) ǫ γ , γ (cid:19) . (A17)We calculate N γ,i ( ǫ ) until the number density of the i times scattered photons is much less than that of thebremsstrahlung in all the range in which we are interested. REFERENCESAartsen, M. G., Abbasi, R., Abdou, Y., et al. 2013a, Science, 342,1242856—. 2013b, Physical Review Letters, 111, 021103Aartsen, M. G., Ackermann, M., Adams, J., et al. 2014, PhysicalReview Letters, 113, 101101—. 2015, Phys. Rev. D, 91, 022001Abdo, A. A., Ackermann, M., Ajello, M., et al. 2009, ApJS, 183,46Abramowicz, M. A., Czerny, B., Lasota, J. P., & Szuszkiewicz, E.1988, ApJ, 332, 646Ackermann, M., Ajello, M., Albert, A., et al. 2014, ArXive-prints, arXiv:1410.3696 Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al.2006, Science, 314, 1424Aharonian, F., Akhperjanian, A. G., Anton, G., et al. 2009a,ApJ, 695, L40—. 2009b, A&A, 503, 817Ahlers, M., & Murase, K. 2014, Phys. Rev. D, 90, 023010Ajello, M., Romani, R. W., Gasparrini, D., et al. 2014, ApJ, 780,73Alvarez-Mu˜niz, J., & M´esz´aros, P. 2004, Phys. Rev. D, 70, 123001Anchordoqui, L. A., Goldberg, H., Halzen, F., & Weiler, T. J.2004, Physics Letters B, 600, 2026 Kimura, Murase, & Toma