ULTRAHIGH-ENERGY COSMIC RAYS FROM THE "EN CAUL" BIRTH OF MAGNETARS^*

The initial rotational energy of the magnetar with moment of inertia I is

$$\begin{eqnarray}&&{E}_{{\rm{rot}}}=\displaystyle \frac{1}{2}I{{\rm{\Omega }}}_{i}^{2}=4.9\times {10}^{51}{I}_{45}{\left(\displaystyle \frac{{P}_{i}}{2{\rm{ms}}}\right)}^{-2}{\rm{erg}},\end{eqnarray} \tag{ 1 }$$

where ${I}_{45}=I/{10}^{45}\,{\rm{g}}\,{\mathrm{cm}}^{2}$. The spin down of the magnetar with magnetic moment $\mu ={R}_{* }^{3}{B}_{* }$ injects energy with a time-dependent rate of

$$\begin{eqnarray}&&{L}_{{\rm{sd}}}(t)=\displaystyle \frac{{\mu }^{2}{\rm{\Omega }}{(t)}^{4}}{{c}^{3}}=\displaystyle \frac{{\mu }^{2}{{\rm{\Omega }}}_{i}^{4}}{{c}^{3}}{\left(1+\displaystyle \frac{t}{{t}_{{\rm{sd}}}}\right)}^{-2},\end{eqnarray} \tag{ 2 }$$

where the spin down timescale is

$$\begin{eqnarray}&&{t}_{{\rm{sd}}}=\displaystyle \frac{{{Ic}}^{3}}{2{\mu }^{2}{{\rm{\Omega }}}_{i}^{2}}=22.8{\mu }_{33}^{-2}{I}_{45}{\left(\displaystyle \frac{{P}_{i}}{2{\rm{ms}}}\right)}^{2}{\rm{minute}}.\end{eqnarray} \tag{ 3 }$$

The ejecta around the magnetar is inflated by this energy input with a time-dependent velocity of

$$\begin{eqnarray}&&{v}_{{\rm{ej}}}(t)=\displaystyle \frac{{{dR}}_{{\rm{ej}}}}{{dt}}={\left[2{M}_{{\rm{ej}}}^{-1}{\int }_{0}^{t}{L}_{{\rm{sd}}}(t^{\prime} ){dt}^{\prime} +{v}_{\mathrm{ej},i}^{2}\right]}^{1/2},\end{eqnarray} \tag{ 4 }$$

which can be simply integrated to find ${R}_{{\rm{ej}}}(t)$. At late times, the velocity of the ejecta can become high,

$$\begin{eqnarray}{v}_{{\rm{ej}}}(t\gg {t}_{{\rm{sd}}}) & \approx & {(2{E}_{{\rm{rot}}}/{M}_{{\rm{ej}}})}^{1/2}\\ & \approx & 0.7{{cI}}_{45}^{1/2}{M}_{-2}^{-1/2}{\left(\displaystyle \frac{{P}_{i}}{2{\rm{ms}}}\right)}^{-1},\end{eqnarray} \tag{ 5 }$$

where ${M}_{-2}={M}_{{\rm{ej}}}/{10}^{-2}\,{M}_{\odot }$. To simplify our analysis, we focus on values of ${E}_{{\rm{rot}}}$ and ${M}_{{\rm{ej}}}$ where ${v}_{{\rm{ej}}}\lesssim c$.

In order to elucidate the UHECR propagation, we divide the perimagnetar environment into two primary regions: (1) the ejecta shell and (2) the nebula (or bubble) of relativistic particles around the magnetar (see the diagram in Metzger & Piro 2014). In each of these regions, both the particle density and radiation content must be followed since these can interact with UHECRs and alter their fate. In this simplified picture, we do not explicitly follow the dynamics of the termination shock between the magnetar wind and the nebula. Nevertheless, this shock is important for converting the initial Poynting-flux dominated outflow of the magnetar into kinetic energy, as has been extensively studied in past work (e.g., Lyubarsky & Kirk 2001; Arons 2008; Lyubarsky 2010; Granot et al. 2011; Porth et al. 2013). Since the presence of the termination shock does not impact the shell dynamics beyond providing a potential site for particle acceleration, we mostly ignore it when solving for the properties of the nebula.

For the ejecta region, the density simply evolves as the ejecta expands as

$$\begin{eqnarray}&&{\rho }_{{\rm{ej}}}(t)=\displaystyle \frac{3{M}_{{\rm{ej}}}}{4\pi {R}_{{\rm{ej}}}{(t)}^{3}}.\end{eqnarray} \tag{ 6 }$$

The composition is thought to be primarily iron-rich elements due to the higher electron fraction produced by neutrino irradiation from the magnetar (Metzger et al. 2009b; Fernández & Metzger 2013). We estimate the number density of nuclei in the ejecta as ${n}_{{\rm{ej}}}\approx {\rho }_{{\rm{ej}}}/56{m}_{p}$, where m_p is the proton mass.

The second region is the relativistic bubble interior to the ejecta. This is inflated by the magnetar spin down luminosity, which can be partitioned into a combination of ${e}^{\pm }$ pairs, magnetic field, radiation, and, importantly for the study here, UHECRs. The high density of radiation and magnetic field within the nebula gives rise to copious ${e}^{\pm }$ pairs, with the strength of the pair creation described by the "compactness parameter":

$$\begin{eqnarray}&&{\ell }=\displaystyle \frac{3{E}_{{\rm{nth}}}{\sigma }_{{\rm{Th}}}}{4\pi {R}_{{\rm{ej}}}^{2}{m}_{e}{c}^{2}},\end{eqnarray} \tag{ 7 }$$

where ${E}_{{\rm{nth}}}$ is the energy in non-thermal radiation within the nebula (which we describe in more detail below) and ${\sigma }_{{\rm{Th}}}$ is the Thomson cross-section. When ${\ell }\gtrsim 1$, the number density of these pairs can be estimated by balancing the rate of pair creation and annihilation (Metzger & Piro 2014; Metzger et al. 2014), resulting in

