© ESO 2018

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 highenergy 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 ∝ t ^−5/3, coherent with fallback 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⁴⁸ erg s⁻¹, and can maintain very high bolometric luminosities (L_bol ~ 10⁴⁷ erg s⁻¹) 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 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 ṅ_tde > 1 Gpc⁻³ yr⁻¹. 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 semianalytic 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 ≡ Q/10 ^x in cgs units unless specified otherwise, and except for particle energies, which are in E_x ≡ 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′ = Γ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={\mathrm{\Gamma }}^{2}c t_{\text{var}}=3\times {10}^{14}\text{cm}{\mathrm{\Gamma }}_1^2{t}_{\text{var,2}}$, where t_var = 10²t_var,2 s and Γ = 10Γ₁ 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 $\text{d}{N}_{\text{CR}}/\text{d}{E}^{\prime }=\mathcal{A}\hspace{0.17em}\mathrm{\Sigma }z fz{E^{\prime }}^{-\alpha }\text{exp}(-E^{\prime }/{E^{\prime }}_{\text{Z,max}})$ with an exponential cutoff at $E_{\text{Z,max}}^{\prime }$ for nuclei of charge Z. Here 𝒜 is a normalization constant and f_Z is the fraction of elements with charge number Z, such that Σ_Z f_Z = 1. The spectral index α can vary (typically between α ~ 1 and α ≳ 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^{\prime }}_{\text{Z,max}}$ is determined by the competition between the acceleration timescale for a nucleus of charge $Z,{t^{\prime }}_{\text{acc}}={\eta }_{\text{acc}}^{-1}{E}^{\prime }/c Z e B^{\prime }$¹, and the energy loss timescales ${t^{\prime }}_{\text{loss}}=\text{min}({t^{\prime }}_{\text{dyn}},{t^{\prime }}_{\text{syn}},{t^{\prime }}_{\text{IC}},{t^{\prime }}_{\text{BH}},{t^{\prime }}_{p\gamma },\dots )$, where ${t^{\prime }}_{\text{dyn}}=R/\beta \mathrm{\Gamma }c$ is the dynamical timescale (see Appendix A for numerical estimates of ${E^{\prime }}_{\text{Z,max}}$). The factor η_acc ≤ 1 describes the efficiency of the acceleration process; for a maximally efficient acceleration, η_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_CR = ξ_CR L_bol, where we define the baryon loading fraction ξ_CR as the fraction of the bolometric luminosity that is injected into cosmic rays of energy E ≥ E_min ≡ 10¹⁶ 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_pk, peak energy ϵ_pk, and width â. The photon energy density then reads

$${{ϵ}^{\prime }}^2{n^{\prime }}_{{ϵ}^{\prime }}=\frac{L_{\text{pk}}}{4\pi {\mathrm{\Gamma }}^2{R}^{2}c}{({ϵ}^{\prime }/{{ϵ}^{\prime }}_{\text{pk}})}^{-\widehat{a}\log ({ϵ}^{\prime }/{{ϵ}^{\prime }}_{\text{pk}})}.$$(1)

The peak luminosity and peak energy set the maximum of the SED. The data can help to constrain the width of the logparabola 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 \text{d}{ϵ}^{\prime }{L}_{\text{pk}}/{ϵ}^{\prime }{({ϵ}^{\prime }/{{ϵ}^{\prime }}_{\text{pk}})}^{-â \log ({ϵ}^{\prime }/{{ϵ}^{\prime }}_{\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_jet,obs ~ L_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_bh > 5 × 10⁷ M_⊙), 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 â and the peak luminosity L_pk of the log-parabola, and to consider a typical peak energy ϵ_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}}^{\prime }$ (Eq. (A.4)). The magnetic field is inferred assuming equipartition between the radiative and magnetic energy densities: ${\xi }_B\int \text{d}{ϵ}^{\prime }{ϵ}^{\prime }{n^{\prime }}_{{ϵ}^{\prime }}=B^2/8\pi$ with ξ_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).

Fig. 1. Photon density for a log-parabola model with fixed peak luminosity Lpk = 1048 erg s−1.

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_dur ~ 10⁵ 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_dur ~ 10⁶ s, with a lower jet efficiency and a similar jet opening angle. For both states, we set a width â = 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 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 cosmic-ray spectra can reach ~40%. The impact of photonuclear models is also strong, whereas more difficult to quantify, especially regarding the channels involving α-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 â = 0.25, â = 0.07, and â = 0.03 for a fixed L_pk = 10⁴⁸ erg s⁻¹. 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 â = 0.25 and â = 0.03. We see that overall, photonuclear interactions dominate over a wide range of particle Lorentz factors γ, 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.

