Point-source and diffuse high-energy neutrino emission from Type IIn supernovae

TL;DR

This paper investigates the potential of Type IIn supernovae to produce high-energy neutrinos, estimating their contribution to IceCube's observed flux and exploring their detectability as point sources.

Contribution

It provides the first detailed Monte Carlo-based estimate of the diffuse neutrino emission from Type IIn supernovae and assesses their role in the astrophysical neutrino flux.

Findings

Type IIn SNe contribute about 10% of IceCube's neutrino flux above 60 TeV.

Detection of individual Type IIn SNe as neutrino sources is feasible within 20 Mpc.

A small fraction of core collapse SNe being Type IIn could dominate the diffuse neutrino background.

Abstract

Type IIn supernovae (SNe), a rare subclass of core collapse SNe, explode in dense circumstellar media that have been modified by the SNe progenitors at their last evolutionary stages. The interaction of the freely expanding SN ejecta with the circumstellar medium gives rise to a shock wave propagating in the dense SN environment, which may accelerate protons to multi-PeV energies. Inelastic proton-proton collisions between the shock-accelerated protons and those of the circumstellar medium lead to multi-messenger signatures. Here, we evaluate the possible neutrino signal of type IIn SNe and compare with IceCube observations. We employ a Monte Carlo method for the calculation of the diffuse neutrino emission from the SN IIn class to account for the spread in their properties. The cumulative neutrino emission is found to be ~ 10 per cent of the observed IceCube neutrino flux above 60 TeV. Type IIn SNe would be the dominant component of the diffuse astrophysical flux, only if 4 per cent of all core collapse SNe were of this type and 20 to 30 per cent of the shock energy was channeled to accelerated protons. Lower values of the acceleration efficiency are accessible by the observation of a single type IIn SN as a neutrino point source with IceCube using up-going muon neutrinos. Such an identification is possible in the first year following the SN shock breakout for sources within 20 Mpc.

Tables

Table 1. Table 1: Parameters used for the two SN scenarios S1 and S2.

Parameter	S1	S2
$r_w$ [cm]	$3.2\times {10}^{16}$	$3.2\times {10}^{17}$
$M_{\mathrm{ej}}$ [M☉]	10
$M_{\mathrm{csm}}$ [M☉]	10
$v_{\mathrm{sh}}$ [km s-1]	$9.5\times {10}^3$
${ϵ}_B$	$3\times {10}^{-4}$
${ϵ}_p$	$3\times {10}^{-2}$
$r_i$ [cm]	$4.4\times {10}^{14}$	$4.5\times {10}^{13}$
$r_o$ [cm]	$3.2\times {10}^{16}$	$3.2\times {10}^{17}$
$n_i$ [cm-3]	${10}^{11}$	${10}^{12}$
$K_w$ [g cm-1]	$4.5\times {10}^{16}$	$4.5\times {10}^{15}$
${\dot{M}}_w$ [$M_{☉}$ yr-1]†	0.1	0.01
$t_i$ [d]	$5.4$	$0.54$
Duration [yr]	$1$	$10.7$

Table 2. Table 2: Expected IceCube ν μ + ν ¯ μ subscript 𝜈 𝜇 subscript ¯ 𝜈 𝜇 \nu_{\mu}+\bar{\nu}_{\mu} events at different energy bands from a SN IIn located at 10 Mpc compared to the number of background events. For the S1 and S2 parameters, see text. The source events for S1 and S2 are calculated by integrating the rates shown in Fig.  8 over a period of 385 days and 10.7 yr, respectively. The number of background events is calculated for a period of 385 days over a 0.75°  circle around the source.

$E_{\nu }$ (TeV)	SN IIn			Background
	S1	S2
		First year	Full duration
0.1-1	0.082	0.003	0.007	4.113
1-10	0.696	0.026	0.057	2.043
10-${10}^2$	1.080	0.039	0.089	0.143
${10}^2$-${10}^3$	0.304	0.010	0.024	0.002
0.1-${10}^3$	2.163	0.079	0.177	6.301

Table 3. Table 3: List of SNe IIn detected within 40 Mpc over a period of eight years (2008-2015).

Name	$T_{\max }$ (MJD)†	Host Galaxy	R.A. (J2000)	Dec. (J2000)	z	$D_L$ (Mpc)	Type	References
SN2008S	54508.5	NGC 6946	20:34:45.4	+60:05:57.8	0.0002	0.7	IIn-pec/LBV	(1)
SN2008X	54502.5	NGC 4141	12:09:48.33	+58:51:01.6	0.006	28.0	IIn	(2)
SN2009ip‡	56206.5	NGC 7259	22:23:08.3	-28:56:52.4	0.006	26.4	IIn	(3)
SN2009kr	55148.5	NGC 1832	05:12:03.3	-15:41:52.2	0.007	16.0	II	(4)
SN2011A	55562.5	NGC 4902	13:01:01.2	-14:31:34.8	0.009	39.7	IIn	(5)
SN2011fh	55797.5	NGC 4806	12:56:14.0	-29:29:54.8	0.008	36.0	IIn	(6)
SN2011ht	55879.5	UGC 5460	10:08:10.58	+51:50:57.1	0.004	16.0	IIn	(7)
SN2013gc	56603.5	ESO 430-G20	08:07:11.9	-28:03:26.3	0.003	15.1	IIn	(8)
PSN J14041297-0938168	56645.5	IC 4363	14:04:13.0	-09:38:16.8	0.003	12.5	IIn	(9)
CSS140111:060437-123740	56734.5	MCG-02-16-02	06:04:36.71	-12:37:40.6	0.007	32.9	IIn	(10)
SN2014G	56674.5	NGC 3448	10:54:34.13	+54:17:56.9	0.005	20.0	II L	(11)
SN2015bh	57163.5	NGC 2770	09:09:34.96	+33:07:20.4	0.006	28.5	Impostor	(12)
SN2015J	57201.5	A073505-6907	07:35:05.2	-69:07:53.1	0.005	24.0	IIn	(13)
PSN J13522411+3941286	57071.5	NGC5337	13:52:24.1	+39:41:28.6	0.007	32.1	IIn	(14)
ASASSN-15lf	57194.5	NGC 4108	12:06:45.56	+67:09:24.00	0.008	37.3	IIn	(15)
SNhunt248	56828.5	NGC 5806	14:59:59.5	+01:54:26.2	0.005	20.1	IIn Pec	(16)

Equations

$$n(r)=n_{\mathrm{i}}{\left(\frac{r_{\mathrm{i}}}{r}\right)}^2=\frac{K_{\mathrm{w}}}{m{r}^2},$$

$$K_{\mathrm{w}}\equiv \frac{{\dot{M}}_{\mathrm{w}}}{4\pi v_{\mathrm{w}}}$$

$$M_{\mathrm{csm}}=4\pi m{\int }_{r_{\ast }}^{r_{\mathrm{w}}}\mathrm{d}r{r}^{2}n(r)\simeq 4\pi n_{\mathrm{i}}m{r}_{\mathrm{i}}^2{r}_{\mathrm{w}},$$

$$r_{\mathrm{i}}=r_{\mathrm{w}}{\left(1+\frac{4\pi mc{r}_{\mathrm{w}}^2}{{\sigma }_{\mathrm{T}}{v}_{\mathrm{sh}}{M}_{\mathrm{csm}}}\right)}^{-1}.$$

$$r_{\mathrm{dec}}=r_{\mathrm{i}}+\frac{M_{\mathrm{ej}}}{4\pi m{n}_{\mathrm{i}}{r}_{\mathrm{i}}^2}$$

$$\frac{\partial N_{\mathrm{p}}(\gamma ,r)}{\partial r}+\frac{N_{\mathrm{p}}(\gamma ,r)}{v_{\mathrm{sh}}{t}_{\mathrm{pp}}(r)}-\frac{\partial }{\partial \gamma }\left[\frac{\gamma }{r}{N}_{\mathrm{p}}(\gamma ,r)\right]$$

$$\frac{\mathrm{d}{N}_{{\nu }_{\mathrm{i}}+{\bar{\nu }}_{\mathrm{i}}}(x_{\nu },r)}{\mathrm{d}r}+\frac{N_{{\nu }_{\mathrm{i}}+{\bar{\nu }}_{\mathrm{i}}}(x_{\nu },r)}{v_{\mathrm{sh}}{t}_{\mathrm{esc}}(r)}$$

$$Q_{\mathrm{p}}(\gamma )\equiv \frac{{\mathrm{d}}^2{N}_{\mathrm{p}}}{\mathrm{d}\gamma \mathrm{d}r}=\frac{9\pi }{8f_{\mathrm{p}}}{ϵ}_{\mathrm{p}}{r}_{\mathrm{i}}^2{n}_{\mathrm{i}}{\left(\frac{v_{\mathrm{sh}}}{c}\right)}^2{\gamma }^{-p}H(\gamma -{\gamma }_{\mathrm{p},\min })H({\gamma }_{\mathrm{p},\max }-\gamma ),$$

