# Ultra-High Energy Cosmic Rays and Neutrinos from Tidal Disruptions by Massive Black Holes

Claire Guépin, Kumiko Kotera, Enrico Barausse, Ke Fang, Kohta Murase

► To cite this version:

Claire Guépin, Kumiko Kotera, Enrico Barausse, Ke Fang, Kohta Murase. Ultra-High Energy Cosmic Rays and Neutrinos from Tidal Disruptions by Massive Black Holes. Astronomy & Astrophysics - A&A, 2018, 616, pp.A179. 10.1051/0004-6361/201732392. hal-01886575

HAL Id: hal-01886575

https://hal.science/hal-01886575v1

Submitted on 23 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Ultra-high-energy cosmic rays and neutrinos from tidal disruptions by massive black holes

Claire Guépin<fnref n="1" />, Kumiko Kotera<fnref n="1" />,<fnref n="2" />, Enrico Barausse<fnref n="1" />, Ke Fang<fnref n="3" />,<fnref n="4" />, and Kohta Murase<fnref n="5" />,<fnref n="6" />,<fnref n="7" />,<fnref n="8" />

Received 30 November 2017 / Accepted 7 May 2018

## ABSTRACT

Tidal disruptions are extremely powerful phenomena that have been designated as candidate sources of ultra-high-energy cosmic rays. The disruption of a star by a black hole can naturally provide protons and heavier nuclei, which can be injected and accelerated to ultra-high energies within a jet. Inside the jet, accelerated nuclei are likely to interact with a dense photon field, leading to a significant production of neutrinos and secondary particles. We model numerically the propagation and interactions of high-energy nuclei in jetted tidal disruption events in order to evaluate consistently their signatures in cosmic rays and neutrinos. We propose a simple model of the light curve of tidal disruption events, consisting of two stages: a high state with bright luminosity and short duration and a medium state, less bright and longer lasting. These two states have different impacts on the production of cosmic rays and neutrinos. In order to calculate the diffuse fluxes of cosmic rays and neutrinos, we model the luminosity function and redshift evolution of jetted tidal disruption events. We find that we can fit the latest ultra-high-energy cosmic-ray spectrum and composition results of the Auger experiment for a range of reasonable parameters. The diffuse neutrino flux associated with this scenario is found to be subdominant, but nearby events can be detected by IceCube or next-generation detectors such as IceCube-Gen2.

**Key words.** astroparticle physics – neutrinos

## 1. Introduction

The detection of ultra-high-energy cosmic rays (UHECRs) with energies > 10²⁰ eV implies the existence of extremely powerful astrophysical accelerators that have not yet been identified. Several UHECR source models have been proposed in the literature, such as radio-loud active galactic nuclei (AGN), cluster accretion shocks for steady objects, gamma-ray bursts (GRBs), fast-rotating neutron stars, or giant AGN flares for transient candidates (see, e.g., Kotera & Olinto 2011 and references therein). Most of these models can successfully fit the observational data of the Auger and Telescope Array experiments for specific choices of astrophysical parameters, and predict associated high-energy neutrino fluxes that could be observed in the next decade by existing and future experiments (see, e.g., Guépin & Kotera 2017 for a systematic study of neutrino signals from transient sources). With the current set of data, however, there is no evidence that allows us to strongly favor one particular scenario over the others.

Many other types of transient powerful events have been discovered lately thanks to unprecedented instrumental performance in terms of time resolution and sensitivity. Among them, tidal disruption events (TDE), and more specifically jetted TDEs observed by the Swift detector (e.g., Cummings et al. 2011; Cenko et al. 2012), appear to be interesting candidate sources of

UHECRs, with their impressive energy reservoirs and estimated occurrence rates.

Tidal disruption events can occur when stars approach massive black holes located at the center of galaxies at distances smaller than the tidal disruption radius. If this radius is larger than the Schwarzschild radius of the black hole, tidal forces can violently disrupt the star and produce luminous and long-lasting flares. After the disruption of the stellar object, which might be a main sequence star or in some extreme cases a white dwarf, part of its material escapes and part is accreted, launching simultaneously a wind or a relativistic outflow.

These transient events were predicted theoretically about 20 yr before their first detections, and TDEs lasting for months (or sometimes years) have been observed in the UV, X-rays and γ rays (e.g., Komossa 2015). The emission mostly shows a fast rising phase and a luminosity decay $L \propto t^{-5/3}$, coherent with fall-back accretion (Phinney 1989). The most luminous events show a higher variability, with sequences of flares of ~1000 s alternating with quiescent periods of ~5 × 10⁴ s. As they can reach luminosities of $L_{max} = 10^{48}$ erg s⁻<fnref n="1" />, and can maintain very high bolometric luminosities ($L_{bol} \sim 10^{47}$ erg s⁻<fnref n="1" />) lasting about 10⁶ s, these powerful emissions are very likely to come from a relativistic jet launched from the central massive black hole (e.g., Bloom et al. 2011; Burrows et al. 2011). To date, it is still not

<footnote n="1">¹ Sorbonne Universités, UPMC Univ. Paris 6 et CNRS, UMR 7095, Institut d'Astrophysique de Paris, 98 bis bd Arago, 75014 Paris, France
e-mail: claire.guepin@iap.fr</footnote>

<footnote n="2">² Laboratoire AIM-Paris-Saclay, CEA/DSM/IRFU, CNRS, Université Paris Diderot, 91191 Gif-sur-Yvette, France</footnote>

<footnote n="3">³ Department of Astronomy, University of Maryland, College Park, MD 20742, USA</footnote>

<footnote n="4">⁴ Joint Space-Science Institute, College Park, MD 20742, USA</footnote>

<footnote n="5">⁵ Department of Physics, The Pennsylvania State University, University Park, Pennsylvania 16802, USA</footnote>

<footnote n="6">⁶ Department of Astronomy & Astrophysics, The Pennsylvania State University, University Park, Pennsylvania 16802, USA</footnote>

<footnote n="7">⁷ Center for Particle and Gravitational Astrophysics, The Pennsylvania State University, University Park, Pennsylvania 16802, USA</footnote>

<footnote n="8">⁸ Yukawa Institute for Theoretical Physics, Kyoto, Kyoto 606-8502, Japan</footnote>
clear if non-jetted and jetted TDEs constitute two distinct populations.

Jetted TDEs could be an ideal site for the production of UHECRs (via the injection of the disrupted stellar material and its acceleration in the jet) and for the production of high-energy neutrinos (produced later by the interaction of the accelerated hadrons with the ambient radiative and/or hadronic backgrounds). Although only a handful of jetted TDEs have been detected so far, these objects have already attracted great interest in the high-energy astroparticle community. High-energy neutrino production in TDE jets was considered before the discovery of IceCube neutrinos (e.g., Murase 2008; Murase & Takami 2009a; Wang et al. 2011), and contributions to IceCube neutrinos have been studied (Senno et al. 2017; Dai & Fang 2017; Lunardini & Winter 2017; Wang & Liu 2016).

The UHECR production in TDE jets was originally suggested by Farrar & Gruzinov (2009), and the external shock scenario was also considered in detail (Farrar & Piran 2014). However, it should be kept in mind that the rate of TDEs is too tight to fit the observed UHECR fluxes, as can be deduced from the constraints derived by Murase & Takami (2009b), who obtained $\dot{n}_{tde} > 1 \text{ Gpc}^{-3} \text{ yr}^{-1}$. Hence a pure proton case is disfavored and the nucleus scenario is required. Recent studies attempted to inject a mixed composition and fit the UHECR flux and composition simultaneously in both the internal and external shock scenarios (Alves Batista & Silk 2017; Zhang et al. 2017).

So far, the existing studies have not attempted to model the production of UHECRs in the inner part of the TDE jet (with acceleration occurring at internal shocks for instance). Modeling this effect requires taking into account the interaction of accelerated nuclei inside the jet in order to calculate consistently the resulting chemical composition. In this work, we study the interaction of accelerated nuclei inside the TDE jet, and the signatures they can produce in UHECRs and neutrinos. For this purpose we developed a new propagation and interaction code that is comprised of modules from CRPropa3 (Alves Batista et al. 2016) and from the code described in Kotera et al. (2009).

In order to calculate the diffuse fluxes of UHECRs and neutrinos, we also introduce a new model for the event rate evolution and luminosity function of TDEs powering jets. The semi-analytic galaxy formation model of Barausse (2012) is used to model the cosmological evolution of massive black holes, which can be related to the jetted TDE comoving event rate density, and thus to the diffuse UHECR and neutrino fluxes.

The properties of TDEs powering jets are subject to large uncertainties. Therefore, we scan the parameter space allowed by TDE observations to model the radiation region (Sect. 2.1) and the typical photon field inside a TDE jet (Sect. 2.2). Inside this region, we consider different interaction processes, detailed in Sect. 2.3. We then calculate mean free path (MFP) tables for the interaction of protons and heavier nuclei with the photon field of the jet. We use these tables in our code to predict UHECR and neutrino signatures (Sect. 3) for single sources. In order to estimate the diffuse particle fluxes from a population of jetted TDEs, we derive the luminosity function and occurrence rate evolution of these events (Sect. 4). We find in Sect. 5 that we can fit the latest UHECR spectrum and composition results of the Auger experiment for a range of reasonable parameters. The diffuse neutrino flux associated with this scenario is found to be detectable with IceCube in the next decade. Transient neutrino signals from single sources would be difficult to detect by IceCube or the upcoming GRAND experiment, except for sources located within ~20 Mpc, associated with a very low event rate.

## 2. Interaction of UHE nuclei inside TDE jets

In the following, all primed quantities are in the comoving frame of the emitting region. Other quantities are in the observer frame. Quantities are labeled $Q_x \equiv Q/10^x$ in cgs units unless specified otherwise, and except for particle energies, which are in $E_x \equiv E/10^x$ eV.

The tidal disruption of a stellar object can occur if it gets close enough to a black hole, and will produce observable flares if it happens outside the black hole event horizon. A part of the stellar material forms a thick accretion disk, and a fraction of this material accretes onto the black hole, most likely in a super-Eddington regime. For most TDEs, the observed radiation comes from the dissipation inside the accretion disk, characterized by a thermal spectrum peaking in extreme ultraviolet or soft X-rays; for a rotating black hole launching a relativistic jet, a non-thermal hard X-ray radiation can be detected, presumably produced through synchro-Compton radiation (e.g., Burrows et al. 2011). The jet radiation should dominate the observed spectrum for black holes with low mass and high spin, jets oriented toward the observer, and large radiative efficiency of the jet.

In the jet comoving frame, using a condition of causality, $R' = \Gamma c t_{var}$ can be considered as the size of the emitting region. In the internal shock model, the distance of the emission region from the black hole is estimated to be $R = \Gamma^2 c t_{var} = 3 \times 10^{14} \text{ cm} \Gamma_1^2 t_{var,2}$, where $t_{var} = 10^2 t_{var,2}$ s and $\Gamma = 10\Gamma_1$ are the typical variability timescale and bulk Lorentz factor for jetted TDEs, respectively. This radius coincides with the radius estimated from high-latitude emission with a duration of ~$t_{var}$ (e.g., Piran 2004). We note that more generally the relationship between $R$ and $R'$ can be modified, for example by subsequent internal shocks caused by merged shells and the existence of multiple emission regions such as subjets (e.g., Murase & Nagataki 2006; Bustamante et al. 2017). However, as long as we consider internal shocks in the jet that expand conically, it is reasonable to consider the expression of $R$ obtained for the one-zone calculation, as has been done in the literature of GRBs.

