REVISITING THE CONTRIBUTIONS OF SUPERNOVA AND HYPERNOVA REMNANTS TO THE DIFFUSE HIGH-ENERGY BACKGROUNDS: CONSTRAINTS ON VERY HIGH REDSHIFT INJECTION

After the early free expansion, the evolution of an SNR enters the Sedov–Taylor self-similar expansion phase (Taylor 1950; Sedov 1959). The ejected remnant shell collides with and sweeps up the ambient medium, forming a shock that can accelerate CRs that subsequently, via pp collisions, produce high-energy neutrinos and γ-rays. We take the beginning of the Sedov–Taylor phase as the onset time for CR acceleration and secondary neutrino and γ-ray production. This occurs at a time

$$\begin{eqnarray}{t}_{0}=\displaystyle \frac{{R}_{{\rm{dec}}}}{{v}_{\mathrm{ej}}} & = & \displaystyle \frac{{(3{M}_{{\rm{ej}}}/4\pi {m}_{p}{n}_{0})}^{1/3}}{{(2{{ \mathcal E }}_{{\rm{SN}}}/{M}_{{\rm{ej}}})}^{1/2}}\\ & & \simeq 1.4\times {10}^{3}{{ \mathcal E }}_{{\rm{SN,51}}}^{-1/2}{n}_{0,0}^{-1/3}{M}_{{\rm{ej,1}}}^{5/6}\,{\rm{yr}}\end{eqnarray} \tag{ 1 }$$

after the initial explosion. The ending time of the Sedov–Taylor phase t_end occurs when the adiabatic expansion approximation is no longer true, and this is determined by equating the cooling time of the post-shock region with the SNR dynamic time. The gas cooling function has been evaluated in detail by Sutherland & Dopita (1993) and has been adopted in our calculation to solve for t_end as a function of different sets of parameters.

We take the usual time dependence of the radius and the shock velocity in the Sedov–Taylor phase $R\,={(25{{ \mathcal E }}_{{\rm{SN}}}/4\pi {m}_{p}{n}_{0})}^{1/5}{t}^{2/5}$ and $v=(2/5){(25{{ \mathcal E }}_{{\rm{SN}}}/4\pi {m}_{p}{n}_{0})}^{1/5}{t}^{-3/5}$ (Taylor 1950; Sedov 1959; Ramirez-Ruiz & MacFadyen 2010), and the commonly used estimate for the post-shock magnetic field $B={(9\pi {\epsilon }_{B}{m}_{p}{n}_{0}{v}^{2})}^{1/2}$, where ${\epsilon }_{B}\leqslant 1$ is the ratio of the post-shock magnetic to thermal energy. In this work, we will assume the CR component to be pure protons. Then, for diffusive shock acceleration, the maximum proton energy is estimated by equating the local acceleration time to the dynamical time (the cooling timescale for pp collisions and synchrotron losses being much longer than the dynamical timescale; e.g., Senno et al. 2015), which gives $(20{\epsilon }_{p,\max }c/3{eB})={Rv}$ (Drury 2011). Thus, in the Sedov–Taylor phase the maximum CR energy also depends on time,

$$\begin{eqnarray}{\epsilon }_{p,\max } & = & 1.58\times {10}^{6}{{ \mathcal E }}_{\mathrm{SN},51}{n}_{0,0}^{1/6}{M}_{{\rm{ej,1}}}^{-2/3}\\ & & \times {\epsilon }_{B,-2}^{1/2}{(t/{t}_{0})}^{-4/5}\,{\rm{GeV}}.\end{eqnarray} \tag{ 2 }$$