$$\displaystyle E_{\nu,\max}=\left\{\begin{array}[]{cc}\frac{eBv_{\rm sh}^{2}}{96\sigma_{\rm pp}nc^{2}}\propto v_{\rm sh}^{3}\epsilon_{\rm B}^{1/2}K_{\rm w}^{-1/2}r&t_{\rm pp}<t_{\rm ad}\\ \\ \frac{eB\kappa_{\rm pp}v_{\rm sh}r}{24c}\propto v_{\rm sh}^{2}\epsilon_{\rm B}^{1/2}K_{\rm w}^{1/2}&t_{\rm pp}>t_{\rm ad}\end{array}\right.$$

$$Q_{{\nu }_{\mu }+{\bar{\nu }}_{\mu }}(x_{\nu },r)$$

$$Q_{{\nu }_{\mathrm{e}}+{\bar{\nu }}_{\mathrm{e}}}(x_{\nu },r)$$

$$E_{\nu }{F}_{\nu }(E_{\nu },r)=\sum _{i=e,\mu }\frac{m_{\mathrm{e}}{c}^2{x}_{\nu }^2{N}_{{\nu }_{\mathrm{i}}+{\bar{\nu }}_{\mathrm{i}}}(x_{\nu },r)}{4\pi D_{\mathrm{L}}^2{t}_{\mathrm{esc}}},$$

$$E_{\nu }{\varphi }_{\nu }(E_{\nu })=(1+z){\int }_{r_{\mathrm{i}}}^{r_{\mathrm{o}}}\mathrm{d}r\frac{E_{\nu }{F}_{\nu }(E_{\nu },r)}{v_{\mathrm{sh}}}.$$

$$E_{\nu }{\Phi }_{\nu }(E_{\nu })$$

$$E_{\nu }{\Phi }_{\nu }(E_{\nu })=\frac{N_{\mathrm{tot}}}{4\pi N_{\mathrm{sim}}}\sum _{j=1}^{10}\sum _{i=1}^{10^4}{E}_{\nu }{\varphi }_{\nu }^{(i,j)}(E_{\nu }),$$

$$F_{\nu }(E_{\nu })\propto {ϵ}_{\mathrm{p}}{v}_{\mathrm{sh}}{K}_{\mathrm{w}}^2{t}^{-1}.$$

$$F_{\nu }(E_{\nu })\propto {ϵ}_{\mathrm{p}}{v}_{\mathrm{sh}}^3{K}_{\mathrm{w}}{t}^0.$$

$$t\lesssim 100\mathrm{d}{K}_{\mathrm{w},16}{\beta }_{\mathrm{sh},-1.5}^{-2},$$

$$N_{\mathrm{IIn}}$$

$$\frac{{K_w\times {ϵ}_{\mathrm{p}}|}_{2011fh}}{{K_w\times {ϵ}_{\mathrm{p}}|}_{S1}}$$

Full text

††pagerange: Point-source and diffuse high-energy neutrino emission from Type IIn supernovae–LABEL:lastpage††pubyear: 2016

1 Introduction

The identification of high-energy ($E_{\nu}>25$ TeV) neutrino sources would provide direct evidence for the acceleration of cosmic rays (CR) up to $\sim$PeV energies. IceCube recently observed astrophysical neutrinos in both starting event and up-going muon neutrino samples (Aartsen et al., 2013a, b, 2016d), but no significant anisotropies are identified in the arrival directions of neutrinos. Additional searches for the origin of these neutrinos have not yet revealed any specific sources (Adrian-Martinez et al., 2014; Aartsen et al., 2017b). No prompt emission of neutrinos was identified coincident with gamma-ray bursts (Aartsen et al., 2016b) and no more than $27\%$ of the sub-PeV neutrino signal may originate from blazars (Aartsen et al., 2017a), a type of radio-loud active galactic nuclei with relativistic jets pointing towards the observer (Blandford & Rees, 1978; Urry & Padovani, 1995). A dominant contribution of blazars beyond PeV energies can now also be excluded (see e.g. Murase & Waxman, 2016), although a $10\%-20\%$ contribution is still viable (Padovani et al., 2015; Aartsen et al., 2016a; Padovani et al., 2016). Star-forming galaxies (i.e., galaxies with vigorous star formation and high gas densities in their central regions) cannot contribute more than $\sim 30\%$ to the diffuse neutrino background between 25 TeV and 2.8 PeV (Bechtol et al., 2015), if recent constraints from the non-blazar extragalactic $\gamma$-ray background (EGB) (Ackermann et al., 2016) are taken into account and the 30 TeV excess of IceCube events is attributed to star-forming galaxies (Aartsen et al., 2015b). Meanwhile, most of the scenarios predicting a dominant Galactic contribution to the high-energy IceCube signal are disfavored (e.g. Ahlers et al., 2016).

There is convincing evidence that CR with energies up to the knee of the CR spectrum are accelerated at Galactic supernova (SN) remnants (for reviews, see Bell, 2013; Blasi, 2013). Acceleration beyond PeV energies may be possible at shocks of interaction-powered SNe, i.e., SNe exploding in dense circumstellar media (CSM) (Katz et al., 2011; Murase et al., 2011; Cardillo et al., 2015; Zirakashvili & Ptuskin, 2016). Due to the presence of (multi-)PeV protons in dense environments, inelastic proton-proton (p-p) collisions with the non-relativistic protons of the shocked CSM may lead to interesting multi-messenger signatures, such as GeV $\gamma$-ray emission, high-energy ($>100$ TeV) neutrino production, and radio emission (Murase et al., 2014; Petropoulou et al., 2016; Zirakashvili & Ptuskin, 2016). In contrast to $\gamma$-rays, which may be attenuated via photon-photon absorption soon after the shock breakout, when their production rate is higher (e.g. Kantzas et al., 2016), neutrinos escape the source unimpeded. In principle, neutrino detection from an individual interaction-powered SN would serve as the smoking gun for CR shock acceleration to PeV energies.

Signs of interaction between the SN ejecta and the CSM are observable on the early-time light curves and spectra (hours to days after the explosion) of type IIn supernovae (SNe IIn), a rare subclass of core collapse (CC) SNe (Schlegel, 1990; Filippenko, 1997). The high CSM densities needed to explain their observational properties can be accounted for, if the SN IIn progenitor has undergone strong mass loss before its explosion (for a review, see Smith, 2014). The inferred mass-loss rates are typically higher than $10^{-3}$ M*☉* yr*-1* (e.g. Salamanca et al., 1998; Chugai et al., 2004; Kiewe et al., 2012; Chandra et al., 2015), but they may be as low as $10^{-5}$ M*☉* yr*-1* (Crowther, 2007). SNe IIn progenitors exhibit wide diversity not only in their estimated mass-loss rates but also in their wind velocities. These typically lie in the range $v_{\rm w}\approx 10-10^{3}$ km s*-1* (e.g. Salamanca et al., 1998; Chevalier & Li, 2000). Giant progenitors have slower and denser winds compared to those those from more compact progenitors which have faster and more tenuous winds.

SNe IIn pose an interesting alternative to existing scenarios for neutrino production, as it has been discussed first by Murase et al. (2011) and more recently by Zirakashvili & Ptuskin (2016) – henceforth ZP16. Murase et al. (2011) demonstrated that multi-TeV neutrinos are detectable by a generic IceCube-like detector for SNe at $\lesssim 20-30$ Mpc, if the cosmic ray acceleration efficiency is 10 per cent. ZP16 calculated the diffuse neutrino emission from the SN IIn population, by solving the hydrodynamical equations for the evolution of the SN shock and taking into account particle acceleration and the CR feedback on the shock structure. The method presented in ZP16 is better suited for a single source, as particle acceleration at both SN shocks (forward and reverse) is treated in more detail, but it is impractical when applied to many sources. The diffuse neutrino emission was, therefore, calculated by adopting the same physical parameters for all SNe (e.g., CSM density, shock velocity, and others). In particular, the CR acceleration efficiency was fixed to be 50 per cent for all SNe shocks, which might be unrealistically high (Caprioli & Spitkovsky, 2014; Park et al., 2015).

In this study, we calculate the neutrino signal from SNe IIn and discuss the possibility of constraining the CR accelerated energy fraction by means of diffuse and point-source neutrino observations with IceCube. We employ a Monte Carlo method for the calculation of the diffuse neutrino emission from the SN IIn class, in an attempt to incorporate the wide spread of their properties into their cumulative emission. The parameter values, which are assigned to the simulated sources, are randomly drawn from distributions that are motivated by observations. For each simulated source, we solve the evolutionary equations for the proton and neutrino distributions, under the assumption that protons are accelerated at the SN forward shock and produce neutrinos via p-p collisions with the non-relativistic protons of the shocked CSM. For the calculation of the diffuse neutrino emission, we adopt the redshift evolution of CC SNe as presented in Hopkins & Beacom (2006).

Interestingly, a SN IIn was recently discovered in a search with the Palomar Transient Factory (PTF) (Law et al., 2009) following an IceCube neutrino doublet trigger alert. However, there was no evidence for a physical connection, as the detection was most likely coincidental (Aartsen et al., 2015c). Motivated by this observation, we complement our analysis by investigating the possibility of detecting an individual SN IIn as a neutrino point source with through-going muons detected with IceCube in the TeV range over a period of seven years (Aartsen et al., 2017b). By comparing the arrival times and directions of IceCube neutrinos with known SNe IIn, we find one starting neutrino event in close spatial and temporal correlation with a close-by SN IIn in the Southern Sky. This finding motivates follow-up searches with up-going neutrino data from the neutrino telescope ANTARES (Adrian-Martinez et al., 2014).

This paper is structured as follows. In Section 2, we describe the theoretical framework and our methods. In Section 3, we compute the diffuse neutrino emission from the SNe IIn class and compare it with IceCube neutrino observations. We calculate the neutrino signal expected from individual sources and discuss the possibility of neutrino detection with IceCube. In Section 4, we discuss our results and model caveats. We also discuss the possible association between SN2011fh and an IceCube event of the starting sample. We finally conclude in Section 5. Here, we adopt a cosmology with $\Omega_{\rm K}=0$, $\Omega_{\rm M}=0.31$, $\Omega_{\Lambda}=0.69$ and $H_{0}=69.6$ km s*-1* Mpc*-1*.

2 Model and Methods

The CSM is modelled as an extended shell with mass density $\propto r^{-2}$ (e.g. Chevalier, 1982) and outer radius $r_{\rm w}=v_{\rm w}t_{\rm w}$, where $v_{\rm w}$ is the expanding velocity of the material that has been ejected from the progenitor star over a period $t_{\rm w}$ (for a recent review, see Smith, 2014). The interaction of the freely expanding SN ejecta with the CSM gives rise to a fast shock wave propagating in the CSM (forward shock) with velocity $v_{\rm sh}$ and a reverse shock that crosses the outer parts of the SN ejecta.111We neglect the contribution of the reverse shock to the neutrino emission, but we discuss our choice in Section 4. As long as the interaction between the SN ejecta and the CSM takes place within a region that is optically thick to Thomson scattering, $\tau_{\rm T}\gg 1$, the SN shock is mediated by radiation, thus prohibiting particle acceleration (e.g. Murase et al., 2011; Katz et al., 2011). The radiation may escape when $\tau_{\rm T}\sim c/v_{\rm sh}$ (Weaver, 1976); this defines the so-called shock breakout time $t_{\rm i}$ and radius $r_{\rm i}\sim v_{\rm sh}t_{\rm i}$ (e.g. Ofek et al., 2014a). For dense CSM environments, as those considered in this work, the shock is expected to break out in the wind, namely $r_{\rm i}\gg r_{*}$, where $r_{*}$ is the typical radius of the stellar envelope.

Particle (electron and ion) acceleration can, in principle, take place at $r\geq r_{\rm i}$ after the SN shock becomes collisionless. Henceforth, we use $r_{\rm i}$ as the normalization radius. The CSM density profile may be written as

[TABLE]

where

[TABLE]

is the wind mass loading parameter, $\dot{M}_{\rm w}$ is the mass-loss rate of the progenitor star, and $m=1.4\,m_{\rm H}$ for a medium with 10 per cent He abundance by number.

The total CSM mass can be then estimated as

[TABLE]

where $r_{\rm w}\gg r_{*}$ was assumed. Combining equation (3) with the condition $\tau_{\rm T}(r_{\rm i})=c/v_{\rm sh}$, where the Thomson optical depth is defined as $\tau_{\rm T}(r)=\int_{r}^{r_{\rm w}}{\rm d}y\,\sigma_{\rm T}n(y)$, we obtain $r_{\rm i}$:

[TABLE]

The CSM density at the shock breakout radius, $n_{\rm i}$, is then derived from equations (3) and (4) knowing $r_{\rm w}$ and $M_{\rm csm}$. Its minimum value, which is obtained when the shock breakout radius becomes maximum for fixed CSM mass and shock velocity, is given by $\left(4\pi^{1/3}m^{1/3}c/\sigma_{\rm T}v_{\rm sh}M_{\rm csm}^{1/3}\right)^{3/2}\simeq 9\times 10^{9}\,{\rm cm}^{-3}\left(M_{\rm csm}/10M_{\sun}\right)^{-1/2}\left(v_{\rm sh}/0.1c\right)^{-3/2}$. We terminate our calculations at a maximum radius beyond which the contribution to the total neutrino fluence is not important. This is set by the deceleration radius

[TABLE]

or by the extent of the CSM, i.e., $r_{\rm o}=\min(r_{\rm dec},r_{\rm w})$. The production rates of accelerated protons and neutrinos (at a fixed energy) are expected to decrease beyond the deceleration radius, since they scale, respectively, as $v_{\rm sh}^{3}$ and $v_{\rm sh}^{2}$ (e.g. Petropoulou et al., 2016). Furthermore, the neutrino production rate is expected to decrease significantly beyond $r_{\rm w}$, where the density of the medium is much lower than $n(r_{\rm w})$. The condition $r_{\rm o}=r_{\rm dec}$ suggests that the ejecta mass is smaller than the total CSM mass, whereas $r_{\rm o}=r_{\rm w}$ implies that $M_{\rm ej}>M_{\rm csm}$.

To determine the temporal evolution of the neutrino emission from a single source we calculate:

the temporal evolution of the proton distribution from the shock breakout $r_{\rm i}$ until the outer radius $r_{\rm o}$. 2. 2.

the maximum proton energy at each radius by comparing the local acceleration and loss timescales ($t_{\rm acc}$ and $t_{\rm loss}$). 3. 3.

the neutrino fluence $\phi_{\nu}$ (in erg cm*-2*) by integrating the neutrino flux from the shock breakout till the time the shock decelerates or reaches the outer radius of the dense CSM.

We solve the evolutionary equation for relativistic protons in the presence of losses and injection. The latter is provided by Fermi acceleration at the shock (e.g. Mastichiadis, 1996; Kirk et al., 1998). Neutrinos are produced at a rate dictated by the relativistic proton distribution and the non-relativistic proton density of the shocked CSM, which is assumed to be uniform and equal to $4n(r)$ (for details, see Petropoulou et al., 2016), and they escape from the shell without any losses. The equations for the proton and neutrino distributions may be written as:

[TABLE]

where $N_{\rm p}(\gamma,r)$ is the total number of protons in the shell with radius $r$ with Lorentz factors between $\gamma$ and $\gamma+{\rm d}\gamma$, $t_{\rm pp}=\left(4\kappa_{\rm pp}\sigma_{\rm pp}cn(r)\right)^{-1}$ is the loss timescale for p-p collisions with inelasticity $\kappa_{\rm pp}\approx 0.5$ and cross-section $\sigma_{\rm pp}\simeq 3\times 10^{-26}$ cm*-2*, $N_{\nu_{\rm i}+\bar{\nu}_{\rm i}}(x_{\nu},r)$ is the total number of neutrinos and anti-neutrinos222For simplicity, we refer to the sum $\nu_{\rm i}+\bar{\nu}_{\rm i}$ as neutrinos. of flavour i=e, $\mu$ with energies (in $m_{\rm e}c^{2}$ units) between $x_{\nu}$ and $x_{\nu}+{\rm d}x_{\nu}$ at a shock radius $r$, $t_{\rm esc}=h/c$, and $h\approx r/4$ is the width of the shocked gas shell. Inelastic p-p collisions can be treated as catastrophic energy losses (e.g. Sturner et al., 1997; Schlickeiser, 2002) in contrast to the energy losses due to adiabatic expansion of the shell (third term in the left hand side of equation 6). Other energy loss processes for protons, such as photopion cooling, are irrelevant (see also Murase et al., 2011). The escape timescale of protons from the shell is assumed to be much longer than all other typical timescales of the system. Any decrease of the relativistic proton number due to the physical escape of the system would lead to lower neutrino fluxes than those presented in the following sections.

For a wind density profile, as adopted here, the proton injection rate is independent of radius and is given by

[TABLE]

where $\epsilon_{\rm p}$ is the fraction of the shock kinetic energy channeled to accelerated protons, $f_{\rm p}=\ln\left(\gamma_{\rm p,\max}/\gamma_{\rm p,\min}\right)$ for power-law index $p=2$ or $f_{\rm p}=(p-2)^{-1}$ for $p>2$, $\gamma_{\rm p,\min}=1$, and $\gamma_{\rm p,\max}$ is the maximum Lorentz factor of the accelerated protons. This is determined by $t_{\rm acc}=\min[t_{\rm pp},t_{\rm ad}]$, where $t_{\rm acc}\sim 6m_{\rm p}\gamma c^{3}/eBv_{\rm sh}^{2}$ assuming Bohm acceleration. Here, $B=\sqrt{9\pi\epsilon_{\rm B}mv_{\rm sh}^{2}n}$, $\epsilon_{\rm B}$ is the fraction of the post-shock magnetic energy, and $t_{\rm ad}\sim r/v_{\rm sh}$ corresponds to the source lifetime (see equations (21) and (29) in Petropoulou et al. 2016). At small radii, where the CSM density is higher, $t_{\rm pp}<t_{\rm ad}$ leading to $\gamma_{\rm p,\max}\propto r$. At larger radii, adiabatic losses dominate and $\gamma_{\rm p,\max}$ becomes independent of radius. The respective maximum neutrino energy is written as:

[TABLE]

The neutrino production rate $Q_{\nu_{\rm i}+\bar{\nu}_{\rm i}}$ for parent proton energies $E_{\rm p}>0.1$ TeV is written as (Kelner et al., 2006):

[TABLE]

where $E_{\rm p}=\gamma m_{\rm p}c^{2}$, $x=x_{\nu}m_{\rm e}c^{2}/E_{\rm p}$, $F_{\nu_{\mu}}^{(1)}(x,E_{\rm p})$ and $F_{\nu_{\mu}}^{(2)}(x,E_{\rm p})\approx F_{\nu_{\rm e}}$ are respectively given by equations (66) and (62) in Kelner et al. (2006). For $E_{\rm p}\lesssim 0.1$ TeV, we adopt the $\delta$-function approximation for the pion production rate as described in Kelner et al. (2006). Equations (13) and (14) result in neutrino energy spectra that can be well described by a power law with index $\sim p$ for $E_{\nu}<E_{\nu,\max}$ and an exponential cutoff at $\sim E_{\nu,\max}$ (see also Fig. 12 in Kelner et al. (2006)).

The all-flavour neutrino energy flux from a source located at a luminosity distance $D_{\rm L}$ is given by:

[TABLE]

where $F_{\nu}$ is the differential neutrino energy flux (i.e., $F_{\nu}\equiv E_{\nu}{\rm d}N_{\nu}/{\rm d}E_{\nu}{\rm d}t$) at the shock radius $r$ and $E_{\nu}=x_{\nu}m_{\rm e}c^{2}/(1+z)$ is the observed neutrino energy. The observed neutrino fluence of a single source is

[TABLE]

The diffuse all-flavour neutrino flux from the SNe IIn class can be estimated as (e.g. Cholis & Hooper, 2013):

[TABLE]

where ${\rm d}V_{\rm c}$ is the comoving volume element (e.g. Hogg, 1999), $z_{\max}=6$, $\dot{n}_{\rm IIn}=\xi\,\dot{n}_{\rm SNII}$ is the volumetric rate of SNe IIn, $\dot{n}_{\rm SNII}$ is the volumetric rate of all CC SNe, and $\xi$ is the fraction of CC SNe that are of the IIn type.

2.1 Physical parameters

The physical parameters of each simulated source which are randomly selected are summarized below:

shock velocity, 2. 2.

ejecta mass, 3. 3.

duration $t_{\rm w}$ of the pre-explosion period of mass loss, 4. 4.

velocity $v_{\rm w}$ of the wind, 5. 5.

total CSM mass, 6. 6.

fraction of the post-shock magnetic energy $\epsilon_{\rm B}$, and 7. 7.

proton accelerated energy fraction $\epsilon_{\rm p}$.

All other physical parameters of the system can be derived from a combination of the above. For example, the extent of the CSM $r_{\rm w}$ is determined by (iii) and (iv). The shock breakout radius $r_{\rm i}$ can be then estimated using equation (4), if $v_{\rm sh}$, $r_{\rm w}$, and $M_{\rm csm}$ are defined. The last two parameters incorporate the details of particle acceleration and magnetic field amplification at the shock, while parameters (i)-(v) are related to the last stages of stellar evolution leading to the supernova explosion. For all the physical parameters mentioned above, we adopt uniform distributions of random variables, as detailed in the following paragraphs. Given that the observations of SNe IIn are still not sufficient for pinpointing the distribution of their physical parameters, our choice ensures that our results will not be biased towards a more probable parameter set for the simulated sources.

Kinetic simulations of particle acceleration at non-relativistic quasi-parallel shocks (i.e., the angle between the shock normal and the magnetic field direction is $\lesssim 45$°) have shown that $\epsilon_{\rm p}=0.05-0.15$, while $\epsilon_{\rm p}\rightarrow 0$ for quasi-perpendicular shocks (e.g. Caprioli & Spitkovsky, 2014). As the proton accelerated energy fraction depends on the pre-existing magnetic field in the unshocked CSM, it is appropriate to assume a range of values in our simulations. Adopting the maximum $\epsilon_{\rm p}$ predicted by theory would imply that the unshocked CSM field in SNe IIn is weak or radial (e.g. Sironi et al., 2013), which is questionable. In our simulations we set $\log\epsilon_{\rm p}=-1-\tilde{r}$, where $\tilde{r}$ denote a uniformly distributed random number in the range (0,1).

The fraction of the post-shock magnetic energy is usually inferred from GHz radio observations of interaction-powered SNe and is typically found to be $10^{-3}-10^{-1}$ (Chevalier, 1998; Weiler et al., 2002; Chandra et al., 2009; Chandra et al., 2015; Kamble et al., 2016). Some of these estimates are obtained assuming that the observed radio emission is produced by shock-accelerated electrons and without taking into account electron cooling. The first assumption may not be valid for the dense environments of SNe IIn, where relativistic electrons produced via the decay of charged pions from p-p collisions, are expected to contribute to the observed radio emission (Murase et al., 2014; Petropoulou et al., 2016). Furthermore, electron cooling cannot be, in general, neglected for the inferred $\epsilon_{\rm B}$ values and the high CSM densities of interaction-powered SNe (see e.g. Martí-Vidal et al., 2011b; Kamble et al., 2016; Petropoulou et al., 2016). It could be therefore possible that $\epsilon_{\rm B}\ll 10^{-3}$, as is usually the case for gamma-ray burst afterglows (e.g. Barniol Duran, 2014). Because of the aforementioned uncertainties, we choose a wide range of $\epsilon_{\rm B}$ values, namely $\log\epsilon_{\rm B}=-5\tilde{r}-1$.

The shock velocity is one of the physical parameters that can be inferred from optical spectroscopy and/or photometry. Ofek et al. (2014b) provided lower limits for the shock velocity at the breakout by fitting the optical light curves of 15 SNe IIn observed with PTF/iPTF (Law et al., 2009). We created a sample of 23 sources by adding the shock velocity estimates of eight more sources that are available in the literature (1986J, Bietenholz et al. 2010; 1994aj, Salamanca et al. 1998; 1994W, Chugai et al. 2004; 1995G, Chugai & Danziger 2003; 1995N, Chandra et al. 2009; 1997eg, Salamanca et al. 2002; 1997ab, Salamanca et al. 1998; 2006jd, Chandra et al. 2012; 2010jl, Chandra et al. 2015). The distribution of observed shock velocities (in logarithmic space) is compatible with a uniform distribution with a median of 4300 km s*-1*, although a Gaussian distribution cannot be excluded. Moreover, most of the values in the observed sample are lower limits (e.g. Ofek et al., 2014b). For these reasons, we have adopted a uniform distribution with a median of 9500 km s*-1* (i.e., $\log(v_{\rm sh}/c)=-1-\tilde{r}$) in our simulations.

All other parameters are less well constrained observationally. The ejecta mass typically ranges from 2-5 M*☉* but it may also be as high as 15 M*☉* depending on the progenitor type and stellar mass at the zero age main sequence (ZAMS) (e.g. Shussman et al., 2016). Here, we adopt a uniform distribution for the ejecta mass between 5 M*☉* and 15 M*☉, noting that $M_{\rm ej}$ does not have a strong effect on the neutrino emission, as it affects only the deceleration radius. There is even larger uncertainty on the distribution of the CSM masses (e.g. Smith, 2014; Margutti et al., 2017). To account for the large spread in the estimated values, in our simulations we assume a uniform distribution (in linear space) between 0.01 M☉* and 20 M*☉*.

The inferred wind velocities also display a diversity, i.e., $v_{\rm w}\simeq 10-10^{3}$ km s*-1* (e.g. Salamanca et al., 1998; Chevalier & Li, 2000). Finally, $t_{\rm w}$ is the least certain parameter, as it can be measured only for isolated cases, namely $\eta$ Carinae, P Cygni and SN 2009ip (e.g. Mauerhan et al., 2013a; Smith, 2014). The duration of the mass-loss events prior to explosion may typically last from a few years up to several decades. Dynamical ages of shells of matter around luminous blue variable stars (LBV) – likely progenitors of SNe IIn – range from hundreds years to several thousand years (see Smith, 2014, and references therein). These are, however, only indicators of the mass-loss duration during the LBV phase. The uncertainty on $t_{\rm w}$ is even larger, as this may depend on the progenitor type and its initial mass. In our simulations we, therefore, assume a wide range of values $\log t_{\rm w}=1+2\tilde{r}$ (in yr), while noting that the diffuse neutrino flux is not sensitive on the maximum adopted duration. The $t_{\rm w}$-distribution in combination with the uniformly distributed $v_{\rm w}$, results in $r_{\rm w}$ values that are distributed around $3\times 10^{16}$ cm.

Figure 1 shows $r_{\rm dec}$ (black points) and $r_{\rm w}$ (red points) versus $r_{\rm i}$ for the simulated sources. No points lie below the solid line, as expected, since this region corresponds to $r_{\rm i}>r_{\rm o}$. The boomerang-like shape of the distribution of points in the middle panel can be understood after inspection of equations (4) and (5). The parameter values for the two scenarios discussed in Section 3.2 are indicated with blue symbols. The mass loss rate, defined as $\dot{M}_{\rm w}=4\pi n_{\rm i}mr_{\rm i}^{2}v_{\rm w}$, is plotted against the CSM density at the shock breakout $n_{\rm i}$ in Fig. 2. Blue symbols have the same meaning as in Fig. 1. The median of density and mass loss rate distributions are respectively $2\times 10^{11}$ cm*-3* and $0.08$ M*☉* yr*-1*. The wind mass-loading parameter, $K_{\rm w}=\dot{M}_{\rm w}/4\pi v_{\rm w}$, is plotted against $t_{\rm i}=r_{\rm i}/v_{\rm sh}$ in Fig. 3 together with estimates of the breakout time and the mass-loading parameter adopted by Ofek et al. (2014b) (red symbols). The two scenarios discussed later in Section 3.2 are also shown as blue symbols. The golf club shape that appears in the $K_{\rm w}-t_{\rm i}$ plot results from the definition of $r_{\rm i}$ and $n_{\rm i}$. For small shock breakout radii, equation (4) reads $r_{\rm i}\approx\sigma_{\rm T}v_{\rm sh}M_{\rm csm}/(4\pi mcr_{\rm w})$ and $n_{\rm i}\approx c/\sigma_{\rm T}v_{\rm sh}r_{\rm i}$. In this limit, $K_{\rm w}\approx 4\pi mct_{\rm i}/\sigma_{\rm T}$ (golf club’s shaft). At larger shock breakout radii where $r_{\rm i}\approx r_{\rm w}$, $n_{\rm i}\approx M_{\rm csm}/4\pi mr_{\rm i}^{3}$ and $K_{\rm w}\approx M_{\rm csm}/v_{\rm sh}t_{\rm i}$ (golf club’s head) with the dispersion arising from the random $M_{\rm csm}$ and $v_{\rm sh}$ values.

2.2 Monte Carlo simulations

We created $10^{4}$ parameter sets as described in the previous section and calculated the respective neutrino emission. A total of $10^{5}$ random redshifts ($0\leq z\leq z_{\max}$) were generated according to the distribution ${\rm d}V_{\rm c}\,\dot{n}_{\rm IIn}(z)\,(1+z)^{-1}$. Henceforth, we adopt the CC SNe volumetric rate of Hopkins & Beacom (2006) and $\xi=0.04$ (see e.g. Table 5 in Cappellaro et al., 2015). Each of the simulated sources was placed at 10 different redshifts and the diffuse neutrino flux was calculated as (see also equation (17)):

[TABLE]

where $N_{\rm tot}=\left(10^{-7}\xi/3.1\right)\int_{0}^{z_{\max}}{\rm d}z\left({\rm d}V_{\rm c}/{\rm d}z\right)\,\dot{n}_{\rm IIn}(z)(1+z)^{-1}$ is the number of SNe IIn per sec, $N_{\rm sim}=10^{5}$ is the number of simulated neutrino fluence spectra, and the indices i,j run over the different parameter sets and redshifts, respectively. Because of the dispersion intrinsic to the simulation process, we repeated the above procedure for 100 different sets of redshifts.

3 Results

3.1 Diffuse neutrino emission

Figure 4 presents our predictions for the diffuse neutrino flux (per flavour) from SNe IIn in comparison with the astrophysical neutrino flux detected by IceCube (The IceCube Collaboration et al., 2015; Aartsen et al., 2016d). The blue solid lines enclose the 68 per cent of the simulated spectra when the accelerated proton energy fraction varies among sources in the range $[0.01,0.1]$. All other curves are obtained under the assumption of a constant acceleration efficiency among sources.

The energy spectrum of the diffuse neutrino emission is best described by a power law with index $\sim 2$ and exponential cutoff at $E_{\nu,\rm c}\sim 50$ TeV, which results from the convolution of individual neutrino spectra with different cutoff energies (equation (12)). The dispersion of the rest frame cutoff energies is mainly driven by the wide range of $\epsilon_{\rm B}$ values used in the simulations (see Section 2.1). $E_{\nu,\rm c}$ would shift towards higher energies, if $\epsilon_{\rm B}\geq 10^{-3}$ in all sources. The location of the exponential cutoff shown in Fig. 4 also depends on the assumption of Bohm acceleration. Were the acceleration process slower or suppressed (e.g., Murase et al., 2011; Metzger et al., 2016), the energy spectrum of the diffuse neutrino flux would steepen at a few TeV energy.

The cumulative neutrino emission is $\sim 10$ per cent of the observed IceCube neutrino flux above 60 TeV, if the accelerated proton energy fraction varies between 0.01 and 0.1 in different sources (solid blue lines). However, if $\epsilon_{\rm p}=0.5$ for all SNe IIn (i.e., higher than most optimistic theoretical predictions), the model-predicted neutrino flux (orange dashed lines) would exceed the observations. The IceCube flux could be still explained by SNe IIn with $\epsilon_{\rm p}=0.5$, if the SNe IIn fraction was one per cent, as shown by ZP16 (green solid line). Regardless, the 4-year HESE IceCube sample already excludes the most extreme value for the proton accelerated energy fraction, namely $\epsilon_{\rm p}<0.5$ for $\xi=0.04$. Similar constraints on $\epsilon_{\rm p}$ can be placed using the 6-year IceCube muon neutrino sample from the Northern Hemisphere (shaded band in Fig. 4). This sample has a higher energy threshold compared to the HESE sample ($\sim 194$ TeV compared to $\sim 60$ TeV, respectively), but higher precision in the derived flux normalization and spectral index. At $E_{\nu}\lesssim 500$ TeV, the model predictions for $\epsilon_{\rm p}=0.3$ and $\xi=0.04$ are in conflict with the $99\%$ contours of IceCube. The tension with IceCube’s 95% confidence limits is resolved, if the acceleration efficiency is as low as $(\epsilon_{\rm p}/0.2)\times(\xi/0.04)\lesssim 1$.

The contribution of SNe IIn at different redshift ranges to the total diffuse neutrino flux is illustrated in Fig. 5. Approximately 50 per cent of the total neutrino flux originates from sources at $z<1$ or $D_{\rm L}<6.6$ Gpc. Sources at $1<z<2$ contribute about 30 per cent to the total flux, while the contribution of all sources beyond $z=2$ drops to 20 per cent, as a result of the decreasing volumetric production rate of CC SNe and decreasing neutrino fluence from a single source, which cannot match the increase of the comoving volume with redshift.

3.2 SNe IIn as neutrino point sources

We present the evolution of the neutrino spectrum from a fiducial SN IIn located at $D_{\rm L}=10$ Mpc. We assume that neutrino production begins after the shock breakout and terminates when the shock decelerates or reaches the outer edge of the CSM extended shell. Neutrino production is expected to be less important after this point because either the shock decelerates or it encounters a sharp negative density gradient at $r\gtrsim r_{\rm w}$.

The neutrino flux is expected to reach its maximum value soon after the shock breakout, with the exact peak time depending on the neutrino energy. This will be exemplified below through numerical calculations. After the peak time, the flux typically decays as a power law with time, i.e., $F_{\nu}(E_{\nu})\propto t^{-2+s}$. Although the exact value of $s$ depends on the model parameters, $s$ is bound between 1 and 2. The two limiting values are obtained in the following regimes:

• •

$s=1$, when proton injection is balanced by adiabatic losses. The neutrino flux then scales as:

[TABLE]

• •

$s=2$, when proton injection is balanced by p-p losses. The neutrino flux is given by (see also Zirakashvili & Ptuskin, 2016):

[TABLE]

By comparing $t_{\rm pp}$ with the shock’s dynamical timescale ($r/v_{\rm sh}$), we find that the second regime is relevant for

[TABLE]

where $\beta_{\rm sh}\equiv v_{\rm sh}/c$.

We consider two scenarios (henceforth, S1 and S2) for the neutrino production. Their parameter values used are listed in Table 1. In S1 the parameters $r_{\rm w},M_{\rm ej},M_{\rm csm},v_{\rm sh},\epsilon_{\rm B},$ and $\epsilon_{\rm p}$ were assigned to the median values of the corresponding distributions (see Section 2.1). Thus, S1 is a representative scenario of the simulated SN IIn population. The impact of different parameters is illustrated via S2, where a ten times larger $r_{\rm w}$ was adopted, while $M_{\rm ej},M_{\rm csm},v_{\rm sh},\epsilon_{\rm B},$ and $\epsilon_{\rm p}$ were kept fixed (see Table 1). Hence, the shock breakout time is $t_{\rm i}=0.54$ d and the neutrino production lasts $\sim 10.7$ yr.

Snapshots of the $\nu_{\mu}+\bar{\nu}_{\mu}$ energy spectrum at different times following the shock breakout are presented in Fig. 6 for S1 (left-hand panel) and S2 (right-hand panel). The results are obtained after taking into account neutrino mixing due to oscillations. The soft energy spectrum of atmospheric muon neutrinos (HKKMS 2007 – Honda et al., 2007; Aartsen et al., 2015d) is overplotted for comparison (thick black line). The neutrino flux increases rapidly at early times, but later decreases with a slower rate, which depends on the specifics of the source. Meanwhile, the maximum neutrino energy is increasing due to the increasing maximum energy of the parent protons. This is evident in both cases during the first year. At late times, where the adiabatic energy losses are more important than those caused by p-p collisions, the maximum energy of protons and, in turn, neutrinos, remains constant. Even in the optimistic scenario, where particle acceleration proceeds at the fastest possible rate, the neutrino spectrum from S1 barely extends beyond 1 PeV, as shown in Fig. 6. However, stronger magnetic fields, faster shocks, and higher mass-loading parameters may result in multi-PeV neutrino cutoff energies (see equation (12)).

The temporal evolution of the $\nu_{\mu}+\bar{\nu}_{\mu}$ neutrino flux at different energies is illustrated in Fig. 7. For S1, the peak flux in the energy range 1-10 TeV is expected within the first 50 days after the shock breakout. On the contrary, the 10-100 TeV flux remains approximately constant for a long period lasting hundreds of days. No significant energy dependance of the temporal evolution of the neutrino flux is found for S2, indicating an approximately constant spectral shape (see also Fig. 6). This, in turn, implies that adiabatic expansion governs the temporal evolution of the parent proton energy distribution. Substitution of the relevant $K_{\rm w}$ and $v_{\rm sh}$ values in equation (21) also shows that p-p collisions are more important than adiabatic expansion only for $t\lesssim 45$ d.

The expected IceCube $\nu_{\mu}+\bar{\nu}_{\mu}$ event rate at different energy bands and as a function of time is shown in Fig. 8 for S1 (thick lines) and S2 (thin lines). The rate has been calculated using the effective area of IceCube for horizontally up-going muon events (Aartsen et al., 2017b, declination -5°  to 30°). Dashed lines indicate the atmospheric neutrino event rate at the same energy bands (see legend). The event rate for S1 increases rapidly at early times in all energy ranges except for the $0.1-1$ PeV band. The late-time increase of the rate in this case ($t>50$ d) reflects the increasing maximum proton and neutrino energies with time. After the peak flux has been reached, the neutrino flux in all energy bands decays very slowly with time. This indicates that the temporal evolution of the parent proton distribution is dictated by p-p collisions (see equations (20) and (21)). On the contrary, the evolution of the proton distribution in S2 is governed by adiabatic losses. This is reflected on the temporal decay of the event rates which scale approximately as $\propto t^{-1}$ (see equation (19)). Furthermore, the event rates in S2 are typically one order of magnitude lower than those otbained in S1, since ${\rm d}N_{\nu}/{\rm d}t\propto K_{\rm w}\epsilon_{\rm p}v_{\rm sh}^{3}$ and $K_{\rm w}$ in S2 is ten times lower than in the first scenario (see Table 1).

Table 2 shows the cumulative event number over a period of 385 days (S1) and 10.7 yr (S2) in different energy bins. For comparison reasons, the atmospheric neutrino background for 385 d of exposure in a location of $0.75^{\circ}$ around the SN location is also shown. IceCube’s sensitivity on a timescale of $\sim 1$ yr (as in S1) translates to a mean of 3.2 events for an unbroken $E_{\nu}^{-2}$ spectrum (Fig. 1 Aartsen et al., 2015a) after taking into account the atmospheric background, while 11 events are required for a 5$\sigma$ discovery. For $E_{\nu}^{-2}$ spectra that have a sharp cutoff between 100 TeV and 1 PeV, as in the scenarios discussed here, the sensitivity worsens by a factor of $2$ in flux $E_{\nu}F_{\nu}$ (Fig. 3 Aartsen et al., 2016c) or by $36\%$ in total event count after taking into account the IceCube effective area. Henceforth, we adopt $N_{\rm s}=4.4$ and $N_{\rm dp}=15$ as the event counts required to reach IceCube’s sensitivity and $5\sigma$ discovery potential, respectively, in a period of 385 d. We note that these event counts are applicable to the Northern Sky only, where IceCube selects a pure sample of muon neutrinos. Due to the cutoff in the neutrino spectrum lying between 100 TeV and 1 PeV, absorption effects in the Earth can be neglected.

Based on the discussion above and on the results presented in Table 2), IceCube is able to limit $\epsilon_{\rm p}<0.06$ for a S1-like SN IIn at 10 Mpc distance. The constraint on the proton acceleration efficiency is more stringent than the one obtained from the analysis of the diffuse neutrino flux, where $\epsilon_{\rm p}\lesssim 0.2$ for a SNe IIn rate of $\xi=0.04$ (see Section 3.1). Were a S1-like SN IIn detected at $\sim 18$ Mpc from Earth, the point source neutrino analysis would result in $\epsilon_{\rm p}\lesssim 0.2$, i.e., it would be as constraining as the diffuse neutrino analysis. The second scenario (S2) is more challenging for neutrino detection due to its lower neutrino rate. Although the duration in S2 is ten times longer than in S1, $\sim 45\%$ of the signal is expected in the first year (Table 2). The background of atmospheric neutrinos is, however, constant in time. Hence, due to the decline in the neutrino rate of S2, the best signal-to-noise ratio is obtained for a period of $200-400$ d, similar to the full duration in S1. As a result, IceCube would be able to constrain $\epsilon_{\rm p}<0.2$, only if a S2-like SN IIn was detected at $\lesssim 3.4$ Mpc.