First, we assume that cosmic rays are injected at the center of a non-evolving radiation region in the comoving frame. The evolution of the region would mainly result in the dilution over time of the radiation and magnetic energy densities, together with adiabatic losses, associated with observable spectral changes. We account for these effects, to a first approximation, by considering two dominant stages for our TDEs: the early stage, when the source is in a high state, at its maximum brightness; and a medium state, reached later, for which the source is typically 1–1.5 orders of magnitude less luminous, but for a longer integrated time. We argue in the following that these two states have different impacts on the production of UHECRs and their associated neutrinos.

### 2.1. UHECR injection and energetics

Cosmic ray nuclei from the stellar material can be accelerated to ultra-high energies inside the TDE jet via one of the various mechanisms advocated for GRBs or AGN jets. We assume that acceleration leads to a rigidity-dependent spectrum in $dN_{CR}/dE' = A \sum_Z f_Z E'^{-\alpha} \exp(-E'/E'_{Z,\text{max}})$ with an exponential cutoff at $E'_{Z,\text{max}}$ for nuclei of charge Z. Here $A$ is a normalization constant and $f_Z$ is the fraction of elements with charge number Z, such that $\sum_Z f_Z = 1$. The spectral index $\alpha$ can vary (typically between $\alpha \sim 1$ and $\alpha \gtrsim 2$) depending on the acceleration
mechanism (e.g., magnetic reconnection or diffusive shock acceleration). The cosmic-ray composition depends on the composition of the disrupted object, but it also strongly depends on what happens to the elements before they get injected and accelerated in the jet. Heavy nuclei could indeed undergo fragmentation during the disruption of the stellar object, or a large fraction of light nuclei could escape as part of the expelled stellar envelope. In this work, the elements injected in the radiation region are protons (p), helium (He), carbon and oxygen (CO), silicium (Si), and iron (Fe).

The maximum injection energy $E'_{Z,\text{max}}$ is determined by the competition between the acceleration timescale for a nucleus of charge Z, $t'_{\text{acc}} = \eta_{\text{acc}}^{-1}E'/cZeB'^{1}$, and the energy loss timescales $t'_{\text{loss}} = \min(t'_{\text{dyn}}, t'_{\text{syn}}, t'_{\text{IC}}, t'_{\text{BH}}, t'_{\text{py}}, \dots)$, where $t'_{\text{dyn}} = R/\beta\Gamma c$ is the dynamical timescale (see Appendix A for numerical estimates of $E'_{Z,\text{max}}$). The factor $\eta_{\text{acc}} \le 1$ describes the efficiency of the acceleration process; for a maximally efficient acceleration, $\eta_{\text{acc}} = 1$. In this study we neglect the re-acceleration of secondary particles, and leave it for future work.

From the energetics point of view, the luminosity injected into cosmic rays is considered related to the bolometric luminosity in photons, such that $L_{\text{CR}} = \xi_{\text{CR}}L_{\text{bol}}$, where we define the baryon loading fraction $\xi_{\text{CR}}$ as the fraction of the bolometric luminosity that is injected into cosmic rays of energy $E \ge E_{\text{min}} \equiv 10^{16}$ eV.

## 2.2. Modeling the TDE spectral energy distribution

As suggested in Senno et al. (2017), we model the spectral energy distribution (SED) inside the TDE jet as a log-parabola with three free parameters: the peak luminosity $L_{\text{pk}}$, peak energy $\epsilon_{\text{pk}}$, and width $\hat{a}$. The photon energy density then reads

$$ \epsilon'^2 n'_{e'} = \frac{L_{\text{pk}}}{4\pi\Gamma^2 R^2 c} (\epsilon'/\epsilon'_{\text{pk}})^{-\hat{a} \log(\epsilon'/\epsilon'_{\text{pk}})}. \quad (1) $$

The peak luminosity and peak energy set the maximum of the SED. The data can help to constrain the width of the log-parabola and a potential high-energy cutoff. However, there are large uncertainties on the observed photon density, due to galaxy absorption, and even more on the photon density inside the jet (see Burrows et al. 2011; Bloom et al. 2011 for the spectrum of Swift J1644+57).

From our SED model, the bolometric luminosity can then be defined as the luminosity integrated over the entire spectrum: $L_{\text{bol}} = \int d\epsilon' L_{\text{pk}}/\epsilon' (\epsilon'/\epsilon'_{\text{pk}})^{-\hat{a} \log(\epsilon'/\epsilon'_{\text{pk}})}$. As we consider a constant photon field, this bolometric luminosity is larger than the peak luminosity. Moreover, as we model the radiation field inside the jet, we should have $L_{\text{jet,obs}} \sim L_{\text{bol}}$. We note that in most cases, the main contribution to the observed luminosity is the jet luminosity, but for high black hole masses ($M_{\text{bh}} > 5 \times 10^7 M_{\odot}$), the thermal luminosity is of the same order of magnitude as the jet luminosity (Krolik & Piran 2012).

In this work, we examine several cases summarized in Table 1 and illustrated in Fig. 1. We choose to only vary the width $\hat{a}$ and the peak luminosity $L_{\text{pk}}$ of the log-parabola, and to consider a typical peak energy $\epsilon_{\text{pk}} = 70$ keV, which is compatible with Swift J1644+57 observations (e.g., Burrows et al. 2011). Each case corresponds to a different magnetic field, and therefore corresponds to a maximum proton energy $E'_{\text{p,max}}$ (Eq. (A.4)). The magnetic field is inferred assuming equipartition between the radiative and magnetic energy densities:

<table><caption>**Table 1.** Properties of the different TDE photon fields considered in this work for a cosmic-ray acceleration efficiency $\eta_{\text{acc}} = 1$.</caption><thead><tr><th>L<sub>bol</sub> [L<sub>pk</sub>] (erg s<sup>-1</sup>)</th><th>â</th><th>B'(G)</th><th>E'<sub>p,max</sub> (eV)</th></tr></thead><tbody><tr><td>3.5 &times; 10<sup>48</sup> [10<sup>48</sup>]</td><td>0.25</td><td>5.1 &times; 10<sup>3</sup></td><td>1.8 &times; 10<sup>18</sup></td></tr><tr><td>6.8 &times; 10<sup>48</sup> [10<sup>48</sup>]</td><td>0.07</td><td>7.0 &times; 10<sup>3</sup></td><td>2.4 &times; 10<sup>18</sup></td></tr><tr><td>1.0 &times; 10<sup>49</sup> [10<sup>48</sup>]</td><td>0.03</td><td>8.7 &times; 10<sup>3</sup></td><td>2.2 &times; 10<sup>18</sup></td></tr><tr><td>6.8 &times; 10<sup>46</sup> [10<sup>46</sup>]</td><td>0.07</td><td>7.0 &times; 10<sup>2</sup></td><td>6.3 &times; 10<sup>18</sup></td></tr><tr><td>6.8 &times; 10<sup>47</sup> [10<sup>47</sup>]</td><td>0.07</td><td>2.2 &times; 10<sup>3</sup></td><td>4.3 &times; 10<sup>18</sup></td></tr><tr><td>6.8 &times; 10<sup>48</sup> [10<sup>48</sup>]</td><td>0.07</td><td>7.0 &times; 10<sup>3</sup></td><td>2.4 &times; 10<sup>18</sup></td></tr></tbody></table>

**Notes.** All the photon fields are modeled by a log-parabola (Eq. (1)), with bolometric luminosity $L_{\text{bol}}$, peak luminosity $L_{\text{pk}}$, peak energy $\epsilon_{\text{pk}}$, and width $\hat{a}$.

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Photon density for a log-parabola model with fixed peak luminosity $L_{\text{pk}} = 10^{48}$ erg s$^{-1}$.</figcaption></figure>

$\xi_B \int d\epsilon' \epsilon' n'_{\epsilon'} = B^2/8\pi$ with $\xi_B = 1$. Rough equipartition is a standard hypothesis for jets that can be argued from measurements of the energy repartition in extragalactic objects, for example blazar jets (Celotti & Ghisellini 2008). It also naturally arises if relativistic reconnection is at play in the outflow, dissipating electromagnetic energy into kinetic energy (Sironi et al. 2015).

The TDE photon spectra evolve in time. As mentioned earlier, although we do not account for proper time evolutions of the SED in this paper, we consider two states of the SED, inferred from the observations of Swift J1644+57, and which are important for our framework: an early state, corresponding to a high state that can typically last $t_{\text{dur}} \sim 10^5$ s with a bright luminosity, a high jet efficiency, and a narrow jet opening angle; and a medium state, 1–1.5 orders of magnitude less bright, but lasting $t_{\text{dur}} \sim 10^6$ s, with a lower jet efficiency and a similar jet opening angle. For both states, we set a width $\hat{a} = 0.07$. These parameters are overall compatible with Swift J1644+57 SED models corrected for galactic absorption (e.g., Burrows et al. 2011).

## 2.3. Interaction processes

All relevant interaction processes for nucleons and heavier nuclei are taken into account in our calculations. Nucleons experience pion production via photohadronic and hadronic interactions, as well as neutron decay. Nuclei undergo photonuclear processes in different regimes (requiring increasing photon energy in the nucleus rest frame): giant dipole resonance, quasi

<footnote n="1">¹ Some papers adopt $t_{\text{acc}} = \eta_{\text{acc}} r'_{L}/c$ due to the historical convention.</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Influence of the log-parabola width $\hat{a}$ on the MFPs and energy loss lengths in the comoving frame for carbon (left panel) and iron (right panel) nuclei of Lorentz factor $\gamma$. The peak luminosity is set to $L_{pk} = 10^{48}$ erg s$^{-1}$ and two examples of widths are presented: $\hat{a} = 0.25$ (blue) and $\hat{a} = 0.03$ (red). The different line styles correspond to different processes: photonuclear (solid), inverse Compton (dashed), Bethe-Heitler (dotted), and synchrotron (double dot-dashed). The black long-dashed line corresponds to the typical comoving size of the region. Wider SED lead to larger MFPs.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Influence of the log-parabola peak luminosity $L_{pk}$. Same as Fig. 2, but for fixed width $\hat{a} = 0.07$ and varying the peak luminosity: $L_{pk} = 10^{46}$ erg s$^{-1}$ (blue) and $L_{pk} = 10^{48}$ erg s$^{-1}$ (red). Higher peak luminosity leads to shorter MFPs.</figcaption></figure>

deuteron, baryon resonance, and photofragmentation. For all particles, including secondary pions and muons, we account for the synchrotron, inverse Compton, and Bethe-Heitler processes.

All the interaction cross sections and products are obtained from analytic formulae (e.g., Rachen 1996; Dermer & Menon 2009) or tabulated using numerical codes like SOPHIA (Mücke et al. 2000) for photopion production, TALYS (Koning et al. 2005) for photonuclear interactions, and EPOS (Werner et al. 2006) for purely hadronic interactions. We assume that the photofragmentation products are similar to the products of hadronic interactions, which is reasonable to a first approximation. In principle, EPOS generates

too many free nucleons in the fragmentation process. The scaling between photofragmentation and hadronic interactions is given as a function of the center of mass energy. However, the discrepancy is less than a factor of 2 compared to the data. We note in any case that we consider here an energy range where uncertainties are very large, and that the photofragmentation process for nuclei is not dominant in our study.

