High-energy cosmic ray nuclei from tidal disruption events: Origin, survival, and implications

TL;DR

This paper investigates whether tidal disruption events can produce and sustain ultrahigh-energy cosmic ray nuclei, analyzing their survival through various shock regions and implications for observed cosmic ray composition.

Contribution

It provides a detailed analysis of UHECR nuclei survival in TDE jets and explores the viability of different stellar progenitors as sources of observed cosmic rays.

Findings

UHECR nuclei can survive in external shocks and winds but are disintegrated in internal shocks.

Main sequence star and CO white dwarf disruptions do not match UHECR observations.

Oxygen-neon-magnesium white dwarfs could explain composition but are rare.

Abstract

Tidal disruption events (TDEs) by supermassive or intermediate mass black holes have been suggested as candidate sources of ultrahigh-energy cosmic rays (UHECRs) and high-energy neutrinos. Motivated by the recent measurements from the Pierre Auger Observatory, which indicates a metal-rich cosmic-ray composition at ultrahigh energies, we investigate the fate of UHECR nuclei loaded in TDE jets. First, we consider the production and survival of UHECR nuclei at internal shocks, external forward and reverse shocks, and nonrelativistic winds. Based on the observations of Swift J1644+57, we show that the UHECRs can survive for external reverse and forward shocks, and disk winds. On the other hand, UHECR nuclei are significantly disintegrated in internal shocks, although they could survive for low-luminosity TDE jets. Assuming that UHECR nuclei can survive, we consider implications of different composition models of TDEs. We find that the tidal disruption of main sequence stars or carbon-oxygen white dwarfs does not successfully reproduce UHECR observations, namely the observed composition or spectrum. The observed mean depth of the shower maximum and its deviation could be explained by oxygen-neon-magnesium white dwarfs, but they may be too rare to be the sources of UHECRs.

Equations

$$t_{A\gamma }^{-1}(E_A)=\frac{c}{2{\gamma }_A^2}{\int }_{{\bar{\epsilon }}_{\mathrm{th}}}^{\infty }d\bar{\epsilon }{\sigma }_{A\gamma }(\bar{\epsilon }){\kappa }_{A\gamma }(\bar{\epsilon })\bar{\epsilon }{\int }_{\bar{\epsilon }/2{\gamma }_A}^{\infty }d\epsilon \frac{1}{{\epsilon }^2}\frac{dn}{d\epsilon },$$

$${\nu }_{m,\mathrm{ob}}^r\simeq 8$$

$${\nu }_{a,\mathrm{ob}}^r\simeq 5.6$$

$$L_{ϵ,\max }^r$$