4 Discussion

SNe surrounded by dense CSM pose an interesting alternative to extragalactic scenarios of high-energy neutrino emission. In this paper, we have calculated the diffuse neutrino emission from this rare class of CC SNe using a Monte Carlo approach that allowed us to incorporate their widely ranging properties, such as mass-loss rates, wind velocities, and proton acceleration efficiencies. We have also demonstrated through two indicative scenarios, the temporal dependence of the neutrino flux from individual SNe and evaluated the possibility of being detected as neutrino point sources with IceCube. In the following, we discuss the main caveats of the model and expand upon several of its aspects, such as possible associations between the IceCube neutrinos and SNe IIn.

4.1 Caveats

The forward shock is expected to be mildly decelerating, unless the SN ejecta have a flat density profile, namely $\rho_{\rm ej}\propto r^{-n_{\rm ej}}$ with $n_{\rm ej}<3$ (e.g. Chevalier, 1982). For example, $v_{\rm sh}\propto r^{-0.1}$ for a wind-like CSM ($n\propto r^{-2}$) and ejecta with $n_{\rm ej}=12$. In contrast, the SN shock would propagate with a constant velocity, at least in the free expansion phase, if both media had uniform densities (e.g. Matzner & McKee, 1999). A slowly decreasing shock velocity would result in decreasing neutrino flux and maximum neutrino energy, as indicated by equations (20)-(19) and (12), respectively. A faster decrease of the neutrino flux with time could also be caused by the propagation of the SN shock to a CSM with steeper density profile than the one considered here (e.g. Chandra et al., 2015). Flatter density profiles, on the other hand, would enhance the neutrino production rate. Such deviations have been inferred for a handful of SNe IIn detected in radio wavelengths (e.g. Fransson et al. (1996); Immler et al. (2001), but see also Martí-Vidal et al. (2011a)).