For extragalactic propagation, photonuclear cross sections and EBL models have a strong influence on the spectrum and composition of cosmic rays, as discussed in Alves Batista et al. (2015). For different EBL models, the discrepancy between
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Cosmic-ray spectra for one source with pure iron injection with spectral index $\alpha = 1.8$, photon field with $\hat{a} = 0.07$, and acceleration efficiency $\eta_{acc} = 0.2$. We show the total spectrum (black) and the composition (other colors), as indicated in the legend, for TDE around a black hole of mass $M_{bh} = 7 \times 10^6 M_{\odot}$, with a corresponding SED in its high state with $L_{pk} = 10^{47.5} \text{ erg s}^{-1}$ and $t_{dur} = 10^5 \text{ s}$ (left) and in its medium state with $L_{pk} = 10^{46} \text{ erg s}^{-1}$ and $t_{dur} = 10^6 \text{ s}$ (right) for a source distance $d_{L,1} = 1.88 \text{ Gpc}$. We assume here $\xi_{CR} = 1$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Neutrino spectra for three flavors for one source with same characteristics as in Fig. 4. We show the total spectra (in eV cm$^{-2}$) for a high state SED with $L_{pk} = 10^{47.5} \text{ erg s}^{-1}$ (left) and a medium state SED with $L_{pk} = 10^{46} \text{ erg s}^{-1}$ (right). We consider two different distances $d_{L,1} = 1.88 \text{ Gpc}$ (dark blue) and $d_{L,2} = 20 \text{ Mpc}$ (light blue). The IceCube and projected GRAND sensitivities (Fang et al. 2017) are also shown (dashed black and green lines). For the IceCube sensitivities, we show two cases depending on the declination: $0^{\circ} < \delta < 30^{\circ}$ (most favorable case, thin line) and $30^{\circ} < \delta < 60^{\circ}$ (Swift J1644+57 case, thick line) (Aartsen et al. 2014).</figcaption></figure>

cosmic-ray spectra can reach ~40%. The impact of photonuclear models is also strong, whereas more difficult to quantify, especially regarding the channels involving $\alpha$-particles.

We show in Fig. 2 the MFPs or energy loss lengths derived for different photon fields. In the top panel, we compare the specific luminosities and photon densities for $\hat{a} = 0.25$, $\hat{a} = 0.07$, and $\hat{a} = 0.03$ for a fixed $L_{pk} = 10^{48} \text{ erg s}^{-1}$. The width of the log-parabola has a strong influence on the MFPs as it substantially changes the radiation energy density. The MFPs for the carbon and iron cases are shown in Fig. 2 for the extreme cases $\hat{a} = 0.25$ and $\hat{a} = 0.03$. We see that overall, photonuclear interactions dominate over a wide range of particle Lorentz factors $\gamma$, up to ultra-high energies where synchrotron losses start taking over. Changing the width of the log-parabola modifies the MFPs by several orders of magnitude, with shorter paths for narrower SED.

In Fig. 3, we compare the specific luminosities and photon densities for $L_{pk} = 10^{46} \text{ erg s}^{-1}$ and $L_{pk} = 10^{48} \text{ erg s}^{-1}$ for a given width $\hat{a} = 0.07$. The influence of peak luminosity on the MFPs is

more moderate; as expected, the MFPs are a power of the peak luminosity, and a higher $L_{pk}$ leads to shorter MFP.

### 3. UHECRs and neutrinos from single TDEs

We calculate the cosmic-ray and neutrino spectra after the propagation of protons or nuclei through the photon field of a jetted TDE. The production of neutrinos should be dominated by the high state when the photon field is brightest (and the opacities greatest), and the UHECR production should be calculated over the longer medium state with a lower luminosity but over longer production timescales. We thus calculate the neutrino and UHECR fluxes at their maximum production states. We first consider one single source and show the outgoing spectra for cosmic rays and neutrinos.

We show in Fig. 4 an example of outgoing cosmic-ray spectrum for a pure iron injection from a single TDE in its high state SED characterized by $L_{pk} = 10^{47.5} \text{ erg s}^{-1}$ and $\hat{a} = 0.07$, and in its medium state SED characterized by $L_{pk} = 10^{46} \text{ erg s}^{-1}$
and $\hat{a} = 0.07$. As is shown below, these two states are associated in our model with a black hole of mass $M_{bh} = 7 \times 10^6 M_{\odot}$. We consider an injection spectral index of $\alpha = 1.8$ and an acceleration efficiency $\eta_{acc} = 0.2$. Here we do not account for the extragalactic propagation of cosmic rays, and the spectrum is normalized by considering the luminosity distance of Swift J1644+57: $d_{L,1} \approx 1.88$ Gpc ($z \approx 0.354$). Two associated neutrino spectra are shown in Fig. 5. One spectrum is normalized by considering the luminosity distance of Swift J1644+57, $d_{L,1} \approx 1.88$ Gpc, and the other by considering a luminosity distance $d_{L,2} = 20$ Mpc. The IceCube sensitivity is characterized by a minimum fluence $S_{IC} = 5 \times 10^{-4} \text{ TeV cm}^{-2}$ over the energy range $10 \text{ TeV} - 10 \text{ PeV}$, which corresponds to a detection limit $s_{IC} \sim 10^{-11} \text{ TeV cm}^{-2} \text{ s}^{-1}$ for a one-year data collection (Aartsen et al. 2015). We give the IceCube sensitivity from the effective area presented in Aartsen et al. (2014) for the optimal declination range $0^{\circ} < \delta < 30^{\circ}$ (thin lines), and for the declination range $30^{\circ} < \delta < 60^{\circ}$ (thick lines) associated with the Swift event J1644+57.

The peak luminosity and width of the photon SED have a strong effect on the cosmic-ray and neutrino spectra as they influence strongly photohadronic and synchrotron losses, which are the two dominant energy loss processes in our framework. For cosmic rays, energy losses due to photohadronic interactions are mainly dominant at low energies, while synchrotron losses dominate at high energies. If the radiation energy density is sufficiently low, the escape time of cosmic rays can be the limiting time at low energies.

Regarding the cosmic-ray spectrum, we see in Fig. 4 that for a medium state SED with $L_{pk} = 10^{46} \text{ erg s}^{-1}$, iron strongly interacts and produces many secondary particles with a large number of nucleons below $E_{cut}/56 \sim 10^{19} \text{ eV}$. Despite these high interaction rates, nuclei can still survive and escape from the region with energies up to $10^{20} \text{ eV}$. For a high state SED with $L_{pk} = 10^{47.5} \text{ erg s}^{-1}$, the iron strongly interacts, as do the secondary cosmic rays produced through iron interactions. No nuclei can survive and escape the region; only protons escape, with a maximum energy around $E_{p,max} \sim 10^{18} \text{ eV}$. The high-energy cutoff for each element with $Z > 1$ results from the competition between the energy loss processes (see Appendix A) or from the maximum injection energy; in Fig. 4, acceleration is the limiting process for $L_{pk} = 10^{46} \text{ erg s}^{-1}$, and photonuclear interactions for $L_{pk} = 10^{47.5} \text{ erg s}^{-1}$.

Figure 5 shows that a nearby medium state TDE at distance 20 Mpc with peak luminosity $L_{pk} = 10^{46} \text{ erg s}^{-1}$ would not be detectable, even with future neutrino detectors such as GRAND (Fang et al. 2017). On the other hand, at early times and in their high states, TDEs would lead to massive production of high-energy neutrinos, and should be marginally detectable with IceCube and with GRAND at the high-energy end for a nearby distance of 20 Mpc. We note that the rate of TDEs at distances smaller than 20 Mpc is between $10^{-6} \text{ yr}^{-1}$ and $3 \times 10^{-5} \text{ yr}^{-1}$, depending on the population model considered, which is extremely low. TDEs in their high states at distances >20 Mpc would not be detectable with IceCube or GRAND because of the flux decrease and the low high-energy cutoff of the neutrino spectrum. We note that our chosen parameter set is consistent with the non-detection of neutrinos from Swift J1644+57 (as was already highlighted in Guépin & Kotera 2017) and allows for baryonic loading at the source $\xi_{CR}$ up to ~ a few 100 to remain consistent with this non-detection.

The presence of a plateau in the neutrino spectrum is due to the contributions of muon and electron neutrinos. The high-energy cutoff is due to pions and muons experiencing energy

losses (mainly synchrotron losses) before they decay. We account for synchrotron and inverse Compton losses, but do not account for the kaon contribution. We note that electron neutrinos have a lower energy cutoff than muon neutrinos. Electron neutrinos are produced through muon decay, and muons are produced through pion decay; therefore, the energy of electron neutrinos is influenced by pion losses and muon losses before they decay. Muon neutrinos, in turn, are produced through pion decay or muon decay. Hence the energy of those produced through pion decay is only influenced by pion energy losses before the decay, which explains the higher energy cutoff for muon neutrinos.

## 4. Modeling the population of TDEs contributing to UHECR and neutrino fluxes

A derivation of the comoving density rate of TDEs can be found in Sun et al. (2015). These authors define the comoving density rate as $\dot{n}(z, L) = \dot{n}_0 \Lambda(L) f(z)$, where $\dot{n}_0$ is the total local event rate density, $f(z)$ the TDE redshift distribution, and $\Lambda(L)$ the TDE luminosity function. The luminosity function is given by a power-law,

$$ \Lambda_{\text{TDE}}(L_{\gamma}) \propto \left( \frac{L_{\gamma, \text{pk}}}{L_{\text{m, pk}}} \right)^{-\alpha_{\text{L}}}, \qquad (2) $$

with $L_{m, pk} = 10^{48} \text{ erg s}^{-1}$ and $\int_{L_{min}}^{L_{max}} dL_{\gamma} \Lambda_{\text{TDE}}(L_{\gamma}) = 1$, with $L_{min} = 10^{45} \text{ erg s}^{-1}$ and $L_{max} = 10^{49} \text{ erg s}^{-1}$, and $\alpha_{L} = 2$.

However, the Sun et al. (2015) model is not well adapted to our framework as their comoving rate density accounts for the entire TDE population and not the subpopulation powering jets. Moreover, the redshift evolution of the luminosity function is neglected, due to the small size of their observational sample.

Thus, in the following we present a prediction for the comoving event rate density of TDEs powering jets by combining the TDE rate per galaxy $\dot{N}_{\text{TDE}}$ and the black hole comoving number density per luminosity $d\eta_{bh}(z, L)/dL$ (i.e., the number of black holes per comoving volume and per bin of jet luminosity).

### 4.1. TDE rate per galaxy

The TDE rate per galaxy $\dot{N}_{\text{TDE}}$ depends on the galaxies considered. Following Wang & Merritt (2004), we consider a lower bound in the case of core galaxies of

$$ \dot{N}_{\text{TDE}} \approx 10^{-5} \text{ yr}^{-1}, \qquad (3) $$

and an upper bound in the case of power-law galaxies of

$$ \dot{N}_{\text{TDE}} \approx 7.1 \times 10^{-4} \text{ yr}^{-1} \left( \frac{\sigma}{70 \text{ km s}^{-1}} \right)^{7/2} M_{\text{bh,6}}^{-1}, \qquad (4) $$

where $M_{bh,6} = M_{bh}/10^6 M_{\odot}$ and $\sigma$ is the stellar velocity dispersion of the bulge. From Kormendy & Ho (2013), the relation between the black hole mass and the bulge velocity dispersion is

$$ \log_{10} M_{\text{bh,9}} = -0.51 + 4.4 \log_{10} \left( \frac{\sigma}{200 \text{ km s}^{-1}} \right). \qquad (5) $$

Using Eqs. (4) and (5), we obtain the TDE rate per galaxy in the case of power-law galaxies, which depends only on the black hole mass:

$$ \dot{N}_{\text{TDE}} \approx 3 \times 10^{-4} \text{ yr}^{-1} M_{\text{bh,6}}^{-0.2}. \qquad (6) $$
## 4.2. Identifying the black hole masses leading to observable TDEs

First, we identify the population of black holes that can lead to observable TDEs. TDEs can occur for stellar objects of mass $M_\star$ at distances $d_\star \le r_t = R_\star(M_{bh}/M_\star)^{1/3}$ (Hills 1975). Following Krolik & Piran (2012), we can estimate the tidal disruption radius $r_t \approx 10R_s M_{bh,6}^{2/3} M_{bh,6}^{-2/3} [(k_\star/f_\star)/0.02]^{1/6}$, where $R_s = 2GM_{bh}/c^2$ is the Schwarzschild radius, $M_{\star,\odot}$ is the mass of the star in solar units, $M_{bh,6} = M_{bh}/10^6 M_\odot$, $k_\star$ is related to the star's radial density profile, and $f_\star$ is its binding energy in units of $GM_\star^2/R_\star$. This radius is obtained for a main sequence star with a typical mass-radius relation $R_\star \approx R_\odot M_{\star,\odot}^{1-\xi}$ with $\xi \approx 0.2$ for $0.1 M_\odot < M_\star \le M_\odot$ or $\xi \approx 0.4$ for $M_\odot < M_\star < 10 M_\odot$ (Kippenhahn & Weigert 1994). Moreover, we consider here and in what follows fully radiative stars, thus $k_\star/f_\star = 0.02$ (Phinney 1989). For white dwarfs, typically $R_\star \sim M_\star^{1/3}$ with $0.5 M_\odot < M_\star \le 0.7 M_\odot$. Their tidal disruption radii are smaller due to the smaller dimensions of these objects; an approximate formula gives $r_t \approx 7.4 \times 10^9 (M_{bh,3.3}/\rho_{\star,7})^{1/3}$ (Luminet & Pichon 1989), where $\rho_\star$ is the white dwarf core density.

With this tidal disruption radius, we can estimate the maximum black hole mass enabling the production of flares. The first-order requirement for flares to be produced reads $r_t \ge R_s$. For a Schwarzschild black hole, this leads to $M_{bh} \le 4 \times 10^7 M_\odot M_{\star,\odot}^{1-3\xi/2} [(k_\star/f_\star)/0.02]^{1/4}$, which ranges from $M_{bh} \le 10^7 M_\odot$ to $M_{bh} \le 10^8 M_\odot$ for $0.1 M_\odot < M_\star \le 10 M_\odot$. However, jetted TDEs are likely to be powered by black holes with moderate to high spin; a general-relativistic treatment accounting for the black hole spin increases the maximum mass of black holes that are able to disrupt a solar-like star: $M_{bh} \sim 7 \times 10^8 M_\odot$ (Kesden 2012).

## 4.3. Relation between black hole mass and jet luminosity

The black hole mass function $dn_{bh}(z, M_{bh})/dM_{bh}$ (i.e., the number of black holes per comoving volume and per mass bin) is obtained with the semi-analytic galaxy formation model review in Sect. 4.4, and we derive $dn_{bh}(z, L)/dL$ by relating the black hole mass and the jet luminosity. Following Krolik & Piran (2012), we consider a TDE which forms a thick accretion disk, powering a jet through the Blandford–Znajek mechanism. We estimate the maximum accretion rate by considering that about 1/3 of the stellar mass is accreted after one orbital period $P_{orb}$ (Lodato et al. 2009). From Krolik & Piran (2012),

$$P_{\text{orb}}(a_{\text{min}}) \approx 5 \times 10^5 \, \text{s} \, M_{\star, \odot}^{(1-3\xi)/2} M_{\text{bh}, 6}^{1/2} \left(\frac{k_{\star}/f_{\star}}{0.02}\right)^{1/2} \beta^3, \quad (7)$$

where $a_{min}$ is the minimum semi-major axis. The parameter $\xi$ comes from the main sequence mass-radius relation $R_\star \approx R_\odot M_{\star,\odot}^{(1-\xi)}$ and $\beta \le 1$ is the penetration factor. We obtain the following accretion rate:

$$\dot{M} \approx 20 M_{\odot} \, \text{yr}^{-1} M_{\star, \odot}^{(1+3\xi)/2} M_{\text{bh}, 6}^{-1/2} \left(\frac{k_{\star}/f_{\star}}{0.02}\right)^{-1/2} \beta^{-3}. \quad (8)$$

The luminosity of a jet powered by a black hole depends on the regime of accretion. In the super-Eddington regime, i.e., for $M_{bh} \le M_{bh,jet}$, where $M_{bh,jet} = 4 \times 10^8 M_\odot (\dot{m}/\dot{m}_0)^{2/3} M_{\star,\odot}^{(1+\xi)/3} [(k_\star/f_\star)/0.02]^{-1/2} \beta^{-3}$, the jet luminosity is

given by Krolik & Piran (2012),

$$L_{\text{jet}} \approx 10^{43} \, \text{erg} \, s^{-1} \frac{f(a)}{\beta_h \alpha_s} M_{\text{bh},6}^{-1/2} \times 8 \times 10^3 q(\dot{m}/\dot{m}_0) \times M_{\star,\odot}^{(1+3\xi)/2} \left(\frac{k_\star/f_\star}{0.02}\right)^{-1/2} \beta^{-3}, \quad (9)$$

where $\alpha_s$ is the ratio of inflow speed to orbital speed of the disk, and $\beta_h$ the ratio of the midplane total pressure near the ISCO to the magnetic pressure in the black hole's stretched horizon, such that $\alpha_s \beta_h \sim 0.1-1$; the function $f(a) \approx a^2$ encodes the dependence of the jet luminosity on the dimensionless spin of the black hole (Piran et al. 2015), $a$, which ranges from $a=0$ (for a Schwarzschild black hole) to $a=1$ (for a maximally spinning black hole); $\dot{m} = \dot{M}c^2/L_{\text{Edd}}$ is the normalized accretion rate in the outer disk (with $L_{\text{Edd}}$ the Eddington luminosity); $\dot{m}_0$ is the peak normalized accretion rate; and $q$ is the fraction of $\dot{m}$ arriving at the black hole, thus accounting for possible outflows. We do not consider the sub-Eddington regime as it involves black holes with higher masses, which should not be able to tidally disrupt main sequence stars. We recall that in the following we assume default values $q=1$, $a=1$, $\alpha_s \beta_h = 1$, $k_\star/f_\star = 0.02$ for the parameters appearing in the jet luminosity. In particular, the choice to set the spin $a=1$ is justified by iron-K$\alpha$ measurements of AGN spins (Brenneman 2013; Reynolds 2013), on which our galaxy formation model is calibrated (Sesana et al. 2014). Clearly, if all black holes had low spins $a \ll 1$ our jetted TDE rates would be significantly decreased, but such a choice seems hard to reconcile with iron-K$\alpha$ observations (see discussion in Sesana et al. 2014).

The total energy release per TDE is given by

$$\mathcal{E}_{\text{jet}} \approx L_{\text{jet}} P_{\text{orb}} \approx 4 \times 10^{52} \, \text{erg} \, \text{erg} \, M_{\star,\odot} \frac{f(a)}{\beta_h \alpha_s} q(\dot{m}/\dot{m}_0), \quad (10)$$

which should be less than $\dot{M}c^2$. We note that a jet luminosity $L_{jet} = \eta_{jet} \dot{M}c^2$ and $\eta_{jet} \sim 1$ are achieved if the disk is magnetically arrested, but the efficiency factor may be smaller. Also, the gravitational binding energy is much lower, so we need to rely on an energy extraction, for example via the Blandford–Znajek mechanism to have powerful jets.

The observable non-thermal luminosity (which, as mentioned before, is identified with the bolometric luminosity of our SED model) is related to the jet luminosity by accounting for the efficiency of energy conversion from Poynting to photon luminosity $\eta_{jet}$ and for the beaming factor $\mathcal{B} = \min(4\pi/\Delta\Omega, 2\Gamma^2)$, where $\Delta\Omega$ is the solid angle occupied by the jet (Krolik & Piran 2012). For a two-sided jet with a jet opening angle $\theta_{jet}$, we have $\mathcal{B} = \min(1/(1 - \cos \theta_{jet}), 2\Gamma^2)$. Therefore, $L_{bol} = L_{jet,\text{obs}} = 2 \eta_{jet, -2} \mathcal{B} L_{jet}$ for $\theta_{jet} \sim 5^\circ$, $\Gamma = 10$, and $\eta_{jet} = 0.01$.

With the observed local rate density, $\dot{n}_{tde,0} = 0.03 \, \text{Gpc}^{-3} \, \text{yr}^{-1}$, the luminosity density is estimated to be

$$\begin{aligned} Q_{\text{TDEjet}} & \approx \eta_{\text{jet}} \mathcal{B} \mathcal{E}_{\text{jet}} \dot{n}_{\text{tde},0} \\ & \approx 2 \times 10^{42} \, \text{erg} \, \text{Mpc}^{-3} \, \text{yr}^{-1} \eta_{\text{jet},-2} \\ & \times M_{\star,\odot} \frac{f(a)}{\beta_h \alpha_s} q(\dot{m}/\dot{m}_0) \frac{\dot{n}_{\text{tde},0}}{0.03 \, \text{Gpc}^{-3} \, \text{yr}^{-1}}. \end{aligned} \quad (11)$$

As discussed in Sect. 5.1, this implies a baryon loading factor $\eta_{CR}$ of order 100 to account for the observed flux of UHECR, which is $\sim 10^{44} \, \text{erg Mpc}^{-3} \, \text{yr}^{-1}$.
### 4.4. Redshift evolution of the black hole mass function

To model the cosmological evolution of massive black holes in their galactic hosts, we utilize the semi-analytic galaxy formation model of Barausse (2012; with incremental improvements described in Sesana et al. 2014; Antonini et al. 2015a,b; Bonetti et al. 2018a,b), adopting the default calibration of Barausse et al. (2017). The model describes the cosmological evolution of baryonic structures on top of Dark Matter merger trees produced with the extended Press-Schechter formalism, modified to more closely reproduce the results of N-body simulations within the ΛCDM model (Press & Schechter 1974; Parkinson et al. 2008). Among the baryonic structures that are evolved along the branches of the merger trees, and which merge at the nodes of the tree, are a diffuse, chemically primordial intergalactic medium, either shock-heated to the Dark Matter halo's virial temperature or streaming into the halo in cold filaments (the former is more common at low redshift and high halo masses, and the latter in small systems at high redshifts; Dekel & Birnboim 2006; Cattaneo et al. 2006; Dekel et al. 2009); a cold interstellar medium (with either disk- or bulge-like geometry), which forms from the cooling of the intergalactic medium or from the above-mentioned cold accretion flows, and which can give rise to star formation in a quiescent fashion or in bursts (Sesana et al. 2014); pc-scale nuclear star clusters, forming from the migration of globular clusters to the galactic center or by in situ star formation (Antonini et al. 2015a,b); a central massive black hole, feeding from a reservoir of cold gas, brought to the galactic center, for example by major mergers and disk bar instabilities. Our semi-analytic model also accounts for feedback processes on the growth or structures (namely from supernovae and from the jets and outflows produced by AGNs), and for the sub-pc evolution of massive black holes, for example the evolution of black hole spins (Barausse 2012; Sesana et al. 2014), migration of binaries due to gas interactions, stellar hardening and triple massive black hole interactions (Bonetti et al. 2018a,b), gravitational-wave emission (Klein et al. 2016).