In the free expansion phase, for an external density ρ = constant, as long as the swept-up mass is small, the velocity v ≃ constant and the radius grows as $R\propto t$. The rate at which mass is swept and shock-heated is $\dot{M}\propto {R}^{2}\rho v\propto {R}^{2}\propto {t}^{2}$, and the bolometric luminosity that can go into CRs grows as $L\propto \dot{M}\propto {t}^{2}$. The equipartition post-shock field behaves as $B\propto v\,\propto$ constant, and the shock strength is also constant. Thus, at the simplest level one expects the CR bolometric luminosity to grow as ${L}_{{\rm{CR}}}\propto {t}^{2}$. In the Sedov phase, on the other hand, $R\propto {t}^{2/5}$ and $v\propto {t}^{-3/5}$, so one expects $\dot{M}\propto {t}^{1/5}$, while $B\propto v\propto {t}^{-3/5}$, and with decreasing velocity the shock strength decreases. Thus, one does not expect a positive proportionality between $\dot{M}$ and ${L}_{{\rm{CR}}}$, but instead a decrease of L_CR can be expected. The SN ejecta begins to slow down at t₀ when it has swept up an amount of external mass comparable to that of the ejecta shell, marking the beginning of the Sedov–Taylor phase, when the CR acceleration is expected to peak. The details of the continued CR injection as the blast wave slows down in the Sedov–Taylor phase are model dependent, and here we shall parameterize it as

$$\begin{eqnarray}&&{L}_{{\rm{CR}}}(t)={ \mathcal A }{(t/{t}_{0})}^{\alpha },\end{eqnarray} \tag{ 3 }$$

where ${ \mathcal A }$ is a normalization factor determined by ${\int }_{{t}_{0}}^{{t}_{{\rm{end}}}}{L}_{{\rm{CR}}}{dt}=\eta {{ \mathcal E }}_{{\rm{SN}}}$. Here $\eta \sim 0.1$ is the CR acceleration efficiency (Caprioli et al. 2015). According to Equation (3), at different stages of the Sedov–Taylor phase, the energy injection rate into CRs will be different. The value α = −1 is a nominal steady injection case (Ohira et al. 2011). Because of the very fast growth ${L}_{{\rm{CR}}}\propto {t}^{2}$ in the initial free expansion phase, for simplicity we neglect contributions from this phase. We shall assume the shock-accelerated proton spectrum to be a power law ${{dN}}_{p}/d{\epsilon }_{p}\propto {\epsilon }_{p}^{-p}$ with index p = 2, so at time t, the proton spectrum can be written as