So far, we have presented results for the high-energy neutrino emission from SNe IIn shocks propagating in dense CSM. Apart from the forward shock, there is a second shock wave that forms in the outer parts of the SN ejecta (reverse shock) and may also contribute to the neutrino signal. Its contribution mainly depends on its velocity, the proton accelerated energy, and the density of the shocked SN ejecta. As the maximum energy of accelerated protons and, in turn, neutrinos scales cubically or quadratically with the shock velocity (see equation (12)), small differences between the forward and reverse shock velocities will lead to large differences in the cutoff energy of the neutrino spectra (see also Fig. 1 in Murase et al. (2011)). Furthermore, the density of the ejecta and the energy of accelerated protons at the reverse shock affects the overall normalization of the neutrino spectrum. A detailed calculation of the neutrino emission from both shocks requires hydrodynamical calculations that are, however, beyond the scope of this paper. Nevertheless, the contribution of the reverse shock to the total non-thermal (i.e., radio synchrotron and neutrino) emission is expected to be small, as long as the density profile of the ejecta is steep (Chevalier & Fransson 2003; see also Fig. 3 in ZP16).

While setting up our Monte Carlo simulations for the diffuse neutrino emission from SNe IIn, we treated the shock velocity and ejecta mass as independent variables and let the SN kinetic energy ($E_{\rm k}$) to be a derived quantity. This choice may, in principle, lead to extreme values of the SN kinetic energy ($E_{\rm k}\gg 10^{51}$ erg). To check this possible caveat, we computed the distribution of $E_{\rm k}$ values using the generated distributions of the shock velocity and ejecta mass (Section 2.1). The median SN kinetic energy in our simulations is found to be $5\times 10^{51}$ erg. Moreover, 68 per cent of the simulated sources had kinetic energies between $10^{51}$ and $2\times 10^{52}$ erg, which is not uncommon for SNe IIn; for example, Smith et al. (2013) showed that $E_{\rm k}>10^{51}$ for 2009ip and 2010mc. The small fraction of sources with $E_{\rm k}\gg 10^{52}$ erg does not affect our main conclusions about the diffuse neutrino emission (Section 3.1), especially when other sources of uncertainty are considered, i.e. production rate of SNe IIn (see following paragraph). The injected energy into accelerated protons can also be estimated as $E_{\rm p}=\epsilon_{\rm p}E_{\rm k}$. The distribution of $E_{\rm p}$ derived by our simulations has a median value of $1.5\times 10^{50}$ erg, while 68 per cent of the simulated sources have $3\times 10^{49}\leq E_{\rm p}\leq 8\times 10^{50}$ erg. These values should be compared with value $5\times 10^{51}$ erg adopted by ZP16 for all sources.