$$\begin{eqnarray}&&{n}_{\pm }(t)={\left[\displaystyle \frac{4{{YL}}_{{\rm{sd}}}(t)}{\pi {R}_{{\rm{ej}}}{(t)}^{3}{\sigma }_{{\rm{Th}}}{m}_{e}{c}^{3}}\right]}^{1/2},\end{eqnarray} \tag{ 8 }$$

where ${\sigma }_{{\rm{Th}}}=6.65\times {10}^{-25}\,{\rm{cm}}$ is the Thomson cross-section and $Y\approx 0.1$ is the pair multiplicity factor for a saturated cascade (essentially the fraction of the magnetar spin down that goes into pairs; Svensson 1987). In our case, ${\ell }\gtrsim 1$ for the times we consider, as we show explicitly below. This density is in contrast to other work on UHECRs from magnetars in a core-collapse scenario (Fang et al. 2012; Lemoine et al. 2015) that instead use a Goldreich–Julian density (Goldreich & Julian 1969). Another difference is that Equation (8) assumes that the ${e}^{\pm }$ are fairly sub-relativistic. This is reasonable because even if the pairs are created with a high Lorentz factor, the ratio of the Compton cooling timescale ${t}_{{\rm{C}}}$ to the time the nebula has been expanding t is

$$\begin{eqnarray}\displaystyle \frac{{t}_{{\rm{C}}}}{t} & = & \displaystyle \frac{\pi {m}_{e}{{cR}}_{{\rm{ej}}}^{3}}{{\sigma }_{{\rm{Th}}}{E}_{{\rm{nth}}}t}\sim \displaystyle \frac{\pi {m}_{e}{{cv}}_{{\rm{ej}}}^{3}t}{{\sigma }_{{\rm{Th}}}{L}_{{\rm{sd}}}}\\ & \sim & 0.1{\mu }_{33}^{2}{I}_{45}^{-1/2}{M}_{-2}^{-3/2}{\left(\displaystyle \frac{t}{1{\rm{day}}}\right)}^{-2},\end{eqnarray} \tag{ 9 }$$

and thus the pairs cool rapidly in the dense nebula.

In addition to ${e}^{\pm }$ pairs, the nebula will be filled with magnetic field. As mentioned above, the general picture is that the magnetar outflow is Poynting-flux dominated inside of the termination shock, but dissipation near the shock results in a much smaller mean field in the rest of the nebula. Since the exact conversion factor is not known, especially for the extreme case of magnetars we consider here, we quantify this uncertainty with a factor ${\eta }_{B}$ that is the fraction of the magnetar spin down that goes into this magnetic field. This then gives

$$\begin{eqnarray}&&{B}_{{\rm{neb}}}(t)\approx {\left(\displaystyle \frac{6{\eta }_{B}}{{R}_{{\rm{ej}}}{(t)}^{3}}{\int }_{0}^{t}{L}_{{\rm{sd}}}(t^{\prime} ){dt}^{\prime} \right)}^{1/2},\end{eqnarray} \tag{ 10 }$$

as an estimate of the time-dependent mean magnetic field in the nebula. We note that this magnetic field was left out of the energy budget of the nebula in the calculations by Metzger & Piro (2014), but this only represents a small correction since ${\eta }_{B}\lesssim 0.1$.

Next, we consider the radiation content in both the nebula and surrounding ejecta. These are highly coupled. The non-thermal radiation in the nebula is dominated by synchrotron and inverse-Compton emission (Metzger et al. 2014; Murase et al. 2015). This is then partially absorbed and thermalized by bound-free reactions of the iron-peak elements in the ejecta. In the Appendix of Metzger & Piro (2014), we showed that this interplay can be accurately followed in each region by simply considering two coupled differential equations. The non-thermal energy of the nebula, ${E}_{{\rm{nth}}}$, is simply given by the energy injection from spin down minus the adiabatic and absorptive losses as a function of time

$$\begin{eqnarray}&&\displaystyle \frac{{{dE}}_{{\rm{nth}}}}{{dt}}={L}_{{\rm{sd}}}-{L}_{{\rm{abs}}}-\displaystyle \frac{{E}_{{\rm{nth}}}}{t+{t}_{i}}.\end{eqnarray} \tag{ 11 }$$

What is absorbed by the ejecta from the non-thermal reservoir is converted to heat, and therefore the evolution of the thermal radiation of the ejecta, ${E}_{{\rm{th}}}$, evolves as

$$\begin{eqnarray}&&\displaystyle \frac{{{dE}}_{{\rm{th}}}}{{dt}}={L}_{{\rm{abs}}}-{L}_{{\rm{rad}}}-\displaystyle \frac{{E}_{{\rm{th}}}}{t+{t}_{i}}.\end{eqnarray} \tag{ 12 }$$

These equations have a number of terms that we now describe in detail. The rate of absorption of non-thermal radiation is

$$\begin{eqnarray}&&{L}_{{\rm{abs}}}=\displaystyle \frac{{f}_{{\rm{abs}}}{E}_{{\rm{nth}}}}{{t}_{{\rm{diff,neb}}}},\end{eqnarray} \tag{ 13 }$$

where ${f}_{{\rm{abs}}}$ is the fraction of non-thermal radiation absorbed by the ejecta and the diffusion time across the nebula is

$$\begin{eqnarray}&&{t}_{{\rm{diff,neb}}}=({\sigma }_{{\rm{Th}}}{R}_{{\rm{ej}}}{n}_{\pm }+1){R}_{{\rm{ej}}}/c.\end{eqnarray} \tag{ 14 }$$

The term ${L}_{{\rm{abs}}}$ appears in both Equations (11) and (12) as it mediates the transport of energy between these two reservoirs of radiation. The rate of radiative energy loss from the ejecta is simply

$$\begin{eqnarray}&&{L}_{{\rm{rad}}}=\displaystyle \frac{{E}_{{\rm{th}}}}{{t}_{{\rm{diff,ej}}}},\end{eqnarray} \tag{ 15 }$$

where

$$\begin{eqnarray}&&{t}_{{\rm{diff,ej}}}=\left(\displaystyle \frac{{\kappa }_{{\rm{es}}}{M}_{{\rm{ej}}}}{4\pi {R}_{{\rm{ej}}}^{2}}+1\right)\displaystyle \frac{{R}_{{\rm{ej}}}}{c},\end{eqnarray} \tag{ 16 }$$