$$\frac{dn^{r}}{d\varepsilon}=\frac{L_{\varepsilon,\rm max}^{r}}{4\pi r_{\times}^{2}c\varepsilon_{m}}\left\{\begin{array}[]{ll}(\frac{\varepsilon_{a}}{\varepsilon_{m}})^{-2/3}(\frac{\varepsilon}{\varepsilon_{a}})^{1}&(\varepsilon<\varepsilon_{a})\\ (\frac{\varepsilon}{\varepsilon_{m}})^{-2/3}&(\varepsilon_{a}\leq\varepsilon<\varepsilon_{m})\\ (\frac{\varepsilon}{\varepsilon_{m}})^{-(p-1)/2-1}&(\varepsilon_{m}\leq\varepsilon<\varepsilon_{c})\\ (\frac{\varepsilon_{c}}{\varepsilon_{m}})^{-(p-1)/2-1}(\frac{\varepsilon}{\varepsilon_{c}})^{-p/2-1}&(\varepsilon\geq\varepsilon_{c})\end{array}\right.,$$

$$E_{A,\max }$$

$${\nu }_{m,\mathrm{ob}}\simeq 1.3\times 10^{11}{\mathcal{E}}_{\mathrm{iso},53.3}^{1/2}{ϵ}_{B,-1.8}^{1/2}{ϵ}_{e,-1.5}^2{f}_{e,-1}^{-2}{t}_6^{-3/2}\mathrm{Hz},$$

$${\nu }_{a,\mathrm{ob}}\simeq 3.1\times 10^9{\mathcal{E}}_{\mathrm{iso},53.3}^{1/5}{ϵ}_{B,-1.8}^{1/5}{ϵ}_{e,-1.5}^{-1}{f}_{e,-1}^{3/5}{\varrho }^{3/5}Hz.$$

$$L_{ϵ,\max }$$

$$\frac{dn^{f}}{d\varepsilon}=\frac{L_{\varepsilon,\rm max}}{4\pi r^{2}c\varepsilon_{m}}\left\{\begin{array}[]{ll}(\frac{\varepsilon_{a}}{\varepsilon_{m}})^{-2/3}(\frac{\varepsilon}{\varepsilon_{a}})^{1}&(\varepsilon<\varepsilon_{a})\\ (\frac{\varepsilon}{\varepsilon_{m}})^{-2/3}&(\varepsilon_{a}\leq\varepsilon<\varepsilon_{m})\\ (\frac{\varepsilon}{\varepsilon_{m}})^{-(p-1)/2-1}&(\varepsilon_{m}\leq\varepsilon<\varepsilon_{c})\\ (\frac{\varepsilon_{c}}{\varepsilon_{m}})^{-(p-1)/2-1}(\frac{\varepsilon}{\varepsilon_{c}})^{-p/2-1}&(\varepsilon\geq\varepsilon_{c})\end{array}\right.,$$

$$E_{A,\max }$$

$${\nu }_m\simeq 1.5\times 10^6{\mathcal{E}}_{k,48}{ϵ}_{B,-1}^{1/2}{ϵ}_{e,-1}^2{f}_{e,-1}^{-2}{\varrho }^{-1/2}(t/t_{\mathrm{dec}}{)}^{-3}\mathrm{Hz},$$

$${\nu }_a\simeq 4.7\times 10^7{\mathcal{E}}_{k,48}^{\frac{p}{p+4}}{ϵ}_{B,-1}^{\frac{p/2+1}{p+4}}{ϵ}_{e,-1}^{\frac{2p-2}{p+4}}{f}_{e,-1}^{\frac{2-2p}{p+4}}{\varrho }^{\frac{3-p/2}{p+4}}(t/t_{\mathrm{dec}}{)}^{\frac{2-3p}{p+4}}\mathrm{Hz},$$

$$\frac{dn}{d\varepsilon}=\frac{L_{\varepsilon,{\rm max}}}{4\pi R_{s}^{2}c\varepsilon_{a}}\left\{\begin{array}[]{ll}(\frac{\varepsilon_{m}}{\varepsilon_{a}})^{5/2}(\frac{\varepsilon}{\varepsilon_{a}})^{-1}(\frac{\varepsilon}{\varepsilon_{m}})^{2}&(\varepsilon<\varepsilon_{m})\\ (\frac{\varepsilon}{\varepsilon_{a}})^{3/2}&(\varepsilon_{m}\leq\varepsilon<\varepsilon_{a})\\ (\frac{\varepsilon}{\varepsilon_{a}})^{-(p-1)/2-1}&(\varepsilon_{a}\leq\varepsilon<\varepsilon_{c})\\ (\frac{\varepsilon_{c}}{\varepsilon_{a}})^{-(p-1)/2-1}(\frac{\varepsilon}{\varepsilon_{c}})^{-p/2-1}&(\varepsilon\geq\varepsilon_{c})\end{array}\right..$$

$$\frac{d{N}_{A^{\prime }}}{d{E}^{\prime }}=f_{A^{\prime }}{N}_0{\left(\frac{E^{\prime }}{Z{E}_0}\right)}^{-s_{\mathrm{esc}}}\exp \left(-\frac{E^{\prime }}{Z{E}_{p,\max }^{\prime }}\right),$$

$${\kappa }_{A\gamma }(\bar{\epsilon })\equiv \frac{\Delta E}{E}=\frac{\Delta N}{N},$$

$${\bar{\kappa }}_{A\gamma }=\frac{1}{{\sigma }_{\mathrm{tot}}}\sum _i{\sigma }_i\frac{[\Delta N{]}_i}{N},$$

$${\Phi }_A=\sum _{A^{\prime }}{\Phi }_{A{A}^{\prime }}=\sum _{A^{\prime }}\frac{c}{4\pi }\frac{d{n}_{A{A}^{\prime }}}{dE},$$

$$d{n}_{A{A}^{\prime }}(E)$$

$$\frac{dt}{dz}=-\frac{1}{H_0(1+z)}\frac{1}{\sqrt{{\Omega }_{\Lambda }+{\Omega }_k(1+z{)}^2+{\Omega }_m(1+z{)}^3}},$$

$$f_{\mathrm{TDE}}(z)={\left[(1+z{)}^{0.2\eta }+{\left(\frac{1+z}{1.43}\right)}^{-3.2\eta }+{\left(\frac{1+z}{2.66}\right)}^{-7.0\eta }\right]}^{\frac{1}{\eta }},$$

$$Q_{\mathrm{CR}}$$

$${\mathcal{E}}_{\mathrm{CR}}^{\mathrm{iso}}=\sum _{A^{\prime }}{\int }_{E_{\min }^{\prime }}^{Z{E}_{p,\max }^{\prime }}d{E}^{\prime }{E}^{\prime }\frac{d{N}_{A^{\prime }}}{d{E}^{\prime }}.$$

$${\rho }_{\mathrm{WD}-\mathrm{IMBH}}\sim 0.06\left(\frac{f_{\mathrm{IMBH}}}{1}\right)\left(\frac{f_{\mathrm{jet}}}{0.1}\right)\left(\frac{f_b}{0.01}\right)\mathrm{Gp}{\mathrm{c}}^{-3}\mathrm{y}{\mathrm{r}}^{-1},$$

Full text

I Introduction

Cosmic rays with energy larger than $10^{18}~{}\rm eV$ are referred to as ultrahigh-energy cosmic rays (UHECRs), and their origin is still largely unknown Hillas (1984); Nagano and Watson (2000); Kotera and Olinto (2011). The observed spectrum of UHECRs has a cutoff at energy around $\sim 4\times 10^{19}\rm~{}eV$ Abraham et al. (2008); Abbasi et al. (2008); Abraham et al. (2010). The flux suppression is consistent with the prediction of the Greisen-Zatsepin-Kuzmin (GZK) effect due to the interaction between UHECRs and cosmic microwave background (CMB) photons Greisen (1966); Zatsepin and Kuzmin (1966). A key clue to the origin of UHECRs is their composition. The primary mass of UHECRs can be inferred from the distributions of the depth of the shower maximum, $X_{\rm max}$. The data from the Telescope Array (TA) and Auger are consistent within systematic and statistical uncertainties Abbasi et al. (2016). The analysis from Auger suggests that the composition of UHECRs is dominated by light nuclei at energy around $10^{18.3}\rm~{}eV$ and becomes gradually heavier with increasing energy up to $10^{19.6}\rm~{}eV$ Aab et al. (2014a, b). The distributions of $X_{\rm max}$ are difficult to explain with a mixture of protons and iron nuclei over the whole energy range if up-to-date hadronic interaction models are correct. Rather, the best fit is reached by including a fraction of intermediate mass nuclei Aab et al. (2014b).

The interpretation of the UHECR composition is still under intense debate, and we avoid discussion that depends on their details. However, considering the larger detection area and lower sampling bias of Auger, it is reasonable to assume that UHECR composition gets gradually heavier at the highest energies. There are not so many candidate sources which satisfy the Hillas criterion to accelerate cosmic rays (CRs) to ultrahigh energies Hillas (1984), and the origin of heavy nuclei has now become an interesting question. In relativistic jets of active galactic nuclei (AGN) (e.g., Norman et al. (1995); Pe’er et al. (2009); Dermer and Razzaque (2010); Murase et al. (2012); Fang and Murase (2017) and references therein), heavy nuclei would be supplied by an accretion disk. For structure formation shocks in galaxy clusters (GCs) Norman et al. (1995); Kang et al. (1996); Inoue et al. (2007), heavy nuclei can be provided from the intracluster medium. However, for these objects, we typically expect a composition similar to the solar composition, unless reacceleration is invoked Kimura et al. (2017); Caprioli (2015). Another possibility is that heavy nuclei can be synthesized in the stellar interiors or outflows from the deaths of massive stars. This scenario is appropriate for gamma-ray bursts (GRBs) Murase et al. (2008); Wang et al. (2008); Horiuchi et al. (2012); Globus et al. (2015); Biehl et al. (2017) (see Waxman (1995); Vietri (1995) for protons), low-luminosity GRBs and transrelativistic supernovae Murase et al. (2006, 2008); Wang et al. (2007); Chakraborty et al. (2011); Liu and Wang (2012), and newborn pulsars and magnetars Fang et al. (2012, 2014); Kotera et al. (2015). However, UHECR nuclei can be depleted before they escape the source environment due to the interaction with background particles or photons, so the “UHECR survival problem” is important to discuss the origin of UHECR nuclei Murase et al. (2008); Wang et al. (2008); Pe’er et al. (2009); Horiuchi et al. (2012); Kotera et al. (2015).

In this work, we revisit tidal disruption events (TDEs) as sources of UHECR nuclei. A TDE is a luminous flare lasting for months to years, which occurs in a galaxy nuclear region Rees (1988). When a star gets very close to the central black hole, it is disrupted if the tidal force is greater than the star’s self-gravity. During the disruption process, nearly half of the stellar debris falls back into the vicinity of the black hole, and the rest becomes unbound from the system. The accretion flow undertakes a fast energy dissipation and circularization process, and then it forms an accretion disk around the central black hole Shen and Matzner (2014); Piran et al. (2015). The TDE may be accompanied by the emergence of a relativistic jet during the super-Eddington accretion phase Giannios and Metzger (2011); Lei and Zhang (2011); Krolik and Piran (2011); De Colle et al. (2012); Shcherbakov et al. (2013); McKinney et al. (2014); MacLeod et al. (2014), and TDEs have been proposed to accelerate particles to ultrahigh energies Farrar and Gruzinov (2009); Farrar and Piran (2014); Pfeffer et al. (2017). Two scenarios are the most popular: disruption of a main-sequence (MS) star by a supermassive black hole (SMBH) and the disruption of a white dwarf (WD) by an intermediate mass black hole (IMBH). The latter could be especially interesting as sources that can inject heavy nuclei.

This paper is organized as follows. In Sec. II, we show that both protons and nuclei can be accelerated to ultrahigh energies, and we study conditions for the survival of UHECR nuclei in jetted TDEs. For the acceleration sites, we consider internal shocks, external reverse and forward shocks, and nonrelativistic winds. In Sec. III, for different composition models, we study the propagation of UHECR nuclei using the public code CRPropa 3 Alves Batista et al. (2016), and compare the results to experimental results by Auger.

Throughout the paper, we use cgs units, and adopt notations such as $Q_{x}\equiv Q/10^{x}$. The cosmological parameters we assume are $H_{0}=67.3\ \rm km\ s^{-1}\ Mpc^{-1}$, $\Omega_{m}=0.315$, and $\Omega_{\Lambda}=0.685$ Olive et al. (2014).

II Survival of cosmic-ray nuclei in TDE Outflows

We discuss possible composition models in the next section. In this section, we consider the fate of UHECR nuclei, following Ref. Murase et al. (2008).

Jetted TDEs (e.g., Swift J1644+57) show clear signatures of nonthermal emission in a wide range of wavelengths from radio to x rays Bloom et al. (2011); Burrows et al. (2011); Zauderer et al. (2011). Diffusive shock acceleration in collisionless shocks is the most popular mechanism of production of nonthermal particles. The first-order Fermi process is predicted to have a power-law distribution of accelerated particles within the “test particle” approximation. The acceleration time scale can be expressed as $t_{\rm acc}=\eta r_{L}/c$, where $r_{L}=E_{A}/ZeB$ is the Larmor radius of a particle with energy $E_{A}$, charge number $Z$, and $B$ the comoving-frame magnetic field strength. The factor $\eta$ depends on details of the turbulence. The minimum value of $\eta\sim 1-10$ can be achieved in the Bohm limit Lemoine and Waxman (2009), and $\eta=1$ is used to demonstrate our results in this section.

The maximum acceleration energy is determined by $t_{\rm acc}\leq{\rm min}(t_{\rm dyn},t_{\rm syn},t_{A\gamma})$, where $t_{\rm dyn}\equiv R/\Gamma\beta c$ is the dynamical time scale; $t_{\rm syn}=3m_{A}^{4}c^{3}\Gamma/(4\sigma_{T}Z^{4}m_{e}^{2}U_{B}E_{A})$ is the synchrotron cooling time scale ($\sigma_{T}$ is the Thompson cross-section, $U_{B}=B^{2}/8\pi$ is the magnetic energy density, $E_{A}$ is the particle energy); and $t_{A\gamma}$ is the energy loss time scale for protons (photomeson interaction) and nuclei (photodisintegration interaction). We can estimate the energy loss time scale using the following formula,

[TABLE]

where $\gamma_{A}$ is the Lorentz factor of UHECRs with mass number $A$, $\bar{\varepsilon}_{\rm th}$ is the threshold energy measured in the rest frame of initial nucleus (NRF), and ${\rm d}n/{\rm d}\varepsilon$ is the differential number density of target photons. Here, $\sigma_{A\gamma}$ is the photohadronic cross section related to the photomeson or photodisintegration process, and $\kappa_{A\gamma}$ is the inelasticity of each process. We show photonuclear and photohadronic cross sections for five typical chemical species as a function of NRF target photon energy ranging from $\sim 10\ \rm MeV$ to $\sim 10^{7}\ \rm GeV$ in Fig. 1; see Appendix A.

Defining the optical depth as $f_{A\gamma}\equiv t_{A\gamma}^{-1}/t_{\rm dyn}^{-1}$, we can expect the survival of UHECR nuclei when $f_{A\gamma}<1$. The value of the interaction time scale $t_{A\gamma-\rm int}$ is equal to the energy loss time scale when $\kappa_{A\gamma}=1$, and the related optical depth is $\tau_{A\gamma}\equiv t_{A\gamma-\rm int}^{-1}/t_{\rm dyn}^{-1}$.

II.1 Internal shock model

Nonthermal hard x-ray emission comes from internal energy dissipation in the inner relativistic jet Burrows et al. (2011). In the internal shock region, fast moving ejecta may catch up with slower ejecta, and a substantial amount of kinetic energy of the relativistic jet may be converted into internal energy. We assume the radius where internal collisions take place is $R_{\rm IS}=\Gamma^{2}c\delta t=3\times 10^{14}\ \Gamma_{1}^{2}\delta t_{2}\ \rm cm$, with $\Gamma=10\Gamma_{1}$ the Lorentz factor of the relativistic jet and $\delta t\sim 100\delta t_{2}\rm\ s$ is the x-ray variability time scale Burrows et al. (2011). The median x-ray luminosity of Swift J1644+57 is $L_{\rm X,iso}=8.5\times 10^{46}\ \rm{erg\ s^{-1}}$, which is well above the Eddington luminosity $L_{\rm Edd}=1.3\times 10^{44}M_{\rm BH,6}\ \rm erg\ s^{-1}$ of a $10^{6}M_{\odot}$ black hole Burrows et al. (2011). In our analysis, the total isotropic luminosity is set to $L_{\rm iso}=10^{48}\ \rm erg\ s^{-1}$. A fraction $\epsilon_{B}$ of the total energy of the outflow is converted into magnetic energy $B^{2}/8\pi=\epsilon_{B}L_{\rm iso}/4\pi R_{\rm IS}^{2}\Gamma^{2}c$. The magnetic field strength in the jet comoving frame is $B=(2\epsilon_{B}L_{\rm iso}/R_{\rm IS}^{2}\Gamma^{2}c)^{1/2}\simeq 860\ L_{\rm iso,48}^{1/2}\epsilon_{B,-1}^{1/2}R_{\rm IS,14.5}^{-1}\Gamma_{1}^{-1}\rm~{}G$. The maximum acceleration energy can be achieved under the condition $t_{\rm acc}\leq t_{\rm dyn}$; we have $E_{A,\rm dyn}\simeq\Gamma\eta^{-1}ZeB(R_{\rm IS}/\Gamma)\sim 1.5\times 10^{20}\ Z\eta^{-1}L_{\rm iso,48}^{1/2}\epsilon_{B,-1}^{1/2}\Gamma_{1}^{-1}\ \rm~{}eV$ measured in the observer frame, where $R_{\rm IS}/\Gamma$ is the comoving frame shell width. The maximum acceleration energy is also limited by the synchrotron energy loss $t_{\rm acc}<t_{\rm syn}$; we have $E_{A,\rm syn}\sim 6.3\times 10^{19}\ A^{2}Z^{-3/2}\eta^{-1/2}L_{\rm iso,48}^{-1/4}\epsilon_{B,-1}^{-1/4}\Gamma_{1}^{3/2}R_{\rm IS,14.5}^{1/2}\rm~{}eV$ measured in the observer frame. To estimate the effect of photonuclear and photohadronic interactions, we use a log-parabola model to fit the high-luminosity state SED of Swift J1644+57 Burrows et al. (2011); Senno et al. (2017).

Our results for the internal shock scenario are shown in Fig. 2. We consider two time scales, one is the interaction time scale, $t_{A\gamma-\rm int}$, and the other is the energy-loss time scale, $t_{A\gamma}$. We see that $t_{A\gamma-\rm int}$ and $t_{A\gamma}$ are shorter than $t_{\rm dyn}$. Our calculation suggests that most CR nuclei in the internal shock region will be disintegrated before they escape, so it is difficult for CR nuclei to survive for luminous TDEs like Swift J1644+57. Note that even partial survival is difficult for this parameter set (i.e., $f_{A\gamma}\gtrsim 1$), implying that the composition becomes light due to efficient photodisintegration.

We also consider the survival of UHECR nuclei in TDEs with a lower luminosity, and we assume an isotropic luminosity of $L_{\rm iso}=10^{44.8}\ \rm erg\ s^{-1}$. The comoving frame magnetic field strength is estimated to be $B\simeq 96\ L_{\rm iso,44.8}^{1/2}\epsilon_{B,-1}^{1/2}R_{\rm IS,14}^{-1}\Gamma_{1}^{-1}\rm~{}G$. The results are shown in Fig. 3, and we find that UHECR nuclei with energy up to $\sim 10^{20}\ \rm eV$ can survive before they escape from the source. Note that this case almost allows the complete survival of nuclei (i.e., $\tau_{A\gamma}\lesssim 1$). For intermediate luminosities corresponding to $\tau_{A\gamma}\gtrsim 1$ and $f_{A\gamma}\lesssim 1$, UHECR nuclei partially survive Murase and Beacom (2010) and secondary nuclei affect the initial composition of UHECRs.

II.2 Reverse shock model

The radio afterglow of Swift J1644+57 was observed a few days after the trigger of the Swift BAT observation Zauderer et al. (2011), and the continued observation extends to $\sim 500$ days Berger et al. (2012); Zauderer et al. (2013). The nonthermal radio emission is consistent with synchrotron radiation from the standard external shock model Metzger et al. (2012); Zauderer et al. (2013). Here, we assume the jet duration time is $t_{j}=10^{6}\ \rm s$ based on the observation of Swift J1644+57 Burrows et al. (2011). The isotropic equivalent kinetic energy of the jet is $L_{\rm iso}\sim 10^{47}\ \rm erg\ s^{-1}$ on average, and it will decrease following $L_{\rm iso}\propto t^{-5/3}$ after the time $t_{j}$. We estimate the total injection energy as ${\mathcal{E}}_{\rm iso}=2L_{\rm iso}t_{j}=2\times 10^{53}\ \rm erg$. The relativistic jet is decelerated when it has swept up a significant amount of circumnuclear medium (CNM). There are two shocks formed in the deceleration radius, one is the reverse shock propagating back into the relativistic jet, and the other is the forward shock propagating into the CNM. For simplicity, we adopt a constant CNM density, $\varrho=1\ \rm cm^{-3}$.

We calculate the evolution of reverse shock following the same method as the one used for GRBs Panaitescu and Kumar (2004); Murase (2007); Murase et al. (2008). The Lorentz factor of the shocked ejecta (relative to unshocked CNM) can be estimated as $\Gamma\approx\Gamma_{0}/(1+2\Gamma_{0}(\varrho/n_{j})^{1/2})^{1/2}$ assuming the pressure equilibrium at the contact discontinuity. We assume $\Gamma_{0}=10$, and the density of the jet $n_{j}={\mathcal{E}}_{\rm iso}/(4\pi m_{p}c^{2}\Gamma_{0}(\Gamma_{0}\Delta)r^{2})$, where ${\mathcal{E}}_{\rm iso}$ is the isotropic equivalent energy of the relativistic outflow, and $\Delta$ is the thickness of the ejecta (relative to the central black hole). The thickness of the ejecta is estimated to be $\Delta_{\times}\approx\Delta_{0}\equiv ct_{j}$. The shock completely crosses the ejecta at $t=t_{\times}$. We write the radius at that time as $r_{\times}\simeq 8.9\times 10^{17}{\mathcal{E}}_{\rm iso,53.3}^{1/4}t_{j,6}^{1/4}\varrho^{-1/4}\ {\rm cm}$ and the Lorentz factor at the crossing time is $\Gamma_{\times}\simeq 3.2{\mathcal{E}}_{\rm iso,53.3}^{1/8}t_{j,6}^{-3/8}\varrho^{-1/8}$. The magnetic field strength is estimated to be $B_{\times}=(32\pi\epsilon_{B}^{r}n_{j}m_{p}c^{2}(\Gamma_{\times}-1)(\Gamma_{\times}+3/4))^{1/2}$, where a fraction $\epsilon_{B}^{r}$ of the post-shock internal energy is converted into magnetic energy. The comoving frame Lorentz factor of the electrons is $\gamma_{e}^{r}=\frac{\epsilon_{e}^{r}}{f_{e}^{r}}\frac{p-2}{p-1}\frac{m_{p}}{m_{e}}(\Gamma_{\times}-1)$, with $p=2.2$, $\epsilon_{e}^{r}=0.02$, and $f_{e}^{r}=0.1$. Here $f_{e}^{r}$ represents the number fraction of accelerated electrons.

Then, the peak synchrotron frequency is written as

[TABLE]

where $g_{1}(\Gamma_{\times})=\Gamma_{\times}(\Gamma_{\times}-1)^{5/2}(\Gamma_{\times}+3/4)^{1/2}$. The synchrotron self-absorption frequency can be estimated when $\tau(\nu_{a})=1$; we have ($\nu_{a}<\nu_{m}$)

[TABLE]

where $g_{2}(\Gamma_{\times})=g_{1}(\Gamma_{\times})(\Gamma_{\times}-1)^{-33/10}(\Gamma_{\times}+3/4)^{-3/10}$. Further, we can express the comoving frame peak luminosity per unit energy from the reverse shock as

[TABLE]

where $g_{3}(\Gamma_{\times})=(\Gamma_{\times}-1)^{1/2}(\Gamma_{\times}+3/4)^{1/2}$ and $N_{e}^{r}={\mathcal{E}}_{\rm iso}/\Gamma_{0}m_{p}c^{2}$ is the number of electrons in the reverse shock. The photon spectrum in the reverse shock is ($\nu_{a}<\nu_{m}<\nu_{c}$)

[TABLE]

where $\varepsilon_{m}=\varepsilon_{m,{\rm ob}}/\Gamma_{\times}$ is the peak photon energy in the reverse shock comoving frame and $\varepsilon_{c}$ is the electron synchrotron cooling energy. The maximum acceleration energy of nuclei can be derived when $t_{\rm acc}=t_{\rm dyn}$,

[TABLE]

We estimate various time scales in the reverse shock model, and our results are shown in Fig. 4. We find that the interaction time scale $t_{A\gamma-\rm int}$ is longer than the dynamical time scale $t_{\rm dyn}$, which means that the optical depth $\tau_{A\gamma}<1$. Our results suggest that UHECR nuclei can survive in the reverse shock model.

II.3 Forward shock model

Once the relativistic jet enters the CNM, the deceleration and transition to the Blandford-Mckee (BM) self-similar regime occurs. The evolution of the Lorentz factor and shock radius are described as $\Gamma(t)\propto t^{-3/8}$ and $r(t)\propto t^{1/4}$, where we assume a constant CNM density. The evolution of the downstream magnetic field follows $B\propto t^{-3/8}$. A fraction of electrons are accelerated in the shock, with a minimum Lorentz factor, $\gamma_{m}$, and the distribution of accelerated electrons is denoted as $dN_{e}/d\gamma_{e}\propto\gamma^{-p}$ with $p=2.2$, and minimum Lorentz factor $\gamma_{m}\propto t^{-3/8}$. In the following, we use the approximation $\Gamma\gg 1$ to derive the simplified expression of typical frequency and luminosity, and we adopt the accurate formula in our calculation due to the mildly relativistic nature of forward shock. The peak frequency from the forward shock can be estimated as

[TABLE]

The self-absorption frequency (assuming $\nu_{a}<\nu_{m}<\nu_{c}$) is

[TABLE]

The comoving frame peak luminosity per unit energy is

[TABLE]

where $N_{e}=(4\pi/3)r^{3}\varrho$ is the total number of swept-up electrons. The comoving frame photon spectrum can be derived in the slow cooling case ($\nu_{a}<\nu_{m}<\nu_{c}$),

[TABLE]

where $\varepsilon_{m}=\varepsilon_{m,{\rm ob}}/\Gamma$ is the peak photon energy in the forward shock comoving frame and $\varepsilon_{c}$ is the electron synchrotron cooling energy.

In the forward shock model, the observed maximum particle energy is achieved when $t_{\rm acc}=t_{\rm dyn}$, with

[TABLE]

where it is assumed that the upstream magnetic field is amplified and comparable to the downstream value. The shock is mildly relativistic and this could be achieved by CR streaming instabilities Lemoine and Revenu (2006); Lemoine et al. (2006). Or one could use the downstream field if particles are accelerated by the second-order Fermi acceleration mechanism Dermer and Humi (2001). We show various time scales in Fig. 5. We expect the survival of UHECR nuclei at the forward shock, because the interaction time scale, $t_{A\gamma}$, is longer than the dynamical time scale, $t_{\rm dyn}$. We also show the evolution of particle maximum acceleration energy and optical depth $f_{A\gamma}$ in Fig. 6. We found that it is easier for UHECR nuclei to survive at later times.

However, it may be difficult to amplify the magnetic field in the upstream region, especially if the shock is ultrarelativistic Murase (2007); Gallant and Achterberg (1999); Sironi et al. (2015). In this case, the magnetic field in the shock region should be similar to the CNM magnetic field $\hat{B}_{\rm cnm}=10^{-5}\hat{B}_{\rm cnm,-5}$. The observed maximum acceleration energy is estimated to be

[TABLE]

It is difficult to accelerate CRs to ultrahigh energies if the magnetic field is not amplified efficiently.

II.4 Nonrelativistic wind model

Along with the formation of a relativistic jet, a nonrelativistic outflow can be driven by the accretion disk. The radio emissions from two TDEs, ASASSIN14li Alexander et al. (2016a) and XMMSL1J0740-85 Alexander et al. (2016b), are consistent with the emission from the interaction between a nonrelativistic outflow and the CNM. We assume that the nonrelativistic outflow has an ejected mass of $M_{\rm ej}\sim 10^{-4}M_{\odot}$ and velocity $v_{\rm ej}\sim 0.1c$. The deceleration radius is estimated to be $r_{\rm dec}\simeq 5.4\times 10^{16}{\mathcal{E}}_{k,48}^{1/3}\varrho^{-1/3}\ \rm cm$ and the deceleration time is $t_{\rm dec}\simeq 1.8\times 10^{7}{\mathcal{E}}_{k,48}^{1/3}\varrho^{-1/3}\ \rm s$. The shocked fluid enters into the Sedov-Taylor evolution phase after it is decelerated by the external medium. In the adiabatic case, the shock velocity and radius follow $V_{\rm s}\propto t^{-3/5}$ and $R_{\rm s}\propto t^{2/5}$, respectively. We assume the accelerated electrons have power-law index $p=3$. The synchrotron peak frequency is

[TABLE]

and the self-absorption frequency ($\nu_{a}>\nu_{m}$) is

[TABLE]

The photon spectrum in this case is($\nu_{m}<\nu_{a}<\nu_{c}$)

[TABLE]

In the nonrelativistic case, the particle acceleration time scale is $t_{\rm acc}\approx(20/3)c^{2}E_{A}/V_{s}^{2}ZeBc$ in the Bohm limit. Our results for various time scales are shown in Fig  7. The dynamical time in the nonrelativistic case is $t_{\rm dyn}=R_{s}/V_{s}$. The maximum acceleration energy is limited by $t_{\rm acc}\leq t_{\rm dyn}$, with $E_{A,\rm max}\approx(3/20)R_{s}ZeBV_{s}/c\sim 1.5\times 10^{15}Z{\mathcal{E}}_{k,48}^{3/5}\epsilon_{B,-1}^{1/2}\varrho^{-1/10}(t/t_{\rm dec})^{-4/5}\ {\rm~{}eV}$. In the wind model, CRs cannot be accelerated to ultrahigh energies.

III Propagation of Escaping Nuclei

In the previous section, we found that UHECR nuclei can survive for external shocks. For internal shocks, the survival is difficult for Swift J1644+57-like TDEs, but lower-luminosity TDEs allow UHECR nuclei to escape without significant disintegration. In this section, for simplicity, we assume that UHECR nuclei survive and escape into intergalactic space and see consequences of different composition models 111This can be a reasonable approximation if the luminosity function has a steep index, i.e., the luminosity density is dominated by lower-luminosity TDEs..

We assume that UHECR nuclei follow a power-law distribution with an exponential cutoff as

[TABLE]

where $f_{A^{\prime}}$ is the number fraction of different particles with mass $A^{\prime}$ at the same rigidity, $N_{0}$ is determined by the total energy per TDE event (see Appendix B), $s_{\rm esc}$ is the spectral index of ejected UHECR nuclei, and $ZE^{\prime}_{p,\rm max}$ is the maximum acceleration energy for particles with charge number $Z$. We assume the minimum particle energy $E^{\prime}_{A,\rm min}=\Gamma^{2}Am_{p}c^{2}$. Strictly speaking, this is justified in the forward shock model but our main conclusion does not change by this assumption.

III.1 Injection spectrum

In order for UHECR nuclei to originate from stellar material (whether it is a MS star or WD), CR injections should occur inside jets or winds, which involve internal shocks and an external reverse shock. Details are highly uncertain but there are various possibilities.

The first possibility is to rely on the shock acceleration at internal shocks. In the diffusive shock acceleration mechanism, accelerated particles have a power-law distribution, $dN/dE\propto E^{-s_{\rm acc}}$, with a typical spectral index $s_{\rm acc}\sim 2.2$ in the ultrarelativistic limit for the small-angle scattering approximation Kirk et al. (2000); Achterberg et al. (2001); Keshet and Waxman (2005). However, the spectral index can be affected by the deflection angle, the ratio of the scattering mean free path to the particle gyroradius and the orientation of the magnetic field to the shock normal. In the large-angle scattering case, where magnetic fluctuations are sufficiently large, particles can gain significant energy in the single scattering and may lead to a harder spectrum Kato and Takahara (2001); Aoi et al. (2008); Baring (2011). Such a possibility has been discussed to explain a hard spectrum of blazars Baring (2009, 2011). In addition, the spectrum of escaping CRs does not have to be the same as that of accelerated CRs (e.g., Ohira et al. (2010); Caprioli et al. (2009); Katz et al. (2010); see also Fang and Murase (2017); Globus et al. (2015); Baerwald et al. (2013)). If CRs could be confined for a long time after the dynamical time, CRs will lose their energies during their diffusive escape, so that the spectrum of escaping CRs is harder. Although details depend on flow dynamics and magnetic field evolution, one can assume that UHECRs escaping from internal shocks may have a small index such as $s_{\rm esc}<2$. For an expanding outflow, one of the most conservative possibilities is to invoke the direct escape of CRs (e.g., Baerwald et al. (2013)), which leads to $s_{\rm esc}=s_{\rm acc}-1$. For $s_{\rm acc}\sim 2$, one may expect $s_{\rm esc}\sim 1$.

The second possibility is to invoke a two-step acceleration mechanism via stochastic acceleration in the downstream of external shocks or possible reverse shock acceleration. The downstream may be highly turbulent and mixed around the contact discontinuity due to Rayleigh-Taylor instabilities Levinson (2009), and a fraction of energized nuclei from jets can be used for further acceleration to ultrahigh energies at the external forward shock Dermer and Humi (2001). If heavy nuclei from jets can be used for CRs and the maximum energy is limited by the amplified magnetic field, accelerated CRs escape from a relativistic decelerating blast wave. The escaping spectrum can be harder in some cases Katz et al. (2010) (see Ref. Ohira et al. (2010) for discussion on supernova remnants and AGN cocoon shocks). The evolution of total energy follows ${\mathcal{E}}_{\rm CR}\propto r^{-\alpha_{\mathcal{E}}}$, the evolution of external medium density is $\varrho\propto r^{-\alpha_{\varrho}}$, and the evolution of shocked fluid bulk Lorentz factor is $\Gamma\propto r^{-\alpha_{\Gamma}}$. The minimum energy of accelerated particles can be estimated as $E_{A,\rm min}\simeq\Gamma^{2}Am_{p}c^{2}\propto r^{-\alpha_{\rm min}}$, and the maximum particle energy is $E_{A,\rm max}\simeq ZeBr\propto r^{-\alpha_{\rm max}}$, as we discussed in the previous section. In the adiabatic expansion case, we have $\alpha_{\rm min}=2\alpha_{\Gamma}$, $\alpha_{\rm max}=\alpha_{\Gamma}+(1/2)\alpha_{\varrho}-1$, and $\alpha_{\mathcal{E}}=2\alpha_{\Gamma}+\alpha_{\varrho}-3$. We assume the accelerated CRs spectrum $\propto E^{-s_{\rm acc}}$, and have spectrum $\propto E^{-s_{\rm esc}}$ after escape from the acceleration site; Ref. Katz et al. (2010) derived a simple analytic relation between $s_{\rm esc}$ and $s_{\rm acc}$ as $s_{\rm esc}=s_{\rm acc}-(\alpha_{\rm min}(s_{\rm acc}-2)+\alpha_{\mathcal{E}})/\alpha_{\rm max}$ with the assumption that the number of ejected UHECRs in an energy interval is the same as that of the CRs at radius $r$, where the maximum acceleration energy is achieved. Let us adopt a typical value of power-law index $s_{\rm acc}\sim 2.2$. In the adiabatic expansion scenario, we have $\alpha_{\mathcal{E}}=0$. If the CNM density is constant, then we have $\alpha_{\varrho}=0$. The escaping particle spectral index can be calculated as $s_{\rm esc}=-5s_{\rm acc}+12=1$. However, if the CNM density decreases with radius as $\varrho\propto r^{-\alpha_{\varrho}}$, we have $s_{\rm esc}=1.4$ for $\alpha_{\varrho}=1$ and $s_{\rm esc}=1.8$ for $\alpha_{\varrho}=2$. The CNM gas density in the galaxy nuclear region is still unclear. It can originate from the stellar winds, and their profile depends on the detailed distribution of stars and star formation history in the galaxy nuclear region Generozov et al. (2015, 2017). We assume that the CNM density is constant in our analysis, which is reasonable for galaxy cores Generozov et al. (2015).

Although the above arguments are rather speculative, it is possible to expect a hard spectral index for escaping CRs, and we use $s_{\rm esc}\sim 1$ in the jetted TDE scenario. In addition, we assume the maximum acceleration energy $E_{A,\rm max}$ according to the rigidity dependence, which is valid in the internal shock model for relatively low-luminosity TDEs and the external forward shock model.

III.2 Composition model

III.2.1 MS-SMBH tidal disruptions

In model I, we consider a solar-type MS star disrupted by a SMBH Rees (1988). The present-day solar composition for H, He and metals is $X=0.793$, $Y=0.2469$, and $Z=0.0141$, and the most abundant heavy nuclei are $\rm O$, $\rm C$, $\rm Ne$, and $\rm Fe$ Lodders (2010). In Fig. 8, we show the CR injection spectra for model I. The mass fractions $X_{A}$ for different chemical species at the same rigidity are similar to the MS star. We assume that the maximum acceleration energy is $ZE_{p,\rm max}=2\times 10^{19}Z\ \rm~{}eV$.

III.2.2 WD-IMBH tidal disruptions

One of the models for Swift J1644+57 is a WD tidally disrupted by a $10^{4}M_{\odot}$ IMBH Burrows et al. (2011); Bloom et al. (2011); Zauderer et al. (2011); Krolik and Piran (2011). The radii of WDs are smaller than those of MS stars, so WDs should be tidally disrupted at smaller radii from the central BHs. We need the factor $\beta_{g}=\frac{R_{t}}{R_{g}}>1$ to avoid WDs being swallowed by BHs, where $R_{t}$ is the tidal disruption radius and $R_{g}$ is the Schwarzschild radius. The upper limit on the BH mass is estimated to be $M_{\rm BH}\lesssim 2.6\times 10^{5}\ {M_{\odot}}\left(\frac{R_{\rm WD}}{10^{9}\ \rm cm}\right)^{3/2}\left(\frac{M_{\rm WD}}{0.6M_{\odot}}\right)^{-\frac{1}{2}}$, where $R_{\rm WD}$ and $M_{\rm WD}$ are the radius and mass of WDs Luminet and Pichon (1989). Most of the WDs are composed of carbon and oxygen as the result of helium burning in the core, and a fraction of WDs may contain neon and magnesium due to the ignition of carbon. The hydrogen rich envelope of WDs could be ejected due to helium thermal pulses in the asymptotic giant branch phase.

In model II-1, we consider that tidally disrupted WDs are carbon-oxygen WDs (CO-WDs). CO-WDs are the end point for low and intermediate mass stars in the mass interval from $0.6M_{\odot}$ to $6M_{\odot}$ D’Antona and Mazzitelli (1990); Weidemann (2000). The ratio of carbon to oxygen depends on the reaction ${}^{12}{\rm C}(\alpha,\gamma){}^{16}{\rm O}$, which is unresolved yet in nuclear physics Horowitz et al. (2010); Fields et al. (2016). We assume that CO-WDs have a roughly equal mass fraction for carbon and oxygen $X_{\rm C}=0.5$, $X_{\rm O}=0.5$. Stars with a mass in the range of $\sim 8M_{\odot}$ to $\sim 12M_{\odot}$ may evolve to form oxygen-neon-magnesium WDs (ONeMg-WDs) Nomoto (1984); Iben and Tutukov (1985). In this case, the burning of carbon will not lead to explosion or collapse. In model II-2, we consider that ONeMg-WDs with a mass fraction of $X_{\rm O}=0.12$, $X_{\rm Ne}=0.76$, and $X_{\rm Mg}=0.12$ are tidally disrupted Isern et al. (1991). In Fig. 9, we show the CR injection spectra for model II-1 (upper panel) and model II-2 (lower panel).

III.2.3 WD-IMBH with ignition

When a WD approaches a massive BH, the tidal compression and relativistic effects can enhance the WD central density and could trigger explosive nuclear burning Luminet and Pichon (1989); Wilson and Mathews (2004); Dearborn et al. (2005); Rosswog et al. (2008, 2009); Haas et al. (2012); Sell et al. (2015); MacLeod et al. (2016); Tanikawa et al. (2017). This kind of ignition has also been suggested as an alternative scenario for type Ia supernovae Wilson and Mathews (2004). In this case, a fraction of nuclear explosive matter can be accreted into the center BH and form an accreting flow Rosswog et al. (2009).

However, ignition and associated nucleosynthesis of heavy nuclei have been questioned by more dedicated simulations, and details are still under debate Tanikawa et al. (2017). Also, the rate of such events is uncertain and may be much smaller. On the other hand, it has been argued that the composition of such TDEs has been considered to explain UHECRs Alves Batista and Silk (2017), so we also consider this scenario for completeness.

In models III-1 and III-2, we adopt the numerical simulation results from Ref. Rosswog et al. (2009). Model III-1 is the case with a $0.2M_{\odot}$ helium WD passing a $10^{3}M_{\odot}$ BH and model III-2 corresponds to a CO-WD ($1.2M_{\odot}$) approaching a $500M_{\odot}$ BH . In Fig. 10, we show the CR injection spectra for model III-1 (upper panel) and model III-2 (lower panel).

III.3 Propagation in intergalactic space

We calculate the propagation of UHECR nuclei, using the public code CRPropa 3 Armengaud et al. (2007); Alves Batista et al. (2016). The main energy-loss process for UHECR nuclei during propagation is photodisintegration due to CMB photons and extragalactic background light (EBL) photons. The EBL photons have a larger effect on the propagation of lower-energy intermediate mass nuclei ($\sim 10^{19}\ \rm~{}eV$) Allard (2012); Aab et al. (2017). In this work, we adopt a semianalytic EBL model derived by Ref. Gilmore et al. (2012). The details of the simulation are shown in Appendix B. The observed UHECR spectrum has been accurately measured by Auger Aab (2015) and TAFukushima (2015). We use an empirical model, which describes the distribution of $X_{\rm max}$ using the generalized Gumbel function ${\mathcal{G}}(X_{\rm max})$ to get the mean value of the depth of shower maximum $\langle X_{\rm max}\rangle$ and $\sigma(X_{\rm max})$ De Domenico et al. (2013). In this work, we adopt the EPOS-LHC hadronic interaction model Werner et al. (2011); Aab et al. (2014a).

In Fig. 11, we show the results derived from model I. This model can fit the observed UHECR spectrum measured by Auger, but failed to fit $\langle X_{\rm max}\rangle$ and $\sigma(X_{\rm max})$; the final spectrum is proton dominated in nearly the entire energy range ($<10^{20}\ \rm~{}eV$). Although the situation is still under debate due to significant uncertainties in hadronic interaction models, this proton dominated scenario seems in strong tension with Auger results Aab et al. (2014b).

In Fig. 12, we show the results from model II-1, where most UHECRs are carbon and oxygen. The energy-loss lengths of C and O nuclei are very short, so most of the observed C and O nuclei should come from nearby sources. There should be a large fraction of secondary protons and helium generated during the propagation. We find that model II-1 can fit the UHECR spectrum measured by Auger reasonably well, except the spectrum becomes softer in the ”ankle” region ($\sim 10^{18.6}\rm~{}eV$). In Fig. 13, we show the results from model II-2. In this model, UHECRs are mostly heavier nuclei, Ne and Mg, which is expected when WDs have higher mass such as ONeMg WDs. We find that the final spectrum can be fitted very well in this scenario. The data of $\langle X_{\rm max}\rangle$ and $\sigma(X_{\rm max})$ in both models are consistent with Auger results. Our results are consistent with their UHECR composition that becomes heavier with increasing energy, and where intermediate mass particles dominate in the high-energy range ($E>10^{19}\ \rm~{}eV$).

For the tidal disruption of lower-mass WDs (helium WDs or CO-WDs) by massive BHs, we also consider the enhancement of heavier nuclei through explosive nuclear reactions. Our results are shown in Fig. 14 for model III-1 and Fig. 15 for model III-2. UHECRs in model III-1 are predominantly He, Si and Fe. This model can fit the UHECR spectrum reasonably in the lower energy range, but the cutoff energy ($\sim 10^{20}\ \rm~{}eV$) is more consistent with TA results and higher than Auger results. Compared to model III-1, model III-2 has a large fraction of iron group nuclei in UHECRs which predict an even higher cutoff energy. We find that model III-1 is consistent with $\langle X_{\rm max}\rangle$ and $\sigma(X_{\rm max})$ measured by Auger, while model III-2 has a large fraction of heavier nuclei compared to the observed composition data in the higher-energy range $(>10^{19.2}\rm~{}eV)$.

IV Discussion and summary

This work consists of two parts. In the first part, we examined the production of UHECRs in TDEs accompanied by relativistic jets. In the internal shock model, we found that CRs can be accelerated to ultrahigh energies. However, it is difficult for UHECR nuclei to survive in luminous TDE jets such as Swift J1644+57, while the survival is allowed for less powerful TDE jets. In the reverse and forward shock model, we can expect the production and survival of UHECRs. We also considered the wind model, where a nonrelativistic outflow interacts with the CNM. We found that in this scenario CRs can be accelerated only up to $\sim 10^{15}~{}Z\rm~{}eV$.

In the second part, we examined different composition models for TDEs. Motivated by the composition data by Auger, we assumed that the injected UHECR nuclei mainly originate from tidally disrupted stars. Although we discussed several possibilities, this should be justified by more detailed work on CR acceleration and escape processes. We consider both MS-SMBH and WD-IMBH tidal disruptions. In the MS-SMBH TDE scenario (model I), the injected UHECRs have a solarlike composition. The UHECR spectrum can be fitted, but a proton dominated composition is expected in the nearly the entire energy range. In the WD-IMBH TDE scenario, model II-1 (CO-WDs) can give a poor fit to the UHECR spectrum but $\langle X_{\rm max}\rangle$ and $\sigma(X_{\rm max})$ can be reasonably accounted for. We found that it is difficult to fit the spectrum and composition simultaneously even if we try a variety of parameter sets, such as higher maximum energy and/or steeper ejection spectra. The main reason is that the attenuation lengths of UHECR carbon or oxygen nuclei are lower than protons or iron nuclei, so most of them will be depleted into secondary protons before reaching Earth Kotera and Olinto (2011). For ONeMg-WDs, we found that the results are more consistent with Auger data. However, the rate density of ONeMg-WDs is expected to be lower than that of CO-WDs ($\sim 1/30$), if we assume the Salpeter initial mass function Livio and Truran (1994); Gil-Pons, P. et al. (2003); Denissenkov et al. (2014) (see also Appendix C). The uncertainties of TDEs redshift evolution do not affect our results very much because most UHECRs mainly come from nearby sources within the GZK horizon ($\sim 100\ \rm Mpc$) Kotera and Olinto (2011).

We also considered special cases, where nucleosynthesis is triggered by an IMBH. Our results show that model III-1 can marginally fit the observed spectrum, and consistent with the composition data, while model III-2 has too many iron group nuclei, making it difficult to reconcile with the Auger data. We also caution that the final composition is sensitive to details of the parameters, such as WD and BH masses Rosswog et al. (2009); MacLeod et al. (2016); Kawana et al. (2017). Reference Tanikawa et al. (2017) performed 3-D smoothed particle hydrodynamics simulations to study the explosive nuclear burning of WD-IMBHs, considering both the adiabatic compression and shock wave generation. They found that the final mass fractions are sensitive to the number of simulation particles, and ignitions occur for low-resolution simulations.

Secondary gamma-ray and neutrino signals are of interest to test the model. TDEs have been discussed as high-energy neutrino sources Murase (2008); Murase and Takami (2009); Wang et al. (2011); Wang and Liu (2016); Dai and Fang (2016); Senno et al. (2017); Lunardini and Winter (2016). Neutrino production is expected to be efficient for luminous TDEs such as Swift J1644+57, while the neutrino production efficiency should be lower for low-luminous TDEs, allowing nucleus survival Murase et al. (2008); Murase and Beacom (2010). We also estimate the cosmogenic neutrino flux according to model II-2 (ONeMg-WDs), and we find $E_{\nu}^{2}\Phi_{\nu}\sim{10}^{-10}\ \rm GeV\ cm^{-2}\ sr^{-1}s^{-1}$. However, this would be too conservative since luminous TDEs also contribute to the neutrino flux. High-energy gamma rays can be produced via neutron pion decay, photodeexcitation, and Bethe-Heitler processes. The survival of UHECR nuclei implies that gamma rays can escape from the sources without efficient two-photon annihilation inside the sources Murase et al. (2008); Murase and Beacom (2010).

Perhaps, one could expect the emission of gravitational waves from WD-IMBH tidal disruptions Kobayashi et al. (2004); Rosswog et al. (2009); Haas et al. (2012), and TDEs can be interesting targets for multimessenger astronomy.

Acknowledgements.

We thank Shigeo S. Kimura and Kumiko Kotera for helpful discussions. The work of K.M. is supported by Alfred P. Sloan Foundation and the U.S. National Science Foundation (NSF) under grants NSF Grant No. PHY-1620777. It is partially supported by the National Natural Science Foundation of China under Grant No. 11273005 and the National Basic Research Program (973 Program) of China under Grant No. 2014CB845800 (Z.L.). B.T.Z. is supported by the China Scholarship Council (CSC) to conduct research at Penn State University.

Appendix A Inelasticity of photodisintegration process

The giant dipole resonance (GDR) is the most important disintegration process at low energies, and the threshold energy is $\sim 8\ \rm MeV$ Rachen (1996); Khan et al. (2005); Allard et al. (2005). The quasideutron (QD) process dominates when the NRF photon energy is larger than $30\ \rm MeV$, and lower than the photopion production threshold energy. In this energy range, the wavelength of projectile photons is comparable with the size of nuclear, leading to the ejection of nucleon pair and additional nucleons. When the NRF photon energy is larger than $\sim 150\ \rm MeV$, baryonic resonances (BR) play a dominant role in the photodisintegration process. At very high energies $\sim 1\ \rm GeV$, a nucleus can be broken into many fragments of much lower energies via photofragmentation (PF).

In this work, the photodisintegration cross section for UHECR nuclei is taken from CRPropa 3 Armengaud et al. (2007); Alves Batista et al. (2016). The cross section for nuclei with mass numbers of $A\geq 12$ and NRF photon energy in the range $0.2-200\ \rm MeV$ can be derived using the nuclear reaction code TALYS Koning et al. (2007). At higher energies, for simplicity, we use the relation that the cross section is proportional to nuclear mass, which can be written as $\sigma_{A\gamma}(\bar{\varepsilon})=A\sigma_{p\gamma}(\bar{\varepsilon})$ (see Fig. 1). The photomeson production cross section for protons is derived using the numerical code Geant4 Agostinelli et al. (2003); Murase et al. (2008).

In the photodisintegration process, a parent nucleus loses energy through the ejection of one or more nucleons. We assume that the Lorentz factor is constant during the interaction; then the relative energy loss can be estimated as Rachen (1996)

[TABLE]

where $N$ is the total number of nucleons in the parent nuclei, and $\Delta N$ is the number of ejected nucleons in each channel. To derive the average energy loss, we consider the contribution from all the dominant channels,

[TABLE]

where $\sigma_{\rm tot}$ is the total cross section, $\sigma_{i}$ is the cross section for the ith channel, and $[\Delta N]_{i}$ is the number of ejected nucleons in the ith channel. In the energy range beyond the photopion production threshold, for simplicity, the inelasticity is assumed to linearly decrease with increasing nuclear mass $\bar{\kappa}_{A\gamma}(\bar{\varepsilon})=\bar{\kappa}_{p\gamma}(\bar{\varepsilon})/A$.

Appendix B UHECRs propagation

The observed CR flux for each component is described by the following formula Unger et al. (2015); Alves Batista et al. (2016)

[TABLE]

where $\Phi_{AA^{\prime}}$ and ${\rm d}n_{AA^{\prime}}(E)$ are the observed cosmic ray flux and number density of particles with mass $A$ generated from parent particles with mass $A^{\prime}$. ${\rm d}n_{AA^{\prime}}(E)$ can be estimated by considering the contribution as a function of redshift

[TABLE]

where

[TABLE]

and we use the redshift evolution of TDEs given by Sun et al. (2015)

[TABLE]

with $\eta\sim-2$. In this work, we adopt $z_{\rm min}=0.0001$ and $z_{\rm max}=2$. Note that the redshift evolution for TDEs is negative. $\rho_{0}$ is the local event rate of jetted TDEs Sun et al. (2015). Also, $\eta_{AA^{\prime}}(E,E^{\prime},z)$ is the fraction of generated CRs with mass $A$ and energy $E$ from parent particles with mass $A^{\prime}$ and energy $E^{\prime}$. ${\rm d}N_{A^{\prime}}/{\rm d}E^{\prime}$ is the UHECR injection spectra per TDE.

The CR luminosity density can be estimated as

[TABLE]

where $\Gamma(2-s_{\rm esc},\frac{E^{\prime}_{\rm min}}{ZE^{\prime}_{p,\rm max}})$ is the incomplete gamma function. We use the following formula to estimate the normalization parameter $N_{0}$:

[TABLE]

The CR energy per TDE (${\mathcal{E}}_{\rm CR}$) can be estimated through the relation ${\mathcal{E}}_{\rm CR}^{\rm iso}=\xi_{\rm CR}{\mathcal{E}}_{\rm rad}^{\rm iso}$, where $\xi_{\rm CR}$ is the cosmic ray loading factor.

Appendix C Energy budget

The event rate of WD-IMBH tidal disruptions is uncertain. IMBHs are believed to exist in dwarf galaxies or globular clusters. In dwarf galaxies, the event rate is estimated to be $R_{\rm TDE-DG}\sim 10f_{\rm IMBH}^{\rm DG}\ \rm Gpc^{-3}\ yr^{-1}$, where $f_{\rm MBH}^{\rm DG}$ is the occupation fraction in dwarf galaxies MacLeod et al. (2014). In globular clusters, the event rate is $R_{\rm TDE-GC}\sim 50f_{\rm MBH}^{\rm GC}\ \rm Gpc^{-3}\ yr^{-1}$, which is slightly higher than the event rate estimated for dwarf galaxies Shcherbakov et al. (2013).

We assume only a fraction of $f_{\rm jet}\sim 10\%$ WD-IMBH tidal disruptions have relativistic jets with the beaming factor $f_{b}\approx\theta_{j}^{2}/2\sim 1/(2\Gamma_{j}^{2})$, with a typical Lorenz factor of $\Gamma_{j}=10$. The apparent rate density of WD-IMBH tidal disruptions is

[TABLE]

with the assumption $f_{\rm IMBH}\simeq f_{\rm IMBH}^{\rm DG}\simeq f_{\rm IMBH}^{\rm GC}\simeq 1$. The event rate can be comparable to the event rate ($\rho_{\rm TDE}\simeq 0.03\ \rm Gpc^{-3}\ yr^{-1}$) obtained from the detection of two jetted TDEs Swift J1655+57 and Swift J2058+05 Sun et al. (2015).

According to our simulation results, we need a total CR energy injection rate of $Q_{\rm CR}\approx 4.2\times 10^{44}\ \rm erg\ Mpc^{-3}\ yr^{-1}$ in model II-2, and the required CR loading factor is $\xi_{\rm CR}\sim 100\ Q_{\rm CR,44.6}{\rho_{0}}_{-10.2}^{-1}\ {\mathcal{E}}_{\rm rad,53}^{-1}$. Note that the luminosity density itself is independent on the beaming factor as in the argument for GRBs. The absolute radiation energy is smaller than the isotropic-equivalent radiation energy by $f_{b}^{-}1$, whereas the true rate density is a factor of $f_{b}^{-}1$ larger than the apparent rate density of TDEs, of which jet points to us.