We have calculated the diffuse neutrino emission from the SN IIn class by adopting the CC volumetric rate of Hopkins & Beacom (2006) (hereafter, HB06) and assuming that 4 per cent of all CC SNe are of the type IIn. The rates predicted by the HB06 model are two times higher than those obtained by the model of Madau & Dickinson (2014)–henceforth MD14, at all redshifts (Cappellaro et al., 2015). Two recent extragalactic surveys performed with the Hubble Space Telescope provided volumetric SN rates in high-z galaxies ($z\lesssim 2.5$), which lend further support to the MD14 model (Strolger et al., 2015). However, the existing observations do not exclude the HB06 model, especially due to the uncertainties that enter in the model, such as the mass range of the progenitor stars (see Section 8.1 in Cappellaro et al. (2015) for details). The rate of SNe with dense and massive CSM is also uncertain. Here, we adopted $\xi=0.04$, but this could range from 1 per cent to $\sim$6 per cent (e.g. Smith, 2014; Strolger et al., 2015; Cappellaro et al., 2015). Adoption of the MD14 model with $\xi\sim 0.08$ would, thus, lead to similar results for the diffuse neutrino emission as those presented in Fig. 4.

The largest uncertainties ($\sim 50$ per cent) entering in the estimation of the diffuse neutrino flux are those described above. If the SNe IIn rate was lower by a factor of two, IceCube HESE observations would limit the proton accelerated energy fraction at $\lesssim 0.45$, or $\epsilon_{\rm p}\lesssim 0.3$ after including the IceCube up-going muon neutrino flux above 195 TeV.