is the photon diffusion time through the ejecta and ${\kappa }_{{\rm{es}}}$ is the electron scattering opacity of the ejecta. The term ${L}_{{\rm{rad}}}$ is important for estimating the thermal radiation observed from the heated ejecta. The timescale ${t}_{i}={R}_{\mathrm{ej},i}/{v}_{\mathrm{ej},i}$ is the initial expansion timescale of the ejecta. It is included to produce well-behaved solutions at early times and the results are largely insensitive to its value. Finally, the last term in each expression accounts for the adiabatic expansion of the relativistic gas, and hence has the familiar scaling of ${E}_{{\rm{(n)th}}}/t$ at late times.

From the above description, given specific properties for the magnetar and ejecta, one can solve for the evolution in ${n}_{{\rm{ej}}}$, ${n}_{\pm }$, ${E}_{{\rm{nth}}}$, and ${E}_{{\rm{th}}}$ to investigate the possible UHECRs evolution. An example solution is presented in Figures 1 and 2 for ${P}_{i}=3\,{\rm{ms}}$, ${B}_{* }=5\times {10}^{14}\,{\rm{G}}$, and ${M}_{{\rm{ej}}}=0.01\,{M}_{\odot }$. In this case, we use f = 1, which generally reproduces the features of more detailed treatments (Metzger & Piro 2014). Starting with Figure 1, the thermal light curve ${L}_{{\rm{rad}}}$ shows the characteristic supernova-like shape, which peaks on a $\sim 1\,{\rm{day}}$ timescale, which is potentially the best electromagnetic signature for identifying such events (note that the non-thermal X-rays have a similar luminosity and timescale as well, Metzger & Piro 2014). The observed ${L}_{{\rm{rad}}}$ is well below ${L}_{{\rm{sd}}}$ because of the large optical depth in pairs hindering the transfer of energy from the nebula to the ejecta, as can be seen from ${\tau }_{\pm }\gtrsim {\tau }_{{\rm{ej}}}$ near the light curve peak. Here, we use ${\tau }_{\pm }={\sigma }_{{\rm{Th}}}{R}_{{\rm{ej}}}{n}_{\pm }$ and ${\tau }_{{\rm{ej}}}={\kappa }_{{\rm{es}}}{M}_{{\rm{ej}}}/4\pi {R}_{{\rm{ej}}}^{2}$. In the bottom panel, we plot ${n}_{{\rm{ej}}}$ and ${n}_{\pm }$, which will be important for understanding UHECR propagation in subsequent sections. We also plot the compactness ℓ in Figure 1, demonstrating that it is greater than unity for most of the evolution, justifying our treatment of the ${e}^{\pm }$ pairs in the saturated regime.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 2, we summarize the more important quantities from the evolution of this model. The top panel shows both the temperature of the ejecta as well as the effective temperature of the thermal optical signature. These become roughly equal once ${\tau }_{{\rm{ej}}}\lesssim 1$. The mean magnetic field in the nebula ${B}_{{\rm{neb}}}$, calculated from Equation (10), will be important for understanding the acceleration of UHECRs as well as for potentially cooling them via synchrotron radiation. In the bottom panel, we summarize the critical densities in each regime, along with the radius of the expanding ejecta ${R}_{{\rm{ej}}}$.

We first summarize the basic energetics to show that the magnetar surface can indeed easily accelerate CRs to ultrahigh energies. Relativistic magnetic rotators have voltage drops across their magnetic field of

$$\begin{eqnarray}&&{\rm{\Phi }}(t)=\displaystyle \frac{\mu {\rm{\Omega }}{(t)}^{2}}{{c}^{2}}.\end{eqnarray} \tag{ 17 }$$

A particle of charge Z that is stripped off the surface of the magnetar (by some combination of strong electric fields, bombardment of particles, and boiling by stellar heat) can be accelerated to an energy of

$$\begin{eqnarray}&&{E}_{{\rm{surf}}}={Ze}{\rm{\Phi }}=3.3\times {10}^{21}Z{\mu }_{33}{\left(\displaystyle \frac{P}{2{\rm{ms}}}\right)}^{-2}\,{\rm{eV}}.\end{eqnarray} \tag{ 18 }$$

This shows that, in practice, it is reasonable to expect nuclei accelerated to sufficient energies to produce UHECRs, especially for large charges ($Z\gt 1$). The problem is that when this acceleration occurs within the light cylinder of the magnetar ${R}_{{\rm{LC}}}={\rm{\Omega }}/c$ (which is much less than ${R}_{{\rm{ej}}}$), then the curvature radiation losses prevent acceleration to such high energies. For this reason, both Arons (2003) and Fang et al. (2012) consider acceleration in the wind at radii $\gg {R}_{{\rm{LC}}}$ to avoid this problem.

Due to this difficulty, we instead consider the energy to which ions can be accelerated just from confinement by magnetic fields within the nebula (Lemoine et al. 2015), which gives

$$\begin{eqnarray}&&{E}_{{\rm{conf}}}(t)\approx {{ZeB}}_{{\rm{neb}}}(t){R}_{{\rm{ej}}}(t).\end{eqnarray} \tag{ 19 }$$

The magnetic field that accelerates particles can also limit the energy of the UHECRs due to synchrotron losses. Equating the synchrotron cooling time to the gyration time of the particles results in

$$\begin{eqnarray}&&{E}_{{\rm{synch}}}(t)=\displaystyle \frac{3}{2}\displaystyle \frac{{({{Am}}_{p}{c}^{2})}^{2}}{{({Ze})}^{3/2}{B}_{{\rm{neb}}}{(t)}^{1/2}}.\end{eqnarray} \tag{ 20 }$$

With this in mind, we consider the highest possible energy for UHECRs at any given time to be the minimum between ${E}_{{\rm{conf}}}$ and ${E}_{{\rm{synch}}}$.