For the purposes of this paper, the crucial input provided by our model is the evolution of the TDE luminosity function. We determine the TDE comoving rate density $\dot{n}_{TDE}(z, L)$ by combining the TDE rate per galaxy $\dot{N}_{TDE}$ and the black hole comoving density $n_{bh}(z, M_{bh})$, and using the black hole mass and jet luminosity relation (Eq. (9)).

For the properties of the jet, we distinguish between the high state, characterized by a high jet efficiency, and the medium state, characterized by a lower jet efficiency. We set these parameters in order to be consistent with the observations of Swift J1644+57, which should be associated with a black hole of mass $M_{bh} \gtrsim 7 \times 10^6 M_{\odot}$ (Seifina et al. 2017) and reaches a bolometric luminosity $L_{bol} \gtrsim 10^{48}$ erg s⁻¹ in the high state. Therefore, we have $\theta_{jet} = 5^{\circ}$ and $\eta_{jet} = 0.35$ in the high state and $\theta_{jet} = 5^{\circ}$ and $\eta_{jet} = 0.01$ in the medium state. For instance, a black hole of mass $M_{bh} \gtrsim 7 \times 10^6 M_{\odot}$ is associated with $L_{bol} \approx 2 \times 10^{48}$ erg s⁻¹ in the high state and $L_{bol} \approx 7 \times 10^{46}$ erg s⁻¹ in the medium state. The other cases that we consider are shown in Table 2. As explained in Sect. 4.2, for high masses $M_{bh} > 10^8 M_{\odot}$, only highly spinning black holes could lead to observable flares. For completeness, we also account for this case in our study. The black hole mass functions at different redshifts and detailed comparisons of the predictions of our model to observational determinations are given in Appendix C.