Fig. 2. Influence of the log-parabola width â on the MFPs and energy loss lengths in the comoving frame for carbon (left panel) and iron (right panel) nuclei of Lorentz factor γ. The peak luminosity is set to Lpk = 1048 erg s−1 and two examples of widths are presented: â = 0.25 (blue) and â = 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.

In Fig. 3, we compare the specific luminosities and photon densities for L_pk = 10⁴⁶ erg s⁻¹ and L_pk = 10⁴⁸ erg s⁻¹ for a given width â = 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.

Fig. 3. Influence of the log-parabola peak luminosity Lpk. Same as Fig. 2, but for fixed width â = 0.07 and varying the peak luminosity: Lpk = 1046 erg s−1 (blue) and Lpk = 1048 erg s−1 (red). Higher peak luminosity leads to shorter MFPs.

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 erg s⁻¹ and â = 0.07, and in its medium state SED characterized by L_pk = 10⁴⁶ erg s⁻¹ and â = 0.07. As is shown below, these two states are associated in our model with a black hole of mass M_bh = 7 × 10⁶ M_⊙. We consider an injection spectral index of α = 1.8 and an acceleration efficiency η_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 ⋍ 1.88 Gpc (z ⋍ 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 ⋍ 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 𝒮_IC = 5 × 10⁻⁴ TeV cm⁻² over the energy range 10 TeV–10 PeV, which corresponds to a detection limit s_IC ~ 10⁻¹¹ TeV cm⁻² s⁻¹ 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° < δ < 30° (thin lines), and for the declination range 30° < δ < 60° (thick lines) associated with the Swift event J1644+57.

Fig. 4. Cosmic-ray spectra for one source with pure iron injection with spectral index α = 1.8, photon field with â = 0.07, and acceleration efficiency η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 Mbh = 7 × 106 M⊙, with a corresponding SED in its high state with Lpk = 1047.5 erg s−1 and tdur = 105 s (left) and in its medium state with Lpk = 1046 erg s−1 and tdur = 106 s (right) for a source distance d L,1 = 1.88 Gpc. We assume here ξCR = 1.

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⁴⁶ erg s⁻¹, iron strongly interacts and produces many secondary particles with a large number of nucleons below E_cut/56 ~ 10¹⁹ eV. Despite these high interaction rates, nuclei can still survive and escape from the region with energies up to 10²⁰ eV. For a high state SED with L_pk = 10^47.5 erg s⁻¹, 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} ~ 10¹⁸ 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⁴⁶ erg s⁻¹, and photonuclear interactions for L_pk = 10^47.5 erg s⁻¹.

Figure 5 shows that a nearby medium state TDE at distance 20 Mpc with peak luminosity L_pk = 10⁴⁶ erg s⁻¹ 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 highenergy 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⁻⁶ yr⁻¹ and 3 × 10⁻⁵ yr⁻¹, 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 ξ_CR up to ~ a few 100 to remain consistent with this nondetection.

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 Lpk = 1047.5 erg s−1 (left) and a medium state SED with Lpk = 1046 erg s−1 (right). We consider two different distances d L,1 = 1.88 Gpc (dark blue) and d L,2 = 20 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° < δ < 30° (most favorable case, thin line) and 30° < δ < 60° (Swift J1644+57 case, thick line) (Aartsen et al. 2014).

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 ṅ(z,L) = ṅ₀Λ (L) f (z), where ṅ₀ is the total local event rate density, f (z) the TDE redshift distribution, and Λ(L) the TDE luminosity function. The luminosity function is given by a power-law,

$${\mathrm{\Lambda }}_{\text{TDE}}(L_{\gamma })\propto {\left(\frac{L_{\gamma ,\text{pk}}}{L_{\text{m},\text{pk}}}\right)}^{-{\alpha }_L},$$(2)

with L_m,pk = 10⁴⁸ erg s⁻¹ and ${\int }_{L_{\text{min}}}^{L_{\text{max}}}\text{d}{L}_{\gamma }{\mathrm{\Lambda }}_{\text{TDE}}\left(L_{\gamma }\right)=1$, with L_min = 10⁴⁵erg s⁻¹ and L_max = 10⁴⁹ erg s⁻¹, and α_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 Ṅ_TDE and the black hole comoving number density per luminosity dn_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 Ṅ_TDE depends on the galaxies considered. Following Wang & Merritt (2004), we consider a lower bound in the case of core galaxies of