This maximum energy is shown in the upper panel of Figure 3 for UHECRs with either proton or Fe composition. At early times, the nebula is dense enough in magnetic fields that synchrotron cooling dominates and the energy is $\approx {E}_{{\rm{sync}}}$. Then, as the nebula expands and the magnetar spins down, synchrotron losses decrease until the particle energy is instead $\approx {E}_{{\rm{conf}}}$. The maximum energy for UHECRs from a given magnetar is where these two energies are roughly equal. Also, note that heavier nuclei generally reach higher energies because their larger charge promotes more acceleration.

Zoom In Zoom Out Reset image size

Even though particles may be accelerated to such high energies, interaction with the environment around the magnetar may prevent the UHECRs from leaving unhindered. We now consider each of these, both in the ejecta and the nebula, in more detail.

For the ejecta, we use a hadron interaction model to estimate the timescale for sapping the UHECRs of their energy. This timescale is given by

$$\begin{eqnarray}&&{t}_{{\rm{ej}}}(E)={[{{cn}}_{{\rm{ej}}}{\sigma }_{{\rm{ej}}}(E)\xi (E)]}^{-1},\end{eqnarray} \tag{ 21 }$$

where ${\sigma }_{{\rm{ej}}}$ is the cross-section for hadronic interactions and ξ is the elasticity. Also, note that all of these terms will have energy dependencies in detail. The cross-section for iron–proton interactions is $1.25\,{\rm{barn}}$, which applies when Z = 1. Since the cross-section should roughly scale as ${Z}^{2/3}$, for iron–iron interactions we estimate the cross-section as $\approx 11\,{\rm{barn}}$. Although ξ can change significantly depending on the number of nucleons, we take it to be of the order of unity given the uncertainties in the cross-sections, which are degenerate with this term. The relevant timescale to compare these to is the crossing timescale for the ejecta:

$$\begin{eqnarray}&&{t}_{{\rm{cross}}}(t)={R}_{{\rm{ej}}}(t)/c.\end{eqnarray} \tag{ 22 }$$

When ${t}_{{\rm{cross}}}\lt {t}_{{\rm{ej}}}$, the UHECRs can escape the ejecta without significant interaction (and vice versa).

In Figure 3, we present an example calculation of the relevant timescales for assessing the interactions of UHECRs with the ejecta along with the maximum possible energy for the generated UHECRs as a function of time. We also compare different compositions of UHECRs: both protons (blue dashed lines) and iron (red dotted lines). This demonstrates that protons satisfy ${t}_{{\rm{cross}}}\lt {t}_{{\rm{ej}}}$ earlier during the evolution and begin escaping earlier (everything to the right of the blue circle escapes). Nevertheless, compared with the energy the protons would have at that time, protons do not reach energies as high as those of iron. Thus, all things being equal, a large charge is favored for producing the highest-energy UHECRs.

With this picture in mind, we next estimate the maximum energy possible for iron-composition UHECRs as a function of the initial spin period P_i and the magnetic field B_*. This is done by considering the highest possible energy for the UHECR past the time when they start escaping the ejecta (past the solid circles from Figure 3). The results are summarized in the contour plot shown in Figure 4, where we use ${M}_{{\rm{ej}}}=0.01\,{M}_{\odot }$. For these calculations, we assume f = 1, which gives a light curve solution that qualitatively matches that of Metzger & Piro (2014) and allows us to cover a wide parameter space of models. In the future, a more detailed treatment of the thermalization should be conducted. This demonstrates that to produce the highest possible energy, UHECRs require a spin period of P_i ≈ 2–4 ms and a rather modest (for a magnetar) magnetic field of B ≈ (1–4) × 10¹⁴ G (note that although we consider ${P}_{i}\lesssim 2\,{\rm{ms}}$ to be unlikely for the reasons mentioned above, we include this parameter space in our plots for completeness). The highest field magnetars are disfavored for producing high-energy particles because they spin down too fast. Even though they produce high-energy particles, these all interact with the ejecta at early times before it can become sufficiently diffuse.

Zoom In Zoom Out Reset image size

When we consider more massive ejecta, the maximum energies for the UHECRs decrease slightly ($\sim 10 \%$ for $0.05\,{M}_{\odot }$) and the "sweet spot" for most efficiently accelerating UHECRs shifts to slightly lower magnetic fields. This is because a lower magnetic field spins down the magnetar more slowly so that the timescale for the energy injection better matches the longer diffusion time through the more massive ejecta. With very high ejecta masses, as is more associated with core-collapse (e.g., Fang et al. 2012), the shift to lower energies can, in some cases, make it difficult to produce UHECRs at all. Overall, however, the ejecta mass impacts the energetics less than B_* and P_i for our model, and thus we focus on these two parameters for most of the remainder of our study.

A fraction of the X-rays coming out of the nebula thermalize with the ejecta. This generates a bath of thermal photons that can lead to photodisintegrations of the iron nuclei when they attempt to pass through. To estimate when this should occur, we largely follow the methodology presented in Kotera et al. (2015).

For ejecta with temperature ${T}_{{\rm{ej}}}$, the photon energy peaks at roughly ${\epsilon }_{{\rm{th}}}\approx 2.8{k}_{B}{T}_{{\rm{ej}}}\approx 2.4{T}_{4}\,{\rm{eV}}$, where ${T}_{4}={T}_{{\rm{ej}}}/{10}^{4}\,{\rm{K}}$ (see Figure 2). In contrast, the threshold energy of the photodisintegration process for iron-like nuclei is ${\epsilon }_{A}\approx 18.3{A}_{56}^{-0.21}\,{\rm{MeV}}$ (Puget et al. 1976; Murase et al. 2009), where ${A}_{56}=A/56$ denotes the atomic mass of the nuclei. The resulting threshold Lorentz factor for photodisintegration is

$$\begin{eqnarray}&&{\gamma }_{{\rm{thres}}}\approx \displaystyle \frac{{\epsilon }_{A}}{{\epsilon }_{{\rm{th}}}}\approx 7.7\times {10}^{6}{A}_{56}^{-0.21}{T}_{4}^{-1},\end{eqnarray} \tag{ 23 }$$