It is interesting to notice that the luminosity function of jetted TDEs is dominated by high luminosities (hence low black

<table><caption>**Table 2.** Observed jet luminosities as a function of black hole mass $M_{bh}(M_{\odot})$ in the medium state $L_{bol,med}$ (erg s⁻¹) for $\theta_{jet} = 5^{\circ}$ and $\eta_{jet} = 0.01$, and in the high state $L_{bol,high}$ (erg s⁻¹) for $\theta_{jet} = 5^{\circ}$ and $\eta_{jet} = 0.35$.</caption><thead><tr><th>M<sub>bh</sub><br>(M<sub>☉</sub>)</th><th>L<sub>bol,med</sub> [L<sub>pk</sub>]<br>(erg s<sup>-1</sup>)</th><th>L<sub>bol,high</sub> [L<sub>pk</sub>]<br>(erg s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>7 × 10<sup>8</sup></td><td>6.8 × 10<sup>45</sup> [10<sup>45</sup>]</td><td>2.1 × 10<sup>47</sup> [10<sup>46.5</sup>]</td></tr><tr><td>7 × 10<sup>7</sup></td><td>2.1 × 10<sup>46</sup> [10<sup>45.5</sup>]</td><td>6.8 × 10<sup>47</sup> [10<sup>47</sup>]</td></tr><tr><td>7 × 10<sup>6</sup></td><td>6.8 × 10<sup>46</sup> [10<sup>46</sup>]</td><td>2.1 × 10<sup>48</sup> [10<sup>47.5</sup>]</td></tr><tr><td>7 × 10<sup>5</sup></td><td>2.1 × 10<sup>47</sup> [10<sup>46.5</sup>]</td><td>6.8 × 10<sup>48</sup> [10<sup>48</sup>]</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Comoving TDE luminosity density in their medium state as a function of redshift, as derived in our model, for $\dot{N}_{TDE} = 10^{-5}$ yr⁻¹. The different luminosities correspond to different black hole masses.</figcaption></figure>

hole masses) in our model, unlike the distribution of Sun et al. (2015). This stems from the flat black hole mass functions at low masses (Fig. C.1) combined with the $L_{jet} \propto M_{bh}^{-1/2}$ relation (Eq. (9)). It implies, quite naturally, which of the observed very bright objects such as Swift J1644+57 are the dominant ones in the population. These objects thus set the maximum bolometric luminosity $L_{max}$ in the luminosity function, which we introduce as a cutoff in our population model. This also implies that the diffuse flux of UHECRs will be dominantly produced by the most luminous objects in their medium state.

Figure 6 shows that the corresponding TDE comoving rate density remains almost constant up to redshift ~3 for luminosities $\gtrsim 10^{45.5}$ erg s⁻¹, which dominate in the production of cosmic-ray and neutrino fluxes in our framework.

### 5. Diffuse UHECR and neutrino fluxes from a TDE population

In the following we calculate the diffuse cosmic-ray and neutrino fluxes, and the composition of cosmic rays by considering a population of jetted TDEs. All primed quantities are in the jet comoving frame, all quantities with superscript c are in the source comoving frame, and all other quantities are in the observer frame. The fluxes depend on the spectra produced by each source, on the comoving rate of TDEs (detailed in the previous subsections), and on the cosmic-ray propagation to the Earth. During the extragalactic propagation, cosmic rays may interact with the cosmic microwave background (CMB) and the extragalactic background light (EBL) through photonuclear
interactions. Because of these processes, they may lose energy and create secondary particles in the case of nuclei. In our work we consider EBL models from Kneiske et al. (2004) and Stecker et al. (2006).

The diffuse cosmic-ray flux is given by

$$ \Phi_{\text{CR}}(E_{\text{CR}}) = \frac{c}{4\pi H_0} \int_{z_{\text{min}}}^{z_{\text{max}}} \int_{L_{\text{min}}}^{L_{\text{max}}} dz \, dL \frac{f_s \xi_{\text{CR}} \dot{N}_{\text{TDE}} dn_{\text{bh}}(z, L)/dL}{\sqrt{\Omega_{\text{M}}(1+z)^3 + \Omega_{\Lambda}}} \times F_{\text{CR,s,p}}^{c}(E_{\text{CR}}^c, z, L) t_{\text{dur}}^c, \quad (12) $$

where $\Omega_{\text{M}} = 0.3$ and $\Omega_{\text{L}} = 0.7$ are our fiducial cosmological parameters, $H_0 = 70 \text{ km s}^{-1} \text{Mpc}^{-1}$ is the Hubble constant, $F_{\text{CR,s,p}}^{c}(E_{\text{CR}}^c, z, L)$ is the spectrum obtained after the propagation of cosmic rays from a source at redshift $z$ (per bin of comoving energy $E_{\text{CR}}^c$ and per unit of comoving time $t^c$), and $t_{\text{dur}}^c$ is the duration of the emission in the source comoving frame. The TDE rate $\dot{N}_{\text{TDE}}$ is the rate of TDEs per galaxy and $dn_{\text{bh}}(z, L)/dL$ is the comoving black hole density per (jet) luminosity bin; $f_s$ is the fraction of the jetted TDE population, calculated in Sect. 4, which contributes to the production of UHECRs.

Similarly, the diffuse neutrino flux reads

$$ \Phi_{\nu}(E_{\nu}) = \frac{c}{4\pi H_0} \int_{z_{\text{min}}}^{z_{\text{max}}} \int_{L_{\text{min}}}^{L_{\text{max}}} dz \, dL \frac{f_s \xi_{\text{CR}} \dot{N}_{\text{TDE}} dn_{\text{bh}}(z, L)/dL}{\sqrt{\Omega_{\text{M}}(1+z)^3 + \Omega_{\Lambda}}} \times F_{\nu,s}^{c}(E_{\nu}^{c}, L) t_{\text{dur}}^{c}, \quad (13) $$

where $F_{\nu,s}^{c}(E_{\nu}^{c}, L)$ is the neutrino flux, per comoving energy and per comoving time, for a source with luminosity $L$.

Due to the flat evolution of the TDE comoving density rates up to $z \sim 3$, we can safely use the jetted TDE luminosity distribution at $z=0$ in the above equations and separate the integrals in $L$ and $z$. We checked in particular that our results were similar when using the distribution function at redshifts $z \le 3$. The impact of the redshift evolution on the cosmic-ray spectrum and on the neutrino flux level is also limited, and close to a uniform evolution as described in Kotera et al. (2010).

## 5.1. Final spectrum and composition of cosmic rays

We show in Fig. 7 the cosmic-ray spectrum obtained for an injection of 70% Si and 30% Fe, a spectral index $\alpha = 1.5$, an acceleration efficiency $\eta_{\text{acc}} = 0.1$, a TDE source evolution, and a fraction $\xi_{\text{CR}}f_s = 1.4\%$ of the maximum local event rate density $\dot{n}_{\text{tde},0} = 155 \text{ Gpc}^{-3} \text{yr}^{-1}$, computed from the TDE rate per galaxy obtained in the case of core galaxies. This heavy composition could be injected for example by the core of disrupted stars. We recall that we consider a production of UHECRs dominated by medium states.

Superimposed are the data from the Auger experiment (The Pierre Auger Collaboration et al. 2015) and from the Telescope Array experiment (Fukushima 2015) shown with their statistical uncertainties. We note that the systematic uncertainty on the energy scale is 14% for the Auger experiment.

We also show in Fig. 8 the corresponding mean and standard deviation for the depth of the maximum of the air showers, $\langle X_{\text{max}} \rangle$ and $\sigma(X_{\text{max}})$. It is represented by a gray band, due to uncertainties related to UHECR-air interaction models EPOS-LHC (Werner et al. 2006; Pierog et al. 2015), SIBYLL 2.1 (Ahn et al. 2009), or QGSJET II-04 (Ostapchenko 2006, 2013). Superimposed are the data from the Auger experiment for the composition of UHECRs, $\langle X_{\text{max}} \rangle$ and $\sigma(X_{\text{max}})$ (Aab et al. 2014),

<figure><img src="image_8.png" /><figcaption>Fig. 7. Diffuse cosmic-ray spectrum from a population of jetted TDEs (calculated in their medium states) obtained for an injection of 70% Si and 30% Fe, with spectral index $\alpha = 1.5$, $\xi_{\text{CR}}f_s = 1.4 \times 10^{-2}$ and source evolution derived in this work, with maximum bolometric luminosity in the population $L_{\text{max}} = 6.8 \times 10^{46} \text{ erg s}^{-1}$. We show the total spectrum (black) and its composition. We superimpose data from the Auger experiment (black dots, The Pierre Auger Collaboration et al. 2015), and from the Telescope Array experiment (gray dots, Fukushima 2015) for which only statistical uncertainties are shown.</figcaption></figure>

which are in good agreement with our model. Only statistical uncertainties are displayed; systematic uncertainties are at most $\pm 10 \text{ g cm}^{-2}$ for $\langle X_{\text{max}} \rangle$ and $\pm 2 \text{ g cm}^{-2}$ for $\sigma(X_{\text{max}})$. Our results are compatible with a light composition at $10^{18.5} \text{ eV}$, shifting toward a heavier composition for increasing energy.

With increased cutoff bolometric luminosities $L_{\text{max}}$, harder injection spectra are needed in order to compensate for the abundant production of nucleons at low energies, which softens the overall spectrum (typically $\alpha = 1$ is required for $L_{\text{max}} = 10^{47} \text{ erg s}^{-1}$). We present the case of the injection of a dominant fraction of heavy elements; the injection of more intermediate elements, such as the CNO group, is possible at the cost of increasing the acceleration efficiency $\eta_{\text{acc}}$, and hardening the injection spectrum further in order to achieve the highest energies.

Our model allows for a fit of both the UHECR spectrum and composition of the Auger observations, as long as the dominant sources supply luminosities $\le 10^{46.5} \text{ erg s}^{-1}$, which is a value that is consistent with the observed Swift event. We note that if the high states were dominant for the UHECR production, we would not be able to fit the data; because of the strong photodisintegration of heavy elements in the very dense radiation background, we would obtain a large production of nucleons below $10^{19} \text{ eV}$ and no survival of heavy elements at the highest energies. However, because of its limited duration ($t_{\text{dur}} \sim 10^5 \text{ s}$ is chosen as an upper bound), the high state is unlikely to be dominant. In a more refined model we should account for the evolution of the luminosity of the jet, which should decrease during the event duration.

The disrupted stellar object provides material (protons and heavier nuclei) that can be injected and accelerated in a jet. As already emphasized, the composition of this material is poorly constrained; it could be similar to the composition of the stellar object or modified during the disruption process. It is interesting that white dwarfs could be commonly disrupted by intermediate mass black holes. These stars could be a source of copious amounts of CNO nuclei, which seem to be observed in the composition of UHECRs measured by Auger, as noted in
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Mean and standard deviation of $X_{\text{max}}$ for the spectrum shown in Fig. 7 (pale gray band). We also show Auger measurements (Aab et al. 2014) with uncertainties (black dots) and simulation results for pure proton injection (red band) and pure iron injection (blue band). The bands are obtained by accounting for hadronic model uncertainties (EPOS-LHC, SIBYLL 2.1 or QGSJET II-04).</figcaption></figure>

Alves Batista & Silk (2017). For completeness, we tested in this study various injection fractions, and we present one case that allows us to fit the Auger data well. We note that a deviation of 5% in the composition does not largely affect the fit to the data within the error bars, given the uncertainties on the other jetted TDE parameters.

Markers of the occurrence of jets associated with TDEs were detected only very recently. Most TDEs should not power jets and only a small fraction of jetted TDEs should point toward the observer, depending on the jet opening angle. Therefore, the properties of these objects are still subject to large uncertainties. From an observational perspective, the jetted TDEs detected recently are very luminous events with a peak isotropic luminosity $L_{\text{pk}} \sim 10^{47}-10^{48}$ erg s$^{-1}$, and a local event rate density is of $\dot{n}_{\text{tde},0} \sim 0.03 \text{ Gpc}^{-3} \text{ yr}^{-1}$ (e.g., Farrar & Piran 2014). On the other hand, normal TDEs are less luminous and are characterized by a higher local event rate density $\dot{n}_{\text{tde},0} = 10^2 \text{ Gpc}^{-3} \text{ yr}^{-1}$ (e.g., Donley et al. 2002). However, the characteristics of this new population, mainly their luminosity distribution and comoving event rate density, are difficult to infer due to the scarcity of observations. From our population model, the maximum local event rate density that we can

<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Diffuse neutrino flux for three flavors from a population of jetted TDEs with the same properties as in Fig. 7 (neutrino production calculated in their high states). We also show the diffuse neutrino flux measured by the IceCube experiment (IceCube Collaboration et al. 2017), and the projected limits for GRAND (Fang et al. 2017), ARA/ARIANNA (Allison et al. 2015; Barwick et al. 2015), and PO-EMMA (Neronov et al. 2017).</figcaption></figure>

expect reaches $\dot{n}_{\text{tde},0} \sim 2 \times 10^2 \text{ Gpc}^{-3} \text{ yr}^{-1}$ for core galaxies and $\dot{n}_{\text{tde},0} \sim 3 \times 10^3 \text{ Gpc}^{-3} \text{ yr}^{-1}$ for power-law galaxies.

The fraction needed to fit the UHECR spectrum of the Auger observations, $\xi_{\text{CR}}f_s = 1.4 \times 10^{-2}$, can account for example for low UHECR injection rates $\xi_{\text{CR}}$, and/or for population constraints, such as the fraction of TDE jets pointing toward the observer. Assuming the low rate inferred from the observations of $\dot{n}_{\text{tde},0} \sim 0.03 \text{ Gpc}^{-3} \text{ yr}^{-1}$ for the jetted events pointed toward the observer, a baryon loading of $\xi_{\text{CR}} \sim 187$ is required. This value is consistent with the non-detection of neutrinos from Swift J1644+57, which implies an upper limit to the baryon density of a few 100 (Senno et al. 2017; Guépin & Kotera 2017).

At the lowest energies, the high state could contribute marginally to the diffuse cosmic-ray flux, as shown in Fig. 5. The strong photodisintegration in the high state leads to a strong production of nucleons below $10^{19}$ eV, which would add a new component to the spectrum and make the composition lighter. This would lead to a better fit of the composition in Fig. 8.

## 5.2. Diffuse neutrino flux

The TDE event density obtained by fitting the Auger data with our UHECR spectrum allows us to calculate the diffuse neutrino flux from a population of TDEs, by considering the fraction $\xi_{\text{CR}}f_s$ calculated above. As shown in Fig. 9, the diffuse neutrino flux from jetted TDEs contributes marginally to the total diffuse neutrino flux observed by IceCube (IceCube Collaboration et al. 2017). As it peaks at high energies, around $10^{16}$ eV, it could be a good target for future generation detectors. However, we note that this flux is too low to be detectable with ARA/ARIANNA, POEMMA, and GRAND at even higher energies. Its high-energy cutoff reduces the flux at higher energies, and therefore it lies just below the GRAND sensitivity limit.

# 6. Discussion and conclusion

We assessed in this study the production of UHECRs and neutrinos by a population of TDEs. In our model, the disruption of a stellar object launches a relativistic jet, where internal shocks
can accelerate a part of the disrupted material, namely light and heavy nuclei. This scenario is connected to recent observations and analytic studies, favoring a jetted model for some very luminous events. In such a case, material from the disrupted object can be injected and accelerated inside the jet, and can experience interactions before escaping and propagating in some cases toward the Earth. However, other scenarios could be contemplated; for instance, a substantial fraction of the accreted material could be ejected as a wind where particles could be linearly accelerated<fnref n="2" />.

The bulk Lorentz factor Γ and the opening angle of the jet $\theta_{jet}$ are two additional important quantities impacting our results. The bulk Lorentz factor of the jet impacts the observed jet isotropic luminosity and the energy of detected cosmic rays and neutrinos. The dynamical time also scales as Γ and the photon energy density as $\Gamma^{-4}$, thus an increase of a factor of a few in Γ could lead to a drastic cut in the photodisintegration rates. Here we use the fiducial value $\Gamma = 10$, but for larger values we expect that the survival of nuclei would be favored, leading to lower nucleon production at lower energies, and thus to a larger parameter space allowing for a good fit for the diffuse UHECR spectrum. Our choice is conservative in this sense. On the other hand, the neutrino production would be consequently reduced. The opening angle of the jet is also not well constrained; therefore, we adopt a small value $\theta_{jet} \sim 5^\circ$ for the high and medium states (e.g., Bloom et al. 2011; Burrows et al. 2011). Like the bulk Lorentz factor, this parameter is also involved in the model that we use to link the black hole mass to the isotropic luminosity of the jet. The jet can be seen only if it is pointing toward the observer. However, we note that the effective opening angle for cosmic rays might be higher than the usual opening angle as cosmic rays can experience small deflections inside the jet; thus, misaligned jetted TDEs characterized by a higher rate than aligned events might also contribute to the diffuse cosmic-ray flux.

While finalizing this paper, we became aware of the independent work of Biehl et al. (2017) on a similar topic. These authors show that the acceleration of nuclei in jets created by the tidal disruption of white dwarfs can lead to a simultaneous fit of the UHECR data and the measured IceCube neutrino flux in the PeV range. One major difference with our study is that we include a detailed jetted TDE population study by modeling the luminosity function and rate evolution in redshift. Our conclusions are also different from theirs, in so far as we cannot fit the observed diffuse IC neutrino flux with our TDE population model. This negative result is consistent with several arguments already highlighted in previous works by Dai & Fang (2017) and Senno et al. (2017). In particular, the absence of observed neutrino multiplets in the IceCube data gives a lower limit of $\ge 100–1000 \text{ Gpc}^{-3} \text{ yr}^{-1}$, which is significantly higher than the rate of jetted TDEs pointing toward us inferred from observations of $\sim 0.03 \text{ Gpc}^{-3} \text{ yr}^{-1}$, and higher than rates with dimmer luminosities also constrained by X-ray observations. In addition, large baryon loadings with $\ge 1000$ are ruled out as such values would imply that Swift J1466+57 would have been observed in neutrinos. Also, a large baryon loading factor implies a total TDE energy of $\ge 10^{34} \text{ erg}$, which violates the energetics argument.

Our model is able to reproduce with a reasonable accuracy and for a reasonable range of parameters the observations from the Auger experiments, and TDE powering jets therefore appear to be good candidates for the production of UHECRs. Our

results are consistent with other TDE studies that also obtain good fits to UHECR data: Zhang et al. (2017) stress that oxygen-neon-magnesium white-dwarf TDE models could provide good fits, but do not account for photodisintegration in the vicinity of the source because they used a steeper luminosity function. Our model can account for these interactions, and allows us to explore the parameter space for the radiation field, the injection and the composition. This is important for our flatter luminosity function, which predicts that the luminosity density is dominated by the highest-power TDEs, i.e., the effective luminosity is $L \sim 2 \times 10^{48} \text{ erg s}^{-1}$ in the high state.

The associated transient HE neutrinos could be detected for single nearby sources (at distances of a few tens of Mpc) with IceCube and upcoming instruments at higher energies such as GRAND or POEMMA. The diffuse flux would be within reach of IceCube in the next decade. Its detection would be more challenging for future generation instruments aiming at the detection of ultra-high-energy neutrinos, due to a high-energy cutoff below $10^{17} \text{ eV}$.

Among the other transient UHECR nuclei models that have been suggested to explain the UHECR data (e.g., fast rotating pulsars, Fang et al. 2012, 2014; Kotera et al. 2015 or GRBs, Wang et al. 2008; Murase et al. 2008; Globus et al. 2015a,b), the jetted TDE model has the interesting property of presenting two different states (low and high) leading to optimal production of both UHECRs and neutrinos. In addition, the jetted TDE scenario appears mildly constrained by photon observations. Within our model, we demonstrated that the observed Swift J1466+57 can be seen as a typical source that would dominate the production of UHECRs and neutrinos. Even under this constraint, the wide range of variation allowed for several free parameters (for example the Lorentz factor of the outflow, as discussed earlier) enables us to correctly fit the cosmic-ray data. A specific signature of this scenario is thus difficult to infer. A direct multi-messenger signal with TDE photons associated with the emission of neutrinos from a single source appears to be the way to validate this scenario.

*Acknowledgements.* We thank Tanguy Pierog for his continuous help on hadronic interactions and EPOS. We are indebted to Marta Volonteri for enlightening conversations on the physics of massive black holes. We also thank Joe Silk, Rafael Alves Batista, and Nicholas Senno for very fruitful discussions. This work is supported by the APACHE grant (ANR-16-CE31-0001) of the French Agence Nationale de la Recherche. CG is supported by a fellowship from the CFM Foundation for Research and by the Labex ILP (reference ANR-10-LABX-63, ANR-11-IDEX-0004-02). EB is supported by the European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No. 690904. For our simulations, we have made use of the Horizon Cluster, hosted by the Institut d'Astrophysique de Paris. We thank Stephane Rouberol for running this cluster smoothly for us. The work of KM is supported by the Alfred P. Sloan Foundation and NSF grant No. PHY-1620777.

## References

Aab, A., Abreu, P., Aglietta, M., et al. 2014, Phys. Rev. D, 90, 122005

Aartsen, M. G., Ackermann, M., Adams, J., et al. 2014, ApJ, 796, 109

Aartsen, M. G., Ackermann, M., Adams, J., et al. 2015, ApJ, 807, 46

Ahn, E.-J., Engel, R., Gaisser, T. K., Lipari, P., & Stanev, T. 2009, Phys. Rev. D, 80, 094003

Aird, J., Nandra, K., Laird, E. S., et al. 2010, MNRAS, 401, 2531

Allison, P., Auffenberg, J., Bard, R., et al. 2015, Astropart. Phys., 70, 62

Alves Batista, R., & Silk, J. 2017, Phys. Rev. D, 96, 103003

Alves Batista, R., Boncioli, D., di Matteo, A., van Vliet, A., & Walz, D. 2015, JCAP, 10, 063

Alves Batista, R., Dundovic, A., Erdmann, M., et al. 2016, JCAP, 5, 038

Antonini, F., Barausse, E., & Silk, J. 2015a, ApJ, 812, 72

Antonini, F., Barausse, E., & Silk, J. 2015b, ApJ, 806, L8

Barausse, E. 2012, MNRAS, 423, 2533

<footnote n="2">² Zhang et al. (2017) show that UHECR acceleration is difficult in the wind model.</footnote>
Barausse, E., Shankar, F., Bernardi, M., Dubois, Y., & Sheth, R. K. 2017, MNRAS, 468, 4782

Barwick, S. W., Berg, E. C., Besson, D. Z., et al. 2015, Astropart. Phys., 70, 12

Biehl, D., Boncioli, D., Lunardini, C., & Winter, W. 2017, ArXiv e-prints [arXiv:1711.03555]

Bloom, J. S., Giannios, D., Metzger, B. D., et al. 2011, Science, 333, 203

Bonetti, M., Haardt, F., Sesana, A., & Barausse, E. 2018a, MNRAS, 477, 3910

Bonetti, M., Sesana, A., Barausse, E., & Haardt, F. 2018b, MNRAS, 477, 2599

Brenneman, L. 2013, Measuring the Angular Momentum of Supermassive Black Holes

Burrows, D. N., Kennea, J. A., Ghisellini, G., et al. 2011, Nature, 476, 421

Bustamante, M., Heinze, J., Murase, K., & Winter, W. 2017, ApJ, 837, 33

Cattaneo, A., Dekel, A., Devriendt, J., Guiderdoni, B., & Blaizot, J. 2006, MNRAS, 370, 1651

Celotti, A., & Ghisellini, G. 2008, MNRAS, 385, 283

Cenko, S. B., Krimm, H. A., Horesh, A., et al. 2012, ApJ, 753, 77

Cummings, J. R., Barthelmy, S. D., Beardmore, A. P., et al. 2011, GRB Coordinates Network, 11823

Dai, L. & Fang, K. 2017, MNRAS, 469, 1354

Dekel, A., & Birnboim, Y. 2006, MNRAS, 368, 2

Dekel, A., Birnboim, Y., Engel, G., et al. 2009, Nature, 457, 451

Dermer, C. D., & Menon, G. 2009, High Energy Radiation from Black Holes: Gamma Rays, Cosmic Rays, and Neutrinos (Princeton University Press)

Donley, J. L., Brandt, W. N., Eracleous, M., & Boller, T. 2002, AJ, 124, 1308

Fang, K., Kotera, K., & Olinto, A. V. 2012, ApJ, 750, 118

Fang, K., Kotera, K., Murase, K., & Olinto, A. V. 2014, Phys. Rev. D, 90, 103005

Fang, K., Álvarez-Muñiz, J., Alves Batista, R., et al. 2017, ArXiv e-prints [arXiv:1708.05128]

Farrar, G. R., & Gruzinov, A. 2009, ApJ, 693, 329

Farrar, G. R., & Piran, T. 2014, Tidal disruption jets as the source of Ultra-High Energy Cosmic Rays

Fukushima, M. 2015, in Eur. Phys. J. Web Conf., 99, 04004

Globus, N., Allard, D., Mochkovitch, R., & Parizot, E. 2015a, MNRAS, 451, 751

Globus, N., Allard, D., & Parizot, E. 2015b, Phys. Rev. D, 92, 021302

Guépin, C., & Kotera, K. 2017, A&A, 603, A76

Hills, J. G. 1975, Nature, 254, 295

Hopkins, P. F., Richards, G. T., & Hernquist, L. 2007, ApJ, 654, 731

IceCube Collaboration (Aartsen, M. G., et al.) 2017, ArXiv e-prints [arXiv:1710.01191]

Kalfountzou, E., Civano, F., Elvis, M., Trichas, M., & Green, P. 2014, MNRAS, 445, 1430

Kesden, M. 2012, Phys. Rev. D, 85, 024037

Kippenhahn, R., & Weigert, A. 1994, Stellar Structure and Evolution, 192

Klein, A., Barausse, E., Sesana, A., et al. 2016, Phys. Rev. D, 93, 024003

Kneiske, T. M., Bretz, T., Mannheim, K., & Hartmann, D. H. 2004, A&A, 413, 807

Komossa, S. 2015, JHEAP, 7, 148

Koning, A. J., Hilaire, S., & Duijvestijn, M. C. 2005, AIP Conf. Proc., 769, 1154

Kormendy, J., & Ho, L. C. 2013, ARA&A, 51, 511

Kotera, K., & Olinto, A. V. 2011, ARA&A, 49, 119

Kotera, K., Allard, D., Murase, K., et al. 2009, ApJ, 707, 370

Kotera, K., Allard, D., & Olinto, A. V. 2010, J. Cosmol. Astropart. Phys., 10, 013

Kotera, K., Amato, E., & Blasi, P. 2015, JCAP, 8, 026

Krolik, J. H., & Piran, T. 2012, ApJ, 749, 92

La Franca, F., Fiore, F., Comastri, A., et al. 2005, ApJ, 635, 864

Lacy, M., Ridgway, S. E., Sajina, A., et al. 2015, ApJ, 802, 102

Lauer, T. R., Tremaine, S., Richstone, D., & Faber, S. M. 2007, ApJ, 670, 249

Lodato, G., King, A. R., & Pringle, J. E. 2009, MNRAS, 392, 332

Luminet, J.-P., & Pichon, B. 1989, A&A, 209, 103

Lunardini, C., & Winter, W. 2017, Phys. Rev. D, 95 123001

Madau, P., & Rees, M. J. 2001, ApJ, 551, L27

Merloni, A., & Heinz, S. 2008, MNRAS, 388, 1011

Mücke, A., Engel, R., Rachen, J. P., Protheroe, R. J., & Stanev, T. 2000, Comput. Phys. Commun., 124, 290

Murase, K. 2008, AIP Conf. Ser., eds. Y.-F. Huang, Z.-G. Dai, & B. Zhang, 1065, 201

Murase, K., & Nagataki, S. 2006, Phys. Rev. D, 73, 063002

Murase, K., & Takami, H. 2009a, in Proceedings, 31st International Cosmic Ray Conference (ICRC 2009): Lodz, Poland, 2009, 1181

Murase, K., & Takami, H. 2009b, ApJ, 690, L14

Murase, K., Ioka, K., Nagataki, S., & Nakamura, T. 2008, Phys. Rev. D, 78, 023005

Neronov, A., Semikoz, D. V., Anchordoqui, L. A., Adams, J. H., & Olinto, A. V. 2017, Phys. Rev. D, 95, 023004

Ostapchenko, S. 2006, Nucl. Phys. B Proc. Suppl., 151, 147

Ostapchenko, S. 2013, Eur. Phys. J. Web Conf., 52, 02001

Parkinson, H., Cole, S., & Helly, J. 2008, MNRAS, 383, 557

Phinney, E. S. 1989, in IAU Symp., The Center of the Galaxy, ed. M. Morris, 136, 543

Pierog, T., Karpenko, I., Katzy, J. M., Yatsenko, E., & Werner, K. 2015, Phys. Rev. C, 92, 034906

Piran, T. 2004, Rev. Mod. Phys., 76, 1143

Piran, T., Sądowski, A., & Tchekhovskoy, A. 2015, MNRAS, 453, 157

Press, W. H., & Schechter, P. 1974, ApJ, 187, 425

Rachen, J. P. 1996, PhD Thesis, Bonn University

Reynolds, C. S. 2013, Class. Quant. Grav., 30, 244004

Schulze, A., Bongiorno, A., Gavignaud, I., et al. 2015, MNRAS, 447, 2085

Seifina, E., Titarchuk, L., & Virgilli, E. 2017, A&A, 607, A38

Senno, N., Murase, K., & Meszaros, P. 2017, ApJ, 838, 3

Sesana, A., Barausse, E., Dotti, M., & Rossi, E. M. 2014, ApJ, 794, 104

Shankar, F. 2013, Class. Quant. Grav., 30, 244001

Sironi, L., Petropoulou, M., & Giannios, D. 2015, MNRAS, 450, 183

Stecker, F. W., Malkan, M. A., & Scully, S. T. 2006, ApJ, 648, 774

Sun, H., Zhang, B., & Li, Z. 2015, ApJ, 812, 33

The Pierre Auger Collaboration (Aab, A., et al.) 2015, ArXiv e-prints [arXiv:1509.03732]

Volonteri, M., Lodato, G., & Natarajan, P. 2008, MNRAS, 383, 1079

Wang, J., & Merritt, D. 2004, ApJ, 600, 149

Wang, X.-Y., & Liu, R.-Y. 2016, Phys. Rev. D, 93, 083005

Wang, X.-Y., Razzaque, S., & Mészáros, P. 2008, ApJ, 677, 432

Wang, X.-Y., Liu, R.-Y., Dai, Z.-G., & Cheng, K. S. 2011, Phys. Rev. D, 84, 081301

Werner, K., Liu, F.-M., & Pierog, T. 2006, Phys. Rev., 74, 044902

Zhang, B. T., Murase, K., Oikonomou, F., & Li, Z. 2017, Phys. Rev. D, 96, 063007
# Appendix A: Cosmic-ray maximum energy

We derive analytic estimates of the effective maximum energies in the comoving frame of the jet as a function of the energy loss processes at play. They depend on the characteristics of the event considered, and its radiation region, for example the bolometric luminosity $L_{bol}$, the comoving mean magnetic field $B'$, the time variability $t_{var}$, or the bulk Lorentz factor $\Gamma$.

## A.1. Maximum injection energy

As explained in Sect. 2.1, the maximum injection energy $E'_{Z,max}$ of a nucleus of charge Z is determined by the competition between the acceleration timescale $t'_{acc} = \eta_{acc}^{-1}E'/cZeB'$ and the energy loss timescales $t'_{loss} = \min(t'_{dyn}, t'_{syn}, t'_{IC}, t'_{BH}, t'_{py}, \dots)$. An upper bound of the maximum energy of accelerated particles is thus given by the competition between the acceleration timescale ($t'_{acc}$) and the dynamical timescale ($t'_{dyn}$):

