Identifying High Energy Neutrino Transients by Neutrino Multiplet-Triggered Followups

The Emission of neutrinos from transient sources can be characterized by the integral luminosity $L_{\nu}$ (defined for the sum of all flavors), the flare duration $\Delta T$ in the source frame, and the neutrino energy spectrum $\phi_{\nu}^{\rm fl}$. The total energy output by a neutrino emission is given by ${\mathcal{E}}_{\nu}^{\rm fl}=L_{\nu}\Delta T$.

The neutrino spectrum $\textstyle{{d\dot{N}_{\nu_{e}+\nu_{\mu}+\nu_{\tau}}}/{d\varepsilon_{\nu}}}$ is assumed to follow a power-law form, with the reference energy $\varepsilon_{0}$ and the flux from a single source at the redshift $z$ is given as

where $\varepsilon_{\nu}$ and $E_{\nu}=\varepsilon_{\nu}{(1+z)}^{-1}$ are the neutrino energies at the time of emission and arrival at Earth’s surface, respectively. In our model, the normalization constant $\kappa$ is bolometrically associated with the source luminosity as

by integrating over the energy range $[0.1\varepsilon_{0},10\varepsilon_{0}]$ to account the neutrino energetics around $\varepsilon_{0}$ reasonably. The proper distance, $d_{z}$, is calculated via

We hereafter adopt $\alpha_{\nu}=2.3$, $\varepsilon_{0}=100\ {\rm TeV}$, and $\Delta T=30\ {\rm days}$ in our generic model construction. These are set to be consistent with current neutrino data ($\alpha_{\nu}$), within the representative energy range in IceCube measurements ($\varepsilon_{0}$) (Abbasi et al., 2022), and within the expected timescale of neutrino flares in the majority of optical transient sources ($\Delta T$).

A population of neutrino sources contribute to the diffuse cosmic background flux. Assuming emission from standard candles (i.e., identical sources over redshifts), the energy flux of diffuse neutrinos from these sources across the universe, $\Phi_{\nu}\equiv dJ_{\nu}/dE_{\nu}$, is calculated by (e.g., Murase et al. 2016)

where $n_{0}\psi(z)$ is the comoving source number density given the local source number density $n_{0}$ and the cosmological evolution factor $\psi(z)$. Sources are distributed between redshift $z_{\rm min}$ and $z_{\rm max}$. We define the burst rate per unit volume as $R_{0}$, and assume the local density is effectively given by $n_{0}=R_{0}\Delta T$. The diffuse cosmic background flux in the present model is, therefore, described by ${\mathcal{E}}_{\nu}^{\rm fl}$, $R_{0}$, $\Delta T$, and the evolution factor $\psi(z)$. The latter is parametrized as $(1+z)^{m}$, the functional form often used in the literature. The evolution factor for the transient neutrino sources is unknown, but many of the proposed sources are related to SNe which approximately trace the star formation rate (SFR). We therefore adopt the following parameterization, which approximately describes SFR (Yoshida & Takami, 2014), as our baseline model

The diffuse cosmic background flux from the sources following an evolution factor other than an SFR-like evolution can be approximately estimated by scaling

where

is the effective evolution term and $\xi_{z}^{\rm SFR}\approx 3$.

The resultant diffuse flux, $\Phi_{\nu}$, limits the range of the source parameters ${\mathcal{E}}_{\nu}^{\rm fl}$ and $R_{0}$ as the flux must be consistent with the IceCube measurements. Figure 2 shows $\textstyle{E_{\nu}^{2}\Phi_{\nu}(3/\xi_{z})}$ as a function of ${\mathcal{E}}_{\nu}^{\rm fl}$ and $R_{0}$. IceCube data suggests $E_{\nu}^{2}\Phi_{\nu}\lesssim 10^{-7}\ {\rm GeV}\ {\rm cm}^{-2}$ s^-1 sr^-1 (e.g. Aartsen et al. 2020a). Any parameter combination that yields $\textstyle{E_{\nu}^{2}\Phi_{\nu}}$ above this limit would overproduce diffuse flux. If $E_{\nu}^{2}\Phi_{\nu}\lesssim 10^{-9}\ {\rm GeV}\ {\rm cm}^{-2}$ s^-1 sr^-1, the contribution to the TeV neutrino sky background would be negligible and so is not considered further. Hence, we limit the source parameter space to meet $10^{-9}\ {\rm GeV}\ {\rm cm}^{-2}{\rm s}^{-1}{\rm sr}^{-1}\leq E_{\nu}^{2}\Phi_{\nu}(3/\xi_{z})\leq 10^{-7}\ {\rm GeV}\ {\rm cm}^{-2}{\rm s}^{-1}{\rm sr}^{-1}$.