or an energy of $E\approx 4\times {10}^{17}\,{\rm{eV}}$ for iron-like nuclei. Thus, this threshold is always easily met in our scenario.

In addition to exceeding this Lorentz factor, the photodisintegration must be sufficiently fast compared to the crossing timescale (Murase et al. 2014). Using a cross-section of ${\sigma }_{A}\approx 8\times {10}^{-26}{A}_{56}\,{\mathrm{cm}}^{2}$, as is appropriate for a Giant Dipole Resonance (Puget et al. 1976; Murase et al. 2008), this timescale is

$$\begin{eqnarray}&&{t}_{{\rm{th}}}\sim {\left(\displaystyle \frac{{\sigma }_{A}{{acT}}_{{\rm{th}}}^{4}}{{\epsilon }_{{\rm{th}}}}\right)}^{-1}\sim 21{A}_{56}^{-1}{T}_{5}^{-3}\,{\rm{s}},\end{eqnarray} \tag{ 24 }$$

and thus strongly increases as the ejecta expands and cools. Nevertheless, during the first $\sim 1\mbox{--}5\,{\rm{days}}$ of evolution, we generally find that ${t}_{{\rm{th}}}\ll {t}_{{\rm{cross}}}$, and thus the UHECRs will be subject to strong photo-hadronic interactions. This does not prevent iron UHECRs, but these estimates imply that lower-mass nuclei will also necessarily be present. We save detailed calculations of the composition of the UHECRs as a function of energy for future work. We do note that once ${t}_{{\rm{th}}}\gtrsim {t}_{{\rm{cross}}}$, it is still possible that $E\gtrsim {10}^{20}\,{\rm{eV}}$, and thus UHECRs from iron-like nuclei will survive.

Another potential location of strong photo-hadronic interactions is within the nebula itself, which is full of non-thermal X-rays. Here, the high energy of the individual photons makes it easier to exceed the Lorentz threshold. Using the total energy in non-thermal X-rays, ${E}_{{\rm{nth}}}$, we thus estimate the timescale for interaction with these as

$$\begin{eqnarray}&&{t}_{{\rm{nth}}}\sim {\left(\displaystyle \frac{{\sigma }_{A}c}{{\epsilon }_{{\rm{nth}}}}\displaystyle \frac{3{E}_{{\rm{nth}}}}{4\pi {R}_{{\rm{ej}}}^{3}}\right)}^{-1},\end{eqnarray} \tag{ 25 }$$

where ${\epsilon }_{{\rm{nth}}}$ is the energy of the photon. Thus, at higher photon energies, this timescale is longer because the power-law distribution of the photons is decreasing. This spectrum continues up to roughly the threshold for pair creation $\sim 2{m}_{e}{c}^{2}\sim {\rm{MeV}}$.

We summarize the main arguments of this section in Figure 5 where we plot the relevant timescales for an example model. This demonstrates that where ${t}_{{\rm{th}}}\lt {t}_{{\rm{cross}}}$, photodisintegration in the ejecta is relevant and this process generally dominates over photodisintegration within the nebula. This likely leads to a mixed composition of UHECRs up to a few days. After this timescale, the UHECRs can leave the magnetar less impacted by the environment.

Zoom In Zoom Out Reset image size

The energy generation rate of UHECRs above ${10}^{19.3}\,{\rm{eV}}$ is estimated to be in the range $\approx {10}^{43.5}\mbox{--}{10}^{44}\,\mathrm{erg}\,{\mathrm{Mpc}}^{-3}\,{\mathrm{yr}}^{-1}$ (Waxman 1995; Berezinsky et al. 2006; Katz et al. 2009; Murase & Takami 2009). Specifically, such estimates depend on the assumed energy spectrum and volume over which the UHECRs can be observed, which may even be impacted by the (uncertain) composition. Nevertheless, such rough rates are helpful for constraining potential mechanisms. In the following discussion, we use ${10}^{43.75}\,\mathrm{erg}\,{\mathrm{Mpc}}^{-3}\,{\mathrm{yr}}^{-1}$ as the characteristic rate to compare to for generating UHECRs.

For the magnetar model, the maximum energy available per event is ${E}_{{\rm{rot}}}$, but we do not expect this entire energy reservoir to be injected into UHECRs. Assuming a fraction ${f}_{{\rm{cr}}}$ does make it into UHECRs, the required rate, estimated from Equation (1) and the volumetric estimate above, is simply

$$\begin{eqnarray}&&R\approx 1.1\times {10}^{-8}{I}_{45}^{-1}{\left(\displaystyle \frac{{P}_{i}}{2{\rm{ms}}}\right)}^{2}{f}_{{\rm{cr}}}^{-1}\,{\mathrm{Mpc}}^{-3}\,{\mathrm{yr}}^{-1}.\end{eqnarray} \tag{ 26 }$$

For AICs, it is helpful to compare to the SN Ia rate since both types of events are produced from WDs. The Lick Observatory Supernova Search found a rate of ${R}_{{\rm{Ia}}}=(3.01\pm 0.062)\,\times {10}^{-5}\,{\mathrm{Mpc}}^{-3}\,{\mathrm{yr}}^{-1}$ (Li et al. 2011), and thus $R/{R}_{{\rm{Ia}}}\,\approx 4\times {10}^{-4}{f}_{{\rm{cr}}}^{-1}$. This can be compared to the AIC rates that have been proposed in the literature, although these are known to be very uncertain. Using population synthesis, Yungelson & Livio (1998) find AIC rates of $8\times {10}^{-7}\mbox{--}8\times {10}^{-5}\,{\mathrm{yr}}^{-1}$ for the Milky Way, depending on assumptions about the common-envelope phase and mass transfer, which correspond to volumetric rates of ${R}_{{\rm{AIC}}}\sim 3\times {10}^{-9}\mbox{--}6\times {10}^{-7}\,{\mathrm{Mpc}}^{-3}\,{\mathrm{yr}}^{-1}$. Similar constraints are also obtained from observed abundances of neutron-rich isotopes (Hartmann et al. 1985; Fryer et al. 1999). This would give a relative rate of ${R}_{{\rm{AIC}}}/{R}_{{\rm{Ia}}}\sim 1\times {10}^{-4}\mbox{--}2\times {10}^{-2}$, which is consistent with our model for values of ${f}_{{\rm{cr}}}\sim 0.01\mbox{--}0.1$.

