# Multimessenger Implications of AT2018cow: High-energy Cosmic-Ray and Neutrino Emissions from Magnetar-powered Superluminous Transients

Ke Fang, Brian D. Metzger, Kohta Murase, Imre Bartos, Kumiko Kotera

► To cite this version:

Ke Fang, Brian D. Metzger, Kohta Murase, Imre Bartos, Kumiko Kotera. Multimessenger Implications of AT2018cow: High-energy Cosmic-Ray and Neutrino Emissions from Magnetar-powered Superluminous Transients. The Astrophysical Journal, 2019, 878 (1), pp.34. 10.3847/1538-4357/ab1b72. hal-01982915

HAL Id: hal-01982915

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

Submitted on 13 Aug 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Multimessenger Implications of AT2018cow: High-energy Cosmic-Ray and Neutrino Emissions from Magnetar-powered Superluminous Transients

Ke Fang¹, Brian D. Metzger², Kohta Murase³,⁴,⁵,⁶, Imre Bartos²,⁷, and Kumiko Kotera⁸,⁹

¹Einstein Fellow, Kavli Institute for Particle Astrophysics and Cosmology (KIPAC), Stanford University, Stanford, CA 94305, USA

²Department of Physics and Columbia Astrophysics Laboratory, Columbia University, New York, NY 10027, USA

³Department of Physics, The Pennsylvania State University, University Park, PA 16802, USA

⁴Department of Astronomy & Astrophysics, The Pennsylvania State University, University Park, PA 16802, USA

⁵Center for Particle and Gravitational Astrophysics, The Pennsylvania State University, University Park, PA 16802, USA

⁶Yukawa Institute for Theoretical Physics, Kyoto, Kyoto 606-8502 Japan

⁷Department of Physics, University of Florida, Gainesville, FL 32611, USA

⁸Sorbonne Universités, UPMC Univ. Paris 6 et CNRS, UMR 7095, Institut d'Astrophysique de Paris, 98 bis bd Arago, F-75014 Paris, France

⁹Laboratoire AIM-Paris-Saclay, CEA/DSM/IRFU, CNRS, Université Paris Diderot, F-91191 Gif-sur-Yvette, France

Received 2019 January 16; revised 2019 April 1; accepted 2019 April 21; published 2019 June 11

## Abstract

Newly born, rapidly spinning magnetars have been invoked as the power sources of superluminous transients, including the class of “fast blue optical transients” (FBOTs). The extensive multiwavelength analysis of AT2018cow, the first FBOT discovered in real time, is consistent with the magnetar scenario and offers an unprecedented opportunity to comprehend the nature of these sources and assess their broader implications. Using AT2018cow as a prototype, we investigate high-energy neutrino and cosmic-ray production from FBOTs and the more general class of superluminous supernovae (SLSNe). By calculating the interaction of cosmic rays and the time-evolving radiation field and baryon background, we find that particles accelerated in the magnetar wind may escape the ejecta at ultrahigh energies. The predicted high-energy neutrino fluence from AT2018cow is below the sensitivity of the IceCube Observatory, and estimates of the cosmically integrated neutrino flux from FBOTs are consistent with the extreme-high-energy upper limits posed by IceCube. High-energy γ rays exceeding GeV energies are obscured for the first months to years by thermal photons in the magnetar nebula, but are potentially observable at later times. Given their potentially higher volumetric rate compared to other engine-powered transients (e.g., SLSNe and gamma-ray bursts), we conclude that FBOTs are favorable targets for current and next-generation multimessenger observatories.

**Key words:** astroparticle physics – supernovae: general

## 1. Introduction

Time-domain optical surveys have in recent years discovered a growing number of luminous, rare, and rapidly evolving extragalactic transients, with visual light curves that rise and decay on a timescale of only a few days (hereafter “fast blue optical transients,” or FBOTs; e.g., Drout et al. 2013, 2014; Arcavi et al. 2016; Rest et al. 2018). Although supernova light curves are often powered by the radioactive decay of $^{56}$Ni, such an energy source is ruled out for FBOTs (even if the entire ejecta mass were composed of nickel, the resulting optical luminosity would be insufficient to explain the observations). Instead, FBOTs must be powered by an additional, centrally concentrated energy source or sources, either due to shock interaction between the ejecta and a massive circumstellar shell (e.g., Kleiser et al. 2018), or from a central compact object, such as a rapidly spinning magnetar or accreting black hole.

Indeed, newly formed magnetars with millisecond rotation periods have widely been invoked as engines responsible for powering the broader class of “superluminous supernovae” (SLSNe; e.g., Kasen & Bildsten 2010; Woosley 2010). Within this framework, the faster evolution of FBOTs requires stellar explosions with lower ejecta masses than of typical SLSNe. Such low ejecta masses $\lesssim 0.1-1 M_\odot$ could plausibly arise in a number of scenarios: explosions of ultra-stripped progenitor stars (e.g., Tauris et al. 2015); electron capture supernovae (e.g., Moriya & Eldridge 2016); initially “failed” explosions of very massive stars such as blue supergiants (in which most of

the stellar mass collapses into the central black hole, with only a tiny ejected fraction; e.g., Fernández et al. 2018); the accretion-induced collapse of a rotating white dwarf (e.g., Metzger et al. 2008; Lyutikov & Toonen 2018); and magnetars created as the long-lived stable remnants of neutron star binary mergers (e.g., Yu et al. 2013; Metzger & Piro 2014; Murase et al. 2018). On the other hand, compared to normal SLSNe or neutron star mergers, FBOTs may be relatively common, occurring at rates of up to $\sim 4\%–7\%$ that of the core-collapse supernova rate (Drout et al. 2014). FBOTs also seem to occur exclusively in star-forming galaxies, supporting an origin associated with the deaths of massive stars (Drout et al. 2014). Rapidly spinning pulsars in ultra-stripped supernovae have also been considered to explain some of the FBOTs (Hotokezaka et al. 2017).

Until recently, all FBOTs were discovered at high redshift ($z>0.1$) and were not identified in real time, thereby precluding key spectroscopic or multiwavelength follow-up observations. This situation changed dramatically with the recent discovery of the optical transient AT2018cow (Prentice et al. 2018; Smartt et al. 2018), which shared many characteristics with FBOTs but took place at a distance of only $\approx 60$ Mpc ($z=0.0141$). AT2018cow reached a peak visual luminosity of $L_{opt} \gtrsim 10^{44}$ erg s$^{-1}$ within only a few days, before fading approximately as a power law $L_{opt} \propto t^{-\alpha}$ with $\alpha \approx 2.5$ (e.g., Perley et al. 2018). The rapid rise time was indicative of an explosion with a low ejecta mass, $M_{ej} \lesssim 0.5 M_\odot$, and a high ejecta velocity, $v_{ej} \gtrsim 0.1c$. The latter
also revealed itself through the self-absorbed radio synchrotron emission from the source (de Ugarte Postigo et al. 2018; Ho et al. 2018; Margutti et al. 2018), which was likely produced by the same fast ejecta undergoing shock interaction with a dense surrounding medium. The early spectra were nearly featureless (also consistent with a high photospheric expansion speed), but after a few weeks narrow helium and hydrogen emission lines appeared (Perley et al. 2018), providing clues to the nature of the progenitor star (e.g., Margutti et al. 2018). The persistent blue colors of the transient (Perley et al. 2018), again unlike the behavior of normal supernovae, provides additional evidence that a central ionizing energy source is responsible for indirectly powering the optical light through absorption and “reprocessing” (Margutti et al. 2018).

Perhaps the most remarkable thing about AT2018cow is the direct visibility of the engine itself. AT2018cow was accompanied by bright X-ray emission (Kuin et al. 2018; Margutti et al. 2018; Rivera Sandoval et al. 2018) with a power-law spectrum and luminosity comparable to that of the optical emission, but showing rapid variability on timescales at least as short as a few days. The fact that X-rays appear to be visible from the same engine responsible for powering the optical luminosity provides evidence that the ejecta is globally asymmetric, such that the observed X-rays preferentially escape along low-density polar channels (see Margutti et al. 2018 for additional evidence and discussion). Prentice et al. (2018) proposed that the engine behind AT2018cow is a millisecond magnetar with a dipole field strength of $B_d \sim 10^{15}$ G and an initial spin period of $P_i \sim 10$ ms. However, other magnetar parameters are consistent with the data and alternative central engine models remain viable, such as the collapse of a blue supergiant star to an accreting black hole, or a deeply embedded shock arising from circumstellar interaction (Margutti et al. 2018).

If AT2018cow is powered by a central compact object, particularly a millisecond magnetar, then it could also be a source of high-energy charged particles (Blasi et al. 2000; Arons 2003; Fang et al. 2012), high-energy neutrinos (Murase et al. 2009, 2014; Fang et al. 2014; Fang 2015), high-energy gamma-rays (Kotera et al. 2013; Murase et al. 2015; Renault-Tinacci et al. 2018), and gravitational waves (Stella et al. 2005; Dall’Osso et al. 2009; Kashiyama et al. 2016). Such a scenario has also been discussed in the context of pulsar-powered supernova, and the similar low ejecta mass context of binary neutron star mergers (Piro & Kollmeier 2016; Fang & Metzger 2017). Indeed, two IceCube neutrino track events were reported in spatial coincidence with AT2018cow during a 3.5 day period following the optical discovery. The events are consistent with the rate of atmospheric neutrino background (Blaufuss 2018). While we will show that it is unlikely these neutrinos are physically associated with the transient, they nevertheless motivate a more thorough study of FBOTs as sources of ultrahigh energy cosmic rays (UHECRs) and neutrinos, particularly given their high volumetric rate as compared to other previously considered UHECR sites, such as gamma-ray bursts and SLSNe.