The number of transient sources that could produce detectable neutrino multiplets for a given neutrino telescope is given by (Murase & Waxman, 2016)

where $\Delta\Omega$ is the solid angle for a given direction of the neutrino multiplet and $P_{\rm p}^{n\geq 2}$ is the Poisson probability of producing multiple neutrinos for the mean number of neutrinos $\mu^{\rm s}$ from a source at redshift $z$, given by

Note that the $1/3$ factor is applied for the conversion of the all-flavor-sum neutrino flux to that of per-flavor assuming the equal neutrino flavor ratio. The muon neutrino detection effective area, $A_{\nu_{\mu}}$, determines the event rate. We use an underground neutrino telescope model with a 1 km³ detection volume  (Gonzalez-Garcia et al., 2009; Murase & Waxman, 2016) to estimate $A_{\nu_{\mu}}$ in our study. Note that $P_{\rm P}^{\rm n=2}\sim\mu_{\rm s}^{2}/2\propto d_{z}^{-4}$ when $\mu^{\rm s}\ll 1$. Therefore, the integrand of $N^{\rm M}_{\Delta\Omega}\propto d_{z}^{-2}$, which indicates that only nearby sources produce multiple events. $N^{\rm M}_{\Delta\Omega}$ is thus sensitive to the minimal value of $d_{z}$, which is determined by $z_{\rm min}$. In our work, it is defined as the average interval of source locations in the local universe, $\textstyle{r=(3/(4\pi n_{0}))^{1/3}}$.

In order to estimate the rate and the resultant sensitivity for detecting the multiple neutrino events, we set the solid angle $\Delta\Omega$ to be comparable to the pointing resolution of neutrino events. We assume $\Delta\Omega=1\ {\rm deg}^{2}$ hereafter, which is comparable to the angular resolution of track-like events recorded by IceCube (Aartsen et al., 2014).

The left panel of Figure 3 shows $N^{\rm M}_{\Delta\Omega}$ as a function of our source characterization parameters, ${\mathcal{E}}_{\nu}^{\rm fl}$ and $R_{0}$. As $N^{\rm M}_{\Delta\Omega}$ represents the number of sources found within a solid angle of $\Delta\Omega$ during the multiplet search time frame $T_{\rm w}$, the expected number of sources to produce multiple neutrino events in the $2\pi$ sky detected in the observation time $T_{\rm obs}$ is given by

Hence, the parameter space of $N^{\rm M}_{\Delta\Omega}\gtrsim 10^{-6}$ is accessible by a 1 km scale neutrino telescope ¹¹1i.e. the absence of neutrino multiples, $N_{\rm all}^{\rm M}\leq 1$, can be used to constrain the population of neutrino source candidates (Murase & Waxman, 2016; Ackermann et al., 2019)..

Note that ${\mathcal{E}}_{\nu}^{\rm fl}=L_{\nu}\Delta T$ represents the total energy of neutrinos if $T_{\rm w}>(1+z)\Delta T$. In general, $L_{\nu}T_{\rm w}$ can be smaller than ${\mathcal{E}}_{\nu}^{\rm fl}$ if $T_{\rm w}<(1+z)\Delta T$.

The number of sources to produce multiplet events $N^{\rm M}_{\Delta\Omega}$ given by Eq. (10) is equivalent to the multiplet signal rate from the viewpoint of a neutrino detector. However, searches for neutrino multiple events from these sources are contaminated by the atmospheric and the unresolved diffuse cosmic neutrino backgrounds. The expected atmospheric background for $T_{\rm w}=30$ days and $\Delta\Omega=1$ deg² reaches $\sim 0.5$ events in the solid angle average for a km³ volume neutrino detector, which would smear out the neutrino signals from the source. Because the atmospheric neutrino spectrum follows $\sim E_{\nu}^{-3.7}$ which is much softer than $\phi_{\nu}^{\rm fl}\propto E_{\nu}^{-\alpha_{\nu}}$, taking into account the energy of each of the multiple neutrino events can adequately remove the background contamination. For this purpose, we construct the extended Poisson likelihood functions for the signal hypothesis $\mathcal{L}^{\sf sig}$ to describe the case when the detected multiple events originate in transient sources, and $\mathcal{L}^{\sf BG}$ for the background hypothesis as