The other scenario for making UHECRs from en caul magnetars which we consider is if some NS mergers produce massive NSs rather than BHs. Indeed, there is now evidence that some NSs could have masses well above $2\,{M}_{\odot }$ (van Kerkwijk et al. 2011; Romani et al. 2012), which favors an equation of state that makes such a scenario more likely. Both population synthesis and observed Galactic NS binaries favor a rate that is comparable to or greater than the AIC rate. For example, Kim et al. (2005) estimate a Galactic rate of $\sim (0.1\mbox{--}3)\times {10}^{-4}\,{\mathrm{yr}}^{-1}$ or a volume rate of $\sim 4\times {10}^{-8}\mbox{--}2\,\times {10}^{-6}\,{\mathrm{Mpc}}^{-3}\,{\mathrm{yr}}^{-1}$. Again, this rate is consistent with what is needed for UHECRs for values of ${f}_{\mathrm{cr}}\sim 0.01$, taxing neither the NS merger rate nor the available energy reservoir of the en caul magnetar population.

The consistency arguments above imply that a fraction of rotational energy ${f}_{{\rm{cr}}}\approx 0.01$ goes into UHECRs. These rates are so uncertain, however, that it is reasonable to ask whether such a value can be understood from independent lines of argument. We now show that, indeed, the order-of-magnitude consistency arguments above agree nicely with more careful calculations.

To understand how representative this value of ${f}_{{\rm{cr}}}$ is requires a slightly more detailed model of the UHECR production rate. To obtain this, we first assume that the spectrum of UHECRs accelerated at any given time follows a ${dN}/{dE}\propto 1/E$ scaling, similar to previous discussions about pulsar magnetospheres (Arons 2003). This assumes a Goldreich–Julian injection (Goldreich & Julian 1969) and has mainly been chosen to illustrate the energy constraints on our model. In general, the spectrum could be different depending on the acceleration type (for example, shock acceleration; Lemoine et al. 2015) and other additional factors, such as the details of the source population and the magnetic field within the nebula which may also soften the spectrum (e.g., Fang et al. 2013). Next, we assume that a fraction of ${\eta }_{{\rm{cr}}}\lesssim 1$ of the magnetar spin down can be used for UHECR acceleration. Energy conservation then requires

$$\begin{eqnarray}&&{\int }_{{{\rm{E}}}_{{\rm{\min }}}}^{{E}_{{\rm{conf}}}}\displaystyle \frac{d\dot{N}}{{dE}}{EdE}\approx {\eta }_{{\rm{cr}}}{L}_{{\rm{sd}}},\end{eqnarray} \tag{ 27 }$$

where ${E}_{{\rm{\min }}}$ is the smallest energy for accelerated particles and can be ignored as long as ${E}_{{\rm{\min }}}\ll {E}_{{\rm{conf}}}$. Equation (27) can be solved to find

$$\begin{eqnarray}&&\displaystyle \frac{d\dot{N}}{{dE}}\approx \displaystyle \frac{{\eta }_{{\rm{cr}}}{L}_{{\rm{sd}}}}{{E}_{{\rm{conf}}}}\displaystyle \frac{1}{E}.\end{eqnarray} \tag{ 28 }$$

This describes the UHECRs produced within the nebula but, as we have discussed previously, they might not all escape the nebula. This is especially true at the highest energies due to synchrotron losses. To capture this physics, we only integrate $d\dot{N}/{dE}$ up to ${E}_{{\rm{synch}}}$ when ${E}_{{\rm{synch}}}\lt {E}_{{\rm{conf}}}$. The time-dependent energy output is then