$$ 
\begin{aligned}
E'_{Z,up,dyn} & \sim cZeB' (1+z)^{-1} \Gamma t_{var} \eta_{acc}, \\
& \approx 6.3 \times 10^{17} \text{ eV } Z_1 B'_{2.8} \Gamma_{10} t_{var,2} \eta_{acc,-1}, \\
& \approx 1.6 \times 10^{19} \text{ eV } Z_{26} B'_{2.8} \Gamma_{10} t_{var,2} \eta_{acc,-1}.
\end{aligned} 
\quad (A.1) $$

The competition between the acceleration timescale ($t'_{acc}$) and the synchrotron timescale ($t'_{syn}$) gives

$$ 
\begin{aligned}
E'_{Z,up,syn} & \sim \left[ \frac{6\pi(m_u c^2)^2 e}{(m_e/m_u)^2 \sigma_T} \right]^{\frac{1}{2}} A^2 Z^{-3/2} B'^{-1/2} \eta_{acc,1}^{1/2}, \\
& \approx 2.4 \times 10^{18} \text{ eV } A_1^2 Z_1^{-3/2} B_{2.8}^{-1/2} \eta_{acc,-1}^{1/2}, \\
& \approx 5.7 \times 10^{19} \text{ eV } A_{56}^2 Z_{26}^{-3/2} B_{2.8}^{-1/2} \eta_{acc,-1}^{1/2}.
\end{aligned} 
\quad (A.2) $$