In this paper, we follow the procedure outlined in Fang & Metzger (2017) to calculate the neutrino emission from millisecond magnetars embedded in low ejecta mass explosions, focusing on parameters motivated by the watershed event, AT2018cow. In Section 2 we describe our model for the radiation and hadronic background. In Section 3 we describe

the acceleration and escape of UHECRs. In Section 4 we discuss neutrino production and address the neutrino detection prospects, both of AT2018cow and the cosmic background of FBOTs. Throughout this paper we adopt the shorthand notation $q_x \equiv q/10^x$ cgs.

## 2. Radiation and Hadron Background

The energetic compact object is embedded within a rapidly expanding ejecta shell filled with baryons and photons. Cosmic rays, which are accelerated close to the compact object, will lose energy and produce neutrinos when traveling through this dense medium. We follow Metzger & Piro (2014) to calculate the density of ejecta baryons, as well as thermal and nonthermal photons. This simplified model assumes the ejecta to be homogeneous and spherically symmetric. While this is reasonable for estimating the environment for cosmic-ray interaction, modeling of AT2018cow indicates that the true ejecta structure is aspherical and thus is more complicated in detail (e.g., Liu et al. 2018; Margutti et al. 2018).

The evolution of nonthermal radiation $E_{\text{nth}}$, thermal radiation $E_{\text{th}}$, magnetic energy $E_B$, ejecta radius $R_{\text{ej}}$, and velocity $v_{\text{ej}}$ are described by a coupled set of differential equations:

$$ \frac{dE_{\text{nth}}}{dt} = L_{\text{sd}} - \frac{E_{\text{nth}}}{R_n} \frac{dR_n}{dt} - \frac{E_{\text{nth}}}{t_d^n}, \quad (1) $$

$$ \frac{dE_{\text{th}}}{dt} = \frac{E_{\text{nth}}}{t_d^n} - \frac{E_{\text{th}}}{R_{\text{ej}}} \frac{dR_{\text{ej}}}{dt} - \frac{E_{\text{th}}}{t_d^{\text{ej}}}, \quad (2) $$

$$ \frac{dE_B}{dt} = \epsilon_B L_{\text{sd}} - \frac{E_B}{R_n} \frac{dR_n}{dt} \quad (3) $$

$$ M_{\text{ej}} v_{\text{ej}} \frac{dv_{\text{ej}}}{dt} = \frac{E_{\text{th}} + E_{\text{nth}} + E_B}{R_{\text{ej}}} \frac{dR_{\text{ej}}}{dt}, \quad (4) $$

$$ \frac{dR_{\text{ej}}}{dt} = v_{\text{ej}}. \quad (5) $$

The system is powered by the dipole spin-down of the magnetar (Ostriker & Gunn 1969), which injects a luminosity given by

$$ L_{\text{sd}} = L_{\text{sd,0}} \left(1 + \frac{t}{t_{\text{sd}}}\right)^{-2}, \quad (6) $$

where $L_{\text{sd,0}} = 4\mu^2\Omega^4/9c^3 = 2.6 \times 10^{45} B_{d,15}^2 P_{i,-2}^{-4} \text{ erg s}^{-1}$ is the initial spin-down power,<fnref n="10" /> $\mu = B_d R_*^3$ is the magnetic moment of a pulsar with dipole field $B_d$ and stellar radius $R_*$, and $P_i = 10^{-2} P_{i,-2}$ is the initial spin period. The spin-down time, $t_{\text{sd}}$, which characterizes the timescale over which most of the initial rotational energy, $E_{\text{rot}} = I\Omega_i^2/2$, is released is given by

$$ t_{\text{sd}} \equiv \frac{E_{\text{rot}}}{L_{\text{sd,0}}} = 1.0 P_{i,-2}^2 B_{d,15}^{-2} \text{ day}. \quad (7) $$

The pulsar moment of inertia is given by $I = 2 M_* R_*^2 / 5$, where $R_* = 10$ km and $M_* = 1.4 M_\odot$ are its radius and mass, respectively. At times $t \gg t_{\text{sd}}$, $L_{\text{sd}} \propto t^{-2}$, corresponding to a pulsar braking index of 3 (Ostriker & Gunn 1969). This decay

<footnote n="10">¹⁰ We use the vacuum dipole formula, following the precedent of Özel et al. (2010). In reality the spin-down is instead best described as a force-free MHD wind (e.g., Spitkovsky 2006), but the results are similar to within a factor of a few (e.g., Kashiyama et al. 2016).</footnote>
rate is moderately shallower than that of the bolometric luminosity of AT2018cow ($L \propto t^{-2.5}$; Margutti et al. 2018; Perley et al. 2018); however, a smaller braking index, or additional mechanisms such as the interaction of fall-back from the explosion with the magnetosphere (Metzger et al. 2018), can result in a steeper decay. We also consider magnetic energy losses due to the expansion, which can be present, e.g., if the magnetic field in the nebula is turbulent.

Finally, internal magnetic fields at the scale of $10^{16}$ G can be generated in the core of newborn magnetars if their initial spin period is on the order of a few milliseconds (Stella et al. 2005; Shore et al. 2009). The anisotropic pressure from the toroidal B-field leads to an ellipticity of $\epsilon_{B,GW} \sim 6.4 \times 10^{-4} B_{i,16.3}^2$ for an average toroidal field strength $B_t$ (Stella et al. 2005). Such a quadrupole moment will cause a spin-down due to gravitational wave emission at the rate $(d\Omega/dt)_{GW} = (32/5)(G/c^5)I\epsilon_{B,GW}^2 \Omega^5$. In comparison, the spin-down rate due to the dipole field is $(d\Omega/dt)_d = \mu^2/(6Ic^3)\Omega^3$. The ratio of the two rates is $(d\Omega/dt)_{GW}/(d\Omega/dt)_d = 5\% P_{i,-3}^{-2} B_{d,15}^{-2} B_{i,16.3}^2$. This fraction could be higher depending on the equation of state of the stellar interior, but should not dominate the energy losses. We thus ignore the quadrupole term in this work, and note that it could impact the spin-down at the level of a few to tens of percent. We note that an associated gravitational wave flux from such a source should not be detectable with current instruments (Arons 2003; Kotera 2011) except that it is extremely close (Kashiyama et al. 2016), but a stochastic gravitational wave signature from a population of magnetar-powered transients could lead to detection with future instruments (Kotera 2011).

A significant portion of the magnetar's rotational energy is ultimately used to accelerate the ejecta, as captured by Equation (4). The mean velocity of the ejecta after time $t$ is thus approximately given by