$$\begin{eqnarray}{L}_{{\rm{cr}}}(t)=\left\{\begin{array}{cc}{\eta }_{{\rm{cr}}}{L}_{{\rm{sd}}}(t), & {E}_{{\rm{synch}}}\gt {E}_{{\rm{conf}}}\\ {\eta }_{{\rm{cr}}}{L}_{{\rm{sd}}}(t)\left(\tfrac{{E}_{{\rm{synch}}}(t)}{{E}_{{\rm{conf}}}(t)}\right), & {E}_{{\rm{synch}}}\lt {E}_{{\rm{conf}}}.\end{array}\right.\end{eqnarray} \tag{ 29 }$$

This is then time integrated over the evolution of the spinning down magnetar to obtain the total energy output from UHECRs. This is summarized in Figure 6 where we plot the contours of the fraction ${f}_{{\rm{cr}}}$ of ${E}_{{\rm{rot}}}$ that goes into UHECRs. This calculation uses ${\eta }_{{\rm{cr}}}=0.1$, which is not dissimilar to the fraction of pulsar power found to experience the full potential drop by Chen & Beloborodov (2014), and these results can simply be scaled to other values of ${\eta }_{{\rm{cr}}}$ because the dependence is simply linear. This shows that ${f}_{{\rm{cr}}}$ is greater than $\approx 0.01$ over a significant fraction of the parameter space.

Zoom In Zoom Out Reset image size

A further comparison is shown in Figure 7. Assuming that a given B_* and P_i are representative of the majority of magnetars that produce UHECRs, we estimate what event rate is needed to produce the observed UHECR rate, in units of the SN Ia rate. As discussed above, typical AIC and NS merger rates are $\lesssim 0.01$ in these units, and thus this scenario favors the lower left quadrant of Figure 7.

Zoom In Zoom Out Reset image size

UHECRs propagating through intergalactic space interact with the CMB or even other sources of extragalactic background light, which limits how far we can detect them. This is known as the GZK effect (Greisen 1966; Zatsepin & Kuz'min 1966). For both protons and Fe-like nuclei, the typical scale is similar and on the order of ${D}_{{\rm{GZK}}}\sim 170\,{\rm{Mpc}}$ (Metzger et al. 2011). Using the required rates estimated in Equation (26), we then estimate

$$\begin{eqnarray}&&{R}_{{\rm{GZK}}}\approx 23{I}_{45}^{-1}{\left(\displaystyle \frac{{P}_{i}}{2{\rm{ms}}}\right)}^{2}{\left(\displaystyle \frac{{f}_{{\rm{cr}}}}{0.01}\right)}^{-1}{\left(\displaystyle \frac{{D}_{{\rm{GZK}}}}{170{\rm{Mpc}}}\right)}^{3}\,{\mathrm{yr}}^{-1},\end{eqnarray} \tag{ 30 }$$

as the rate of events within the GZK volume. Such a distance scale is in fact interesting for other astrophysical interests. For example, this is similar to the horizon expected for detected NS mergers from ground-based interferometric detectors such as Advanced LIGO (Abadie et al. 2010). This will either constrain or provide a direct measurement of NS merging events within the GZK volume.

In addition, the apparent bolometric magnitude of the electromagnetic transient of such a UHECR-producing event would be $\approx 15.5$ (Metzger & Piro 2014) at $200\,{\rm{Mpc}}$. In principle, these events could have been detected by wide-field, transient surveys like the Palomar Transient Factory (Rau et al. 2009) and the Panoramic Survey Telescope and Rapid Response System (Kaiser et al. 2002). Unfortunately, given their particularly short timescales of $\sim 1\,{\rm{day}}$ and their blue color, it is possible that they have so far been missed. Also, note that the events best suited for producing UHECRs will have smaller magnetic fields and correspondingly fainter electromagnetic counterparts.

In the future, this will change with surveys that have particularly rapid cadences, such as the Zwicky Transient Facility (Law et al. 2009), the All-sky Automated Survey for Supernovae (Shappee et al. 2014), or the Large Synoptic Survey Telescope (LSST Science Collaboration et al. 2009, depending on the final cadence). In fact, such an event may have already been observed in GRB 080503 (Perley et al. 2009), which showed a rebrightening in both the optical and X-rays. Thus, perhaps another way to constrain the rate of these events is through GRBs (although the X-ray/optical signature is expected to be more isotropic than the beamed gamma-rays). Associated radio, either from the plerion associated with the magnetar (Piro & Kulkarni 2013) or from the interaction with the interstellar medium (Nakar & Piran 2011), may also provide constraints. Either way, it appears that these events may be observable by independent means on multiple different fronts if they are indeed making UHECRs (for example, see Metzger et al. 2015; Horesh et al. 2016). Within the next ∼3–4 years, we will have a much better empirical handle on the rates of such events and their role in UHECR acceleration.

UHECRs with large charges can be significantly deflected by the intergalactic magnetic field, and so it is a good check to consider this effect in the context of our model. As a particle with charge Ze travels over a magnetic field with correlation length λ, it typically is deflected by an angle $\alpha \sim d/{r}_{{\rm{L}}}$, where ${r}_{{\rm{L}}}=E/{{ZeB}}_{{\rm{IGM}}}$ is the Larmor radius and ${B}_{{\rm{IGM}}}$ is the intergalactic medium magnetic field. Traveling over a distance D, the typical total deflection will be

$$\begin{eqnarray}&&\theta \sim {\left(\displaystyle \frac{D}{\lambda }\right)}^{1/2}\displaystyle \frac{\lambda }{{r}_{{\rm{L}}}}.\end{eqnarray} \tag{ 31 }$$

This deflection needs to be less than $\sim 1\,{\rm{rad}}$ within the GZK volume; otherwise, the path length of a UHECR will increase too much. This would decrease the effective GZK volume to the point that it may cause problems obtaining a sufficient rate of events. This limits the magnetic field to be

$$\begin{eqnarray}{B}_{{\rm{IGM}}} & \lesssim & 8\times {10}^{-9}{Z}^{-1}{E}_{20}{\left(\displaystyle \frac{{D}_{{\rm{GZK}}}}{170{\rm{Mpc}}}\right)}^{-1/2}\\ & & \times {\left(\displaystyle \frac{\lambda }{1{\rm{Mpc}}}\right)}^{-1/2}\,{\rm{G}},\end{eqnarray} \tag{ 32 }$$

where ${E}_{20}=E/{10}^{20}\,{\rm{eV}}$.

Another consideration is that the intergalactic magnetic field should provide enough smoothing such that at least one source within the GZK volume is providing UHECRs at any given time, otherwise there would be large angular anisotropies observed in the UHECR rate. For this to be satisfied, the diffusion time, which is roughly $\tau \sim {\theta }^{2}D/c$, should be longer than ${R}_{{\rm{GZK}}}^{-1}$, which we derived in Equation (30) above. This yields

$$\begin{eqnarray}{B}_{{\rm{IGM}}} & \gtrsim & 7\times {10}^{-14}{Z}^{-1}{E}_{20}{I}_{45}^{1/2}{\left(\displaystyle \frac{{P}_{i}}{2{\rm{ms}}}\right)}^{-1}{\left(\displaystyle \frac{{f}_{{\rm{cr}}}}{0.01}\right)}^{1/2}\\ & & \times {\left(\displaystyle \frac{{D}_{{\rm{GZK}}}}{170{\rm{Mpc}}}\right)}^{-5/2}{\left(\displaystyle \frac{\lambda }{1{\rm{Mpc}}}\right)}^{-1/2}\,{\rm{G}},\end{eqnarray} \tag{ 33 }$$

as a lower limit for the intergalactic magnetic field. Interestingly, these upper and lower limits show that, in principle, measuring the properties of the intergalactic magnetic field could help to constrain the properties of the magnetars generating the UHECRs, although the wide range of constraints would probably make this difficult in practice.

One is finally left to wonder, "What if an AIC occurred in the Milky Way?" Similar to the estimates made above, the propagation time across the disk of the Galaxy is

$$\begin{eqnarray}\tau & \sim & 58{Z}^{2}{E}_{20}^{-2}{\left(\displaystyle \frac{D}{15{\rm{kpc}}}\right)}^{2}\left(\displaystyle \frac{\lambda }{0.1\,{\rm{kpc}}}\right)\\ & & \times {\left(\displaystyle \frac{{B}_{{\rm{ISM}}}}{3\times {10}^{-6}{\rm{G}}}\right)}^{2}{\rm{year}}.\end{eqnarray} \tag{ 34 }$$

For Z = 26, this implies $\tau \sim 4\times {10}^{4}\,{\rm{year}}$, and thus the diffusion time for heavy nuclei is longer than the timescale between AIC or NS merger events in the Galaxy. Therefore, a non-zero fraction of UHECRs that we measure today could have originated from an event in the Milky Way. The deflection angle is

$$\begin{eqnarray}\theta & \sim & 5\times {10}^{-2}{{ZE}}_{20}^{-1}{\left(\displaystyle \frac{D}{15{\rm{kpc}}}\right)}^{1/2}\\ & & \times {\left(\displaystyle \frac{\lambda }{0.1{\rm{kpc}}}\right)}^{1/2}\left(\displaystyle \frac{{B}_{{\rm{ISM}}}}{3\times {10}^{-6}\,{\rm{G}}}\right){\rm{rad}},\end{eqnarray} \tag{ 35 }$$

which, for Z = 26, implies deflection angles of $\sim 50^\circ$, which is obviously not in the small angle approximation from which this formula was derived! The conclusion is that either $Z\sim 1$ with very little contribution from the Galaxy, or Z = 26 and, despite some contribution from the Galaxy, associating them with locations in the Milky Way will be difficult due to their large deflections from their origins.

Hadronic interactions between the UHECRs and the ejecta surrounding the magnetar will lead to charged pion production, which subsequently decay into high-energy neutrinos via ${\pi }^{\pm }\to {e}^{\pm }+{\nu }_{e}({\bar{\nu }}_{e})+{\nu }_{\mu }+{\bar{\nu }}_{\mu }$. Although a detailed study of this process is outside the scope of this paper, we make a simple estimate of the expected neutrino background to check for consistency with constraints from the IceCube Observatory (The IceCube Collaboration 2011; Aartsen et al. 2013).

The optical depth of UHECRs to hadron interactions, which we estimated in Section 3.2 in the context of preventing UHECRs from escaping, is $\sim {t}_{{\rm{ej}}}/{t}_{{\rm{cross}}}$. Thus, we estimate the efficiency of pion production as

$$\begin{eqnarray}&&{f}_{\pi }\sim \min ({t}_{{\rm{ej}}}/{t}_{{\rm{cross}}},1).\end{eqnarray} \tag{ 36 }$$

Furthermore, we assume that 0.2 of the energy of a given UHECR can be converted to the energy of a pion, and 0.25 of the pion energy then goes into the neutrino (as was done in Fang et al. 2016). This provides us with an overall conversion factor of ${f}_{\nu }\sim 0.05$ from a given UHECR to a neutrino.

Integrating over the entire UHECR energy spectrum, we can derive the total energy of neutrinos that are generated. Note that in this particular case, f_π can be taken outside of the energy integral since we are using a simplistic model for the hadronic interactions, which may not be true in a more detailed treatment. This greatly simplifies the evaluation of the integral, and the total neutrino luminosity comes out to

$$\begin{eqnarray}&&{L}_{\nu }(t)\sim {f}_{\nu }{f}_{\pi }(t){\eta }_{{\rm{cr}}}{L}_{{\rm{sd}}}(t).\end{eqnarray} \tag{ 37 }$$

Note that here we are using the entire luminosity that goes into UHECRs, which is represented by the term ${\eta }_{{\rm{cr}}}{L}_{{\rm{sd}}}$, rather than just the fraction that escapes the magnetar environment, which is given in Equation (29). This is because it is this luminosity of UHECRs that is important for generating the high-energy neutrinos. To relate this to the observed high-energy neutrino flux, we integrate L_ν over the entire time for a single event and over a distance D, resulting in

$$\begin{eqnarray}&&{F}_{\nu }=\displaystyle \frac{{DR}}{4\pi }{\int }_{0}^{\infty }{L}_{\nu }(t){dt},\end{eqnarray} \tag{ 38 }$$

where R is the volumetric rate of events. We do not include additional neutrinos produced by secondary nuclei and pions since our ejecta are so much less massive than in supernova-related scenarios (Murase et al. 2009).

The calculation of this high-energy neutrino background is summarized in Figure 8. At each location in this plot, we assume that B_* and P_i are representative of all UHECR-producing magnetars. Furthermore, we set R to be the rate needed to match the observed UHECR rate (as summarized in Figure 7). Finally, we use $D=4000\,{\rm{Mpc}}$, essentially assuming a constant rate of UHECR production over the entire history of the universe. A more detailed treatment that includes the time-dependent cosmic star formation rate and connects this to the rate of en caul magnetar births over time is outside the scope of this work.

Zoom In Zoom Out Reset image size

The IceCube upper limits are resolved as a function of energy and reach a minimum of $\sim {10}^{-8}\,\mathrm{GeV}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1}\,{\mathrm{sr}}^{-1}$ at a neutrino energy of $\sim {10}^{14}\,{\rm{eV}}$. Furthermore, going up to energies of $\sim {10}^{18}\mbox{--}{10}^{20}\,{\rm{GeV}}$, the upper limits are $\sim 3\times {10}^{-8}\mbox{--}3\times {10}^{-7}\,\mathrm{GeV}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1}\,{\mathrm{sr}}^{-1}$. Compared to Figure 8, this is above our expected background rate over most of the parameter space, and we note that we are integrating over neutrinos of all energies, and thus likely overestimating the potential background. The reason we are able to meet the neutrino constraints is that the low amount of ejecta naturally gives ${t}_{{\rm{ej}}}/{t}_{{\rm{cross}}}\ll 1$ over most of the time when UHECRs are produced in the greatest amounts. Thus, for the same reason that we are able to efficiently produce UHECRs, we are also able to not overproduce high-energy neutrinos. We note that there may also be photo-hadronic interactions within the nebula that may also produce neutrinos (Lemoine et al. 2015; Fang et al. 2016). Since this process is stronger for protons than heavy nuclei, and since our main focus here is on the heavy nuclei as UHECRs, we save such a calculation for future work where the spectra of UHECRs and the non-thermal radiation is computed in more detail.