High-energy cosmic ray nuclei from tidal disruption events: Origin, survival, and implications
aa r X i v : . [ a s t r o - ph . H E ] S e p High-energy cosmic ray nuclei from tidal disruption events:Origin, survival, and implications
B. Theodore Zhang , , , Kohta Murase , , , , Foteini Oikonomou , , and Zhuo Li , Department of Astronomy, School of Physics, Peking University, Beijing 100871, China Kavli Institute for Astronomy and Astrophysics, Peking University, Beijing 100871, China Department of Physics, The Pennsylvania State University, University Park, Pennsylvania 16802, USA Department of Astronomy & Astrophysics, The Pennsylvania State University, University Park, Pennsylvania 16802, USA Center for Particle and Gravitational Astrophysics,The Pennsylvania State University, University Park, Pennsylvania 16802, USA and Yukawa Institute for Theoretical Physics, Kyoto, Kyoto 606-8502 Japan (Dated: September 27, 2017)Tidal disruption events (TDEs) by supermassive or intermediate mass black holes have been sug-gested as candidate sources of ultrahigh-energy cosmic rays (UHECRs) and high-energy neutrinos.Motivated by the recent measurements from the Pierre Auger Observatory, which indicates a metal-rich cosmic-ray composition at ultrahigh energies, we investigate the fate of UHECR nuclei loadedin TDE jets. First, we consider the production and survival of UHECR nuclei at internal shocks,external forward and reverse shocks, and nonrelativistic winds. Based on the observations of SwiftJ1644+57, we show that the UHECRs can survive for external reverse and forward shocks, anddisk winds. On the other hand, UHECR nuclei are significantly disintegrated in internal shocks, al-though they could survive for low-luminosity TDE jets. Assuming that UHECR nuclei can survive,we consider implications of different composition models of TDEs. We find that the tidal disruptionof main sequence stars or carbon-oxygen white dwarfs does not successfully reproduce UHECR ob-servations, namely the observed composition or spectrum. The observed mean depth of the showermaximum and its deviation could be explained by oxygen-neon-magnesium white dwarfs, althoughthey may be too rare to be the sources of UHECRs.
I. INTRODUCTION
Cosmic rays with energy larger than 10 eV are re-ferred to as ultrahigh-energy cosmic rays (UHECRs), andtheir origin is still largely unknown [1–3]. The observedspectrum of UHECRs has a cutoff at energy around ∼ × eV [4–6]. The flux suppression is consis-tent with the prediction of the Greisen-Zatsepin-Kuzmin(GZK) effect due to the interaction between UHECRsand cosmic microwave background (CMB) photons [7, 8].A key clue to the origin of UHECRs is their composition.The primary mass of UHECRs can be inferred from thedistributions of the depth of the shower maximum, X max .The data from the Telescope Array (TA) and Auger areconsistent within systematic and statistical uncertainties[9]. The analysis from Auger suggests that the composi-tion of UHECRs is dominated by light nuclei at energyaround 10 . eV and becomes gradually heavier with in-creasing energy up to 10 . eV [10, 11]. The distributionsof X max are difficult to explain with a mixture of protonsand iron nuclei over the whole energy range if up-to-datehadronic interaction models are correct. Rather, the bestfit is reached by including a fraction of intermediate massnuclei [11].The interpretation of the UHECR composition is stillunder intense debate, and we avoid discussion that de-pends on their details. However, considering the largerdetection area and lower sampling bias of Auger, it is rea-sonable to assume that UHECR composition gets grad-ually heavier at the highest energies. There are not somany candidate sources which satisfy the Hillas crite- rion to accelerate cosmic rays (CRs) to ultrahigh en-ergies [1], and the origin of heavy nuclei has now be-come an interesting question. In relativistic jets of ac-tive galactic nuclei (AGN) (e.g., [12–16] and referencestherein), heavy nuclei would be supplied by an accre-tion disk. For structure formation shocks in galaxyclusters (GCs) [12, 17, 18], heavy nuclei can be pro-vided from the intracluster medium. However, for theseobjects, we typically expect a composition similar tothe solar composition, unless reacceleration is invoked[19, 20]. Another possibility is that heavy nuclei canbe synthesized in the stellar interiors or outflows fromthe deaths of massive stars. This scenario is appropri-ate for gamma-ray bursts (GRBs) [21–25] (see [26, 27]for protons), low-luminosity GRBs and transrelativisticsupernovae [21, 28–31], and newborn pulsars and mag-netars [32–34]. However, UHECR nuclei can be depletedbefore they escape the source environment due to the in-teraction with background particles or photons, so the“UHECR survival problem” is important to discuss theorigin of UHECR nuclei [13, 21–23, 34].In this work, we revisit tidal disruption events (TDEs)as sources of UHECR nuclei. A TDE is a luminous flarelasting for months to years, which occurs in a galaxy nu-clear region [35]. When a star gets very close to thecentral black hole, it is disrupted if the tidal force isgreater than the star’s self-gravity. During the disrup-tion process, nearly half of the stellar debris falls backinto the vicinity of the black hole, and the rest becomesunbound from the system. The accretion flow under-takes a fast energy dissipation and circularization pro-cess, and then it forms an accretion disk around the cen-tral black hole [36, 37]. The TDE may be accompaniedby the emergence of a relativistic jet during the super-Eddington accretion phase [38–44], and TDEs have beenproposed to accelerate particles to ultrahigh energies [45–47]. Two scenarios are the most popular: disruption ofa main-sequence (MS) star by a supermassive black hole(SMBH) and the disruption of a white dwarf (WD) by anintermediate mass black hole (IMBH). The latter couldbe especially interesting as sources that can inject heavynuclei.This paper is organized as follows. In Sec. II, we showthat both protons and nuclei can be accelerated to ultra-high energies, and we study conditions for the survivalof UHECR nuclei in jetted TDEs. For the accelerationsites, we consider internal shocks, external reverse andforward shocks, and nonrelativistic winds. In Sec. III,for different composition models, we study the propaga-tion of UHECR nuclei using the public code CRPropa 3[48], and compare the results to experimental results byAuger.Throughout the paper, we use cgs units, and adoptnotations such as Q x ≡ Q/ x . The cosmological pa-rameters we assume are H = 67 . − Mpc − ,Ω m = 0 . Λ = 0 .
685 [49].
II. SURVIVAL OF COSMIC-RAY NUCLEI INTDE OUTFLOWS
We discuss possible composition models in the nextsection. In this section, we consider the fate of UHECRnuclei, following Ref. [21].Jetted TDEs (e.g., Swift J1644+57) show clear signa-tures of nonthermal emission in a wide range of wave-lengths from radio to x rays [50–52]. Diffusive shockacceleration in collisionless shocks is the most popularmechanism of production of nonthermal particles. Thefirst-order Fermi process is predicted to have a power-lawdistribution of accelerated particles within the “test par-ticle” approximation. The acceleration time scale canbe expressed as t acc = ηr L /c , where r L = E A /ZeB isthe Larmor radius of a particle with energy E A , chargenumber Z , and B the comoving-frame magnetic fieldstrength. The factor η depends on details of the turbu-lence. The minimum value of η ∼ −
10 can be achievedin the Bohm limit [53], and η = 1 is used to demonstrateour results in this section.The maximum acceleration energy is determined by t acc ≤ min( t dyn , t syn , t Aγ ), where t dyn ≡ R/ Γ βc is the dy-namical time scale; t syn = 3 m A c Γ / (4 σ T Z m e U B E A ) isthe synchrotron cooling time scale ( σ T is the Thompsoncross-section, U B = B / π is the magnetic energy den-sity, E A is the particle energy); and t Aγ is the energy losstime scale for protons (photomeson interaction) and nu-clei (photodisintegration interaction). We can estimate -2 -1 ¯ε[GeV] -2 -1 σ A γ [ m b ] FeSiO Hep
FIG. 1. Photonuclear and photomeson production cross sec-tions for five chemical species, which are used in this work:Fe, Si, O, He, and proton, as a function of NRF target photonenergy [48, 54–56]. For simplicity, the superposition model isassumed for the photomeson production. the energy loss time scale using the following formula, t − Aγ ( E A ) = c γ A Z ∞ ¯ ε th d ¯ εσ Aγ (¯ ε ) κ Aγ (¯ ε )¯ ε Z ∞ ¯ ε/ γ A dε ε dndε , (1)where γ A is the Lorentz factor of UHECRs with massnumber A , ¯ ε th is the threshold energy measured in therest frame of initial nucleus (NRF), and d n/ d ε is the dif-ferential number density of target photons. Here, σ Aγ isthe photohadronic cross section related to the photome-son or photodisintegration process, and κ Aγ is the inelas-ticity of each process. We show photonuclear and pho-tohadronic cross sections for five typical chemical speciesas a function of NRF target photon energy ranging from ∼
10 MeV to ∼ GeV in Fig. 1; see Appendix A.Defining the optical depth as f Aγ ≡ t − Aγ /t − , we canexpect the survival of UHECR nuclei when f Aγ <
1. Thevalue of the interaction time scale t Aγ − int is equal to theenergy loss time scale when κ Aγ = 1, and the relatedoptical depth is τ Aγ ≡ t − Aγ − int /t − . A. Internal shock model
Nonthermal hard x-ray emission comes from inter-nal energy dissipation in the inner relativistic jet [51].In the internal shock region, fast moving ejecta maycatch up with slower ejecta, and a substantial amountof kinetic energy of the relativistic jet may be con-verted into internal energy. We assume the radiuswhere internal collisions take place is R IS = Γ cδt =3 × Γ δt cm, with Γ = 10Γ the Lorentz factorof the relativistic jet and δt ∼ δt s is the x-rayvariability time scale [51]. The median x-ray luminos-ity of Swift J1644+57 is L X , iso = 8 . × erg s − ,which is well above the Eddington luminosity L Edd =1 . × M BH , erg s − of a 10 M ⊙ black hole [51].In our analysis, the total isotropic luminosity is set to L iso = 10 erg s − . A fraction ǫ B of the total en-ergy of the outflow is converted into magnetic energy B / π = ǫ B L iso / πR Γ c . The magnetic field strengthin the jet comoving frame is B = (2 ǫ B L iso /R Γ c ) / ≃ L / , ǫ / B, − R − , . Γ − G. The maximum accelera-tion energy can be achieved under the condition t acc ≤ t dyn ; we have E A, dyn ≃ Γ η − ZeB ( R IS / Γ) ∼ . × Zη − L / , ǫ / B, − Γ − eV measured in the observerframe, where R IS / Γ is the comoving frame shell width.The maximum acceleration energy is also limited by thesynchrotron energy loss t acc < t syn ; we have E A, syn ∼ . × A Z − / η − / L − / , ǫ − / B, − Γ / R / , . eV mea-sured in the observer frame. To estimate the effect ofphotonuclear and photohadronic interactions, we use alog-parabola model to fit the high-luminosity state SEDof Swift J1644+57 [51, 57].Our results for the internal shock scenario are shownin Fig. 2. We consider two time scales, one is the interac-tion time scale, t Aγ − int , and the other is the energy-losstime scale, t Aγ . We see that t Aγ − int and t Aγ are shorterthan t dyn . Our calculation suggests that most CR nucleiin the internal shock region will be disintegrated beforethey escape, so it is difficult for CR nuclei to survivefor luminous TDEs like Swift J1644+57. Note that evenpartial survival is difficult for this parameter set (i.e., f Aγ & L iso = 10 . erg s − . The co-moving frame magnetic field strength is estimated to be B ≃ L / , . ǫ / B, − R − , Γ − G. The results are shownin Fig. 3, and we find that UHECR nuclei with energyup to ∼ eV can survive before they escape from thesource. Note that this case almost allows the completesurvival of nuclei (i.e., τ Aγ . τ Aγ & f Aγ .
1, UHECRnuclei partially survive [58] and secondary nuclei affectthe initial composition of UHECRs.
B. Reverse shock model
The radio afterglow of Swift J1644+57 was observeda few days after the trigger of the Swift BAT obser-vation [52], and the continued observation extends to ∼
500 days [59, 60]. The nonthermal radio emissionis consistent with synchrotron radiation from the stan-dard external shock model [60, 61]. Here, we assumethe jet duration time is t j = 10 s based on the obser-vation of Swift J1644+57 [51]. The isotropic equivalentkinetic energy of the jet is L iso ∼ erg s − on av-erage, and it will decrease following L iso ∝ t − / afterthe time t j . We estimate the total injection energy as
10 12 14 16 18 20−8−6−4−2024 l og (t − [ s − ] ) t −1dyn t −1acc,Fe t −1acc,p t −1syn,p FeSiOHep10 12 14 16 18 20 log(E A [eV]) −8−6−4−2024 l og (t − [ s − ] ) t −1dyn t −1acc,Fe t −1acc,p t −1syn,p FIG. 2. Various time scales for five typical chemical species(Fe, Si, O, He, and proton) in the internal shock regionas a function of particle energy (measured in the observerframe). We show interaction time scales in the upper paneland energy-loss time scales in the lower panel. The thin(thick) black line represents the acceleration time scale forthe proton (Fe). We show the synchrotron cooling time scalefor protons as the dotted black line. The parameters are L iso = 10 erg s − , Γ = 10, ǫ B = 0 .
1, and r = 3 × cm. E iso = 2 L iso t j = 2 × erg. The relativistic jet isdecelerated when it has swept up a significant amountof circumnuclear medium (CNM). There are two shocksformed in the deceleration radius, one is the reverse shockpropagating back into the relativistic jet, and the otheris the forward shock propagating into the CNM. For sim-plicity, we adopt a constant CNM density, ̺ = 1 cm − .We calculate the evolution of reverse shock followingthe same method as the one used for GRBs [21, 62, 63].The Lorentz factor of the shocked ejecta (relative tounshocked CNM) can be estimated as Γ ≈ Γ / (1 +2Γ ( ̺/n j ) / ) / assuming the pressure equilibrium atthe contact discontinuity. We assume Γ = 10, and thedensity of the jet n j = E iso / (4 πm p c Γ (Γ ∆) r ), where E iso is the isotropic equivalent energy of the relativisticoutflow, and ∆ is the thickness of the ejecta (relativeto the central black hole). The thickness of the ejectais estimated to be ∆ × ≈ ∆ ≡ ct j . The shock com-pletely crosses the ejecta at t = t × . We write the ra-
10 12 14 16 18 20−8−6−4−2024 l og (t − [ s − ] ) t −1dyn t −1acc,Fe t −1acc,p t −1syn,p FeSiOHep10 12 14 16 18 20 log(E A [eV]) −8−6−4−2024 l og (t − [ s − ] ) t −1dyn t −1acc,Fe t −1acc,p t −1syn,p FIG. 3. Same as Fig. 2, but the parameters are L iso =10 . erg s − , Γ = 10, ǫ B = 0 .
1, and r = 10 cm. dius at that time as r × ≃ . × E / , . t / j, ̺ − / cmand the Lorentz factor at the crossing time is Γ × ≃ . E / , . t − / j, ̺ − / . The magnetic field strength is esti-mated to be B × = (32 πǫ rB n j m p c (Γ × − × + 3 / / ,where a fraction ǫ rB of the post-shock internal energy isconverted into magnetic energy. The comoving frameLorentz factor of the electrons is γ re = ǫ re f re p − p − m p m e (Γ × − p = 2 . ǫ re = 0 .
02, and f re = 0 .
1. Here f re representsthe number fraction of accelerated electrons.Then, the peak synchrotron frequency is written as ν rm, ob ≃ × g (Γ × ) g (3 . E / , . × ǫ rB, − / ǫ re, − . f re, − − ̺ / t − / j, Hz , (2)where g (Γ × ) = Γ × (Γ × − / (Γ × + 3 / / . Thesynchrotron self-absorption frequency can be estimatedwhen τ ( ν a ) = 1; we have ( ν a < ν m ) ν ra, ob ≃ . × g (Γ × ) g (3 . E / , . ǫ rB, − / × ǫ re, − . − f re, − / ̺ / t − / j, Hz , (3)where g (Γ × ) = g (Γ × )(Γ × − − / (Γ × + 3 / − / .Further, we can express the comoving frame peak lumi- nosity per unit energy from the reverse shock as L rǫ, max = 12 π ~ N re f re √ e B × m e c ≃ . × g (Γ × ) g (3 . E / , . × ǫ rB, − / f re, − ̺ / t − / j, s − , (4)where g (Γ × ) = (Γ × − / (Γ × + 3 / / and N re = E iso / Γ m p c is the number of electrons in the reverseshock. The photon spectrum in the reverse shock is( ν a < ν m < ν c ) dn r dε = L rε, max πr × cε m ( ε a ε m ) − / ( εε a ) ( ε < ε a )( εε m ) − / ( ε a ≤ ε < ε m )( εε m ) − ( p − / − ( ε m ≤ ε < ε c )( ε c ε m ) − ( p − / − ( εε c ) − p/ − ( ε ≥ ε c ) , (5)where ε m = ε m, ob / Γ × is the peak photon energy in thereverse shock comoving frame and ε c is the electron syn-chrotron cooling energy. The maximum acceleration en-ergy of nuclei can be derived when t acc = t dyn , E A, max = Γ × η − ZeB × ( r × / Γ × ) ≃ . × Zη − g (Γ × ) g (3 . × E / , . ǫ rB, − / t − / j, eV . (6)We estimate various time scales in the reverse shockmodel, and our results are shown in Fig. 4. We findthat the interaction time scale t Aγ − int is longer than thedynamical time scale t dyn , which means that the opticaldepth τ Aγ <
1. Our results suggest that UHECR nucleican survive in the reverse shock model.
C. Forward shock model
Once the relativistic jet enters the CNM, the deceler-ation and transition to the Blandford-Mckee (BM) self-similar regime occurs. The evolution of the Lorentz fac-tor and shock radius are described as Γ( t ) ∝ t − / and r ( t ) ∝ t / , where we assume a constant CNM den-sity. The evolution of the downstream magnetic fieldfollows B ∝ t − / . A fraction of electrons are acceler-ated in the shock, with a minimum Lorentz factor, γ m ,and the distribution of accelerated electrons is denotedas dN e /dγ e ∝ γ − p with p = 2 .
2, and minimum Lorentzfactor γ m ∝ t − / . In the following, we use the approxi-mation Γ ≫ ν m, ob ≃ . × E / , . ǫ / B, − . ǫ e, − . f − e, − t − / Hz , (7)
17 18 19 20 21 22−12−10−8−6−4−20 l og (t − [ s − ] ) t −1dyn t −1acc,Fe t −1acc,p t −1syn,p FeSiOHep17 18 19 20 21 22 log(E A [eV]) −12−10−8−6−4−20 l og (t − [ s − ] ) t −1dyn t −1acc,Fe t −1acc,p t −1syn,p FIG. 4. Various time scales for five typical chemical species(Fe, Si, O, He, and proton) in the reverse shock region as afunction of particle energy (measured in the observer frame).We show the interaction time scales in the upper panel andenergy-loss time scales in the lower panel. The thin (thick)black line represents the acceleration time scale for the proton(Fe). We show the synchrotron cooling time scale for theproton as the dotted black line. The parameters we use are L iso = 10 erg s − , t j = 10 s, Γ = 10, ǫ rB = 0 . ǫ re = 0 . f re = 0 . The self-absorption frequency (assuming ν a < ν m < ν c )is ν a, ob ≃ . × E / , . ǫ / B, − . ǫ − e, − . f / e, − ̺ / Hz . (8)The comoving frame peak luminosity per unit energyis L ǫ, max = 12 π ~ N e f e √ e Bm e c ≃ . × E / , . ǫ / B, − . f e, − ̺ / t / s − , (9)where N e = (4 π/ r ̺ is the total number of swept-upelectrons. The comoving frame photon spectrum can be derived in the slow cooling case ( ν a < ν m < ν c ), dn f dε = L ε, max πr cε m ( ε a ε m ) − / ( εε a ) ( ε < ε a )( εε m ) − / ( ε a ≤ ε < ε m )( εε m ) − ( p − / − ( ε m ≤ ε < ε c )( ε c ε m ) − ( p − / − ( εε c ) − p/ − ( ε ≥ ε c ) , (10)where ε m = ε m, ob / Γ is the peak photon energy in theforward shock comoving frame and ε c is the electron syn-chrotron cooling energy.In the forward shock model, the observed maximumparticle energy is achieved when t acc = t dyn , with E A, max = Γ η − ZeB ( r/ Γ) (11) ≃ . × eV Zη − E / , . ǫ / B, − . ̺ / t − / , where it is assumed that the upstream magnetic field isamplified and comparable to the downstream value. Theshock is mildly relativistic and this could be achievedby CR streaming instabilities [64, 65]. Or one coulduse the downstream field if particles are accelerated bythe second-order Fermi acceleration mechanism [66]. Weshow various time scales in Fig. 5. We expect the sur-vival of UHECR nuclei at the forward shock, because theinteraction time scale, t Aγ , is longer than the dynamicaltime scale, t dyn . We also show the evolution of particlemaximum acceleration energy and optical depth f Aγ inFig. 6. We found that it is easier for UHECR nuclei tosurvive at later times.However, it may be difficult to amplify the magneticfield in the upstream region, especially if the shock isultrarelativistic [63, 67, 68]. In this case, the magneticfield in the shock region should be similar to the CNMmagnetic field ˆ B cnm = 10 − ˆ B cnm , − . The observed max-imum acceleration energy is estimated to be E A, max = Γ η − Ze Γ ˆ B cnm ( r/ Γ)= 4 . × eV Z E / , . ̺ − / ˆ B cnm , − t − / . (12)It is difficult to accelerate CRs to ultrahigh energies ifthe magnetic field is not amplified efficiently. D. Nonrelativistic wind model
Along with the formation of a relativistic jet, a non-relativistic outflow can be driven by the accretion disk.The radio emissions from two TDEs, ASASSIN14li [69]and XMMSL1J0740-85 [70], are consistent with the emis-sion from the interaction between a nonrelativistic out-flow and the CNM. We assume that the nonrelativisticoutflow has an ejected mass of M ej ∼ − M ⊙ and ve-locity v ej ∼ . c . The deceleration radius is estimated tobe r dec ≃ . × E / k, ̺ − / cm and the decelerationtime is t dec ≃ . × E / k, ̺ − / s. The shocked fluidenters into the Sedov-Taylor evolution phase after it isdecelerated by the external medium. In the adiabatic
17 18 19 20 21 22−12−10−8−6−4−20 l og (t − [ s − ] ) t −1dyn t −1acc,Fe t −1acc,p t −1syn,p FeSiOHeP17 18 19 20 21 22 log(E A [eV]) −12−10−8−6−4−20 l og (t − [ s − ] ) t −1dyn t −1acc,Fe t −1acc,p t −1syn,p FIG. 5. Various time scales for five typical chemical species(Fe, Si, O, He, and proton) in the forward shock region as afunction of particle energy (measured in the observer frame).We show the interaction time scales in the upper panel andenergy-loss time scales in the lower panel. The thin (thick)black line represents the acceleration time scale for the proton(Fe). We show the synchrotron cooling time scale for theproton as the dotted black line. The parameters we use are L iso = 10 erg s − , t = 10 s, Γ = 10, ǫ B = 0 . ǫ e = 0 . f e = 0 . case, the shock velocity and radius follow V s ∝ t − / and R s ∝ t / , respectively. We assume the accelerated elec-trons have power-law index p = 3. The synchrotron peakfrequency is ν m ≃ . × E k, ǫ / B, − ǫ e, − f − e, − ̺ − / ( t/t dec ) − Hz , (13)and the self-absorption frequency ( ν a > ν m ) is ν a ≃ . × E pp +4 k, ǫ p/ p +4 B, − ǫ p − p +4 e, − f − pp +4 e, − ̺ − p/ p +4 ( t/t dec ) − pp +4 Hz , (14)The photon spectrum in this case is( ν m < ν a < ν c ) dndε = L ε, max πR s cε a ( ε m ε a ) / ( εε a ) − ( εε m ) ( ε < ε m )( εε a ) / ( ε m ≤ ε < ε a )( εε a ) − ( p − / − ( ε a ≤ ε < ε c )( ε c ε a ) − ( p − / − ( εε c ) − p/ − ( ε ≥ ε c ) . (15) l og ( f A γ ) FeSiOHep log(t[s]) l og ( E A , m a x [ e V ] ) FIG. 6. The time evolution of maximum energy (lower panel)and optical depth f Aγ (upper panel) for five typical chemicalspecies: Fe, Si, O, He, and proton in the forward shock model.The optical depth is calculated when UHECR nuclei haveenergy 10 eV (measured in the observer frame). In the nonrelativistic case, the particle accelerationtime scale is t acc ≈ (20 / c E A /V s ZeBc in the Bohmlimit. Our results for various time scales are shown inFig 7. The dynamical time in the nonrelativistic case is t dyn = R s /V s . The maximum acceleration energy is lim-ited by t acc ≤ t dyn , with E A, max ≈ (3 / R s ZeBV s /c ∼ . × Z E / k, ǫ / B, − ̺ − / ( t/t dec ) − / eV. In the windmodel, CRs cannot be accelerated to ultrahigh energies. III. PROPAGATION OF ESCAPING NUCLEI
In the previous section, we found that UHECR nucleican survive for external shocks. For internal shocks, thesurvival is difficult for Swift J1644+57-like TDEs, butlower-luminosity TDEs allow UHECR nuclei to escapewithout significant disintegration. In this section, forsimplicity, we assume that UHECR nuclei survive andescape into intergalactic space and see consequences ofdifferent composition models [71].We assume that UHECR nuclei follow a power-law dis-
13 14 15 16 17 18 19 20 21 22−14−12−10−8−6−4−2 l og (t − [ s − ] ) t −1dyn t −1acc,Fe t −1acc,p t −1syn,p FeSiOHep13 14 15 16 17 18 19 20 21 22 log(E A [eV]) −14−12−10−8−6−4−2 l og (t − [ s − ] ) t −1dyn t −1acc,Fe t −1acc,p t −1syn,p FIG. 7. Various time scales for five typical chemical species(Fe, Si, O, He, and proton) in the wind model as a functionof particle energy. We show the interaction time scales in theupper panel and energy-loss time scales in the lower panel.The thin (thick) black line represents the acceleration timescale for the proton (Fe). We show the synchrotron coolingtime scale for protons as the dotted black line. The parame-ters we use are E k = 10 erg, V s = 0 . c , ǫ B = 0 . ǫ e = 0 . f e = 0 . tribution with an exponential cutoff as dN A ′ dE ′ = f A ′ N (cid:18) E ′ ZE (cid:19) − s esc exp (cid:18) − E ′ ZE ′ p, max (cid:19) , (16)where f A ′ is the number fraction of different particleswith mass A ′ at the same rigidity, N is determined bythe total energy per TDE event (see Appendix B), s esc isthe spectral index of ejected UHECR nuclei, and ZE ′ p, max is the maximum acceleration energy for particles withcharge number Z . We assume the minimum particle en-ergy E ′ A, min = Γ Am p c . Strictly speaking, this is justi-fied in the forward shock model but our main conclusiondoes not change by this assumption. A. Injection spectrum
In order for UHECR nuclei to originate from stellarmaterial (whether it is a MS star or WD), CR injections should occur inside jets or winds, which involve internalshocks and an external reverse shock. Details are highlyuncertain but there are various possibilities.The first possibility is to rely on the shock accelera-tion at internal shocks. In the diffusive shock accelerationmechanism, accelerated particles have a power-law distri-bution, dN/dE ∝ E − s acc , with a typical spectral index s acc ∼ . s esc <
2. For an expandingoutflow, one of the most conservative possibilities is toinvoke the direct escape of CRs (e.g., [82]), which leadsto s esc = s acc −
1. For s acc ∼
2, one may expect s esc ∼ E CR ∝ r − α E , the evolution of ex-ternal medium density is ̺ ∝ r − α ̺ , and the evolution ofshocked fluid bulk Lorentz factor is Γ ∝ r − α Γ . The min-imum energy of accelerated particles can be estimatedas E A, min ≃ Γ Am p c ∝ r − α min , and the maximum par-ticle energy is E A, max ≃ ZeBr ∝ r − α max , as we dis-cussed in the previous section. In the adiabatic expan-sion case, we have α min = 2 α Γ , α max = α Γ + (1 / α ̺ − α E = 2 α Γ + α ̺ −
3. We assume the acceleratedCRs spectrum ∝ E − s acc , and have spectrum ∝ E − s esc after escape from the acceleration site; Ref. [81] de-rived a simple analytic relation between s esc and s acc as s esc = s acc − ( α min ( s acc −
2) + α E ) /α max with the as-sumption that the number of ejected UHECRs in an en-ergy interval is the same as that of the CRs at radius r ,where the maximum acceleration energy is achieved. Letus adopt a typical value of power-law index s acc ∼ . α E = 0.If the CNM density is constant, then we have α ̺ = 0.The escaping particle spectral index can be calculated as s esc = − s acc + 12 = 1. However, if the CNM densitydecreases with radius as ̺ ∝ r − α ̺ , we have s esc = 1 . α ̺ = 1 and s esc = 1 . α ̺ = 2. The CNM gas densityin the galaxy nuclear region is still unclear. It can orig-inate from the stellar winds, and their profile dependson the detailed distribution of stars and star formationhistory in the galaxy nuclear region [84, 85]. We assumethat the CNM density is constant in our analysis, whichis reasonable for galaxy cores [84].Although the above arguments are rather speculative,it is possible to expect a hard spectral index for escapingCRs, and we use s esc ∼ E A, max according to the rigidity dependence, whichis valid in the internal shock model for relatively low-luminosity TDEs and the external forward shock model. B. Composition model
1. MS-SMBH tidal disruptions
In model I, we consider a solar-type MS star disruptedby a SMBH [35]. The present-day solar composition forH, He and metals is X = 0 . Y = 0 . Z =0 . X A for differentchemical species at the same rigidity are similar to theMS star. We assume that the maximum accelerationenergy is ZE p, max = 2 × Z eV.
2. WD-IMBH tidal disruptions
One of the models for Swift J1644+57 is a WD tidallydisrupted by a 10 M ⊙ IMBH [40, 50–52]. The radiiof WDs are smaller than those of MS stars, so WDsshould be tidally disrupted at smaller radii from thecentral BHs. We need the factor β g = R t R g > R t is thetidal disruption radius and R g is the Schwarzschild ra-dius. The upper limit on the BH mass is estimated tobe M BH . . × M ⊙ (cid:0) R WD cm (cid:1) / (cid:16) M WD . M ⊙ (cid:17) − , where R WD and M WD are the radius and mass of WDs [87].Most of the WDs are composed of carbon and oxygenas the result of helium burning in the core, and a frac-tion of WDs may contain neon and magnesium due tothe ignition of carbon. The hydrogen rich envelope ofWDs could be ejected due to helium thermal pulses inthe asymptotic giant branch phase. log(E A [eV]) E d N / d E [ e r g ] pHeC ONeMg SiFe FIG. 8. CR injection spectra for model I: MS stars tidallydisrupted by SMBHs. The mass fraction is X H = 73 . X He = 24 . X C = 0 . X N = 0 . X O = 0 . X Ne = 0 . X Mg = 0 . X Si = 0 . X Fe = 0 . ZE p, max = 2 × Z eV, and the spectral index is s esc = 1.The total CR injection energy is E CR = 10 erg. In model II-1, we consider that tidally disrupted WDsare carbon-oxygen WDs (CO-WDs). CO-WDs are theend point for low and intermediate mass stars in the massinterval from 0 . M ⊙ to 6 M ⊙ [88, 89]. The ratio of carbonto oxygen depends on the reaction C( α, γ ) O, whichis unresolved yet in nuclear physics [90, 91]. We assumethat CO-WDs have a roughly equal mass fraction forcarbon and oxygen X C = 0 . X O = 0 .
5. Stars with amass in the range of ∼ M ⊙ to ∼ M ⊙ may evolve toform oxygen-neon-magnesium WDs (ONeMg-WDs) [92,93]. In this case, the burning of carbon will not lead toexplosion or collapse. In model II-2, we consider thatONeMg-WDs with a mass fraction of X O = 0 . X Ne =0 .
76, and X Mg = 0 .
12 are tidally disrupted [94]. In Fig. 9,we show the CR injection spectra for model II-1 (upperpanel) and model II-2 (lower panel).
3. WD-IMBH with ignition
When a WD approaches a massive BH, the tidal com-pression and relativistic effects can enhance the WD cen-tral density and could trigger explosive nuclear burning[87, 95–102]. This kind of ignition has also been sug-gested as an alternative scenario for type Ia supernovae[95]. In this case, a fraction of nuclear explosive mattercan be accreted into the center BH and form an accretingflow [98].However, ignition and associated nucleosynthesis ofheavy nuclei have been questioned by more dedicated E d N / d E [ e r g ] CO NeMg SiFe log(E A [eV]) E d N / d E [ e r g ] ONe MgSi Fe
FIG. 9. CR injection spectra for WDs disrupted byIMBHs. Upper panel: Model II-1, CO-WDs with mass frac-tion X C = 0 . X O = 0 .
5. Lower panel: Model II-2, ONeMg-WDs with mass fraction X O = 0 . X Ne = 0 . X Mg = 0 . ZE p, max = 6 . × Z eV, and thespectral index is s esc = 1. The total CR injection energy is E CR = 10 erg. simulations, and details are still under debate [102]. Also,the rate of such events is uncertain and may be muchsmaller. On the other hand, it has been argued that thecomposition of such TDEs has been considered to ex-plain UHECRs [103], so we also consider this scenariofor completeness.In models III-1 and III-2, we adopt the numerical simu-lation results from Ref. [98]. Model III-1 is the case witha 0 . M ⊙ helium WD passing a 10 M ⊙ BH and modelIII-2 corresponds to a CO-WD (1 . M ⊙ ) approaching a500 M ⊙ BH . In Fig. 10, we show the CR injection spec-tra for model III-1 (upper panel) and model III-2 (lowerpanel). E d N / d E [ e r g ] HeCO NeMg SiFe log(E A [eV]) E d N / d E [ e r g ] HeCO NeMg SiFe
FIG. 10. CR injection spectra for WDs disrupted by anIMBH. Upper panel: Model III-1, 0 . M ⊙ helium WDs dis-rupted by 10 M ⊙ IMBHs. The mass fraction is X He = 77 . X C = 0 . X Si = 7 . X Fe = 0 . . M ⊙ CO-WDs disrupted by 500 M ⊙ IMBHs[98]. The mass fraction is X He = 15 . X C = 3 . X O = 10 . X Ne = 0 . X Mg = 0 . X Si = 21 . X Fe = 66 . ZE p, max =6 . × Z eV, and the spectral index is s esc = 1. Thetotal CR injection energy is E CR = 10 erg. C. Propagation in intergalactic space
We calculate the propagation of UHECR nuclei, usingthe public code CRPropa 3 [48, 104]. The main energy-loss process for UHECR nuclei during propagation is pho-todisintegration due to CMB photons and extragalacticbackground light (EBL) photons. The EBL photons havea larger effect on the propagation of lower-energy inter-mediate mass nuclei ( ∼ eV) [105, 106]. In thiswork, we adopt a semianalytic EBL model derived byRef. [107]. The details of the simulation are shown inAppendix B. The observed UHECR spectrum has been0accurately measured by Auger [108] and TA[109]. Weuse an empirical model, which describes the distributionof X max using the generalized Gumbel function G ( X max )to get the mean value of the depth of shower maximum h X max i and σ ( X max ) [110]. In this work, we adopt theEPOS-LHC hadronic interaction model [10, 111].In Fig. 11, we show the results derived from model I.This model can fit the observed UHECR spectrum mea-sured by Auger, but failed to fit h X max i and σ ( X max ); thefinal spectrum is proton dominated in nearly the entireenergy range ( < eV). Although the situation is stillunder debate due to significant uncertainties in hadronicinteraction models, this proton dominated scenario seemsin strong tension with Auger results [11].In Fig. 12, we show the results from model II-1, wheremost UHECRs are carbon and oxygen. The energy-losslengths of C and O nuclei are very short, so most ofthe observed C and O nuclei should come from nearbysources. There should be a large fraction of secondaryprotons and helium generated during the propagation.We find that model II-1 can fit the UHECR spectrummeasured by Auger reasonably well, except the spectrumbecomes softer in the ”ankle” region ( ∼ . eV). InFig. 13, we show the results from model II-2. In thismodel, UHECRs are mostly heavier nuclei, Ne and Mg,which is expected when WDs have higher mass such asONeMg WDs. We find that the final spectrum can befitted very well in this scenario. The data of h X max i and σ ( X max ) in both models are consistent with Augerresults. Our results are consistent with their UHECRcomposition that becomes heavier with increasing energy,and where intermediate mass particles dominate in thehigh-energy range ( E > eV).For the tidal disruption of lower-mass WDs (heliumWDs or CO-WDs) by massive BHs, we also consider theenhancement of heavier nuclei through explosive nuclearreactions. Our results are shown in Fig. 14 for model III-1 and Fig. 15 for model III-2. UHECRs in model III-1are predominantly He, Si and Fe. This model can fit theUHECR spectrum reasonably in the lower energy range,but the cutoff energy ( ∼ eV) is more consistent withTA results and higher than Auger results. Comparedto model III-1, model III-2 has a large fraction of irongroup nuclei in UHECRs which predict an even highercutoff energy. We find that model III-1 is consistent with h X max i and σ ( X max ) measured by Auger, while modelIII-2 has a large fraction of heavier nuclei compared tothe observed composition data in the higher-energy range( > . eV). IV. DISCUSSION AND SUMMARY
This work consists of two parts. In the first part, weexamined the production of UHECRs in TDEs accompa-nied by relativistic jets. In the internal shock model, wefound that CRs can be accelerated to ultrahigh energies.However, it is difficult for UHECR nuclei to survive in log(E[eV]) E Φ ( E ) [ e V m − s r − s − ] TA (2015, energy scale - 13%)Auger(ICRC 2015)ExGalZ = 1Z = 2 3 <= Z <= 89 <= Z <= 1314 <= Z <= 1920 <= Z <= 25Z = 26 < X m a x > [ g c m − ] protoniron Epos-LHCQGSJet II-04Sibyll 2.1Auger(Aab et al, 2014) log(E[eV]) σ ( X m a x ) [ g c m − ] protoniron FIG. 11. Model I: MS stars with the solar composition. Weuse a maximum proton energy of E p, max = 2 × eV andspectral index of s esc = 1. luminous TDE jets such as Swift J1644+57, while thesurvival is allowed for less powerful TDE jets. In thereverse and forward shock model, we can expect the pro-duction and survival of UHECRs. We also consideredthe wind model, where a nonrelativistic outflow interactswith the CNM. We found that in this scenario CRs canbe accelerated only up to ∼ Z eV.In the second part, we examined different composi-tion models for TDEs. Motivated by the compositiondata by Auger, we assumed that the injected UHECR1 log(E[eV]) E Φ ( E ) [ e V m − s r − s − ] TA (2015, energy scale - 13%)Auger(ICRC 2015)ExGalZ = 1Z = 2 3 <= Z <= 89 <= Z <= 1314 <= Z <= 1920 <= Z <= 25Z = 26 < X m a x > [ g c m − ] protoniron Epos-LHCQGSJet II-04Sibyll 2.1Auger(Aab et al, 2014) log(E[eV]) σ ( X m a x ) [ g c m − ] protoniron FIG. 12. Model II-1: CO-WDs with an initial mass composi-tion, X C = 0 . X O = 0 .
5. We use a maximum proton energyof E p, max = 6 . × eV and spectral index of s esc = 1. nuclei mainly originate from tidally disrupted stars. Al-though we discussed several possibilities, this should bejustified by more detailed work on CR acceleration andescape processes. We consider both MS-SMBH and WD-IMBH tidal disruptions. In the MS-SMBH TDE scenario(model I), the injected UHECRs have a solarlike com-position. The UHECR spectrum can be fitted, but aproton dominated composition is expected in the nearlythe entire energy range. In the WD-IMBH TDE sce-nario, model II-1 (CO-WDs) can give a poor fit to the log(E[eV]) E Φ ( E ) [ e V m − s r − s − ] TA (2015, energy scale - 13%)Auger(ICRC 2015)ExGalZ = 1Z = 2 3 <= Z <= 89 <= Z <= 1314 <= Z <= 1920 <= Z <= 25Z = 26 < X m a x > [ g c m − ] protoniron Epos-LHCQGSJet II-04Sibyll 2.1Auger(Aab et al, 2014) log(E[eV]) σ ( X m a x ) [ g c m − ] protoniron FIG. 13. Model II-2: ONeMg-WDs with an initial masscomposition X O = 0 . X Ne = 0 . X Mg = 0 .
12. We usea maximum proton energy of E p, max = 6 . × eV andspectral index of s esc = 1. UHECR spectrum but h X max i and σ ( X max ) can be rea-sonably accounted for. We found that it is difficult tofit the spectrum and composition simultaneously even ifwe try a variety of parameter sets, such as higher maxi-mum energy and/or steeper ejection spectra. The mainreason is that the attenuation lengths of UHECR carbonor oxygen nuclei are lower than protons or iron nuclei,so most of them will be depleted into secondary protonsbefore reaching Earth [3]. For ONeMg-WDs, we found2 log(E[eV]) E Φ ( E ) [ e V m − s r − s − ] TA (2015, energy scale - 13%)Auger(ICRC 2015)ExGalZ = 1Z = 2 3 <= Z <= 89 <= Z <= 1314 <= Z <= 1920 <= Z <= 25Z = 26 < X m a x > [ g c m − ] protoniron Epos-LHCQGSJet II-04Sibyll 2.1Auger(Aab et al, 2014) log(E[eV]) σ ( X m a x ) [ g c m − ] protoniron FIG. 14. Model III-1: 0 . M ⊙ (ignited) helium WDs dis-rupted by 10 M ⊙ IMBHs. We use a maximum proton energyof E p, max = 6 . × eV and spectral index of s esc = 1. that the results are more consistent with Auger data.However, the rate density of ONeMg-WDs is expected tobe lower than that of CO-WDs ( ∼ / ∼
100 Mpc) [3].We also considered special cases, where nucleosynthe-sis is triggered by an IMBH. Our results show that model log(E[eV]) E Φ ( E ) [ e V m − s r − s − ] TA (2015, energy scale - 13%)Auger(ICRC 2015)ExGalZ = 1Z = 2 3 <= Z <= 89 <= Z <= 1314 <= Z <= 1920 <= Z <= 25Z = 26 < X m a x > [ g c m − ] protoniron Epos-LHCQGSJet II-04Sibyll 2.1Auger(Aab et al, 2014) log(E[eV]) σ ( X m a x ) [ g c m − ] protoniron FIG. 15. Model III-2: 1 . M ⊙ (ignited) CO-WDs disruptedby 500 M ⊙ IMBHs. We use a maximum proton energy of E p, max = 6 . × eV and spectral index of s esc = 1. III-1 can marginally fit the observed spectrum, and con-sistent with the composition data, while model III-2 hastoo many iron group nuclei, making it difficult to recon-cile with the Auger data. We also caution that the finalcomposition is sensitive to details of the parameters, suchas WD and BH masses [98, 101, 115]. Reference [102]performed 3-D smoothed particle hydrodynamics simu-lations to study the explosive nuclear burning of WD-IMBHs, considering both the adiabatic compression andshock wave generation. They found that the final mass3fractions are sensitive to the number of simulation parti-cles, and ignitions occur for low-resolution simulations.Secondary gamma-ray and neutrino signals are of in-terest to test the model. TDEs have been discussedas high-energy neutrino sources [57, 116–121]. Neutrinoproduction is expected to be efficient for luminous TDEssuch as Swift J1644+57, while the neutrino productionefficiency should be lower for low-luminous TDEs, allow-ing nucleus survival [21, 58]. We also estimate the cos-mogenic neutrino flux according to model II-2 (ONeMg-WDs), and we find E ν Φ ν ∼ − GeV cm − sr − s − .However, this would be too conservative since luminousTDEs also contribute to the neutrino flux. High-energygamma rays can be produced via neutron pion decay,photodeexcitation, and Bethe-Heitler processes. The sur-vival of UHECR nuclei implies that gamma rays can es-cape from the sources without efficient two-photon anni-hilation inside the sources [21, 58].Perhaps, one could expect the emission of gravitationalwaves from WD-IMBH tidal disruptions [98, 99, 122],and TDEs can be interesting targets for multimessengerastronomy. ACKNOWLEDGMENTS
We thank Shigeo S. Kimura and Kumiko Kotera forhelpful discussions. The work of K.M. is supported byAlfred P. Sloan Foundation and the U.S. National Sci-ence Foundation (NSF) under grants NSF Grant No.PHY-1620777. It is partially supported by the NationalNatural Science Foundation of China under Grant No.11273005 and the National Basic Research Program (973Program) of China under Grant No. 2014CB845800(Z.L.). B.T.Z. is supported by the China ScholarshipCouncil (CSC) to conduct research at Penn State Uni-versity.
Appendix A: Inelasticity of photodisintegrationprocess
The giant dipole resonance (GDR) is the most im-portant disintegration process at low energies, and thethreshold energy is ∼ ∼
150 MeV, baryonic reso-nances (BR) play a dominant role in the photodisintegra-tion process. At very high energies ∼ A ≥ . −
200 MeVcan be derived using the nuclear reaction code TALYS[56]. At higher energies, for simplicity, we use the relationthat the cross section is proportional to nuclear mass,which can be written as σ Aγ (¯ ε ) = Aσ pγ (¯ ε ) (see Fig. 1).The photomeson production cross section for protons isderived using the numerical code Geant4 [21, 55].In the photodisintegration process, a parent nucleusloses energy through the ejection of one or more nucleons.We assume that the Lorentz factor is constant during theinteraction; then the relative energy loss can be estimatedas [54] κ Aγ (¯ ε ) ≡ ∆ EE = ∆ NN , (A1)where N is the total number of nucleons in the parentnuclei, and ∆ N is the number of ejected nucleons in eachchannel. To derive the average energy loss, we considerthe contribution from all the dominant channels,¯ κ Aγ = 1 σ tot X i σ i [∆ N ] i N , (A2)where σ tot is the total cross section, σ i is the cross sectionfor the ith channel, and [∆ N ] i is the number of ejectednucleons in the ith channel. In the energy range beyondthe photopion production threshold, for simplicity, theinelasticity is assumed to linearly decrease with increas-ing nuclear mass ¯ κ Aγ (¯ ε ) = ¯ κ pγ (¯ ε ) /A . Appendix B: UHECRs propagation
The observed CR flux for each component is describedby the following formula [48, 125]Φ A = X A ′ Φ AA ′ = X A ′ c π dn AA ′ dE , (B1)where Φ AA ′ and d n AA ′ ( E ) are the observed cosmic rayflux and number density of particles with mass A gener-ated from parent particles with mass A ′ . d n AA ′ ( E ) canbe estimated by considering the contribution as a func-tion of redshift dn AA ′ ( E ) = Z z max z min dz (cid:12)(cid:12)(cid:12)(cid:12) dtdz (cid:12)(cid:12)(cid:12)(cid:12) f TDE ( z ) ρ × Z E ′ max E ′ min dE ′ dN A ′ dE ′ η AA ′ ( E, E ′ , z ) , (B2)where dtdz = − H (1 + z ) 1 p Ω Λ + Ω k (1 + z ) + Ω m (1 + z ) , (B3)4and we use the redshift evolution of TDEs given by [126] f TDE ( z ) = " (1 + z ) . η + (cid:18) z . (cid:19) − . η + (cid:18) z . (cid:19) − . η η , (B4)with η ∼ −
2. In this work, we adopt z min = 0 . z max = 2. Note that the redshift evolution for TDEsis negative. ρ is the local event rate of jetted TDEs[126]. Also, η AA ′ ( E, E ′ , z ) is the fraction of generatedCRs with mass A and energy E from parent particleswith mass A ′ and energy E ′ . d N A ′ / d E ′ is the UHECRinjection spectra per TDE.The CR luminosity density can be estimated as Q CR = X A ′ Z ZE ′ p, max E ′ min dE ′ E ′ dN A ′ dE ′ ρ = X A ′ f A ′ N ρ E Z (cid:18) E ′ p, max E (cid:19) − s esc × Γ(2 − s esc , E ′ min ZE ′ p, max ) , (B5)where Γ(2 − s esc , E ′ min ZE ′ p, max ) is the incomplete gamma func-tion. We use the following formula to estimate the nor-malization parameter N : E isoCR = X A ′ Z ZE ′ p, max E ′ min dE ′ E ′ dN A ′ dE ′ . (B6)The CR energy per TDE ( E CR ) can be estimated throughthe relation E isoCR = ξ CR E isorad , where ξ CR is the cosmic rayloading factor. Appendix C: Energy budget
The event rate of WD-IMBH tidal disruptions is un-certain. IMBHs are believed to exist in dwarf galaxiesor globular clusters. In dwarf galaxies, the event rateis estimated to be R TDE − DG ∼ f DGIMBH
Gpc − yr − ,where f DGMBH is the occupation fraction in dwarf galaxies[44]. In globular clusters, the event rate is R TDE − GC ∼ f GCMBH
Gpc − yr − , which is slightly higher than theevent rate estimated for dwarf galaxies [42].We assume only a fraction of f jet ∼
10% WD-IMBHtidal disruptions have relativistic jets with the beamingfactor f b ≈ θ j / ∼ / (2Γ j ), with a typical Lorenz factorof Γ j = 10. The apparent rate density of WD-IMBHtidal disruptions is ρ WD − IMBH ∼ . (cid:18) f IMBH (cid:19) (cid:18) f jet . (cid:19) (cid:18) f b . (cid:19) Gpc − yr − , (C1)with the assumption f IMBH ≃ f DGIMBH ≃ f GCIMBH ≃
1. Theevent rate can be comparable to the event rate ( ρ TDE ≃ .
03 Gpc − yr − ) obtained from the detection of twojetted TDEs Swift J1655+57 and Swift J2058+05 [126].According to our simulation results, we need a total CRenergy injection rate of Q CR ≈ . × erg Mpc − yr − in model II-2, and the required CR loading factor is ξ CR ∼ Q CR , . ρ − − . E − , . Note that the lumi-nosity density itself is independent on the beaming fac-tor as in the argument for GRBs. The absolute radiationenergy is smaller than the isotropic-equivalent radiationenergy by f − b
1, whereas the true rate density is a factorof f − b [1] A. M. Hillas, Ann. Rev. Astron. Astrophys. , 425(1984).[2] M. Nagano and A. A. Watson, Rev. Mod. Phys. , 689(2000).[3] K. Kotera and A. V. Olinto, Ann. Rev. Astron. Astro-phys. , 119 (2011), 1101.4256.[4] J. Abraham et al. (Pierre Auger Collaboration), Phys.Rev. Lett. , 061101 (2008), 0806.4302.[5] R. U. Abbasi et al. (HiRes Collaboration), Phys. Rev.Lett. , 101101 (2008), astro-ph/0703099.[6] J. Abraham et al. (Pierre Auger Collaboration), Phys.Lett. B685 , 239 (2010), 1002.1975.[7] K. Greisen, Phys. Rev. Lett. , 748 (1966).[8] G. T. Zatsepin and V. A. Kuzmin, JETP Lett. , 78(1966), [Pisma Zh. Eksp. Teor. Fiz.4,114(1966)].[9] R. Abbasi et al. (Pierre Auger and Telescope ArrayCollaborations), JPS Conf. Proc. , 010016 (2016),1503.07540.[10] A. Aab et al. (Pierre Auger Collaboration), Phys. Rev. D90 , 122005 (2014), 1409.4809.[11] A. Aab et al. (Pierre Auger Collaboration), Phys. Rev.
D90 , 122006 (2014), 1409.5083. [12] C. A. Norman, D. B. Melrose, and A. Achterberg, As-trophys. J. , 60 (1995).[13] A. Pe’er, K. Murase, and P. Meszaros, Phys. Rev.
D80 ,123018 (2009), 0911.1776.[14] C. D. Dermer and S. Razzaque, Astrophys. J. , 1366(2010), 1004.4249.[15] K. Murase, C. D. Dermer, H. Takami, and G. Migliori,Astrophys. J. , 63 (2012), 1107.5576.[16] K. Fang and K. Murase (2017), 1704.00015.[17] H. Kang, D. Ryu, and T. W. Jones, Astrophys. J. ,422 (1996), astro-ph/9507113.[18] S. Inoue, G. Sigl, F. Miniati, and E. Armengaud (2007),astro-ph/0701167.[19] S. S. Kimura, K. Murase, and B. T. Zhang (2017),1705.05027.[20] D. Caprioli, Astrophys. J. , L38 (2015), 1505.06739.[21] K. Murase, K. Ioka, S. Nagataki, and T. Nakamura,Phys. Rev.
D78 , 023005 (2008), 0801.2861.[22] X.-Y. Wang, S. Razzaque, and P. Meszaros, Astrophys.J. , 432 (2008), 0711.2065.[23] S. Horiuchi, K. Murase, K. Ioka, and P. Meszaros, As-trophys. J. , 69 (2012), 1203.0296. [24] N. Globus, D. Allard, R. Mochkovitch, and E. Parizot,Mon. Not. Roy. Astron. Soc. , 751 (2015), 1409.1271.[25] D. Biehl, D. Boncioli, A. Fedynitch, and W. Winter(2017), 1705.08909.[26] E. Waxman, Phys. Rev. Lett. , 386 (1995), astro-ph/9505082.[27] M. Vietri, Astrophys. J. , 883 (1995), astro-ph/9506081.[28] K. Murase, K. Ioka, S. Nagataki, and T. Nakamura,Astrophys. J. , L5 (2006), astro-ph/0607104.[29] X.-Y. Wang, S. Razzaque, P. Meszaros, and Z.-G. Dai,Phys. Rev. D76 , 083009 (2007), 0705.0027.[30] S. Chakraborty, A. Ray, A. M. Soderberg, A. Loeb, andP. Chandra, Nature Commun. , 175 (2011), 1012.0850.[31] R.-Y. Liu and X.-Y. Wang, Astrophys. J. , 40(2012), 1111.6256.[32] K. Fang, K. Kotera, and A. V. Olinto, Astrophys. J. , 118 (2012), 1201.5197.[33] K. Fang, K. Kotera, K. Murase, and A. V.Olinto, Phys. Rev. D90 , 103005 (2014), [Phys.Rev.D90,103005(2014)], 1311.2044.[34] K. Kotera, E. Amato, and P. Blasi, JCAP , 026(2015), 1503.07907.[35] M. J. Rees, Nature , 523 (1988).[36] R.-F. Shen and C. D. Matzner, Astrophys. J. , 87(2014), 1305.5570.[37] T. Piran, G. Svirski, J. Krolik, R. M. Cheng, and H. Sh-iokawa, Astrophys. J. , 164 (2015), 1502.05792.[38] D. Giannios and B. D. Metzger, Mon. Not. Roy. Astron.Soc. , 2102 (2011), 1102.1429.[39] W.-H. Lei and B. Zhang, Astrophys. J. , L27 (2011),1108.3115.[40] J. H. Krolik and T. Piran, Astrophys. J. , 134(2011), 1106.0923.[41] F. De Colle, J. Guillochon, J. Naiman, and E. Ramirez-Ruiz, Astrophys. J. , 103 (2012), 1205.1507.[42] R. V. Shcherbakov, A. Pe’er, C. S. Reynolds, R. Haas,T. Bode, and P. Laguna, Astrophys. J. , 85 (2013),1212.4837.[43] J. C. McKinney, A. Tchekhovskoy, A. Sadowski, andR. Narayan, Mon. Not. Roy. Astron. Soc. , 3177(2014), 1312.6127.[44] M. MacLeod, J. Goldstein, E. Ramirez-Ruiz, J. Guil-lochon, and J. Samsing, Astrophys. J. , 9 (2014),1405.1426.[45] G. R. Farrar and A. Gruzinov, Astrophys. J. , 329(2009), 0802.1074.[46] G. R. Farrar and T. Piran (2014), 1411.0704.[47] D. N. Pfeffer, E. D. Kovetz, and M. Kamionkowski,Mon. Not. Roy. Astron. Soc. , 2922 (2017),1512.04959.[48] R. Alves Batista, A. Dundovic, M. Erdmann, K.-H.Kampert, D. Kuempel, G. Mller, G. Sigl, A. van Vliet,D. Walz, and T. Winchen, JCAP , 038 (2016),1603.07142.[49] K. A. Olive et al. (Particle Data Group), Chin. Phys.
C38 , 090001 (2014).[50] J. S. Bloom et al., Science , 203 (2011), 1104.3257.[51] D. N. Burrows et al., Nature , 421 (2011), 1104.4787.[52] B. A. Zauderer et al., Nature , 425 (2011),1106.3568.[53] M. Lemoine and E. Waxman, J. Cosmol. Astropart.Phys. , 009 (2009).[54] J. Rachen, Ph.D. thesis, Universitt Bonn. (1996). [55] S. Agostinelli et al. (GEANT4), Nucl. Instrum. Meth.
A506 , 250 (2003).[56] A. Koning, S. Hilaire, and M. Duijvestijn,J.Nucl.Sci.Tech. p. 211 (2007).[57] N. Senno, K. Murase, and P. Meszaros, Astrophys. J. , 3 (2017), 1612.00918.[58] K. Murase and J. F. Beacom, Phys. Rev.
D82 , 043008(2010), 1002.3980.[59] E. Berger, A. Zauderer, G. G. Pooley, A. M. Soderberg,R. Sari, A. Brunthaler, and M. F. Bietenholz, Astro-phys. J. , 36 (2012), 1112.1697.[60] B. A. Zauderer, E. Berger, R. Margutti, G. G. Pooley,R. Sari, A. M. Soderberg, A. Brunthaler, and M. F.Bietenholz, Astrophys. J. , 152 (2013), 1212.1173.[61] B. D. Metzger, D. Giannios, and P. Mimica, Mon. Not.Roy. Astron. Soc. , 3528 (2012), 1110.1111.[62] A. Panaitescu and P. Kumar, Mon. Not. Roy. Astron.Soc. , 511 (2004), astro-ph/0406027.[63] K. Murase, Phys. Rev.
D76 , 123001 (2007), 0707.1140.[64] M. Lemoine and B. Revenu, Mon. Not. Roy. Astron.Soc. , 635 (2006), astro-ph/0510522.[65] M. Lemoine, G. Pelletier, and B. Revenu, Astrophys. J. , L129 (2006), astro-ph/0606005.[66] C. D. Dermer and M. Humi, Astrophys. J. , 479(2001), astro-ph/0012272.[67] Y. A. Gallant and A. Achterberg, Mon. Not. Roy. As-tron. Soc. , 6 (1999), astro-ph/9812316.[68] L. Sironi, U. Keshet, and M. Lemoine, Space Sci. Rev. , 519 (2015), 1506.02034.[69] K. D. Alexander, E. Berger, J. Guillochon, B. A. Za-uderer, and P. K. G. Williams, Astrophys. J. , L25(2016), 1510.01226.[70] K. D. Alexander, M. H. Wieringa, E. Berger, R. D. Sax-ton, and S. Komossa (2016), 1610.03861.[71] Note1, this can be a reasonable approximation if the lu-minosity function has a steep index, i.e., the luminositydensity is dominated by lower-luminosity TDEs.[72] J. G. Kirk, A. W. Guthmann, Y. A. Gallant, andA. Achterberg, Astrophys. J. , 235 (2000), astro-ph/0005222.[73] A. Achterberg, Y. A. Gallant, J. G. Kirk, and A. W.Guthmann, Mon. Not. Roy. Astron. Soc. , 393(2001), astro-ph/0107530.[74] U. Keshet and E. Waxman, Phys. Rev. Lett. , 111102(2005), astro-ph/0408489.[75] T. N. Kato and F. Takahara, Mon. Not. Roy. Astron.Soc. , 642 (2001), astro-ph/0012514.[76] J. Aoi, K. Murase, and S. Nagataki, Mon. Not. Roy.Astron. Soc. , 1431 (2008), 0711.2772.[77] M. G. Baring, Adv. Space Res. , 1427 (2011),1002.3848.[78] M. G. Baring, AIP Conf. Proc. , 294 (2009),0901.2535.[79] Y. Ohira, K. Murase, and R. Yamazaki, Astron. Astro-phys. , A17 (2010), 0910.3449.[80] D. Caprioli, P. Blasi, and E. Amato, Mon. Not. Roy.Astron. Soc. , 2065 (2009), 0807.4259.[81] B. Katz, P. Meszaros, and E. Waxman, JCAP ,012 (2010), 1001.0134.[82] P. Baerwald, M. Bustamante, and W. Winter, Astro-phys. J. , 186 (2013), 1301.6163.[83] A. Levinson, Astrophys. J. , L213 (2009), 0906.4483.[84] A. Generozov, N. C. Stone, and B. D. Metzger, Mon.Not. Roy. Astron. Soc. , 775 (2015), 1505.00268. [85] A. Generozov, P. Mimica, B. D. Metzger, N. C. Stone,D. Giannios, and M. A. Aloy, Mon. Not. Roy. Astron.Soc. , 2481 (2017), 1605.08437.[86] K. Lodders (2010), pp. 379–417, 1010.2746.[87] J.-P. Luminet and B. Pichon, Astron. Astrophys ,103 (1989).[88] F. D’Antona and I. Mazzitelli, Annu. Rev. Astron. andAstrophys. , 139 (1990).[89] V. Weidemann, Astron. Astrophys. , 647 (2000).[90] C. J. Horowitz, A. S. Schneider, and D. K. Berry, Phys.Rev. Lett. , 231101 (2010).[91] C. E. Fields, R. Farmer, I. Petermann, C. Iliadis, andF. X. Timmes, Astrophys. J. , 46 (2016).[92] K. Nomoto, Astrophys. J. , 791 (1984).[93] I. Iben, Jr. and A. V. Tutukov, Astrophys. J. Suppl. Ser. , 661 (1985).[94] J. Isern, R. Canal, and J. Labay, Astrophys. J. ,L83 (1991).[95] J. R. Wilson and G. J. Mathews, Astrophys. J. , 638(2004), [Astrophys. J.610,368(2004)], astro-ph/0307337.[96] D. S. P. Dearborn, J. R. Wilson, and G. J. Mathews,Astrophys. J. , 309 (2005).[97] S. Rosswog, E. Ramirez-Ruiz, and W. R. Hix, Astro-phys. J. , 1385 (2008), 0712.2513.[98] S. Rosswog, E. Ramirez-Ruiz, and R. Hix, Astrophys.J. , 404 (2009), 0808.2143.[99] R. Haas, R. V. Shcherbakov, T. Bode, and P. Laguna,Astrophys. J. , 117 (2012), 1201.4389.[100] P. H. Sell, T. J. Maccarone, R. Kotak, C. Knigge, andD. J. Sand, Mon. Not. Roy. Astron. Soc. , 4198(2015), 1504.05584.[101] M. MacLeod, J. Guillochon, E. Ramirez-Ruiz,D. Kasen, and S. Rosswog, Astrophys. J. , 3 (2016),1508.02399.[102] A. Tanikawa, Y. Sato, K. Nomoto, K. Maeda,N. Nakasato, and I. Hachisu (2017), 1703.08278.[103] R. Alves Batista and J. Silk (2017), 1702.06978.[104] E. Armengaud, G. Sigl, T. Beau, and F. Miniati, As-tropart. Phys. , 463 (2007), astro-ph/0603675.[105] D. Allard, Astropart. Phys. , 33 (2012),1111.3290.[106] A. Aab et al. (Pierre Auger Collaboration), JCAP ,038 (2017), 1612.07155.[107] R. C. Gilmore, R. S. Somerville, J. R. Primack, andA. Dom´ınguez, Mon. Not. R. Astron. Soc. , 3189(2012). [108] The Pierre Auger Observatory: Contributions to the34th International Cosmic Ray Conference (ICRC2015) (2015), 1509.03732.[109] M. Fukushima (Telescope Array Collaboration), EPJWeb Conf. , 04004 (2015), 1503.06961.[110] M. De Domenico, M. Settimo, S. Riggi, and E. Bertin,JCAP , 050 (2013), 1305.2331.[111] K. Werner, I. Karpenko, and T. Pierog, Phys. Rev. Lett. , 122004 (2011), 1011.0375.[112] M. Livio and J. W. Truran, Astrophys. J. , 797(1994).[113] Gil-Pons, P., Garca-Berro, E., , J. Jo, Hernanz, M., andTruran, J. W., Astron. Astrophys , 1021 (2003).[114] P. A. Denissenkov, J. W. Truran, M. Pignatari, R. Trap-pitsch, C. Ritter, F. Herwig, U. Battino, K. Setood-ehnia, and B. Paxton, Mon. Not. R. Astron. Soc. ,2058 (2014).[115] K. Kawana, A. Tanikawa, and N. Yoshida (2017),1705.05526.[116] K. Murase, in American Institute of Physics ConferenceSeries , edited by Y.-F. Huang, Z.-G. Dai, and B. Zhang(2008), vol. 1065 of
American Institute of Physics Con-ference Series , pp. 201–206.[117] K. Murase and H. Takami, in
Proceedings, 31th Inter-national Cosmic Ray Conference (ICRC 2009): Lodz,Poland, 2009 (2009), p. 1181.[118] X.-Y. Wang, R.-Y. Liu, Z.-G. Dai, and K. S. Cheng,Phys. Rev.
D84 , 081301 (2011), 1106.2426.[119] X.-Y. Wang and R.-Y. Liu, Phys. Rev.
D93 , 083005(2016), 1512.08596.[120] L. Dai and K. Fang (2016), 1612.00011.[121] C. Lunardini and W. Winter (2016), 1612.03160.[122] S. Kobayashi, P. Laguna, E. S. Phinney, andP. Meszaros, Astrophys. J. , 855 (2004), astro-ph/0404173.[123] E. Khan, S. Goriely, D. Allard, E. Parizot, T. Suomi-jarvi, A. J. Koning, S. Hilaire, and M. C. Duijvestijn,Astropart. Phys. , 191 (2005), astro-ph/0412109.[124] D. Allard, E. Parizot, E. Khan, S. Goriely, and A. V.Olinto, Astron. Astrophys. , L29 (2005), astro-ph/0505566.[125] M. Unger, G. R. Farrar, and L. A. Anchordoqui, Phys.Rev. D92 , 123001 (2015), 1505.02153.[126] H. Sun, B. Zhang, and Z. Li, Astrophys. J.812