$$v_{ej} \approx \left( \frac{2 \int_{0}^{t} L_{sd} dt'}{M_{ej}} + v_{0}^{2} \right)^{1/2}, \quad (8)$$

where $v_0 \lesssim 0.05c$ is the initial ejecta velocity from the supernova explosion (Prentice et al. 2018) and $M_{ej} = 0.3 M_{ej,-0.5} M_\odot$ is the ejecta mass. The nebular radius $R_n$ is taken to be a fixed fraction of the mean ejecta radius $R_n = R_{ej}/5$, where $R_{ej} \approx v_{ej}t$, consistent with observations indicating that the photosphere velocity of AT2018cow declines from $\gtrsim 0.05c$ at 4 days to $\sim 0.01-0.02c$ over its first two weeks of evolution (Prentice et al. 2018).

At the wind termination shock, the majority of the spin-down power is converted into nonthermal radiation (e.g., Tanaka & Takahara 2010, 2013). A portion of this is then converted to thermal radiation via absorption by the ejecta walls, as described by Equations (1) and (2), respectively. The time-scales $t_d^n$ and $t_d^{ej}$ are those required for photons to diffuse radially through the nebula or ejecta shell, respectively. Specifically,

$$t_d^n \approx \frac{R_n}{c}(1 + \tau_{es}^n), \quad (9)$$

where $\tau_{es}^n = (4 Y L_{sd} \sigma_T / \pi m_e c^3 R_n)^{1/2} \approx 17 B_{15} \beta_{ej,-1.3}^{-1/2} t_5^{-1/2} P_{i,-2}^{-2}$ is the optical depth due to Thomson scattering by

electron positron pairs. The latter are generated by a pair cascade process in the nebula due to its high compactness (Metzger & Piro 2014; Fang & Metzger 2017), where $Y \approx 0.1$ is the pair multiplicity (fraction of the spin-down power converted into pair rest mass). Once the photon diffusion time becomes less than the nebula expansion time $t_d^n = R_n / \beta_n c$, photons are able to travel freely through the nebula; this occurs after a characteristic time

$$t_d^n \approx 0.03 B_{15}^2 P_{i,-2}^{-4} \beta_{ej,-1.3} \text{ day.} \quad (10)$$

Similarly, photons diffuse through the ejecta on a timescale

$$t_d^{ej} \approx \frac{R_{ej}}{c}(1 + \tau_{es}^{ej}), \quad (11)$$

where $\tau_{es}^{ej} = 3 M_{ej} \kappa / 4 \pi R_{ej}^2$ is the optical depth and $\kappa \sim 0.01-0.4 \text{ cm}^2 \text{ g}^{-1}$ is the ejecta opacity at UV/visual frequencies. The free escape of photons takes place once $t_d^{ej} < R_{ej} / \beta_{ej} c$, as occurs after the time

$$t_{d,0}^{ej} = 2.1 M_{ej,-0.5}^{1/2} \left( \frac{\kappa}{0.01 \text{ cm}^2 \text{ g}^{-1}} \right)^{1/2} \beta_{-1.3}^{-1/2} \text{ days.} \quad (12)$$

Equations (1) and (2) were derived under the approximation that all of the nonthermal radiation is reprocessed into thermal radiation following its absorption by the ejecta (Metzger et al. 2014). This is justified by the observational fact that the X-ray luminosity of AT2018cow (that escaping from the engine) is less than the optical luminosity (that thermalized) over the first several weeks of evolution (e.g., Margutti et al. 2018).

The thermal luminosity (supernova light curve) is given by the loss term in Equation (2),

$$L_{\text{th}} = \frac{E_{\text{th}}}{t_d^{\text{ej}}}. \quad (13)$$

At early times ($t \ll t_d^n$), the kinetic term (PdV work) dominates the energy loss, such that $E_{\text{nth}} \sim L_{sd,0}t$ and $E_{\text{th}} \sim L_{sd,0}t / \tau_{es}^n \propto t^{3/2}$. The luminosity thus scales $L_{\text{th}} \sim E_{\text{th}}/t_d^{\text{ej}} \propto t^{5/2}$ at $t \ll t_{d,0}^{\text{ej}}$. At late times ($t \gg t_{d,0}^{\text{ej}}$), $E_{\text{nth}} \sim L_{sd,0}t_{sd}/(t/t_{sd})^{-1}$ and $E_{\text{th}} \sim E_{\text{nth}} t_d^{\text{ej}}/t_d^n \propto t^{-1}$. As a result, $L_{\text{th}}$ declines $\propto t^{-2}$, i.e., following the magnetar spin-down luminosity.

Figure 1 compares the time-dependent thermal luminosity calculated from Equation (13) to the observed optical light curve of AT2018cow (as integrated from the UV to near-IR wavelength bands). We consider two models for the properties of the central magnetar, referred to hereafter as Case I and Case II, respectively (see Table 1). Case I corresponds to a magnetar with dipole field $B_d = 10^{15}$ G and initial spin period $P_i = 10$ ms. This model reproduces well the observed optical light curve, as shown in Figure 1, and is consistent with previous fits to the magnetar model (Margutti et al. 2018; Prentice et al. 2018). We use parameters from Case I in our analytic estimates below. Case II corresponds to a magnetar with a somewhat stronger field $B_d = 2 \times 10^{15}$ G and larger rotational energy $P_i = 2$ ms. The light curve in this case, as shown by a dashed line in Figure 1, is also roughly consistent with that of AT2018cow (though decaying a bit too quickly at late times).
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Luminosity of thermal emission from the magnetar-powered transient (as calculated from Equation (13)) in comparison to the UVONIR light curve of AT2018cow (Prentice et al. 2018). Model results are shown separately for a magnetar with an initial rotation period $P_i = 10$ ms, surface magnetic field $B_d = 10^{15}$ G, ejecta opacity $\kappa = 0.01$ cm$^2$ g$^{-1}$ (blue solid line), and $P_i = 2$ ms, $B_d = 2 \times 10^{15}$ G, $\kappa = 0.1$ cm$^2$ g$^{-1}$ (orange-dashed line). In both cases, the ejecta mass is $M_{ej} = 0.3 M_\odot$ and the velocity of the magnetar-inflated nebula is taken to be five times lower than the mean ejecta velocity (to match the range of photosphere velocities observed in AT2018cow).</figcaption></figure>

<table><caption>**Table 1**

Magnetar Parameters of Fiducial Models</caption><thead><tr><th></th><th>P<sub>i</sub> (ms)</th><th>B<sub>d</sub> (G)</th></tr></thead><tbody><tr><td>Case I</td><td>10</td><td>10<sup>15</sup></td></tr><tr><td>Case II</td><td>2</td><td>2 &times; 10<sup>15</sup></td></tr></tbody></table>

The magnetic field in the magnetar-wind nebula can be estimated as

$$B_n \approx \left( \frac{6\epsilon_B L_{sd} t}{R_n^3} \right)^{1/2} \approx 720 \epsilon_{B,-4}^{1/2} B_{15} P_{i,-2}^{-2} \beta_{n,-2}^{-3/2} t_5^{-2} \text{ G}, \quad (14)$$

where a fraction $\epsilon_B$ of the spin-down energy $\sim L_{sd}t$ is assumed to be placed in the magnetic energy of the nebula, $E_B = (B_n^2/8\pi)V_n$ (e.g., at the wind termination shock), where $V_n = 4\pi R_n^3/3$ is the nebula volume. Typical values $\epsilon_B \sim 10^{-4}-10^{-2}$ are obtained by modeling pulsar wind nebulae (e.g., Kennel & Coroniti 1984; Torres et al. 2014).

The temperature of the ejecta, and thus that of the thermal radiation field, is approximately given by $T_{th} = (3E_{th}/4\pi a R_{ej}^3)^{1/4}$. The number density of thermal photons is then given by

$$\begin{aligned} n_{\text{th}} &\approx 16\pi\zeta(3)(kT_{\text{th}}/hc)^3 \\ &\approx 2.3 \times 10^{16} B_{15}^{3/2} P_{i,-2}^{-3} t_5^{-3/2} \beta_{ej,-1.3}^{-3/2} \text{ cm}^{-3}. \end{aligned} \quad (15)$$

The thermal emission of the hot magnetar as considered in Kotera et al. (2015) is subdominant here as the acceleration and interaction sites are distant from the star. As detailed in Section 3, the acceleration should happen beyond the light cylinder. The density of thermal photons from the magnetar is a factor of $(R_{\text{star}}/R_{\text{acc}})^2$ times lower than that near the star, where $R_{\text{star}}$ is the stellar radius and $R_{\text{acc}}$ is the radius of the acceleration site. The contribution of these photons to neutrino production is thus negligible.

The number density of nonthermal photons can likewise be estimated to be

$$\begin{aligned} n_{\text{nth}} &\sim \frac{L_{\text{nth}}}{4\pi R_n^2 c \varepsilon_{\text{min}} \ln(\varepsilon_{\text{max}}/\varepsilon_{\text{min}})} \\ &\approx 3.6 \times 10^{15} B_{15}^{3/2} P_{i,-2}^{-3} t_5^{-3/2} \beta_{n,-2}^{-3/2} \text{ cm}^{-3}, \end{aligned} \quad (16)$$

where we have assumed a flat power-law spectrum $n_{\text{nth}}(\varepsilon) \propto \varepsilon^{-2}$, extending from the energy of the thermal radiation $\varepsilon_{\text{min}} \sim 3 k_B T_{\text{th}}$ to the pair creation threshold $\varepsilon_{\text{max}} \sim 2 m_e c^2$ (Svensson 1987). The observed X-ray spectrum of AT2018cow ($n(\varepsilon) \sim \varepsilon^{-1.5}$; Margutti et al. 2018) is somewhat harder than assumed in the model. However, this difference does not critically affect our conclusions because the density of the nonthermal radiation is lower than that of the thermal photons (as well as higher in energy), and therefore are generally less important targets for neutrino production.

The baryon density of the ejecta is given by

$$n_p \approx \frac{3M_{ej}}{4\pi R_{ej}^3 m_p} \approx 2.5 \times 10^{13} M_{ej, -0.5} t_5^{-3} \beta_{ej, -1.3}^{-3} \text{ cm}^{-3}. \quad (17)$$

While $n_p$ scales with time as $\propto t^{-3}$, the evolution of $n_{\text{th}}$ changes at $t_{d,0}^n$ and $t_{d,0}^{ej}$, introducing features to the cosmic-ray interaction, as discussed below.

# 3. Acceleration and Escape of UHECRs

## 3.1. Cosmic-Ray Injection

Pulsars and magnetars offer promising sites for particle acceleration. In general, charged particles may tap the open field voltage and gain an energy (Arons 2003)

$$\begin{aligned} E_{\text{CR}} &\approx \eta Z e \Phi_{\text{mag}} \\ &\approx 1.3 \times 10^{19} Z \eta_{-1} B_{15} P_{i,-2}^{-2} \left(1 + \frac{t}{t_{\text{sd}}}\right)^{-1} \text{eV}, \end{aligned} \quad (18)$$

where $\eta = 0.1 \eta_{-1}$ is the acceleration efficiency and Z is the particle charge. For simplicity, we assume cosmic rays of proton composition. The effects of a heavier composition on neutrino production are discussed in Section 5.

Assuming that ions follow the Goldreich-Julian charge density, the rate of cosmic-ray injection from the magnetic polar cap is given by $\dot{N} = \mu\Omega^2/Zec$ (Arons 2003). Equation (18) then implies the energy spectrum of accelerated particles is given by

$$\frac{dN}{dE} = \frac{9 c^2 I}{8 e \mu E} = 2 \times 10^{42} E^{-1} B_{15}^{-1}. \quad (19)$$

The specific acceleration mechanism of high-energy particles in the pulsar magnetosphere is still debated (e.g., Cerutti & Beloborodov 2017). For UHECRs, Arons (2003) hypothesized that particles achieve their energy in the relativistic wind through surf-riding acceleration at a radius $10^3-10^4$ times the light cylinder. The possibility of wake-field acceleration has also been discussed (e.g., Murase et al. 2009; Iwamoto et al. 2017). Recent particle-in-cell simulations (Philippov & Spitkovsky 2018) support a significant fraction of the accelerated particle energy flux being carried by ions extracted from the stellar surface. The ion energy may reach up to 10%–30% of the polar cap voltage. As the misalignment of the rotational and magnetic axes impact the current distribution, the spectrum also
depends on the inclination angle of the pulsar. Alternatively, if the acceleration happens close to the star, curvature radiation will be nonnegligible and limit the maximum energy. However, particles could be picked up and reaccelerated by the pulsar wind. If the pair multiplicity is low, a significant wind power would go into ions and still allow a UHECR production (Kotera et al. 2015).

## 3.2. Interaction and Cooling Rates

Accelerated cosmic rays interact with the background of baryons or photons, producing charged pions ($p + p(\gamma) \rightarrow \pi^{\pm}$), which decay into muons and neutrinos ($\pi^{\pm} \rightarrow \mu^{\pm} + \nu_{\mu}(\bar{\nu}_{\mu})$). The muons further decay into electrons and neutrinos ($\mu^{\pm} \rightarrow e^{\pm} + \nu_{e} + \bar{\nu}_{e} + \nu_{e}(\bar{\nu}_{e})$). Each step of this process may be suppressed by particle cooling.

A proton with Lorentz factor $\gamma_p$ interacts with the photon field of spectrum $n(\varepsilon) = dn/d\varepsilon$ on a characteristic timescale given by

$$t_{p\gamma,\text{int}}^{-1} = \frac{c}{2\gamma_p^2} \int_0^\infty d\varepsilon \frac{n(\varepsilon)}{\varepsilon^2} \int_0^{2\gamma_p\varepsilon} d\varepsilon' \varepsilon' \sigma_{p\gamma}(\varepsilon'), \quad (20)$$

where $\sigma_{p\gamma}$ is the cross section of photopion production. The cooling time is $t_{p\gamma} = t_{p\gamma,\text{int}}/\kappa_{p\gamma} \sim (n_{\gamma} \sigma_{p\gamma} \kappa_{p\gamma} c)^{-1}$, where $\sigma_{p\gamma} \kappa_{p\gamma} \sim 10^{-28} \text{ cm}^{-2}$ is the inelastic component of the p\gamma interaction cross section. The rate for hadronuclear interaction is likewise given by

$$t_{pp}^{-1} = n_p \sigma_{pp} \kappa_{pp} c, \quad (21)$$

where $\sigma_{pp} \sim 10^{-25} \text{ cm}^2$ (at energies $\sim 10^{18}$ eV) and $\kappa_{pp} \sim 0.5$ (e.g., Eidelman et al. 2004). The total interaction rate of protons, due to both pp and p\gamma processes, is then $t_p^{-1} = t_{pp}^{-1} + t_{p\gamma}^{-1}$.

The gyroradii of protons in the magnetic field of the nebula, $r_L$, is comparable to or larger than the nebula size,

$$\frac{r_L}{R_n} \approx 1.7 \eta_{-1} \epsilon_{B,-4}^{-1/2} \beta_{n,-2}^{1/2}, \quad t \gg t_{\text{sd}}. \quad (22)$$

Synchrotron or adiabatic losses of cosmic rays crossing the nebula are thus generally only marginally important. Inverse Compton cooling at such high energies is furthermore suppressed and negligible due to the Klein–Nishina effect.

Protons travel freely when their crossing time is shorter than the interaction time, $t_{\text{cross,ej}} < t_p$, as occurs after a characteristic timescale

$$t_{p,0}^{\text{ej}} \approx 16 M_{\text{ej},-0.5}^{1/2} \beta_{\text{ej},-1.3}^{-1}, \quad \text{days}, \quad (23)$$

where we have assumed that at late times pp interactions dominate over p\gamma interactions (however, note that p\gamma interactions are included in our full numerical calculations).

Figure 2 shows the proton interaction timescales as a function of time since explosion. In Case I, the dominant process for protons is hadronuclear interaction. In Case II, the photopion production with the thermal photon background becomes important at late times, $t \gtrsim t_{d,0}^n$. In both cases, the ejecta becomes optically thin at roughly one week, after which time the accelerated protons freely escape to infinity.

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

<figure><img src="image_5.png" /><figcaption>Figure 2. Rates of interaction and cooling processes of protons as a function of time since explosion. The photopion production rates are computed using Equation (20) for the thermal (blue solid) and nonthermal background (light blue solid) density computed in Section 2. The hadron interaction rates are from Equations (21) (red solid). The proton energy is determined by Equation (18). At all times, the particle’s Lamor radius (black dashed) is comparable to the size of the nebula (black-dotted), meaning that both synchrotron and adiabatic losses of the protons are negligible. The magnetar parameters are $Pi = 10$ ms and $B_d = 10^{15}$ G in Case I (top) and $Pi = 2$ ms and $B_d = 2 \times 10^{15}$ G in Case II (bottom).</figcaption></figure>

## 3.3. Ultrahigh Energy Cosmic-Ray Production

Figure 3 shows the maximum energy of the accelerated UHECR protons at the time when they escape the ejecta freely. The shaded area corresponds to an event with a low ejecta mass $M_{ej} = 0.3 M_\odot$, similar to that inferred for FBOTs such as AT2018cow. Depending on the magnetic field of the magnetar, an AT2018cow-like event with $Pi \lesssim 10$ ms is a promising source for cosmic rays of energy $\gtrsim 10^{18}$ eV (see also Piro & Kollmeier 2016). If the magnetar wind is composed of nuclei with charge Z instead of protons ($Z=1$), then their maximum energy is a factor of Z times higher than for protons and thus in the range of the highest energy cosmic rays observed by Auger (The Pierre Auger Collaboration et al. 2017) and TA (Matthews 2018), as detailed in Fang et al. (2012).

For comparison, millisecond magnetars born in the normally considered class of SLSNe, for which the ejecta mass is
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Maximum energy of cosmic rays at the time when the ejecta first become optically thin to pp and pγ interactions as a function of the surface magnetic field B<sub>d</sub> and birth period Pi of the neutron star. Magnetars with Pi < 10 ms are promising sources of UHE protons in transients with low ejecta masses M<sub>ej</sub> ≲ 0.3 M<sub>☉</sub> such as AT2018-cow (the parameters of our fiducial models, Case I and II, are denoted as pink crosses). By contrast, millisecond magnetars in standard SLSNe with larger ejecta mass M<sub>ej</sub> = 10 M<sub>☉</sub> require values of (B<sub>d</sub>, Pi) to the left of the blue dashed curve to allow the escape of protons with E<sub>p</sub> ≥ 10<sup>18</sup> eV. White points show magnetar parameters fit to a sample of SLSNe light curves by Nicholl et al. (2017). All calculations assume values of η = 0.1 and κ = 0.1 cm<sup>2</sup> g<sup>-1</sup>.</figcaption></figure>

typically much higher (M<sub>ej</sub> ~ 10 M<sub>☉</sub>), only allow the escape of UHE protons for magnetar parameters to the left of the blue dashed curve. This is because massive ejecta shells require longer to become optically thin, delaying the time of cosmic-ray escape and thus reducing the pulsar voltage at this epoch. White data points show magnetar parameters fit to the sample of SLSNe in Nicholl et al. (2017), roughly one-third of which appear to be promising UHECR sources.

We conclude this section by estimating the UHECR energy budget of FBOTs similar to AT2018cow. Assuming Pi = 2–10 ms and that a fraction η = 0.1 of the magnetar’s rotational energy goes into cosmic rays, the UHECR energy yield per event is (0.2–5) × 10<sup>50</sup> erg. The total rate of FBOTs is ~4%–7% of the core-collapse supernovae rate (Drout et al. 2014). However, the rate of the most luminous members of this class like AT2018cow is likely lower (M. Drout 2019, private communication), perhaps only ~0.5%–1% of the core-collapse supernova (CCSN) rate, corresponding to an estimated volumetric rate of AT2018cow-like events of 600–1200 Gpc<sup>-3</sup> yr<sup>-1</sup>. Using an FBOT rate with ~1% of the CCSN rate, the integrated cosmic-ray luminosity density from AT2018cow-like FBOTs is thus roughly estimated to be ~ (0.3–6) × 10<sup>44</sup> erg Mpc<sup>-3</sup> yr<sup>-1</sup>. This is comparable to the UHECR luminosity density, which is in the order of 10<sup>44</sup> erg Mpc<sup>-3</sup> yr<sup>-1</sup> (e.g., Murase & Fukugita 2018).

# 4. High-energy Neutrino Production

## 4.1. IceCube Observation of AT2018cow

Two IceCube neutrino track events with ~2° angular resolution were found in spatial coincidence with AT2018cow during a 3.5 day period between the last nondetection and the discovery, corresponding to a 1.8σ chance coincidence. Assuming an E<sup>-2</sup> spectrum, a time-integrated ν<sub>μ</sub> flux upper

limit (E<sub>ν</sub><sup>2</sup> J<sub>ν</sub>)<sub>UL</sub><sup>IC</sup> = 6.1 × 10<sup>-2</sup> GeV cm<sup>-2</sup> is found for this observation period at 90% CL.<sup>11</sup>

Neutrinos with TeV–PeV energies lie in the best sensitivity window for the IceCube Observatory, and are predicted to arrive within the first few hours in our model. Using the neutrino effective area of the IceCube Observatory with its complete configuration of 86 string detectors<sup>12</sup> (Aartsen et al. 2014), we compute the detector sensitivity at the decl. of AT2018cow as a function of energy. It is shown by the gray dashed curve in Figure 6. Point sources with fluxes above the curve, assuming that they follow an E<sup>-2</sup> spectrum over a half decade in energy, are excluded at 90% C.L.

## 4.2. Competitive Cooling and Decay of Pions and Muons

Charged pions and muons decay into neutrinos following their lab-frame lifetimes, γ<sub>x</sub>τ<sub>x</sub>, where x denotes either π or μ, and τ<sub>π±</sub> = 2.6 × 10<sup>-8</sup> s and τ<sub>μ</sub> = 2.2 × 10<sup>-6</sup> s are the rest-frame lifetimes. These lifetimes can be long enough to allow cooling processes to occur. Mesons and muons interact with the hadron background at a rate given by

$$t_{\text{xp}}^{-1} = n_p \sigma_{\text{xp}} \kappa_{\text{xp}} c. \qquad (24)$$

Mesons and muons can also interact with the photon background, but as in the case of protons the smaller πγ cross section renders the πp interaction more important, especially at early times.

Figure 4 compares the lifetime of a pion (muon), which is assumed to carry away 20% (15%) of the energy of the parent proton, to the cooling time due to πp interactions. This shows that most pions interact with the ejecta baryons, rather than decay, in the first day or longer.

Pion–proton interactions produce lower-energy pions and protons, which then undergo further interactions with background particles. Eventually, higher-order pions reach sufficiently low energies to enable their decay into neutrinos. This happens once t<sub>πp</sub> > γ<sub>π</sub>τ<sub>π</sub>, as occurs after a time

$$t_{\pi p,0} \approx 2.2 \eta_{-1}^{1/4} B_{15}^{-1/4} \beta_{ej,-1.3}^{-3/4} M_{ej,-0.5}^{1/4} \text{ days.} \qquad (25)$$

This effect introduces a break at time t in the neutrino spectrum at a characteristic energy (Murase et al. 2009)

$$E_{\nu,b} \approx 59.7 M_{ej,-0.5}^{-1} t_5^3 \beta_{ej,-1.3}^3 \text{ PeV.} \qquad (26)$$

In addition to the higher-order products, multiple pions and other types of mesons such as kaons may be produced from each pp and pγ interaction. We take these into account in our numerical simulations presented in Section 4.3.

Finally, similar to protons, the gyroradii of pions and muons are comparable to the nebula size. Synchrotron cooling thus becomes important if r<sub>L</sub>/R<sub>n</sub> ≪ 1, as occurs for high nebular magnetization ϵ<sub>B</sub> ≳ 10<sup>-2</sup>.

## 4.3. Numerical Procedure

We compute the neutrino flux from magnetar-powered supernovae with low ejecta masses according to the following numerical procedure. At each time t following the explosion, the ejecta radius R<sub>ej</sub>, temperature of the thermal background T<sub>th</sub>, number density of ejecta baryons n<sub>b</sub>, and strength of the

<footnote n="11"><sup>11</sup> http://www.astronomerstelegram.org/?read=11785</footnote>

<footnote n="12"><sup>12</sup> https://icecube.wisc.edu/science/data/PS-IC86-2011</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Interaction and decay rate of pions and muons as a function of time since explosion. The pion and muon are assumed to carry 20% and 15%, respectively, of the proton energy determined by Equation (18). Pions undergo severe $\pi p$ losses at early times, but are free to decay after a day or so. Due to their longer rest-frame lifetimes, muons experience strong cooling for a longer duration than pions. Like protons, the gyroradii of pions and muons in the magnetic field of the nebula are comparable to the nebula size and thus these particles are not subject to strong synchrotron losses unless the nebula magnetization is high, $\epsilon_B \gtrsim 10^{-2}$.</figcaption></figure>

nebular magnetic field $B_n$, are calculated from Equations (1)–(5). Protons of energy $E_p$ determined by Equation (18) are injected into the ejecta. Their interaction with the baryon ejecta is calculated using a Monte Carlo approach employing the hadronic interaction model EPOS-LHC (Pierog et al. 2015) as in Fang et al. (2012). Photomeson interactions between cosmic-ray protons and the thermal background is computed based on SOPHIA (Mücke et al. 2000) through CRPropa 3 (Alves Batista et al. 2016).

As we are particularly interested in the earliest phases of the transient (near the time of the observed neutrino coincidence), we take a different approach in treating the pion–proton ($\pi p$) interaction from previous works (e.g., Murase et al. 2009; Fang et al. 2016). The suppression of neutrino production due to $\pi p$ interaction is usually described by a suppression factor $f_{sup} \equiv \min(t_{cool}/\gamma_\pi \tau_\pi, 1)$, with $\gamma_\pi$ and $\tau_\pi$ being the Lorentz factor and the rest-frame lifetime of the charged pion. This analytical approach, however, misses the secondary and higher-

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Distribution of products from the interaction of a pion of initial energy $10^{19}$ eV and a proton at rest, including charged pions $\pi^\pm$, gamma-rays $\gamma$ from $\pi^0$ decays, kaons $K^\pm$, protons $p$, and neutrons $n$. The fraction of the energy carried by each group is indicated in the legend.</figcaption></figure>

order pions produced by the interaction, which have lower energy and thus decay into neutrinos more easily than their parents. This is a secondary effect, but can play an important role in neutrino production in dense environments.

For reference, Figure 5 shows the distribution of outgoing particles from the interaction of a pion of energy $E = 10^{19}$ eV with a proton at rest as calculated using the EPOS-LHC model. Roughly half of the energy of the primary pion is carried away by charged mesons, while 15% is in the form of baryons. These products will ultimately generate neutrinos.

To account for neutrinos produced by secondary pions, we record all charged meson and baryon products from the pp and $p\gamma$ interaction. At any time, each of these intermediate products may undergo one of four processes, depending on the background densities at the current step: (1) cooling by $\pi p$ interaction; (2) cooling by synchrotron radiation; (3) decay into a muon and a muon neutrino (with a 100% branching ratio for $\pi^\pm$ and a 63% branching ratio for $K^\pm$); and (4) free propagation. The products of $\pi p$ interaction are computed using the EPOS-LHC model (Pierog et al. 2015). All secondaries and their higher-order products are tracked until either their energy falls below TeV (where the atmospheric background dominates over astrophysical sources), or they escape the source without further interaction.

### 4.4. Neutrino Fluence of AT2018cow

Figure 6 presents the neutrino fluence in the two fiducial models. In Case I, we show results in time intervals normalized to the magnetar spin-down time $t_{sd}$ (roughly one day). Neutrino production is low at both early times $t \ll t_{sd}$, when the bulk of the cosmic rays have yet to be injected, and at late times $t_{pp} > t_{cross,ej}$, when pp interactions become inefficient. Most neutrinos are generated around the time $t_{sd}$, when the magnetar releases most of its rotational energy and a sufficiently dense baryon background still exists for pion production. The turnover in the neutrino spectrum is determined at early times by the break energy $E_{\nu,b}$ (Equation (26)), while at late times the break is determined by the maximum cosmic-ray energy $E_{CR}$ (Equation (18)). The spectral index before the break is a convolution of the energy distribution of charged pions from
<figure><img src="image_9.png" /><figcaption>**Figure 6.** All-flavor fluence of high-energy neutrinos from a magnetar-powered low ejecta mass supernova on timescales ranging from one hour to ~50 days after the explosion. For comparison, the time-integrated upper limit on $\nu_{\mu}$ flux ($E_{\nu}^{\text{UL}}/E_{\nu}^{\text{C}}$) placed by the IceCube Observatory is shown (dotted line), which corresponds to a 1.8σ chance coincidence of neutrino events with the transient AT2018cow in 3.5 days. In addition, the IceCube sensitivity (90% C.L.) for point sources at $\delta = 22^\circ$ is shown as a dashed curve, assuming an $E^{-2}$ spectrum over a half decade in energy (Aartsen et al. 2014). The parameters of the magnetar model are initial spin period $Pi = 10$ ms, surface magnetic field $B_s = 10^{15}$ G (Case I; top), and $Pi = 2$ ms, $B_s = 2 \times 10^{15}$ G (Case II; bottom).</figcaption></figure>

the pp and pγ interaction (similar to that from $\pi p$ as shown in Figure 5) and the history of particle injection, as described by Equation (19).

The bottom panel of Figure 6 shows our results for Case II ($Pi = 2$ ms, $B_s = 2 \times 10^{15}$ G). The spin-down time in this case is much shorter $t_{sd} = 0.01$ days, whereas the neutrino break energy exceeds 1 TeV only after times $t = 0.03 M_{\text{ej}}^{-1/3} \beta_{\text{ej}}^{-1} = 1.3$ day. As a result, most cosmic rays are injected too early to generate neutrinos in the energy range of interest. As the magnetar releases most of its energy before the environment becomes optically thin, significant TeV–PeV neutrinos are produced in the first 0.5 days due to the $\pi p$ interaction. The neutrino spectrum at the earliest epoch features two peaks; the low-energy bump is from pion decay, while the tail at high energies arises from the decay of short-lived mesons other than $\pi^{\pm}$. After about one day, the ejecta becomes sufficiently dilute that meson cooling is no longer severe and the neutrino flux becomes maximal when the suppression is not

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Luminosity of neutrinos with energies >1 TeV as a function of time since explosion, shown separately for Case I (blue solid line) and Case II (red dashed line). As indicated by black lines, the light curves approximately obey $L_{\nu} \propto t^{-2}$ at times $t_{\pi p,0} \lesssim t \lesssim t_{p,0}^{\text{ej}}$, and $L_{\nu} \propto t^{-4}$ at times $t > t_{p,0}^{\text{ej}}$.</figcaption></figure>

important. Then, the neutrino spectrum comes to resemble that in Case I. This kind of time evolution of neutrino spectra owing to meson cooling in magnetars was first found in Murase et al. (2009) and Fang et al. (2014).

Figure 7 shows the total neutrino luminosity, $L_{\nu} = \int_{E_{\nu}>\text{TeV}} E_{\nu}(dN_{\nu}/dE_{\nu})dE_{\nu}dt$, at energies >1 TeV as a function of time since explosion. In both Cases I and II the light curve obeys $L_{\nu} \propto t^{-2}$ at times $t_{\pi p,0} \lesssim t \lesssim t_{p,0}^{\text{ej}}$, similar to the optical and X-ray light curves of AT2018cow (Figure 1; e.g., Margutti et al. 2018). In the time window when the ejecta is still optically thick to cosmic rays, yet after the pion cooling no longer suppresses neutrino production, the neutrino luminosity tracks the spin-down power of the magnetar, $L_{sd} \propto t^{-2}$ (Equation (6)).

The origin of the time dependence of the neutrino luminosity is more complex at earlier and later times. At early times $t \ll t_{\pi p,0}$, both the cosmic-ray injection rate $dN/dt$ and maximum energy $E_p$ decrease in time as $(1+t/t_{sd})^{-1}$, while the suppression factor due to $\pi$ cooling rises as $f_{\text{sup}} \propto t^3(1+t/t_{sd})$. If, as in Case II, the majority of the cosmic-ray energy is injected at early times when the system is still optically thick to protons and pions, then the time evolution is also influenced by neutrinos released from pγ interaction, $\pi p$ interaction, and kaon decay. At late times $t \gg t_{p,0}^{\text{ej}}$, the effective optical depth decreases as $\tau_{\nu} \propto t^{-2}$, and thus the neutrino light curve obeys a steeper decay, $L_{\nu} \sim L_{sd}\tau_{\nu} \propto t^{-4}$ (see also Murase et al. 2009; Fang 2015).

Our predictions for the neutrino fluence from AT2018cow is well below the upper limits placed by the IceCube Observatory for both models, supporting a conclusion that the two detected events are from the background rather than of astrophysical origin. A future event otherwise similar to AT2018cow but occurring ~5 times closer (at a distance of $\lesssim 15$ Mpc), would be a promising IceCube source. Prospects may be better for future neutrino telescopes with sensitivity in the EeV energy range, enabling a direct test of FBOTs as cosmic particle accelerators.
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Peak neutrino fluence as estimated by Equation (28) for each member of the sample of SLSNe with known distances and best-fit magnetar parameters from Nicholl et al. (2017). The color of each symbol indicates the initial spin period while the symbol size indicates the strength of the magnetar’s dipole magnetic field. For reference, the dashed line shows the 5σ discovery potential of a time-dependent search of IceCube at decl. $\delta = 16^{\circ}$ with an $E^{-2}$ spectrum (Aartsen et al. 2015).</figcaption></figure>

### 4.5. Neutrino Flux of the SLSNe Population

The neutrino flux in our models peak around the time $t_{\pi p, 0}$, after which pions have sufficient time to decay into neutrinos before cooling. We therefore can estimate the peak flux as

$$ (E^2 F_{\nu})_{\text{pk}} \approx \frac{1}{2} E_p(t_{\pi p, 0}) \frac{dN_p}{dt}(t_{\pi p, 0}) \frac{1}{4\pi D^2}. \quad (27) $$

For a magnetar with $t_{sd} \ll t_{\pi p, 0}$ (as in our Case II), the peak flux is somewhat higher than this estimate due to additional contributions from secondaries.

Under the assumption that the peak neutrino flux lasts for a duration $\Delta t \sim t_{\pi p, 0}$, the corresponding fluence can be estimated by

$$ (E^2 J_{\nu})_{\text{pk}} \sim (E^2 F_{\nu})_{\text{pk}} t_{\pi p, 0}. \quad (28) $$

Figure 8 shows the peak fluence for a sample of SLSNe (Nicholl et al. 2017) estimated by Equation (28) using the magnetar parameters best-fit to the optical light curves and the distances of the sources. For comparison, we show the discovery potential of a time-dependent search by IceCube for a source at the decl. $\delta = 16^{\circ}$ (Aartsen et al. 2015) for different source time integrations $\sim t_{\pi p, 0}$. The IceCube point-source sensitivity also depends on the decl. of the source and is greatest for $\delta \approx 0$ (Aartsen et al. 2017).

### 4.6. Integrated Neutrino Flux of FBOT Population

The integrated neutrino flux over cosmological distances of magnetar-powered FBOTs with properties similar to AT2018cow is given by

$$ \Phi_{\nu} = \frac{c}{4\pi} \int_{0}^{z_{\text{max}}} dz \rho(z) \left. \frac{dN}{dE'} \right|_{E'=(1+z)E} (1+z) \left| \frac{dt}{dz} \right|, \quad (29) $$

where $z$ is the source redshift, $|dz/dt| = H_0(1+z)\sqrt{\Omega_M(1+z)^3 + \Omega_{\Lambda}}$ taking a flat $\Lambda$CDM with $\Omega_M = 0.308$ and $H_0 = 67.8 \text{ km s}^{-1} \text{ Mpc}^{-1}$ (Planck

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Integrated neutrino flux over a cosmological source population of magnetar-powered FBOTs similar to AT2018cow, assuming an event rate of 1% of the core-collapse supernovae rate in the local universe that evolves with redshift according the cosmic star formation rate (Hopkins & Beacom 2006). Results are shown separately for magnetar parameters corresponding to Cases I and II, respectively. Flux constraints from the IceCube 6 yr high-energy starting events (HESE; Kopper 2018), 6 yr muon neutrino events (Aartsen et al. 2016), and 9 yr extreme-high-energy (EHE) 90% upper limit (IceCube Collaboration et al. 2018) are shown for comparison as labeled.</figcaption></figure>

Collaboration et al. 2016). $\rho(z) = \rho_0 g(z)$ is the source birth rate at given redshift. We assume that AT2018cow-like FBOT events occur at $\sim 1\%$ of the core-collapse supernova rate ($\rho_0 = 10^{-6} \text{ Mpc}^{-3} \text{yr}^{-1}$), and track the cosmological star formation history (Hopkins & Beacom 2006) with $g(z) \propto (1+z)^{3.4}$ at $0 < z < 1$, then $g(z) \propto (1+z)^{-0.3}$ up to $1 < z < 4$, and $(1+z)^{-3.5}$ at $z > 4$. Figure 9 shows the integrated neutrino flux, calculated separately using neutrino fluences based on Case I and Case II, respectively. The peak flux is similar to that computed in Murase et al. (2009) after taking into account the difference in the source rates. This is because the peak fluence is determined by the time when pions start to decay rather than interact, and is thus insensitive to the factor $\beta_{ej}^{3/4} M_{ej}^{-1/4}$ (see Equation (25) of this work and Equation (4) of Murase et al. 2009). The flux is consistent with the extreme-high-energy upper limits from IceCube at 90% C.L. (IceCube Collaboration et al. 2018). Interestingly, if all FBOTs were as luminous as AT2018cow (i.e., if their rate was 4%–7% of the CCSN rate; Drout et al. 2014), then the integrated neutrino flux would have overproduced the IceCube limit.

### 4.7. Gamma-Ray Emission from FBOTs and SLSNe

Neutral pions created by UHE proton interactions produce high-energy $\gamma$ rays via $\pi^0 \to 2 \gamma$. In pp and $p\gamma$ channels, the energy passed by protons to $\gamma$ rays is about 2/3 and 4/3 times that given to neutrinos, respectively. High-energy gamma-rays will quickly undergo pair production with low-energy photons in the nebula ($\gamma\gamma \to e^{\pm}$). For GeV–TeV $\gamma$ rays, the peak of the $\gamma\gamma$ pair production cross section lies in the optical to soft X-ray range. The optical depth of the background photon to high-
<figure><img src="image_13.png" /><figcaption>Figure 10. Optical depth to γ rays at day 10. The dominant loss process for MeV, sub-GeV, and high-energy γ rays are Compton scattering, Bethe-Heitler pair production, and γγ pair production respectively. At the plotting time, pair production of γ rays from GeV to 0.1 TeV is dominated by nonthermal photons from the nebula, and that above 0.1 TeV is dominated by thermal photons from the ejecta.</figcaption></figure>

energy γ rays is approximately given by

$$ \tau_{\gamma\gamma} \sim n_{\text{th}} \sigma_{\gamma\gamma} R_{\text{ej}} \\ \approx 2.0 \times 10^5 P_{i,-2}^{-3/2} \beta_{\text{ej},-1.3}^{-1/2} \left(\frac{t}{\text{1 day}}\right)^{-2}, \quad (30) $$

where $\sigma_{\gamma\gamma} = 3/16 \sigma_T$ is the γγ pair production cross section and $u_{\text{th}} \sim L_{\text{sd}}/(4\pi R_{\text{ej}}^2 c)$ is the thermal photon density at late times. Thus high-energy γ rays cannot escape from the ejecta freely until after a time

$$ t_{\gamma,0} \approx 450 P_{i,-2}^{-1} \beta_{\text{ej},-1.3}^{-1/2} \text{ days.} \quad (31) $$

The pair production and inverse Compton scattering of the resulting electrons leads to an electromagnetic cascade (e.g., Metzger & Piro 2014; Metzger et al. 2014; Murase et al. 2015). Figure 10 presents the optical depth of the nebula and ejecta to γ rays at different energies around day 10. For MeV, sub-GeV, and >GeV γ rays, the dominate energy loss process is Compton scattering ($e\gamma \rightarrow e\gamma$), Bethe-Heitler pair production ($p\gamma \rightarrow p\gamma$) and γγ pair production respectively. The effective optical depth is defined as $\tau_{\text{eff}} = Y_b \sigma \kappa$, where $Y_b = 3 M_{\text{ej}}/(4\pi R_{\text{ej}}^2 \mu_e)$ is the column density of the ejecta, and $\mu_e \approx 2$ is the mean molecular weight. The cross section $\sigma$ and inelasticity $\kappa$ of the interactions depend on the γ-ray energy (Equations (40), (46), and (48) of Murase et al. 2015, also see Dermer & Menon 2009 and references therein). Figure 10 shows that comparing to high-energy γ rays, MeV to sub-GeV γ rays have a better chance to be observed at early times. Note that the optical depths due to photon-matter interactions are estimated assuming the spherical geometry. If X-rays from the engine are observed, a more complicated geometry seems necessary (Margutti et al. 2018), and then the escape of gamma-rays is likely to be easier in such more realistic setups. Also, we remark that the optical depth to the two-photon annihilation process decreases as energy, so UHE photons can escape from the system.

Murase et al. (2015) computed the γ-ray and hard X-ray emission from a magnetar-powered CCSN, taking into account details of electromagnetic cascades. Applying the numerical code to AT2018cow-like events to a low-mass ejecta with a magnetar, the GeV gamma-ray flux is estimated to be at the level of $10^{-12} \text{ erg cm}^{-2} \text{ s}^{-1}$ at the distance of AT2018cow. If the ejecta is asymmetric, as suggested by X-ray observations, MeV-GeV γ rays are likely to leak out from low-density regions in the similar geometry and be observed.

Renault-Tinacci et al. (2018) searched for GeV γ rays in the directions of an SLSNe sample with the Fermi-LAT data and found no signals. Assuming that the SLSNe population is equally luminous and has an $E^{-2}$ spectrum, Renault-Tinacci et al. (2018) concluded an upper limit at 95% C.L. to the 1–10 GeV luminosity $L_{\gamma,1-10 \text{ GeV}} < 9.1 \times 10^{41} \text{ erg s}^{-1}$. Our scenario is consistent with this limit, as the total spin-down power is below $3 \times 10^{41} \text{ erg s}^{-1}$ after 100 days.

For AT2018cow, a gamma-ray nondetection was reported by Fermi-LAT.<fnref n="13" /> during an one-week interval from day 3 to day 10 since the epoch of detection. The flux upper limit at 10 days is approximately $\sim 10^{-11} \text{ erg cm}^{-2} \text{ s}^{-1}$ (see Murase et al. 2018; Renault-Tinacci et al. 2018). HESS reported a 95% C.L. upper limit of $5 \times 10^{-13} \text{ cm}^{-2} \text{ s}^{-1}$ based on their observation in the third week.<fnref n="14" /> Both are consistent with our model.

## 5. Discussion and Conclusions

Rapidly spinning magnetars have been proposed as the engines responsible for powering the optical light curves of SLSNe. The recent discovery of AT2018cow, the first local example of a superluminous transient with a low ejecta mass (a so-called FBOT), enabled the discovery of coincident time-variable X-ray emission consistent with the presence of a central engine, given that X-rays escape due to the asymmetry of the ejecta (e.g., Ho et al. 2018; Margutti et al. 2018; Perley et al. 2018; Prentice et al. 2018). The engine behind AT2018cow released more than $\sim 10^{50} \text{ erg}$ over a timescale of less than a few days, behind a low-mass ejecta shell $M_{\text{ej}} \lesssim 0.3 M_{\odot}$. If the engine is a millisecond magnetar spinning down in isolation (i.e., neglecting fall-back accretion), then fits to the light curve require an initial spin period of $\sim 2–10 \text{ ms}$ and a dipole magnetic field of $\sim 1–2 \times 10^{15} \text{ G}$.

For the same reasons they are capable of powering large luminosities, the extreme magnetars responsible for SLSNe and FBOTs also possess strong electromagnetic potentials, making them potential sites for the acceleration of relativistic particles and even UHECRs. Events with particularly low ejecta masses, such as AT2018cow, provide a way to directly view X-rays from the magnetar nebula (through prompt photoionization of the ejecta shell; Metzger et al. 2014) and allow for the timely escape of accelerated cosmic rays relative to the bulk of the normal SLSNe population.

Motivated by the discovery of AT2018cow as a local example of an FBOT, which shows direct evidence for a central engine, we have calculated the interaction of high-energy cosmic rays accelerated near the magnetar with background particles of the nebula and ejecta. We take into account the time-evolving thermal and nonthermal radiation field of the nebula, and track all primary and higher-order interaction products down to TeV energies.

<footnote n="13">¹³ http://www.astronomerstelegram.org/?read=11808</footnote>

<footnote n="14">¹⁴ http://www.astronomerstelegram.org/?read=11956</footnote>
Our results are largely insensitive to the modeling of the background radiation fields. The photopion interaction is dominated by the optical/UV emission, which is directly observed, and the proton–proton interaction depends on the density of baryons, which is well constrained by the rise time of the explosion and the spectroscopically observed ejecta velocity. Our results do, however, depend on whether the engine behind AT2018cow is truly a magnetar (e.g., as opposed to an internal shock from CSM interaction or an accreting black hole), and whether millisecond magnetars are indeed efficient particle accelerators. Nevertheless, the direct detection of time-variable X-rays from AT2018cow provides greater confidence in the central engine scenario than was available from previous SLSNe samples (which generally show no coincident X-ray emission, likely due to photoelectric attenuation by the larger ejecta shells; Margalit et al. 2018; Margutti et al. 2018).

UHECR protons accelerated by millisecond magnetars in FBOTs and SLSNe escape the source with characteristic energies ~EeV. In addition, heavy elements may be synthesized efficiently in magnetar winds (Metzger et al. 2011) and destroyed inside nebulae (Horiuchi et al. 2012; Murase et al. 2014). Nuclei heavier than protons could tap Z times more energy from the same electric potential. Their energy losses are dominated by photodisintegration and hadronnuclear interaction. Fang et al. (2012) shows that they may escape from a massive ejecta of CCSNe with $E > 10^{20}$ eV if the pulsar has a millisecond spin period and surface dipole field in the range $10^{12}-10^{13}$ G. Neutrino production by heavy nuclei depends on the interaction channel. It could be comparable to that by protons when hadronnuclear interactions dominate (Fang 2015), but lower when photodisintegration dominates. We will leave a more detailed study to a future work. Engine-driven supernovae, including jet-driven ones, have been discussed as the sources of UHECRs, where particle acceleration sites have been attributed to internal shocks in outflows, as well as external forward and reverse shocks (e.g., Wang et al. 2007; Murase et al. 2008; Chakraborti et al. 2011; Zhang et al. 2018). Our physical model is different in the sense that UHECR acceleration occurs inside pulsar wind nebulae.

For AT2018cow, the cosmic-ray energies that can be achieved are lower than for the bulk of the SLSNe population considered previously. On the other hand, if relatively common events similar to AT2018cow provide the dominant UHECR source just above the ankle energy, a tail of rarer more powerful FBOTs or SLSNe (e.g., those born with initial spin periods close to their minimum break-up value $P_0 \sim 1$ ms instead of $P_0 \gtrsim 2$ ms) could dominate the UHECR budget at the highest energies.

Magnetars with birth properties similar to those required to power AT2018cow may be in some ways optimal for neutrino production. This is because of the apparent coincidence that the magnetar spin-down timescale is comparable to the narrow time window within which the optical depth for cosmic-ray interaction is still high, but the secondary pions are no longer efficiently cooled before decaying. FBOTs similar to AT2018cow can therefore be in principle ideal targets for neutrino telescopes.

Depending on the uncertain volumetric rate of FBOTs giving birth to millisecond magnetars with properties similar to those required to explain AT2018cow, we find that such a population could explain ~10%–20% of the IceCube astrophysical

neutrino background. However, given the sensitivity of current-generation facilities, the detection of individual FBOTs will be challenging without the fortuitous discovery of a source located several times closer than AT2018cow (which, however, was itself already much nearer than the previous cosmological populations of FBOTs and SLSNe). Nevertheless, IceCube neutrino events can be used to optimize the strategies of transient telescopes such as Zwicky Transient Facility for observing FBOTs.

The same interactions giving rise to neutrinos in FBOTs and SLSNe also inevitably give rise to high-energy gamma-rays. However, because of the high pair creation optical depth created by the thermal photons of the transient, high-energy gamma-rays require months to years to escape from the magnetar nebula and thus to become detectable from Earth. Nevertheless, well-studied optical light curves of some SLSNe show evidence for “leakage” of engine energy at late times $\gtrsim$months (e.g., Nicholl et al. 2018), which is not observed to be escaping in the soft X-ray band (Bhirombhakdi et al. 2018; Margutti et al. 2018) and therefore may emerge first in the hard X-ray or gamma-ray band.

Our work points to FBOTs and SLSNe as potentially promising sources for future neutrino and $\gamma$-ray telescopes. The predicted neutrino spectrum peaks at 100 PeV–EeV, which can be favorably probed by both EeV neutrino detectors such as Askaryan Radio Array (Ara Collaboration et al. 2012), ARIANNA (Barwick et al. 2017), GRAND (GRAND Collaboration et al. 2018), and POEMMA (Olinto et al. 2017), as well as TeV–PeV detectors such as KM3NeT (Margiotta 2014) and IceCube-Gen2 (IceCube-Gen2 Collaboration et al. 2014). Depending on the structure of the ejecta and the distance of the source, the $\gamma$ rays may be observed by existing wide-field telescopes such as Fermi and HAWC, and the next-generation Cerenkov Telescopes Array (The CTA Consortium 2017).

The authors thank Tanguy Pierog for helpful discussions on the pion–proton interaction calculations, Erik Blaufuss and the IceCube Collaboration for useful feedback. B.D.M. is supported by NASA (grants HST-AR-15041.001-A, 80NSSC18K1708, a80NSSC17K0501). K.M. is supported by the Alfred P. Sloan Foundation and NSF grant No. PHY-1620777. I.B. is grateful for the generous support of the University of Florida. K.K. is supported by the APACHE grant (ANR-16-CE31-0001) of the French Agence Nationale de la Recherche.

## ORCID iDs

Ke Fang https://orcid.org/0000-0002-5387-8138

Brian D. Metzger https://orcid.org/0000-0002-4670-7509

Kohta Murase https://orcid.org/0000-0002-5358-5642

Imre Bartos https://orcid.org/0000-0001-5607-3637

Kumiko Kotera https://orcid.org/0000-0002-6337-0696

## References

* Aartsen, M. G., Abraham, K., Ackermann, M., et al. 2016, ApJ, 833, 3

* Aartsen, M. G., Abraham, K., Ackermann, M., et al. 2017, ApJ, 835, 151

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

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

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

* Ara Collaboration, Allison, P., Auffenberg, J., et al. 2012, APh, 35, 457

* Arcavi, I., Wolf, W. M., Howell, D. A., et al. 2016, ApJ, 819, 35

* Arons, J. 2003, ApJ, 589, 871
Barwick, S. W., Besson, D. Z., Burgman, A., et al. 2017, *Aph*, 90, 50

Bhirombhakdi, K., Chornock, R., Margutti, R., et al. 2018, *ApJL*, 868, L32

Blasi, P., Epstein, R. I., & Olinto, A. V. 2000, *ApJL*, 533, L123

Blaufuss, E. 2018, *ATel*, 11785, 1

Cerutti, B., & Beloborodov, A. M. 2017, *SSRv*, 207, 111

Chakraborti, S., Ray, A., Soderberg, A. M., Loeb, A., & Chandra, P. 2011, *NatCo*, 2, 175

Dall'Osso, S., Shore, S. N., & Stella, L. 2009, *MNRAS*, 398, 1869

de Ugarte Postigo, A., Bremer, M., Kann, D. A., et al. 2018, *ATel*, 11749, 1

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

Drout, M. R., Chornock, R., Soderberg, A. M., et al. 2014, *ApJ*, 794, 23

Drout, M. R., Soderberg, A. M., Mazzali, P. A., et al. 2013, *ApJ*, 774, 58

Eidelman, S., Hayes, K. G., Olive, K. A., et al. 2004, *PhLB*, 592, 1

Fang, K. 2015, *JCAP*, 6, 004

Fang, K., Kotera, K., Murase, K., & Olinto, A. V. 2014, *PhRvD*, 90, 103005

Fang, K., Kotera, K., Murase, K., & Olinto, A. V. 2016, *JCAP*, 4, 010

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

Fang, K., & Metzger, B. D. 2017, *ApJ*, 849, 153

Fernández, R., Quataert, E., Kashiyama, K., & Coughlin, E. R. 2018, *MNRAS*, 476, 2366

GRAND Collaboration, Alvarez-Muniz, J., Alves Batista, R., et al. 2018, *arXiv:1810.09994*

Ho, A. Y. Q., Phinney, E. S., Ravi, V., et al. 2019, *ApJ*, 871, 73

Hopkins, A. M., & Beacom, J. F. 2006, *ApJ*, 651, 142

Horiuchi, S., Murase, K., Ioka, K., & Mészáros, P. 2012, *ApJ*, 753, 69

Hotokezaka, K., Kashiyama, K., & Murase, K. 2017, *ApJ*, 850, 18

IceCube Collaboration, Aartsen, M. G., Ackermann, M., et al. 2018, *PhRvD*, 98, 062003

IceCube-Gen2 Collaboration, Aartsen, M. G., et al. 2014, *arXiv:1412.5106*

Iwamoto, M., Amano, T., Hoshino, M., & Matsumoto, Y. 2017, *ApJ*, 840, 52

Kasen, D., & Bildsten, L. 2010, *ApJ*, 717, 245

Kashiyama, K., Murase, K., Bartos, I., Kiuchi, K., & Margutti, R. 2016, *ApJ*, 818, 94

Kennel, C. F., & Coroniti, F. V. 1984, *ApJ*, 283, 710

Kleiser, I. K. W., Kasen, D., & Duffell, P. C. 2018, *MNRAS*, 475, 3152

Kopper, C. 2018, ICRC (Busan), 35, 981

Kotera, K. 2011, *PhRvD*, 84, 023002

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

Kotera, K., Phinney, E. S., & Olinto, A. V. 2013, *MNRAS*, 432, 3228

Kuin, N. P. M., Wu, K., Oates, S., et al. 2019, *MNRAS*, tmp stz053

Liu, L.-D., Zhang, B., Wang, L.-J., & Dai, Z.-G. 2018, *ApJL*, 868, L24

Lyutikov, M., & Toonen, S. 2018, *arXiv:1812.07569*

Margalit, B., Metzger, B. D., Berger, E., et al. 2018, *MNRAS*, 481, 2407

Margiotta, A. 2014, *NIMPA*, 766, 83

Margutti, R., Metzger, B. D., Chornock, R., et al. 2018, *MNRAS*, 481, 2407

Matthews, J. 2018, ICRC (Busan), 35, 1096

Metzger, B. D., Beniamini, P., & Giannios, D. 2018, *ApJ*, 857, 95

Metzger, B. D., Giannios, D., & Horiuchi, S. 2011, *MNRAS*, 415, 2495

Metzger, B. D., & Piro, A. L. 2014, *MNRAS*, 439, 3916

Metzger, B. D., Quataert, E., & Thompson, T. A. 2008, *MNRAS*, 385, 1455

Metzger, B. D., Vurm, I., Hascoët, R., & Beloborodov, A. M. 2014, *MNRAS*, 437, 703

Moriya, T. J., & Eldridge, J. J. 2016, *MNRAS*, 461, 2155

Mücke, A., Engel, R., Rachen, J. P., Protheroe, R. J., & Stanev, T. 2000, *CoPhC*, 124, 290

Murase, K., Dasgupta, B., & Thompson, T. A. 2014, *PhRvD*, 89, 043012

Murase, K., & Fukugita, M. 2019, *PhRvD*, 99, 063012

Murase, K., Ioka, K., Nagataki, S., & Nakamura, T. 2008, *PhRvD*, 78, 023005

Murase, K., Kashiyama, K., Kiuchi, K., & Bartos, I. 2015, *ApJ*, 805, 82

Murase, K., Mészáros, P., & Zhang, B. 2009, *PhRvD*, 79, 103001

Murase, K., Toomey, M. W., Fang, K., et al. 2018, *ApJ*, 854, 60

Nicholl, M., Blanchard, P. K., Berger, E., et al. 2018, *ApJL*, 866, L24

Nicholl, M., Guillochon, J., & Berger, E. 2017, *ApJ*, 850, 55

Olinto, A. V., Adams, J. H., Aloisio, R., et al. 2017, *ICRC*, 35, 542

Ostriker, J. P., & Gunn, J. E. 1969, *ApJ*, 157, 1395

Özel, F., Psaltis, D., Ransom, S., Demorest, P., & Alford, M. 2010, *ApJL*, 724, L199

Perley, D. A., Mazzali, P.A., Yan, L., et al. 2019, *MNRAS*, 484, 1031

Philippov, A. A., & Spitkovsky, A. 2018, *ApJ*, 855, 94

Pierog, T., Karpenko, I., Katzy, J. M., Yatsenko, E., & Werner, K. 2015, *PhRvC*, 92, 034906

Piro, A. L., & Kollmeier, J. A. 2016, *ApJ*, 826, 97

Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, *A&A*, 594, A13

Prentice, S. J., Maguire, K., Smartt, S. J., et al. 2018, *ApJL*, 865, L3

Renault-Tinacci, N., Kotera, K., Neronov, A., & Ando, S. 2018, *A&A*, 611, A45

Rest, A., Garnavich, P. M., Khatami, D., et al. 2018, *NatAs*, 2, 307

Rivera Sandoval, L. E., Maccarone, T. J., Corsi, A., et al. 2018, *MNRAS*, 480, 146

Shore, S. N., Dall'Osso, S., & Stella, L. 2009, *MNRAS*, 398, 1869

Smartt, S. J., Clark, P., Smith, K. W., et al. 2018, *ATel*, 11727, 1

Spitkovsky, A. 2006, *ApJL*, 648, L51

Stella, L., Dall'Osso, S., Israel, G. L., & Vecchio, A. 2005, *ApJL*, 634, L165

Stella, L., Dall'Osso, S., Israel, G. L., & Vecchio, A. 2005, *ApJL*, 634, L165

Svensson, R. 1987, *MNRAS*, 227, 403

Tanaka, S. J., & Takahara, F. 2010, *ApJ*, 715, 1248

Tanaka, S. J., & Takahara, F. 2013, *MNRAS*, 429, 2945

Tauris, T. M., Langer, N., & Podsiadlowski, P. 2015, *MNRAS*, 451, 2123

The CTA Consortium 2017, *arXiv:1709.05434*

The Pierre Auger Collaboration, Aab, A., Abreu, P., et al. 2017, *arXiv:1708.06592*

Torres, D., Cillis, A., Martín, J., & de Oña Wilhelmi, E. 2014, *JHEAp*, 1–2, 31

Wang, X.-Y., Razzaque, S., Mészáros, P., & Dai, Z.-G. 2007, *PhRvD*, 76, 083009

Woosley, S. E. 2010, *ApJL*, 719, L204

Yu, Y.-W., Zhang, B., & Gao, H. 2013, *ApJL*, 776, L40

Zhang, B. T., Murase, K., Kimura, S. S., Horiuchi, S., & Mészáros, P. 2018, *PhRvD*, 97, 083010