$$\begin{eqnarray}{\epsilon }_{p}{L}_{{\epsilon }_{p}}(t)=\left\{\begin{array}{ll}\tfrac{{L}_{{\rm{CR}}}(t)}{{ \mathcal C }} & {\rm{if}}\ {\epsilon }_{p}\leqslant {\epsilon }_{p,\max }(t)\\ \tfrac{{L}_{{\rm{CR}}}(t)}{{ \mathcal C }}{e}^{-{\epsilon }_{p}/{\epsilon }_{p,\max }} & {\rm{if}}\ {\epsilon }_{p}\gt {\epsilon }_{p,\max }(t)\end{array}\right.,\end{eqnarray} \tag{ 4 }$$

where ${ \mathcal C }=\mathrm{ln}({\epsilon }_{p,\max }/{\epsilon }_{p,\min })$.

Hypernovae (HNe) are a subtype of SNe that are more energetic but relatively rarer. Recent studies indicate that the total rate of all core-collapse SNe is ${{ \mathcal R }}_{{\rm{CCSNe}}}\,$ $=\,(1.06\ \pm 0.11(\mathrm{stat})\pm 0.15(\mathrm{sys}))\,$ $\times \,{10}^{-4}$ ${({h}_{0}/0.7)}^{3}\,$ ${\mathrm{Mpc}}^{-3}{\mathrm{yr}}^{-1}$ (Taylor et al. 2014b). For typical HNe, the kinetic energy is of order ${10}^{52\mbox{--}53}\,{\rm{erg}}$, with a typical rate ${{ \mathcal R }}_{{\rm{HNe}}}\leqslant 4 \% {{ \mathcal R }}_{{\rm{CCSNe}}}$ (Guetta & Della Valle 2007; Arcavi et al. 2010; Senno et al. 2015), but with substantial uncertainties. HNe have been discussed as possible sources for high-energy CRs (Dermer 2002; Sveshnikova 2003; Wang et al. 2007; Ioka & Mészáros 2010). Due to larger kinetic energy, HNe can lead to larger maximum proton energies and will dominate the neutrino flux at PeV energies.

The total local CR energy budget is

$$\begin{eqnarray}{\epsilon }_{p}{Q}_{{\epsilon }_{p}}(z=0) & = & {{ \mathcal R }}_{{\rm{SNe}}}{\displaystyle \int }_{{t}_{{\rm{0,SNe}}}}^{{t}_{{\rm{end,SNe}}}}{\epsilon }_{p}{L}_{{\epsilon }_{p}}^{{\rm{SNe}}}(t){dt}\\ & & +{{ \mathcal R }}_{{\rm{HNe}}}{\displaystyle \int }_{{t}_{{\rm{0,HNe}}}}^{{t}_{{\rm{end,HNe}}}}{\epsilon }_{p}{L}_{{\epsilon }_{p}}^{{\rm{HNe}}}(t){dt},\end{eqnarray} \tag{ 5 }$$

and the cosmological evolution of sources can be expressed as ${\epsilon }_{p}{Q}_{{\epsilon }_{p}}(z)={\epsilon }_{p}{Q}_{{\epsilon }_{p}}(z=0)S(z)$, in which the evolution term is $S(z)={\left[{(1+z)}^{-34}+{\left(\tfrac{1+z}{5000}\right)}^{3}+{\left(\tfrac{1+z}{9}\right)}^{35}\right]}^{-0.1}$ (Hopkins & Beacom 2006; Yüksel et al. 2008). In this work, the typical neutrino energy is approximated to be ${E}_{\nu }\sim 0.05{\epsilon }_{p}/(1+z)$.

The diffuse neutrino flux per flavor is calculated as (e.g., Murase et al. 2013, 2016; Bechtol et al. 2015; Senno et al. 2015)

$$\begin{eqnarray}{E}_{\nu }^{2}{{\rm{\Phi }}}_{{\nu }_{i}} & = & \displaystyle \frac{c}{4\pi {H}_{0}}{\displaystyle \int }_{0}^{{z}_{{\rm{\max }}}}\displaystyle \frac{\min [1,{f}_{{pp}}]}{6}\\ & & \times \displaystyle \frac{{\epsilon }_{p}{Q}_{{\epsilon }_{p}}(z=0)S(z)}{{(1+z)}^{2}\sqrt{{{\rm{\Omega }}}_{M}{(1+z)}^{3}+{{\rm{\Omega }}}_{{\rm{\Lambda }}}}}{dz}\end{eqnarray} \tag{ 6 }$$

where z_max = 4, the cosmology parameters are ${H}_{0}=67.8\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{Mpc}}^{-1},{{\rm{\Omega }}}_{M}=0.308$ (Planck Collaboration 2015), and ${f}_{{pp}}=\bar{n}\kappa {\sigma }_{{pp}}c\cdot \min [{t}_{{\rm{diff}}},{t}_{{\rm{adv}}}]$. Here the diffusive escape timescale of CRs from the host galaxy is ${t}_{{\rm{diff}}}={h}^{2}/4D$, where h is the galaxy scale height, the diffusion coefficient is $D({\epsilon }_{p})={D}_{c}[{({\epsilon }_{p}/{\epsilon }_{p,c})}^{1/3}+{({\epsilon }_{p}/{\epsilon }_{p,c})}^{2}]$, and the timescale of advective escape is ${t}_{{\rm{adv}}}=h/{V}_{w}$. We use the values of $h={10}^{21}\,{\rm{cm}}$, ${D}_{c}\sim 3.1\times {10}^{29}\,{\mathrm{cm}}^{2}\,{{\rm{s}}}^{-1}$, ${\epsilon }_{p,c}\sim 9.3\times {10}^{9}{\rm{GeV}}$, ${V}_{w}=1500\,\mathrm{km}\,{{\rm{s}}}^{-1}$, $B=4\,\mathrm{mG}$, $\bar{n}=100\,{\mathrm{cm}}^{-3}$ for SBGs (Strickland & Heckman 2009) and h = 1 kpc, ${D}_{c}\sim 3.1\times {10}^{29}\,{\mathrm{cm}}^{2}\,{{\rm{s}}}^{-1}$, ${\epsilon }_{p,c}\sim 9.3\times {10}^{6}\,{\rm{GeV}}$, ${V}_{w}=500\,\mathrm{km}\,{{\rm{s}}}^{-1}$, $B=1\,\mu {\rm{G}}$, $\bar{n}=1\,{\mathrm{cm}}^{-3}$ for SFGs (Keeney et al. 2006; Crocker 2012). We assume a starburst fraction of ${\xi }_{{\rm{SBG}}}=0.1$ (Tamborra et al. 2014). Note that the upstream magnetic field of SBGs can be stronger than in the usual ISM due to the superbubble structure. Also, the upstream field amplification by CRs may be expected (Beresnyak et al. 2009). A further contribution to f_pp from diffusion in the host galaxy cluster is considered in Senno et al. (2015).