$${\stackrel{}{Ṅ}}_{\text{TDE}}\approx {10}^{-5}{\text{yr}}^{-1},$$(3)

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

$${\stackrel{}{Ṅ}}_{\text{TDE}}\approx 7.1\times {10}^{-4}{\text{yr}}^{-1}{\left(\frac{\sigma }{70\text{km}\hspace{0.17em}{\text{s}}^{-1}}\right)}^{7/2}{M}_{\text{bh,6}}^{-1},$$(4)

where M_bh,6 = M_bh/10⁶M_⊙ and σ 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{kms}}^{-1}}\right).$$(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:

$${\stackrel{}{Ṅ}}_{\text{TDE}}\approx 3\times {10}^{-4}{\text{yr}}^{-1}{M}_{\text{bh},6}^{-0.2}·$$(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_⋆ at distances d_⋆ ≤ r_t = R_⋆(M_bh/M_⋆)^1/3 (Hills 1975). Following Krolik & Piran (2012), we can estimate the tidal disruption radius $r_{\text{t}}\simeq 10R_{\text{s}}{M}_{\star ,\odot }^{2/3-\xi }{M}_{\text{bh,6}}^{-2/3}{\left[(k_{\star }/f_{\star })/0.02\right]}^{1/6}$, where R_s = 2GM_bh/c ² is the Schwarzschild radius, M_⋆,⊙ is the mass of the star in solar units, M_bh,6 = M_bh/10⁶ M_⊙,k_⋆ is related to the star’s radial density profile, and f_⋆ 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 ξ ⋍ 0.2 for 0.1 M_⊙ < M_⋆ ≤ M_⊙ or ξ ⋍ 0.4 for M_⊙ < M_⋆ < 10 M_⊙ (Kippenhahn & Weigert 1994). Moreover, we consider here and in what follows fully radiative stars, thus k_⋆/ f_⋆ = 0.02 (Phinney 1989). For white dwarfs, typically $R_{\star }\sim M_{\star }^{-1/3}$ with 0.5 M_⊙ < M_⋆ ≤ 0.7 M_⊙. Their tidal disruption radii are smaller due to the smaller dimensions of these objects; an approximate formula gives r_t ⋍ 7.4 × 10⁹ (M_bh,3.3/ρ_⋆,7)^1/3 (Luminet & Pichon 1989), where ρ_⋆ 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 ≳ R_s. For a Schwarzschild black hole, this leads to $M_{\text{bh}}\lesssim 4\times {10}^7{M}_{\odot }{M}_{\star ,\odot }^{1-3\xi /2}{\left[(k_{\star }/f_{\star })/0.02\right]}^{1/4}$, which ranges from M_bh ≲ 10⁷ M_⊙ to M_bh ≲ 10⁸ M_⊙ for 0.1 M_⊙ < M_⋆ ≤ 10 M_⊙. 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 ~ 7 × 10⁸ M_⊙ (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\hspace{0.17em}\text{s}\hspace{0.17em}{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,$$(7)

where a_min is the minimum semi-major axis. The parameter ξ comes from the main sequence mass–radius relation $R_{\star }\approx R_{\odot }{M}_{\star ,\odot }^{(1-\xi )}$ and β ≲ 1 is the penetration factor. We obtain the following accretion rate:

$$\dot{M}\approx 20\hspace{0.17em}{M}_{\odot }\hspace{0.17em}{\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}.$$(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 ≲ M_bh,jet, where M_bh,jet = 4 × 10⁸ M_⊙ $(\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}\hspace{0.17em}{\text{s}}^{-1}\frac{f(a)}{{\beta }_h{\alpha }_{\text{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},$$(9)

where α_s is the ratio of inflow speed to orbital speed of the disk, and β_h the ratio of the midplane total pressure near the ISCO to the magnetic pressure in the black hole’s stretched horizon, such that α_s β_h ~ 0.1–1; the function f (a) ≈ a ² 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); ṁ = Ṁc ²/L_Edd is the normalized accretion rate in the outer disk (with L_Edd the Eddington luminosity); ṁ₀ is the peak normalized accretion rate; and q is the fraction of ṁ 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, α_s β_h = 1, k_⋆/ f_⋆ = 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α 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 ≪ 1 our jetted TDE rates would be significantly decreased, but such a choice seems hard to reconcile with iron-Kα observations (see discussion in Sesana et al. 2014).

The total energy release per TDE is given by

$${Ɛ}_{\text{jet}}\approx L_{\text{jet}}{P}_{\text{orb}}\simeq 4\times {10}^{52}\hspace{0.17em}\text{erg}\hspace{0.17em}{M}_{\star ,\odot }\frac{f(a)}{{\beta }_h{\alpha }_{\text{s}}}q(\dot{m}/{\dot{m}}_0),$$(10)

which should be less than Ṁc ². We note that a jet luminosity L_jet = η_jet Ṁc ² and η_jet ~ 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 η_jet and for the beaming factor ℬ = min(4π/ΔΩ, 2Γ²), where ΔΩ is the solid angle occupied by the jet (Krolik & Piran 2012). For a two-sided jet with a jet opening angle θ_jet, we have ℬ = min(1/(1 − cos θ_jet), 2Γ²). Therefore, L_bol = L_jet,obs = 2 η_jet,−2ℬL_jet for θ_jet ~ 5°, Γ = 10, and η_jet = 0.01.

With the observed local rate density, ṅ_tde,0 = 0.03 Gpc⁻³ yr⁻¹, the luminosity density is estimated to be

$$\begin{array}{ll}{Q}_{\text{TDEjet}}\hfill & \approx {\eta }_{\text{jet}}\hspace{0.17em}ℬ{Ɛ}_{\text{jet}}\hspace{0.17em}{\dot{n}}_{\text{tde},0}\hfill \\ \hfill & \simeq 2\times {10}^{42}\text{erg}\hspace{0.17em}{\text{Mpc}}^{-3}\hspace{0.17em}{\text{yr}}^{-1}\hspace{0.17em}{\eta }_{\text{jet},-2}\hfill \\ \hfill & \times M_{\star ,\odot }\frac{f(a)}{{\beta }_h{\alpha }_{\text{s}}}q(\dot{m}/{\dot{m}}_0)\frac{{\dot{n}}_{\text{tde},0}}{0.03\hspace{0.17em}{\text{Gpc}}^{-3}\hspace{0.17em}{\text{yr}}^{-1}}\cdot \hfill \end{array}$$(11)

As discussed in Sect. 5.1, this implies a baryon loading factor η_CR of order 100 to account for the observed flux of UHECR, which is ~10⁴⁴ erg Mpc⁻³ yr⁻¹.

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), gravitationalwave 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 ṅ_TDE(z,L) by combining the TDE rate per galaxy Ṅ_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 ≳ 7 × 10⁶ M_⊙ (Seifina et al. 2017) and reaches a bolometric luminosity L_bol ≳ 10⁴⁸ erg s⁻¹ in the high state. Therefore, we have θ_jet = 5° and η_jet = 0.35 in the high state and θ_jet = 5° and η_jet = 0.01 in the medium state. For instance, a black hole of mass M_bh ≳ 7 × 10⁶ M_⊙ is associated with L_bol ⋍ 2 × 10⁴⁸ erg s⁻¹ in the high state and L_bol ⋍ 7 × 10⁴⁶ 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⁸ M_⊙, 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 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_{\text{jet}}\propto M_{\text{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 ≥ 10^45.5 erg s ⁻¹, which dominate in the production of cosmicray and neutrino fluxes in our framework.

Fig. 6. Comoving TDE luminosity density in their medium state as a function of redshift, as derived in our model, for ṄTDE = 10−5 yr−1. The different luminosities correspond to different black hole masses.

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

$${\mathrm{\Phi }}_{\text{CR}}(E_{\text{CR}})=\frac{c}{4\pi H_0}{\displaystyle \underset{z_{\text{min}}}{\overset{z_{\text{max}}}{\int }}{\displaystyle \underset{L_{\text{min}}}{\overset{L_{\text{max}}}{\int }}\text{d}z\hspace{0.17em}\text{d}L\frac{f_{\text{s}}\hspace{0.17em}{\xi }_{\text{CR}}\hspace{0.17em}{\dot{N}}_{\text{TDE}}\hspace{0.17em}\text{d}{n}_{\text{bh}}(z,L)/\text{d}L}{\sqrt{{\mathrm{\Omega }}_{\text{M}}{(1+z)}^3+{\mathrm{\Omega }}_{\mathrm{\Lambda }}}}}}\times F_{\text{CR},\text{s},\mathrm{p}}^c(E_{\text{CR}}^c,z,L)t_{\text{dur}}^c,$$(12)

where Ω_M = 0.3 and Ω_L = 0.7 are our fiducial cosmological parameters, H₀ = 70 km s⁻¹ Mpc⁻¹ 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 enegy $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 Ṅ_TDE is the rate of TDEs per galaxy and dn_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

$${\mathrm{\Phi }}_{\nu }(E_{\nu })=\frac{c}{4\pi H_0}{\displaystyle \underset{z_{\text{min}}}{\overset{z_{\text{max}}}{\int }}{\displaystyle \underset{L_{\text{min}}}{\overset{L_{\text{max}}}{\int }}\text{d}z\hspace{0.17em}\text{d}L\frac{f_{\text{s}}\hspace{0.17em}{\xi }_{\text{CR}}{\dot{N}}_{\text{TDE}}\hspace{0.17em}\text{d}{n}_{\text{bh}}(z,L)/\text{d}L}{\sqrt{{\mathrm{\Omega }}_{\text{M}}{(1+z)}^3+{\mathrm{\Omega }}_{\mathrm{\Lambda }}}}\times F_{\nu ,\text{s}}^c(E_{\nu }^c,L)t_{\text{dur}}^c,}}$$(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 ~ 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 ≲ 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 α = 1.5, an acceleration efficiency η_acc = 0.1, a TDE source evolution, and a fraction ξ_CR f_s = 1.4% of the maximum local event rate density ṅ_tde,0 = 155 Gpc⁻³ yr⁻¹, 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.

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 α = 1.5, ξCR fs = 1.4 × 10−2 and source evolution derived in this work, with maximum bolometric luminosity in the population Lmax = 6.8 × 1046 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.

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, 〈X_max〉 and σ(X_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, 〈X_max〉 and σ(X_max) (Aab et al. 2014), which are in good agreement with our model. Only statistical uncertainties are displayed; systematic uncertainties are at most ±10 g cm⁻² for 〈X_max〉 and ±2 g cm⁻² for σ(X_max). Our results are compatible with a light composition at 10^18.5 eV, shifting toward a heavier composition for increasing energy.

Fig. 8. Mean and standard deviation of Xmax 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).

With increased cutoff bolometric luminosities L_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 α = 1 is required for L_max = 10⁴⁷ erg s⁻¹). 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 η_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 ≲ 10^46.5 erg s⁻¹, 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¹⁹ eV and no survival of heavy elements at the highest energies. However, because of its limited duration (t_dur ~ 10⁵ 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 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_pk ~ 10⁴⁷–10⁴⁸ erg s⁻¹, and a local event rate density is of ṅ_tde,0 ~ 0.03 Gpc⁻³ yr⁻¹ (e.g., Farrar & Piran 2014). On the other hand, normal TDEs are less luminous and are characterized by a higher local event rate density ṅ_tde,0 = 10² Gpc⁻³ yr⁻¹ (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 expect reaches ṅ_tde,0 ~ 2 × 10² Gpc⁻³ yr⁻¹ for core galaxies and ṅ_tde,0 ~ 3 × 10³ Gpc⁻³ yr⁻¹ for power-law galaxies.

The fraction needed to fit the UHECR spectrum of the Auger observations, ξ_CR f_s = 1.4 × 10⁻², can account for example for low UHECR injection rates ξ_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 ṅ_tde,0 ~ 0.03 Gpc⁻³ yr⁻¹ for the jetted events pointed toward the observer, a baryon loading of ξ_CR ~ 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¹⁹ 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 ξ_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¹⁶ 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 highenergy cutoff reduces the flux at higher energies, and therefore it lies just below the GRAND sensitivity limit.

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 POEMMA (Neronov et al. 2017).

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².

The bulk Lorentz factor Γ and the opening angle of the jet θ_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 Γ⁻⁴, 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 Γ = 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 θ_jet ~ 5° 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 ≳100–1000 Gpc⁻³ yr⁻¹, which is significantly higher than the rate of jetted TDEs pointing toward us inferred from observations of ~0.03 Gpc⁻³ yr⁻¹, and higher than rates with dimmer luminosities also constrained by X-ray observations. In addition, large baryon loadings with ≳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 ≳10⁵⁴ 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 ~ 2 × 10⁴⁸ erg s⁻¹ 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¹⁷ 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 multimessenger signal with TDE photons associated with the emission of neutrinos from a single source appears to be the way to validate this scenario.

Acknowledgments

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 Sklodowska-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

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 Γ.

A.1. Maximum injection energy

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

$$\begin{array}{ll}{E^{\prime }}_{Z,\text{up},\text{dyn}}\hfill & ~c Z e B^{\prime }{(1+z)}^{-1}\mathrm{\Gamma }{t}_{\text{var}}{\eta }_{\text{acc}},\hfill \\ \hfill & \simeq \hspace{0.17em}6.3\times {10}^{17}\text{eV} Z_1{B^{\prime }}_{2.8}{\mathrm{\Gamma }}_{10}{t}_{\text{var},2}{\eta }_{\text{acc},-1},\hfill \\ \hfill & \simeq \hspace{0.17em}1.6\times {10}^{19}\text{eV} Z_{26}{B^{\prime }}_{2.8}{\mathrm{\Gamma }}_{10}{t}_{\text{var},2}{\eta }_{\text{acc},-1}.\hfill \end{array}$$(A.1)

The competition between the acceleration timescale ($t_{\text{acc}}^{\prime }$) and the synchrotron timescale ($t_{\text{syn}}^{\prime }$) gives

$$\begin{array}{ll}{E^{\prime }}_{Z,\text{up},\text{syn}}\hfill & ~\hspace{0.17em}{\left[\frac{6\pi {(m_{\text{u}}{c}^2)}^{2}e}{{(m_e/m_{\text{u}})}^2{\sigma }_{\text{T}}}\right]}^{\frac{1}{2}}{A}^2{Z}^{-3/2}{B^{\prime }}^{-1/2}{\eta }_{\text{acc}}^{1/2},\hfill \\ \hfill & \simeq 2.4\times {10}^{18}\text{eV}{A}_1^2{Z}_1^{-3/2}{B^{\prime }}_{2.8}^{-1/2}{\eta }_{\text{acc},-1}^{1/2},\hfill \\ \hfill & \simeq 5.7\times {10}^{19}\text{eV}{A}_{56}^2{Z}_{26}^{-3/2}{B^{\prime }}_{2.8}^{-1/2}{\eta }_{\text{acc},-1}^{1/2}.\hfill \end{array}$$(A.2)

Here, the mean magnetic field $B^{\prime }=\sqrt{8\pi \int \text{d}{ϵ}^{\prime }{ϵ}^{\prime }{n^{\prime }}_{{ϵ}^{\prime }}}\simeq {10}^{2.85}\text{G}$ is obtained for a log-parabola SED with peak luminosity L _pk = 10⁴⁶ erg s⁻¹ and a width â = 0.07.

The upper bound given by the competition between the acceleration timescale ($t_{\text{acc}}^{\prime }$) and the photohadronic timescale ($t_{\mathrm{p}\gamma }^{\prime }$) is computed numerically. For the parameters considered above, we obtain for protons $E_{\text{p,up,}\mathrm{p}\gamma }^{\prime }\simeq 2.5\times {10}^{18}\text{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_{\text{syn}}^{\prime }$) and escape ($t_{\text{esc}}^{\prime }$) for protons,

$$\frac{6\pi {(m_{\mathrm{p}}{c}^2)}^2}{{(m_{\text{e}}/m_{\mathrm{p}})}^2{\sigma }_{\text{T}}c{E^{\prime }}_{\mathrm{p}}{B^{\prime }}^2}=\mathrm{\Gamma }{t}_{\text{var}},$$(A.3)

with $E_{\mathrm{p}}^{\prime }$ the proton energy in the comoving frame, σ _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{array}{ll}{E^{\prime }}_{\mathrm{p},\text{max}}\hfill & =\frac{6\pi {(m_{\mathrm{p}}{c}^2)}^2}{{(m_{\text{e}}/m_{\mathrm{p}})}^2{\sigma }_{\text{T}}c{B^{\prime }}^2\mathrm{\Gamma }{t}_{\text{var}}},\hfill \\ \hfill & \simeq 9.0\times {10}^{18}\text{eV}{B^{\prime }}_{2.8}^{-2}{\mathrm{\Gamma }}_1^{-1}{t}_{\text{var},2}^{-1}.\hfill \end{array}$$(A.4)

The competition between synchrotron losses ($t_{\text{syn}}^{\prime }$) and pion production ($t_{\mathrm{p}\gamma }^{\prime }$) for a transient event characterized by a hard spectrum gives

$$\frac{6\pi {(m_{\mathrm{p}}{c}^2)}^2}{{(m_{\text{e}}/m_{\mathrm{p}})}^2{\sigma }_{\text{T}}c{E^{\prime }}_{\mathrm{p}}{B^{\prime }}^2}=\frac{4\pi {\mathrm{\Gamma }}^5{c}^2{t}_{\text{var}}^2{ϵ}_{\text{pk}}}{\langle {\sigma }_{\mathrm{p}\gamma }{\kappa }_{\mathrm{p}\gamma }\rangle L_{\text{pk}}},$$(A.5)

with σ_pγ and _pγ the photopion cross section and inelasticity, ϵ_pk the peak energy, and L_pk the peak luminosity. For ϵ_pk = 70 keV, we obtain a high-energy cutoff:

$$\begin{array}{ll}{E^{\prime }}_{\mathrm{p},\text{max}}\hfill & =\frac{3{(m_{\mathrm{p}}{c}^2)}^2\langle {\sigma }_{\mathrm{p}\gamma }{\kappa }_{\mathrm{p}\gamma }\rangle L_{\text{pk}}}{2{(m_e/m_p)}^2{\sigma }_{\text{T}}{c}^3{B^{\prime }}^2{\mathrm{\Gamma }}^5{t}_{\text{var}}^2{ϵ}_{\text{pk}}},\hfill \\ \hfill & \simeq 5.0\times {10}^{21}\text{eV}{B^{\prime }}_{2.8}^{-2}{\mathrm{\Gamma }}_1^{-5}{L}_{\text{pk},46}{t}_{\text{var},2}^{-2}.\hfill \end{array}$$(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 }_{\text{N}}=E_{\text{N}}^{\prime }/A{m}_{\text{u}}{c}^2$, where m _u is the atomic mass unit, for the competition between synchrotron losses ($t_{\text{syn}}^{\prime }$) and escape ($t_{\text{esc}}^{\prime }$) we obtain

$$\begin{array}{ll}{E^{\prime }}_{\text{N},\text{max}}\hfill & =\frac{6\pi {(m_{\text{u}}{c}^2)}^2{(A/Z)}^4}{{(m_e/m_{\text{u}})}^2{\sigma }_{\text{T}}c{B^{\prime }}^2\mathrm{\Gamma }{t}_{\text{var}}},\hfill \\ \hfill & \simeq 1.4\times {10}^{20}\text{eV}{B^{\prime }}_{2.8}^{-2}{\mathrm{\Gamma }}_1^{-1}{t}_{\text{var},2}^{-1}{A}_{56}^4{Z}_{26}^{-4},\hfill \end{array}$$(A.7)

where A ₅₆ = A/56 and Z ₂₆ = 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_{\nu }^c=(1+z)E_{\nu }$ and the bolometric luminosity of the source: $N_{\nu ,\text{s}}(E_{\nu }^c,L)$.

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

$$\frac{\text{d}V(z)}{\text{d}z}=\frac{c}{H_0}\frac{4\pi D_{\text{c}}^2}{\sqrt{{\mathrm{\Omega }}_{\text{M}}{(1+z)}^3+{\mathrm{\Omega }}_{\mathrm{\Lambda }}}},$$(B.1)

where Ω_M = 0.3 and Ω_L = 0.7 are our fiducial cosmological parameters, and H₀ = 70 km s⁻¹ Mpc⁻¹ is the Hubble constant. The comoving event rate density dṅ(z,L)/dL = Ṅ_TDE dn_bh(z,L)/dL depends on the TDE rate per galaxy Ṅ_TDE and the black hole comoving density per luminosity dn_bh(z,L)/dL. Therefore, the diffuse neutrino flux is given by

$$\begin{array}{ll}\frac{\text{d}{N}_{\nu }}{\text{d}{E}_{\nu }}(E_{\nu })\hfill & =\frac{1}{4\pi }{\displaystyle \underset{z_{\text{min}}}{\overset{z_{\text{max}}}{\int }}{\displaystyle \underset{L_{\text{min}}}{\overset{L_{\text{max}}}{\int }}[\text{d}z\text{d}L{f}_{\text{s}}{\xi }_{\text{CR}}\frac{\text{d}\dot{n}(z,L)}{\text{d}L}}} \hfill \\ \hfill & \times \frac{1}{1+z}\frac{\text{d}V(z)}{\text{d}z}\frac{1}{4\pi D_{\text{c}}^2}\frac{\text{d}{N}_{\nu ,\text{s}}(E_{\nu }^c,L)}{\text{d}{E}_{\nu }^c}\frac{\text{d}{E}_{\nu }^c}{\text{d}{E}_{\nu }}].\hfill \end{array}$$(B.2)

As $E_{\nu }^c=(1+z)E_{\nu }$, we obtain a diffuse neutrino flux

$${\mathrm{\Phi }}_{\nu }(E_{\nu })=\frac{c}{4\pi H_0}{\displaystyle \underset{z_{\text{min}}}{\overset{z_{\text{max}}}{\int }}{\displaystyle \underset{L_{\text{min}}}{\overset{L_{\text{max}}}{\int }}\text{d}z \text{d}L\frac{f_{\text{s}}{\xi }_{\text{CR}}\text{d}\dot{n}(z,L)/\text{d}L}{\sqrt{{\mathrm{\Omega }}_{\text{M}}{(1+z)}^3+{\mathrm{\Omega }}_{\mathrm{\Lambda }}}}}}\times F_{\nu ,\text{s}}^c(E_{\nu }^c,L)t_{\text{dur}}^c,$$(B.3)

where Φ_ν(E_ν) = d²N_ν(E_ν)/dE_νdt is the diffuse neutrino flux per unit time (observer frame), $F_{\nu ,\text{s}}^c(E_{\nu }^c,L)={\text{d}}^2{N}_{\nu ,\text{s}}(E_{\nu }^c,L)/\text{d}{E}_{\nu }^c{\text{d}t}^c$ 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

$$\begin{array}{ll}{\mathrm{\Phi }}_{\text{CR}}(E_{\text{CR}})\hfill & =\frac{c}{4\pi H_0}{\displaystyle \underset{z_{\text{min}}}{\overset{z_{\text{max}}}{\int }}{\displaystyle \underset{L_{\text{min}}}{\overset{L_{\text{max}}}{\int }}\text{d}z\text{d}L\frac{f_{\text{s}}{\xi }_{\text{CR}}\text{d}\dot{n}(z,L)/\text{d}L}{\sqrt{{\mathrm{\Omega }}_{\text{M}}{(1+z)}^3+{\mathrm{\Omega }}_{\mathrm{\Lambda }}}}}} \hfill \\ \hfill & \times F_{\text{CR},\text{s},\mathrm{p}}^c(E_{\text{CR}}^c,z,L)t_{\text{dur}}^c,\hfill \end{array}$$(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 semianalytic 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 _⊙ (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⁵ M _⊙) 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⁸ M _⊙). 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).

Fig. C.1. Predictions of our model in two different scenarios, light-seed (“Pop III”) and heavy-seed (“Qc = 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).

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.

All Tables

All Figures

	Fig. 1. Photon density for a log-parabola model with fixed peak luminosity Lpk = 1048 erg s−1.
In the text

Fig. 2. Influence of the log-parabola width â on the MFPs and energy loss lengths in the comoving frame for carbon (left panel) and iron (right panel) nuclei of Lorentz factor γ. The peak luminosity is set to Lpk = 1048 erg s−1 and two examples of widths are presented: â = 0.25 (blue) and â = 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.

In the text

	Fig. 3. Influence of the log-parabola peak luminosity Lpk. Same as Fig. 2, but for fixed width â = 0.07 and varying the peak luminosity: Lpk = 1046 erg s−1 (blue) and Lpk = 1048 erg s−1 (red). Higher peak luminosity leads to shorter MFPs.
In the text

Fig. 4. Cosmic-ray spectra for one source with pure iron injection with spectral index α = 1.8, photon field with â = 0.07, and acceleration efficiency η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 Mbh = 7 × 106 M⊙, with a corresponding SED in its high state with Lpk = 1047.5 erg s−1 and tdur = 105 s (left) and in its medium state with Lpk = 1046 erg s−1 and tdur = 106 s (right) for a source distance d L,1 = 1.88 Gpc. We assume here ξCR = 1.

In the text

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 Lpk = 1047.5 erg s−1 (left) and a medium state SED with Lpk = 1046 erg s−1 (right). We consider two different distances d L,1 = 1.88 Gpc (dark blue) and d L,2 = 20 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° < δ < 30° (most favorable case, thin line) and 30° < δ < 60° (Swift J1644+57 case, thick line) (Aartsen et al. 2014).

In the text

	Fig. 6. Comoving TDE luminosity density in their medium state as a function of redshift, as derived in our model, for ṄTDE = 10−5 yr−1. The different luminosities correspond to different black hole masses.
In the text

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 α = 1.5, ξCR fs = 1.4 × 10−2 and source evolution derived in this work, with maximum bolometric luminosity in the population Lmax = 6.8 × 1046 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.

In the text

	Fig. 8. Mean and standard deviation of Xmax 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).
In the text

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 POEMMA (Neronov et al. 2017).

In the text

Fig. C.1. Predictions of our model in two different scenarios, light-seed (“Pop III”) and heavy-seed (“Qc = 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).

In the text

	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.
In the text