4.2 Constraints on $\epsilon_{\rm p}$

Numerical simulations of particle acceleration at non-relativistic shocks are the main means of deriving the proton acceleration efficiency $\epsilon_{\rm p}$ from first principles (Caprioli & Spitkovsky, 2014; Park et al., 2015). Additional constraints on the acceleration efficiency can be placed through the detection (or non-detection) of multi-messenger signatures associated with relativistic protons. In this study, we have demonstrated how high-energy neutrino observations may be used to assess $\epsilon_{\rm p}$ at SN shocks propagating in dense CSM. In particular, the IceCube’s measurement of the diffuse neutrino flux already constrains $\epsilon_{\rm p}\lesssim 0.3\times(0.04/\xi)$ or $\epsilon_{\rm p}\lesssim 0.2\times(0.04/\xi)$ when including the up-going muon neutrino sample from the Northern Sky (Section 3.1). Even more stringent constraints on the acceleration efficiency may be placed by direct observations of close-by SNe IIn with neutrino telescopes (Section 3.2).

The neutrino flux from a SN IIn depends on its distance, shock velocity, mass-loading parameter, and proton acceleration efficiency (Section 3.2). All parameters, except for the latter, can be, in principle, inferred from radio and optical observations of SNe IIn. For example, $v_{\rm sh}$ is routinely inferred from the width of optical emission lines, while $K_{\rm w}$ can be derived by the fitting of radio observations (Chevalier, 1998; Chevalier et al., 2006; Chandra et al., 2015). Besides the intrinsic source parameters that affect the neutrino luminosity, as exemplified through scenarios S1 and S2, the actual distance of the source has the strongest impact on the IceCube expected neutrino rate (see equation (15)).