Here, the mean magnetic field $B' = \sqrt{8\pi \int d\epsilon' \epsilon' n'_{\epsilon'} \approx 10^{2.85}}$ G is obtained for a log-parabola SED with peak luminosity $L_{pk} = 10^{46}$ erg s$^{-1}$ and a width $\hat{a} = 0.07$.

The upper bound given by the competition between the acceleration timescale ($t'_{acc}$) and the photohadronic timescale ($t'_{py}$) is computed numerically. For the parameters considered above, we obtain for protons $E'_{p,up,py} \approx 2.5 \times 10^{18}$ eV. The comparison between the different energy loss timescales allows us to determine the limiting energy loss process: for the previous example the dynamical timescale is the limiting timescale.

## A.2. Competition between energy loss processes in the radiation region

The competition between the energy loss processes in the radiation region influences the outgoing cosmic-ray spectrum, and in particular the high-energy cutoffs. By considering the competition between synchrotron losses ($t'_{syn}$) and escape ($t'_{esc}$) for protons,

$$ \frac{6\pi(m_p c^2)^2}{(m_e/m_p)^2 \sigma_T c E'_p B'^2} = \Gamma t_{var}, \quad (A.3) $$

with $E'_p$ the proton energy in the comoving frame, $\sigma_T$ the Thomson cross section, $m_p$ the proton mass, $m_e$ the electron mass, and c the speed of light, we can derive the corresponding high-energy cutoff:

$$ 
\begin{aligned}
E'_{p,max} &= \frac{6\pi(m_p c^2)^2}{(m_e/m_p)^2 \sigma_T c E'_p B'^2} \Gamma t_{var}, \\
&\approx 9.0 \times 10^{18} \text{ eV } B'^{-2.8} \Gamma^{-1} t^{-1}_{var,2}.
\end{aligned} 
\quad (A.4) $$

The competition between synchrotron losses ($t'_{syn}$) and pion production ($t'_{py}$) for a transient event characterized by a hard spectrum gives

$$ \frac{6\pi(m_p c^2)^2}{(m_e/m_p)^2 \sigma_T c E'_p B'^2} = \frac{4\pi\Gamma^5 c^2 t_{var}^2 \epsilon_{pk}}{\langle \sigma_{py} \kappa_{py} \rangle L_{pk}}, \quad (A.5) $$

with $\sigma_{py}$ and $\kappa_{py}$ the photopion cross section and inelasticity, $\epsilon_{pk}$ the peak energy, and $L_{pk}$ the peak luminosity. For $\epsilon_{pk} = 70$ keV, we obtain a high-energy cutoff:

$$ 
\begin{aligned}
E'_{p,max} &= \frac{3(m_p c^2)^2 \langle \sigma_{py} \kappa_{py} \rangle L_{pk}}{2(m_e/m_p)^2 \sigma_T c^3 B'^2 \Gamma^5 t_{var}^2 \epsilon_{pk}}, \\
&\approx 5.0 \times 10^{21} \text{ eV } B'^{-2.8} \Gamma^{-1} L_{pk,46} t^{-2}_{var,2}.
\end{aligned} 
\quad (A.6) $$

Our numerical estimates are evaluated for typical parameters of jetted TDEs (e.g., the characteristics of Swift J1644+57).

For nuclei, the high-energy cutoffs depend on the mass and atomic numbers. As $\gamma_N = E'_{N}/Am_u c^2$, where $m_u$ is the atomic mass unit, for the competition between synchrotron losses ($t'_{syn}$) and escape ($t'_{esc}$) we obtain

$$ 
\begin{aligned}
E'_{N,max} &= \frac{6\pi(m_u c^2)^2 (A/Z)^4}{(m_e/m_u)^2 \sigma_T c B'^2 \Gamma t_{var}}, \\
&\approx 1.4 \times 10^{20} \text{ eV } B'^{-2.8} \Gamma^{-1} t^{-1}_{var,2} A_{56}^4 Z_{26}^{-4},
\end{aligned} 
\quad (A.7) $$

where $A_{56} = A/56$ and $Z_{26} = Z/26$ for iron nuclei.

# Appendix B: Derivation of diffuse neutrino and cosmic-ray fluxes

To calculate the neutrino diffuse flux, we integrate the neutrino flux of a single source over the TDE population. Primed quantities are in the jet comoving frame, quantities with superscript c are in the source comoving frame, and other quantities are in the observer frame. We account for the total number of neutrinos produced by one single source, which depends on the neutrino energy in the source comoving frame $E'_v = (1+z)E_v$ and the bolometric luminosity of the source: $N_{v,s}(E'_v, L)$.

Moreover, for a redshift $z$, we can observe a population of sources characterized by a comoving event rate density $\dot{n}(z, L)$ during an observation time $T_{obs} = (1+z)T_{obs}^c$ (where $T_{obs}$ is the time in the observer frame and $T_{obs}^c$ is the time in the source comoving frame), in a comoving volume,

$$ \frac{dV(z)}{dz} = \frac{c}{H_0} \frac{4\pi D_c^2}{\sqrt{\Omega_M (1+z)^3 + \Omega_\Lambda}}, \quad (B.1) $$

where $\Omega_M = 0.3$ and $\Omega_L = 0.7$ are our fiducial cosmological parameters, and $H_0 = 70$ km s$^{-1}$ Mpc$^{-1}$ is the Hubble constant. The comoving event rate density $\dot{n}(z, L)/dL = \dot{N}_{TDE} dn_{bh}(z, L)/dL$ depends on the TDE rate per galaxy $\dot{N}_{TDE}$ and the black hole comoving density per luminosity $dn_{bh}(z, L)/dL$. Therefore, the diffuse neutrino flux is given by

$$ \frac{dN_{\nu}}{dE_{\nu}}(E_{\nu}) = \frac{1}{4\pi} \int_{z_{\min}}^{z_{\max}} \int_{L_{\min}}^{L_{\max}} \left[ dz \, dL \, f_s \xi_{CR} \frac{d\dot{n}(z, L)}{dL} \right. \\ \left. \times \frac{1}{1+z} \frac{dV(z)}{dz} \frac{1}{4\pi D_c^2} \frac{dN_{v,s}(E_{\nu}^c, L)}{dE_{\nu}^c} \frac{dE_{\nu}^c}{dE_{\nu}} \right] \quad (B.2) $$
As $E_v^c = (1 + z)E_v$, we obtain a diffuse neutrino flux

$$ \Phi_v(E_v) = \frac{c}{4\pi H_0} \int_{z_{\min}}^{z_{\max}} \int_{L_{\min}}^{L_{\max}} dz \, dL \frac{f_s \xi_{\text{CR}} d\dot{h}(z, L)/dL}{\sqrt{\Omega_{\text{M}}(1+z)^3 + \Omega_{\Lambda}}} \times F_{v,s}^c(E_v^c, L)t_{\text{dur}}^c \quad (B.3) $$

where $\Phi_v(E_v) = d^2N_v(E_v)/dE_v dt$ is the diffuse neutrino flux per unit time (observer frame), $F_{v,s}^c(E_v^c, L) = d^2N_{v,s}(E_v^c, L)/dE_v^c dt$ is the number of neutrinos emitted by one single source per bin of comoving energy and per unit of comoving time, and $t_{\text{dur}}^c$ is the duration of the emission in the source comoving frame.

The cosmic-ray diffuse flux is calculated in a similar manner, but in this case we also need to account for the large-scale propagation of cosmic rays

$$ \Phi_{\text{CR}}(E_{\text{CR}}) = \frac{c}{4\pi H_0} \int_{z_{\min}}^{z_{\max}} \int_{L_{\min}}^{L_{\max}} dz \, dL \frac{f_s \xi_{\text{CR}} d\dot{h}(z, L)/dL}{\sqrt{\Omega_{\text{M}}(1+z)^3 + \Omega_{\Lambda}}} \times F_{\text{CR,s,p}}^c(E_{\text{CR}}^c, z, L)t_{\text{dur}}^c \quad (B.4) $$

where $F_{\text{CR,s,p}}^c(E_{\text{CR}}^c, z, L)$ is the spectrum obtained after the propagation of cosmic rays from a source at redshift $z$.

## Appendix C: Evolution of the black hole mass function in redshift

We compare the black hole mass function predicted by our semi-analytic galaxy formation model with the observational determinations of Shankar (2013) and Lauer et al. (2007) (at $z=0$), and with those of Merloni & Heinz (2008) and Schulze et al. (2015) (at $z>0$) in Fig. C.1. The model’s predictions are shown as red bars or blue dots, the first referring to a scenario in which black holes form from low-mass seeds of a few hundred $M_{\odot}$ (e.g., the remnants of Pop III stars; Madau & Rees 2001), and the second representing a model in which black holes descend from “heavy” (~$10^5 M_{\odot}$) seeds arising, for example, from instabilities of protogalactic disks. In more detail, for the latter case we use the model of Volonteri et al. (2008), setting their critical Toomre parameter, which regulates the onset of the instability, to their preferred value $Q_c = 2.5$. The error bars of the model’s points are Poissonian. As can be seen, the agreement with the data is rather good, especially in the mass range relevant for our purposes ($M_{bh} < 10^8 M_{\odot}$). As a further test, we also compared the predictions of our model for the AGN (bolometric) luminosity function with the observations of Hopkins et al. (2007), Lacy et al. (2015), La Franca et al. (2005), and Aird et al. (2010), whose envelope we show in Fig. C.2 as a shaded orange area. We note that we only consider the luminosity function of Aird et al. (2010) at $z < 3$ as it may be underestimated at larger redshifts (Kalfountzou et al. 2014).
<figure><img src="image_12.png" /><figcaption>**Fig. C.1.** Predictions of our model in two different scenarios, light-seed ("Pop III") and heavy-seed ("$Q_c = 2.5$"), for the mass function of black holes as a function of redshift. The model's error bars are Poissonian. For comparison, also shown are the observational determinations of Shankar (2013) and Lauer et al. (2007) (at $z=0$), and those of Merloni & Heinz (2008) and Schulze et al. (2015) (at $z>0$).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. C.2.** Predictions of our model (with the same seed scenarios as in Fig. C.1) for the bolometric AGN luminosity function, compared with observational determinations – Hopkins et al. (2007), Lacy et al. (2015), La Franca et al. (2005), and Aird et al. (2010), the last only considered at $z>3$ –, whose envelope is shown by a shaded orange area.</figcaption></figure>