Based on the standard branching ratio between charged and neutral pion production of pp interaction, the gamma-ray flux can initially be related to the neutrino flux through ${E}_{\gamma }^{2}{{\rm{\Phi }}}_{\gamma }=2{E}_{\nu }^{2}{{\rm{\Phi }}}_{\nu }{| }_{{\epsilon }_{\nu }=0.5{\epsilon }_{\gamma }}$. However, once they are produced, high-energy gamma rays undergo $\gamma \gamma$ interactions. The main target soft photons are from the extragalactic background light (EBL; e.g., Salamon & Stecker 1998; Dai & Lu 2002; Dai et al. 2002; Murase et al. 2012). The optical depth for gamma rays depends on their energy ${\epsilon }_{\gamma }$, as well as on redshift, because the EBL intensity varies at different z (Finke et al. 2010). We can set ${\tau }_{\gamma \gamma }({E}_{\gamma },z)=1$ to get the cutoff energy ${E}_{\gamma }^{{\rm{cut}}}$. Beyond ${E}_{\gamma }^{{\rm{cut}}}$, electromagnetic cascade develops, and the resulting spectrum has a universal form (Berezinskii & Smirnov 1975; Murase et al. 2012),

$$\begin{eqnarray}{E}_{\gamma }\tfrac{{{dN}}_{\gamma }}{{{dE}}_{\gamma }}\propto \left\{\begin{array}{lc}{\left(\tfrac{{E}_{\gamma }}{{E}_{\gamma }^{{\rm{br}}}}\right)}^{-1/2} & {\rm{if}}\ {E}_{\gamma }\leqslant {E}_{\gamma }^{{\rm{br}}}\\ {\left(\displaystyle \frac{{E}_{\gamma }}{{E}_{\gamma }^{{\rm{br}}}}\right)}^{-1} & {\rm{if}}\ {E}_{\gamma }^{{\rm{br}}}\lt {E}_{\gamma }\leqslant {E}_{\gamma }^{{\rm{cut}}}\end{array}\right.\end{eqnarray} \tag{ 7 }$$

where ${E}_{\gamma }^{{\rm{br}}}=0.0085{(1+z)}^{2}{\left(\tfrac{{E}_{\gamma }^{{\rm{cut}}}}{100\mathrm{GeV}}\right)}^{2}$ (Senno et al. 2015). The normalization is determined by integrating the gamma-ray flux over ${E}_{\gamma }^{{\rm{cut}}}$.

Below ${E}_{\gamma }^{{\rm{cut}}}$, the gamma-ray intensity is attenuated by a factor ${e}^{-{\tau }_{\gamma \gamma }}$. This is calculated similarly to the diffuse neutrino flux, with a typical energy ${E}_{\gamma }\sim 0.1{\epsilon }_{p}/(1+z)$,

$$\begin{eqnarray}{E}_{\gamma }^{2}{{\rm{\Phi }}}_{\gamma } & = & \displaystyle \frac{c}{4\pi {H}_{0}}{\displaystyle \int }_{0}^{{z}_{{\rm{\max }}}}\displaystyle \frac{\min [1,{f}_{{pp}}]}{3}\\ & & \times \displaystyle \frac{{\epsilon }_{p}{Q}_{{\epsilon }_{p}}(z=0)S(z)}{{(1+z)}^{2}\sqrt{{{\rm{\Omega }}}_{M}{(1+z)}^{3}+{{\rm{\Omega }}}_{{\rm{\Lambda }}}}}{e}^{-{\tau }_{\gamma \gamma }({E}_{\gamma },z)}{dz}.\end{eqnarray} \tag{ 8 }$$

The final diffuse gamma-ray flux is the sum of these two components.

Here we use the above-discussed values of the luminosity functions, the EBL density, and other parameters that are commonly used for redshifts $z\lesssim 4$. Figure 1 shows a combined fit of the diffuse gamma-ray and neutrino flux for the nominal case of $\alpha =-1$ in Equation (3). We use this nominal value to delimit the physically expected range of $\alpha \lesssim 0.5$. The IceCube neutrino and Fermi-LAT EGB observations are shown by the blue and red data points, respectively. The black solid line is the predicted total diffuse neutrino flux, and the red solid line represents the gamma-ray flux. The black dashed and dotted lines represent the contribution to the neutrino flux from SNe and HNe, respectively. We can see that HNe dominate the PeV neutrino flux, while SNe only contribute up to sub-PeV neutrino energies. However, this case based on the conventional assumptions cannot explain the IceCube data, suggesting that the HN rate would need to be enhanced to explain the neutrino data. This difficulty in getting a good combined $\nu ,\gamma$ fit for this conventional case is compatible with the independent generic arguments of Bechtol et al. (2015) against SBGs being the dominant sources of IceCube neutrinos.

Zoom In Zoom Out Reset image size

Obviously, we would expect to increase the fraction of HN contribution because they can achieve higher CR energies. Figure 2 shows what happens if we increase the HN kinetic energy to ${10}^{52.5}$ erg. The blue solid line is the diffuse neutrino flux, which satisfies the observations better now, although it is still below the 1σ range of the last data point. The SN contribution is the blue dashed line, and the HN contribution is the blue dotted line. Except for ${{ \mathcal E }}_{{\rm{HNe}}}$, all other parameters are the same for the black lines in this figure. The corresponding gamma-ray flux is the green solid line, which is still in the allowed range. In this case, the ratio between the HN and SN contribution is $\zeta \equiv \tfrac{\eta {{ \mathcal E }}_{{\rm{HN}}}{{ \mathcal R }}_{{\rm{HN}}}}{\eta {{ \mathcal E }}_{{\rm{SN}}}{{ \mathcal R }}_{{\rm{SN}}}}=\tfrac{{10}^{52.5}\times 3 \% }{5\times {10}^{50}\times 97 \% }\sim 2$, compared to $\zeta \sim 0.6$ for the ${{ \mathcal E }}_{{\rm{HN}}}={10}^{52}$ erg case. We can draw a rough conclusion that $\zeta \gtrsim 1$ is needed for a good fit.

Zoom In Zoom Out Reset image size

In addition, we consider more realistic cases where ${L}_{{\rm{CR}}}\propto {t}^{\alpha }$, for various α. We plot the combined fits for different values of α in Figure 3. The neutrino spectrum shape is seen to depend on α. Looking at the results for $\alpha =0$ and for positive values $\alpha =1,2$ (although the positive values are unlikely to occur), we see that all three neutrino spectral curves decrease sharply at the high-energy end, undershooting the data. One way to address this might be to increase the rate of HNe (or to choose a very large SBG fraction ${\xi }_{{\rm{SBG}}}$), thus increasing the energy input of HNe so as to reach the IceCube data, but this inevitably gives rise to an increased gamma-ray flux, which would overshoot the Fermi data. However, for the physically more reasonable negative α values, the neutrino spectrum is "flatter" (since at later times in the HN evolution both the maximum energy and the flux decrease) and more energy is evident in the flux at the high-energy end, so we do not need such a high fraction of HNe as before. We can see that the last IceCube data point is reachable for $\alpha =-2$ and $\alpha =-3$ using a quite conservative HN fraction ${{ \mathcal R }}_{{\rm{HNe}}}=3 \% {{ \mathcal R }}_{{\rm{CCSNe}}}$. Clearly, negative values of α are more favorable from an observational point of view.

Zoom In Zoom Out Reset image size

Generally, the parameter space gets larger for smaller α. The major parameters that we can adjust in our fits include the SBG upstream magnetic field; the kinetic energies ${{ \mathcal E }}_{{\rm{SN}}},\,{{ \mathcal E }}_{{\rm{HN}}};$ the pre-shock medium number density n₀; the HN fraction ${{ \mathcal R }}_{{\rm{HN}}};$ and the CR acceleration efficiency η. The nominal values used here are in the observationally reasonable range, e.g., $B=4\,\mathrm{mG}$ (e.g., Thompson et al. 2006; Robishaw et al. 2008; Beck 2016), ${{ \mathcal E }}_{{\rm{SN}}}\sim {10}^{50\mbox{--}51}\,\mathrm{erg},{{ \mathcal E }}_{{\rm{HN}}}\sim {10}^{52\mbox{--}53}\,\mathrm{erg}$ (e.g., Iwamoto et al. 1998; Mazzali et al. 2003), ${n}_{0}\sim 1\,{\mathrm{cm}}^{-3}$(Draine & Woods 1991; Chevalier & Claes 2001; Koo et al. 2007; Fox et al. 2011), ${{ \mathcal R }}_{{\rm{HNe}}}\leqslant 0.04{{ \mathcal R }}_{{\rm{CCSNe}}}$ (Guetta & Della Valle 2007; Arcavi et al. 2010), and $\eta \lt 1$ (physically required). For illustration, Figure 4 shows what happens if we change the SBG upstream magnetic field. A stronger field helps to explain the observations. The key issue here is that we wish to achieve higher maximum proton energies in order to account for PeV neutrino flux.

Zoom In Zoom Out Reset image size

Alternatively, given the contribution of other sources such as blazars to the Fermi gamma-ray background (Bechtol et al. 2015; Fermi-LAT Collaboration 2016), we can provide limits on the SNR contributions to the neutrino flux. If we normalize the gamma-ray flux to the best-fit 14% of the total EGB, we find that for our nominal parameters the corresponding diffuse neutrino flux is about 50% of the observed value, provided that we choose not to satisfy the ∼30 TeV neutrino data point, which, however, has small error bars. If this data point has a different physical origin than the higher-energy data points, SBGs may still contribute a potentially interesting fraction of the $\gtrsim 100$ TeV IceCube neutrino flux.

One may wonder whether very high energy gamma rays may be attenuated inside SFGs and SBGs. This could happen above ∼10 TeV energies (Lacki & Thompson 2013; Murase et al. 2013; Chang & Wang 2014). However, this would not change the results essentially. The synchrotron losses of the cascade electron–positron pairs lead to a reduction in the diffuse gamma-ray flux. However, gamma rays in the 0.1–10 TeV range remain unattenuated in the sources and will be injected into intergalactic cascades. Also, electron–positron pairs including the contribution from muon decays, as well as primary electrons, produce some gamma rays, which can make the constraints tighter (Lacki et al. 2014).

An interesting but more uncertain contribution to the diffuse neutrino and gamma-ray background can be expected from the Pop-III star formation occurring at high redshifts (e.g., Umeda & Nomoto 2003; Iwamoto et al. 2005). For these epochs, besides the increased redshift of the received photons, also (other things being equal) the ${(1+z)}^{3}$ increase of the proper densities affects the $\gamma \gamma$ interaction and the gamma-ray attenuation. Pop-III stars are extremely metal-poor and can be very massive (e.g., Abel et al. 2002; Hosokawa et al. 2011; Whalen et al. 2013; Yang et al. 2015). They may end their life as HN explosions, ${{ \mathcal E }}_{{\rm{kin}}}={10}^{51\mbox{--}53}\,\mathrm{erg}$ (e.g., Chen 2014; Tominaga et al. 2014; Toma et al. 2016), and such Pop-III SNe, including jet-driven SNe and GRBs, are also expected to accelerate CRs and produce neutrinos (Schneider et al. 2002; Iocco et al. 2008; Gao et al. 2011; Berezinsky & Blasi 2012).

The theoretical uncertainties for this Pop-III epoch are significant, but we can assume reasonable HN parameters used in the literature for illustrative purposes. For instance, the Pop-III star formation rate is under debate (Wise & Abel 2005; Bromm & Loeb 2006; Trenti & Stiavelli 2009; de Souza et al. 2011; Wise et al. 2012) and can be 10–100 times smaller (Trenti & Stiavelli 2009; Wise et al. 2012). Here we use the HN explosion parameters of Chen (2014), Tominaga et al. (2014), and the EBL model and the upper Pop-III model given in Inoue et al. (2013), with a typical HN kinetic energy of ${{ \mathcal E }}_{{\rm{kin}}}={10}^{52.5}$ erg and an external gas density of ${n}_{0}={10}^{3}$ cm⁻³, the latter value being compatible with gas densities inferred from the $z\sim 6.3$ GRB afterglow analysis of Gou et al. (2007). Neglecting at first entirely the contribution of SNe and HNe at $z\leqslant 4$, we plot in Figure 5 the contribution of only these Pop-III HNe to the neutrino and gamma-ray flux received from the redshift range $4\leqslant z\leqslant 10$, using the $\alpha =-1$ case as an example.

Zoom In Zoom Out Reset image size

From Figure 5 it is seen that such Pop-III HNe would likely make only a minor contribution to the diffuse neutrino flux. However, if the explosion energy is larger, ${{ \mathcal E }}_{\mathrm{POP} \mbox{-} \mathrm{III}}\sim {10}^{53}$ erg, these Pop-III HNe by themselves may be able to fit the neutrino data well, without violating the Fermi-LAT 14% residual background. Actually, there is still a larger parameter space for this fitting, as evidenced by the upper bound shown in Figure 5, making Pop-III stars more interesting sources for the present purposes. Alternatively, our results suggest that their explosion energy is limited to ${{ \mathcal E }}_{\mathrm{POP} \mbox{-} \mathrm{III}}\lesssim 7\times {10}^{53}$ erg.

Based on this Pop-III model, we also explore the influence that the initial CR spectral index p has on the results. Figure 6 shows the case p = 2.2. Solid lines are for $\eta =0.1$, and dashed lines are for an enhancement by a factor of 20, while other parameters remain unchanged from Figure 5. We see that at high energies the upper bound is lower than that of the p = 2 case. This means that for p = 2.2, under the same constraint given by Fermi-LAT, less neutrino flux can be attributed to the Pop-III stars. From this point of view, the smaller p are more favorable. We can give a rough constraint of $p\leqslant 2.2$.

Zoom In Zoom Out Reset image size

Of course, there is no reason to assume that only Pop-III HNRs contribute neutrinos in the IceCube energy range, so a more realistic scenario would involve at least two components contributing to the neutrino flux, one being the standard SNe and HNe in the usual redshift range $z\leqslant 4$, and the other being the Pop-III HNe in the $4\leqslant z\leqslant 10$ range. The occurrence rates in the lower redshift range are observationally better constrained than for the Pop-III component, but the ratio between the two has large uncertainties. Since ${E}_{\nu }^{2}{{\rm{\Phi }}}_{{\nu }_{i}}\propto \eta { \mathcal E }{ \mathcal R }$, the ratio between the Pop-III and Pop-I/II CR contribution can be written as

$$\begin{eqnarray}\displaystyle \frac{{({E}_{\nu }^{2}{{\rm{\Phi }}}_{{\nu }_{i}})}_{{\rm{III}}}}{{({E}_{\nu }^{2}{{\rm{\Phi }}}_{{\nu }_{i}})}_{{\rm{I/II}}}} & = & \displaystyle \frac{\eta {{ \mathcal E }}_{{\rm{HN}}}{{ \mathcal R }}_{{\rm{III}}}({\bar{z}}_{{\rm{III}}})}{\eta {{ \mathcal E }}_{{\rm{SN}}}0.97{{ \mathcal R }}_{{\rm{I/II}}}({\bar{z}}_{{\rm{I/II}}})+\eta {{ \mathcal E }}_{{\rm{HN}}}0.03{{ \mathcal R }}_{{\rm{I/II}}}({\bar{z}}_{{\rm{I/II}}})}\\ & & \times \displaystyle \frac{{(1+{\bar{z}}_{{\rm{I/II}}})}^{2}\sqrt{{{\rm{\Omega }}}_{M}{(1+{\bar{z}}_{{\rm{I/II}}})}^{3}+{{\rm{\Omega }}}_{{\rm{\Lambda }}}}}{{(1+{\bar{z}}_{{\rm{III}}})}^{2}\sqrt{{{\rm{\Omega }}}_{M}{(1+{\bar{z}}_{{\rm{III}}})}^{3}+{{\rm{\Omega }}}_{{\rm{\Lambda }}}}}\displaystyle \frac{{f}_{{pp}}^{{\rm{III}}}}{{f}_{{pp}}^{{\rm{I/II}}}}.\end{eqnarray} \tag{ 9 }$$

Here (as seen also below in Figure 8) the first term is ∼0.3, where the typical redshifts of Pop-I/II and Pop-III remnants are ${\bar{z}}_{{\rm{I/II}}}\sim 1$ and ${\bar{z}}_{{\rm{III}}}\sim 5$, respectively. The second term gives ∼0.02, and the last term is ∼20, leading to a ratio of ∼0.1. The star formation rate will be lower for Pop-III than it is at low redshifts. On the other hand, $\eta { \mathcal E }{ \mathcal R }(z)$ could be substantially larger due to a plausibly larger kinetic energy. Thus, introducing this Pop-III contribution potentially allows one to produce a fraction approaching unity of the low- and high-redshift IceCube neutrino flux, without violating the Fermi residual gamma-ray background. An example with these two components is shown in Figure 7, where the Pop-III/Pop-I/II efficiency is set to the same $\eta =0.1$. We can see that Pop-III HNe can contribute an interesting fraction to the diffuse neutrino flux, providing a good argument in favor of the SBG scenario.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 8 shows the CR input power ${\epsilon }_{p}{Q}_{{\epsilon }_{p}}$ versus redshift for this fit, the Pop-III being subdominant (red solid line). However, with a higher density $\bar{n}$, an ${f}_{{pp}}=1$ (at the low-energy end) is obtained for Pop-III, while ${f}_{{pp}}\sim 0.3$ for SBGs and ${f}_{{pp}}\sim 0.03$ for SFGs, which compensates for the lower star formation rates of Pop-III stars. (Note that the assumed energy injection might be optimistic because of the uncertainty in f_pp and ${\xi }_{{\rm{SBG}}}$.)