The IceCube $\nu_{\mu}+\bar{\nu}_{\mu}$ event number expected within the first year in both scenarios discussed in Section 3.2, is plotted in Fig. 9 (bottom panel) as a function of the source distance. Different types of lines correspond to three values of the proton acceleration efficiency marked on the plot. The curves are obtained after scaling the total number of muon neutrinos obtained in S1 and S2 (see Table 2) with the source distance and $\epsilon_{\rm p}$. The top panel in Fig. 9 shows the all-sky number of SNe IIn expected over 8 years up to a distance $D_{\rm L}$ (black thick line); the number is derived using the CC SN rate of HB06 and $\xi=0.04$. For $\epsilon_{\rm p}$ values allowed by the diffuse flux measurements, i.e., $\epsilon_{\rm p}\lesssim 0.2\times(0.04/\xi)$, only S1-like SNe IIn at distances $\lesssim 18$ Mpc are strong enough neutrino emitters to constrain less efficient acceleration scenarios, in agreement with the predictions by Murase et al. (2011). If $\epsilon_{\rm p}=0.2$, IceCube would be capable of claiming a neutrino detection from an S1-like SN IIn exploding at $\sim 10$ Mpc. For SNe IIn with lower mass-loading parameter $K_{\rm w}$, as discussed in S2, the neutrino rate decreases significantly, even for $\epsilon_{\rm p}=0.2$. This restricts the accessible distance to no more than 4 Mpc.

S1-like sources with faster shocks are more promising candidates for neutrino detection. For higher shock velocities, the neutrino rate increases as $\propto v_{\rm sh}^{3}$ (see equation (20)), while the duration of the neutrino production decreases with $v_{\rm sh}^{-1}$. The latter leads to an accordingly lower number of background events, thus increasing the signal-to-noise. Furthermore, the higher shock velocity results in an increased cutoff value of the neutrino energy spectrum, as this scales with $v_{\rm sh}^{2}$ (Section 3.2). For example, a source with a three times higher shock velocity than in S1 ($v_{\rm sh}=3\times 10^{4}$ km s*-1*) and the same $K_{\rm w},\epsilon_{\rm p}$ values, would yield $\sim 38$ events above $10$ TeV. Hence, IceCube would be sensitive to such a source up to 40 Mpc even for $\epsilon_{\rm p}=0.03$. A positive detection would also be possible for a source located at $\lesssim 20$ Mpc. Given that IceCube’s sensitivity improves when the neutrino spectrum extends above 100 TeV (Fig. 3 in Aartsen et al., 2016c), these estimates are rather conservative.

Another important aspect in the search of neutrino point sources is the rate of SNe IIn explosions at the relevant distances (i.e., $\lesssim 40$ Mpc). The expected number of SNe IIn at distances $\lesssim 22$ Mpc (40 Mpc) is $0.7^{+0.9}_{-0.5}$ ($4.4^{+2.1}_{-1.4}$)333The errors correspond to 1 $\sigma$ errors of Gaussian statistics. within 8 years in the Northern Sky (top panel in Fig. 9). During the period of 2010-2016 three SNe IIn were detected in the Northern Sky at a mean distance of $\sim 29$ Mpc (Table 3). Stacking of their neutrino signal could place stronger constraints on $\epsilon_{\rm p}$ than the diffuse neutrino flux by a factor of $\sqrt{N_{\rm IIn}}$.

Neutrino production in p-p collisions is accompanied by the injection of relativistic electron-positron pairs in the post-shock region and the production of GeV $\gamma$-ray photons via the decay of neutral pions. The first systematic search for $\gamma$-ray emission in Fermi-LAT  data from the ensemble of 147 SNe IIn exploding in a dense CSM was recently presented by Ackermann et al. (2015). No significant excess above the background was found leading to model-independent $\gamma$-ray flux upper limits. These $\gamma$-ray non-detections constrained the ratio of the $\gamma$-ray to optical luminosity in the range 0.01-1. These can, in principle, be translated to limits on $\epsilon_{\rm p}$. However, due to the uncertainty in the escaping fraction of $\gamma$-rays, no stringent limits could be placed on the proton acceleration efficiency.

4.3 IceCube neutrinos in coincidence with SNe IIn?