Here $\mu_{\rm atm}\simeq 0.5$, and $\mu_{\rm dif}\lesssim 10^{-3}$ are the expected mean number events from the atmospheric and the unresolved diffuse neutrino backgrounds, respectively. $P_{\rm s}^{\rm E},P_{\rm atm}^{\rm E}$, and $P_{\rm dif}^{\rm E}$ are the probability density functions (pdf) of the energies of neutrinos from the transient sources, the atmospheric background, and the diffuse flux. These are obtained by multiplying $A_{\nu_{\mu}}$ with each of the neutrino fluxes. In the limit $N_{\rm PS}^{\rm M}\ll 1$ and $\mu_{\rm atm}+\mu_{\rm dif}\ll 1$ where we consider only the doublet case, this likelihood can be described by a more intuitive formulas

The test statistic for the background-only hypothesis is constructed with the log likelihood ratio

where the hat notation represents the value that maximizes the likelihood. The test statistic for the signal hypothesis with $({\mathcal{E}}_{\nu}^{\rm fl},R_{0})$ is

Figure 4 shows the calculated p-values as a function of neutrino energy in a doublet, $E_{\nu}=E_{\nu}^{1}=E_{\nu}^{2}$. We set $N^{\rm M}_{\Delta\Omega}$ = $3\times 10^{-6}$ for illustrative purpose. The p-values that support the signal hypothesis reach a plateau $\textstyle{e^{-(\mu_{\rm atm}+\mu_{\rm dif})}N_{\Delta\Omega}^{\rm M}\simeq 1.8\times 10^{-6}}$ when the doublet neutrino energy gets higher, as expected. The $2\pi$ sky converted annual false alarm rate (FAR) is $\sim 0.25~{}\mathrm{yr}^{-1}$ when the doublet energy is higher than $\varepsilon_{0}=100$ TeV. Note that the intensity of the prompt atmospheric neutrino component (Bhattacharya et al., 2015) produced from the decay of heavy charmed hadrons was found to be subdominant compared to that of the astrophysical neutrinos (Abbasi et al., 2022), and the statistical test discussed here is robust against the uncertainties originating in the heavy hadron physics.

The sensitivity for a given $N_{\rm PS}^{\rm M}({\mathcal{E}}_{\nu}^{\rm fl},R_{0})$ is evaluated by convolution of p-values for the multiplet source hypothesis, given by the test statistic Eq. (16), with the probability of detecting $E_{\nu}$ energy neutrinos

where $p(E_{\nu}^{1},E_{\nu}^{2})$ represents the p-value for doublet $(E_{\nu}^{1},E_{\nu}^{2})$. This indicates the effective rate (or p-value) of finding a source that produces multiple neutrino events, considering the atmospheric background contamination. The middle panel of Figure 3 shows the effective p-values, $P^{\rm eff}_{\rm m}$. For a given neutrino source parameter set $({\mathcal{E}}_{\nu}^{\rm fl},R_{0})$, this shows how often we can see multiple neutrinos that are inconsistent with the atmospheric neutrino background for $\Delta\Omega=1$ deg² and $T_{\rm w}=$ 30 days. The domain of $P^{\rm eff}_{\rm m}\gtrsim 10^{-6}$ is reachable by a five year observation with $2\pi$ sky coverage.

A null detection of any multiple neutrino events by a 1 km³ neutrino telescope constrains the neutrino source parameters $({\mathcal{E}}_{\nu}^{\rm fl},R_{0})$. Hence, criteria for rejecting the neutrino background-only hypothesis globally across the $2\pi$ sky must be adequately introduced. As an example, in the right panel of Figure 3, $N^{\rm M}_{\Delta\Omega}$ requires the criterion that with the local p-value of a neutrino multiplet for the background-only hypothesis must be less than $10^{-6}$. This corresponds to an annual global FAR of $\sim 0.25$ in the $2\pi$ sky. The parameter space where $N^{\rm M}_{\Delta\Omega}\gtrsim 10^{-6}$ can be constrained by a few years of observations when the rate of multiplet detection under this condition is consistent with the global FAR.

Neutrino sources that yield multiple neutrino events are likely to be optical transient objects. Searching for optical/NIR counterparts in followup observations to identify the neutrino emitter can be contaminated by the detection of more dominant SNe. We anticipate finding $\gtrsim 100$ SNe in a field of view of $\Delta\Omega=1\ {\rm deg}^{2}$ during a time window of $T_{\rm w}=30$ days. Their redshift distribution is, however, quite different from that expected for the neutrino multiplet sources. As we have already shown in Figure 1, most of the sources that produce multiple neutrino events are confined within the local universe. The resultant difference in the probability distribution of redshift between the neutrino multiplet sources and the unrelated SNe allows for a statistical test to indicate which hypothesis is more consistent with the observational data.

Among the candidate optical transient counterparts, the object with the minimum redshift, $\scalebox{1.2}{\it z}_{\scalebox{0.6}{ min}}^{\scalebox{0.6}{ trans}}$, is the most likely source of the neutrino multiplet event. The pdf of $\scalebox{1.2}{\it z}_{\scalebox{0.6}{ min}}^{\scalebox{0.6}{ trans}}$, which is defined by our signal hypothesis, is obtained by normalizing $N^{\rm M}_{\Delta\Omega}$ from Eq. (10):

and the likelihood is constructed by $\mathcal{L}_{\rm S}^{\rm M}(\scalebox{1.2}{\it z}_{\scalebox{0.6}{ min}}^{\scalebox{0.6}{ trans}})=\rho_{z_{\rm min}}^{\rm M}({\mathcal{E}}_{\nu}^{\rm fl},R_{0},z_{\rm min}^{\rm tran})$

In the background hypothesis case, the closest counterpart object belongs to the population of SNe not associated with the neutrino detection. The number of unrelated SNe within redshift $z$ in a field of view of $\Delta\Omega$ is given by

Here, $\psi_{\rm SN}(z)$ is the cosmological evolution term, which is assumed to follow an SFR-like distribution, represented by Eq. (7), and $n^{\rm SN}_{0}$ is the average number density of SNe in the present epoch, which is effectively obtained from the SNe rate, $R^{\rm SN}_{0}$, as

$R^{\rm SN}_{0}=1.3\times 10^{-4}$ Mpc^-3 yr^-1 is assumed as the nominal value in our work (Madau & Dickinson 2014 and references therein). The pdf of its redshift is given by

The pdf of $\scalebox{1.2}{\it z}_{\scalebox{0.6}{ min}}^{\scalebox{0.6}{ trans}}$ in the background hypothesis is given by

The likelihood for the background hypothesis is constructed by $\mathcal{L}_{\rm BG}^{\rm SN}(\scalebox{1.2}{\it z}_{\scalebox{0.6}{ min}}^{\scalebox{0.6}{ trans}})=\rho_{z_{\rm min}}^{\rm SN}(R^{\rm SN}_{0},\scalebox{1.2}{\it z}_{\scalebox{0.6}{ min}}^{\scalebox{0.6}{ trans}})$.

Figure 5 shows the probability distribution of $\scalebox{1.2}{\it z}_{\scalebox{0.6}{ min}}^{\scalebox{0.6}{ trans}}$ at a given redshift for the signal and the background hypothesis, respectively. Note that the lower bound of the allowed redshift for the signal hypothesis is determined by the average source distance interval $\textstyle{r=(3/(4\pi n_{0}))^{1/3}}$. The substantial difference between the two distributions provides a statistical power to calculate the p-values to reject the unrelated SN hypothesis for a given $\scalebox{1.2}{\it z}_{\scalebox{0.6}{ min}}^{\scalebox{0.6}{ trans}}$. The test statistic to reject the background hypothesis is constructed by

Table 1 lists the p-values for the background hypothesis. If we find the closest counterpart at redshift of 0.05, the statistical significance against the incorrect coincident SN detection hypothesis is $\sim 2.4\sigma$. Any distant counterpart with $z\gtrsim 0.15$ would exhibit a $\sim 1\sigma$ significance at most, and its association with the neutrino multiplet cannot be claimed.

The SNe rate may not be exactly a volume-averaged value. The inhomogeneous distribution of galaxies in our local universe impacts the expected value of the nearby SNe rate. For example, the local SFR suggests that the CC SNe rate within 10 Mpc is larger than the volume-averaged value, which enhances the detectability of high-energy neutrinos from nearby SNe (Kheirandish & Murase, 2022). The local anisotropic and inhomogeneous structure appears on a distance scale of $\lesssim$ 100 Mpc, governed by the local clusters of galaxies. In order to evaluate the robustness of the present statistical approach, we build an empirical local universe model to calculate the local SNe rate within 100 Mpc of our Galaxy.