Table 3 shows the detected SNe IIn in the local Universe (i.e., within a distance $D_{\rm L}=40$ Mpc). Out of the 29 high-energy muon neutrinos detected in the Northern Sky above 200 TeV (Aartsen et al., 2016d), none is found to be in spatial or temporal coincidence with these SNe IIn. However, one cascade-like event444ID 16 with deposited energy of 30.6 TeV and 19.4 degrees median angular uncertainty. from the 4-year HESE IceCube sample (The IceCube Collaboration et al., 2015) was detected on MJD 55798.63, i.e. $\delta t=1.13$ days after the peak of the SN2011fh optical light curve, with an angular offset of 7.12 degrees.555This neutrino event has another possible connection to SN2011A which was detected 236.13 days prior to the neutrino at an angular separation of 8.61 degrees.

The coincident rate of SNe IIn in the local Universe ($z\ll 0.1$) during the time window $\delta t$ and within the positional uncertainty $\sigma_{\nu}$ of a neutrino event can be estimated as:

[TABLE]

For a IIn fraction of $\xi=0.04$ and a local rate of CC SNe $10^{-4}$ Mpc*-3* yr*-1*, the observation of the cascade event ID 16 deviates from the background expectation (i.e., $N_{\rm IIn}=0$ in $\delta t$) at $2.76\,\sigma$. A similar result was obtained after scrambling the neutrinos of the 4-year HESE sample in right ascension and time or the SNe IIn (Table 3, in direction and time if detected within IceCube livetime) and counting the number of random associations of one SN IIn at an angular distance smaller than the neutrino uncertainty and at a time $\delta t<1.13$ days. If the cascade-like event was indeed physically associated to SN2011fh, additional muon neutrino events should have been detected by the ANTARES neutrino telescope, as this has three to four times larger effective area at 30 TeV than IceCube for up-going muon neutrino events (Adrian-Martinez et al., 2014).

Not much information is available for SN2011fh besides the velocity of its ejecta $v_{\rm ej}=7400$ km s*-1*, which is a good indicator of the shock velocity ($v_{\rm sh}\gtrsim v_{\rm ej}$). Inspection of Fig. 9 shows that neither IceCube nor ANTARES would be sensitive to neutrinos from SN2011fh located at 36 Mpc, if this had the same properties as the source in scenario S1. In addition, ID 16 is an event produced within the detector (member of the HESE sample). The rate of such events is much lower (by a factor $\sim 180$ at 30 TeV) than the rate of up-going muon neutrinos (Table 2), which was used in producing Fig. 9. Regardless, a physical association of SN2011fh with the cascade event ID 16 would require higher $K_{\rm w}$ or/and $\epsilon_{\rm p}$ than those used in S1.

An estimate of the mass-loading parameter and proton acceleration efficiency in SN2011fh can be derived assuming a physical connection with the cascade event ID 16 and using our estimates of the expected neutrino rate. In the energy range of 10–100 TeV, we expect $N_{\rm S1}\sim 1$ up-going muon neutrino events (see Table 2). One neutrino (ID 16) was observed in the starting event channel, which has a $\sim 200$ times smaller signal expectation in the energy range 15–45 TeV. The $90\%$ limits for the detection of one starting event with no background is $N_{\rm 2011fh}=0.11\dots 4.36$ (Feldman & Cousins, 1998). Assuming a year-long duration, the product $K_{\rm w}\times\epsilon_{\rm p}$ for SN2011fh can be estimated as (see also equation (20)):

[TABLE]

Assuming that in SN2011fh $\epsilon_{\rm p}=0.2$ (i.e., the largest allowed value from the diffuse neutrino flux measurements), its mass-loading parameter has to be at least $K_{\rm w}^{\rm 2011fh}\sim 87\,K_{w}^{\rm S1}\sim 4\times 10^{18}$ g cm*-1*. Radio observations of SN2011fh could place strong constraints on $K_{\rm w}^{\rm 2011fh}$ and, ultimately, exclude (or not) a physical connection to neutrino ID 16.

4.4 Diffuse $\gamma$-ray emission

The neutrino emission resulting from the SN IIn class is accompanied by a diffuse $\gamma$-ray component, which may contribute to the isotropic $\gamma$-ray background (IGRB) at GeV energies due to cascades initiated by the absorption of multi-TeV photons on the extragalactic background light (e.g. Murase et al., 2012). Assuming that the sources are optically thin to photon-photon absorption, Murase et al. (2013) showed that generic p-p scenarios of neutrino production with spectra ${\rm d}N_{\nu}/{\rm d}E_{\nu}\propto E_{\nu}^{-2}$, for $E_{\nu}\leq E_{\nu,\rm b}\sim 1$ PeV and $\propto E_{\nu}^{-p}$ otherwise, could explain IceCube observations above 100 TeV, if $p\lesssim 2.1-2.2$. In general, tighter constraints on the power-law index of the relativistic proton distribution can be obtained (i.e., $p\sim 2$), if $E_{\nu,\rm b}\lesssim 30$ TeV. The induced $\gamma$-ray emission from generic p-p scenarios of neutrino emission was shown to be marginally consistent with the IGRB (Bechtol et al., 2015; Murase et al., 2016), if $p=2.5$ and the break energy was fixed to the lowest energy bin of the combined neutrino data between 25 TeV and 2.8 PeV (Aartsen et al., 2015b). The tension with the IGRB can be relaxed if a fraction of the multi-TeV $\gamma$-rays is absorbed in the source (Murase et al., 2016), which is not unexpected for the early-time evolution of SNe with dense CSM (Kantzas et al., 2016). An alternative way of reconciling the p-p scenario with the IGRB measurements would be $E_{\nu,\rm b}>25$ TeV, which is still consistent with the neutrino data.

We showed that the diffuse neutrino flux (per flavour) from SNe IIn is $\sim 4\times 10^{-9}$ GeV cm*-2* s*-1* sr*-1* or $\sim 10$ per cent of the observed IceCube flux above 100 TeV, if the proton acceleration efficiency varies between 0.01 and 0.1 among the sources (Section 3.1). Assuming no attenuation of multi-TeV $\gamma$-rays in the sources, the diffuse $\gamma$-ray flux from SNe IIn at $\sim 50$ GeV is expected to be $\sim 10^{-8}$ GeV cm*-2* s*-1* sr*-1* (see also Fig. 4 in Bechtol et al., 2015). Only if all shocks in SNe IIn had the same acceleration efficiency $\epsilon_{\rm p}\sim 0.2\times(0.04/\xi)$, could the SN IIn class explain the $>100$ TeV neutrino observations. In this case, the accompanying diffuse $\gamma$-ray emission would be marginally consistent with the IGRB, but it would exceed the non-blazar contribution to the EGB (Bechtol et al., 2015; Ackermann et al., 2016).

5 Conclusions

We have evaluated the possible neutrino signal of SNe IIn and placed constraints on the proton accelerated energy fraction, $\epsilon_{\rm p}$, by means of diffuse and point-source neutrino observations with IceCube. By employing a Monte Carlo method that takes into account the wide spread in the physical properties of SNe IIn, we showed that the diffuse neutrino emission from SNe IIn can account for $\sim 10$ per cent of the observed IceCube neutrino flux above 100 TeV. In the less realistic scenario, where the proton accelerated energy fraction is the same for all SN shocks, we showed that the observed diffuse astrophysical neutrino spectrum could be explained totally by SNe IIn, if $\epsilon_{\rm p}\lesssim 0.2$ and 4 per cent of the CC SNe were of the IIn type. However, the identification of a single SN IIn as a neutrino point source with IceCube using up-going muon neutrinos could place stronger constraints on $\epsilon_{\rm p}$. We concluded that such an identification is possible in the first year following the SN shock breakout for sources within $\lesssim$18 Mpc and $\epsilon_{\rm p}\lesssim 0.2$ or $\lesssim 7$ Mpc and $\epsilon_{\rm p}\lesssim 0.03$. Interestingly, one cascade-like event (ID 16) from the 4-year HESE sample of IceCube was found to be in spatial agreement with SN2011fh ($D_{\rm L}=$36 Mpc) and was detected only 27.12 hours after the maximum optical SN light. The probability that this association was not a chance coincidence was found to be $2.76\,\sigma$. In case of a positive connection, additional muon neutrinos should be detected by ANTARES in coincidence with SN2011fh, which should have a very high mass-loading parameter ($>3.9\times 10^{18}$ g cm*-1*). Analysis of propitiatory ANTARES data and of SN2011fh radio observations are strongly encouraged to resolve the nature of this association.

Acknowledgments

We thank the anonymous referee, E. Resconi, and K. Kotera for insightful comments on the manuscript. We also thank K. Murase and A. Mastichiadis for useful discussions. We thank A. Hopkins for providing the theoretical curves of the volumetric production rate of Type II supernovae. SC is supported by the cluster of excellence “Origin and Structure of the Universe” of the Deutsche Forschungsgemeinschaft. GV acknowledges support from the BMWi/DLR grant FKZ 50 OR 1208.