The local SFR $R^{\rm SFR}_{0}$ is associated with the average B-band luminosity $\bar{L}_{\rm B}$ (Kennicutt, 1998). We assume $\textstyle{R^{\rm SN}_{0}\propto R^{\rm SFR}_{0}\propto\bar{L}_{\rm B}}$. The ratio of $R^{\rm SFR}_{0}$ and $\bar{L}_{\rm B}$ is approximately constant for various galaxy types (James et al., 2008):

$R^{\rm SN}_{0}$ in the local universe, $R^{\rm SN,local}_{0}$, is related to $R^{\rm SFR}_{0}$ via

Here $\rho_{\rm B}$ is the B-band luminosity density of the local universe, and $1/M_{\rm SN}$ corresponds to the probability of SNe per mass, which is given by

where $M_{{\rm min}}$ and $M_{{\rm min}}^{\rm SN}$ are the minimum masses for the stellar population and for SNe, respectively. We assume $M_{\rm min}=0.1M_{\odot}$ and $M_{\rm min}^{\rm SN}=8M_{\odot}$ by following Madau & Dickinson (2014). For the initial mass function, $dN_{s}/dM_{s}\propto M_{s}^{-2.35}$ is assumed.

We consider the B-band luminosity density $\rho_{\rm B}$ for a given patch of the local universe as a measure to quantify a departure from the cosmologically averaged star formation activity. It is calculated by

The galaxy number distribution $dn_{\rm gal}/dL_{\rm B}$ can be estimated by the actual observations of galaxies. By using the GLADE galaxy catalog (Dálya et al., 2018), we estimated $\rho_{\rm B}$ within 100 Mpc for various directions with $\Delta\Omega=1$ deg² by counting galaxies brighter than -18 mag which results in $\sim 60$ % catalog completeness in terms of the integrated B-band luminosity. Figure 6 shows the estimated $\rho_{\rm B}$ skymap in equatorial coordinates. The mean number was found to be $\rho_{\rm B}\simeq 10^{8}L_{\odot}$ Mpc^-3 while $\rho_{\rm B}\lesssim 10^{9}L_{\odot}$ Mpc^-3 in 98% of the sky patches. We take the latter value as the upper bound for conservative estimates of the p-values for the background SNe hypothesis.

When $10^{8}L_{\odot}\ {\rm Mpc}^{-3}\lesssim\rho_{\rm B}\lesssim 10^{9}L_{\odot}$ Mpc^-3, $R^{\rm SN,local}_{0}$ calculated by Eq. (25) ranges from $7\times 10^{-5}$ Mpc^-3 yr^-1 to $7\times 10^{-4}$ Mpc^-3 yr^-1. Table 1 shows the p-value for these two cases for the local universe, for the average and upper bounds of the galaxy density.

As described in Section 3.1, the multiplet source is most likely located at a relatively small redshift ($z<0.15$ with an 88% probability). Given that most astrophysical neutrino multiplets will therefore originate in objects with small redshifts, we can construct a follow-up search strategy for identifying the optical counterpart. Figure 7 shows a flow chart to illustrate our proposed procedure to find a source candidate. The first step in this procedure for the followup observations is to see if the optical transient counterparts are close enough. We propose $z=0.15$ as the threshold for the first level selection.

To detect the optical emission from transients at $z\lesssim 0.15$, the required sensitivity for optical follow-up observations is about 23 mag. Figure 8 shows the typical peak optical magnitude of different transients as a function of redshift. The sensitivity of 23 mag is sufficient to detect hypernovae and broad-lined, energetic Type Ic SNe at $z=0.15$. To perform such an optical survey for 1 deg² area, wide-field optical telescopes with a diameter of $>$ 4 m, such as DECam on the 4 m Blonco telescope (Flaugher et al., 2015), HSC on the 8.2 m Subaru telescope (Miyazaki et al., 2018), and the 8.4 m Rubin observatory telescope (Ivezić et al., 2019), are required.

Such a deep survey can identify not only the true multiplet source but also unrelated transient objects (i.e., contaminants). In the following sections, we first estimate the number of contamination sources that may be found with such a survey, and then discuss the strategy to identify the neutrino multiplet source.