{'fraction_non_alphanumeric': 0.03823461803693204, 'fraction_numerical': 0.017853038081373077, 'mean_word_length': 4.457764148313354, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
High-energy neutrinos from choked-jet supernovae: searches and implications

The origin of the high-energy astrophysical neutrinos discovered by IceCube remains largely unknown. Multi-messenger studies have indicated that the majority of these neutrinos come from gamma-ray-dark sources. Choked-jet supernovae (cjSNe), which are supernovae powered by relativistic jets stalled in stellar materials, may lead to neutrino emission via photohadronic interactions while the coproduced gamma rays are absorbed. In this paper, we perform an unbinned-maximum-likelihood analysis to search for correlations between IceCube's 10-year muon-track events and our SN Ib/c sample, collected from publicly available catalogs. In addition to the conventional power-law models, we also consider the impacts of more realistic neutrino emission models for the first time, and study the effects of the jet beaming factor in the analyses. Our results show no significant correlation. Even so, the conservative upper limits we set to the contribution of cjSNe to the diffuse astrophysical neutrino flux still allow SNe Ib/c to be the dominant source of astrophysical neutrinos observed by IceCube. We discuss implications to the cjSNe scenario from our results and the power of future neutrino and supernova observations.

The origin of the high-energy astrophysical neutrinos discovered by IceCube remains largely unknown. Multi-messenger studies have indicated that the majority of these neutrinos come from gamma-ray-dark sources. Choked-jet supernovae (cjSNe), which are supernovae powered by relativistic jets stalled in stellar materials, may lead to neutrino emission via photohadronic interactions while the coproduced gamma rays are absorbed. In this paper, we perform an unbinned-maximumlikelihood analysis to search for correlations between IceCube's 10-year muon-track events and our SN Ib/c sample, collected from publicly available catalogs. In addition to the conventional powerlaw models, we also consider the impacts of more realistic neutrino emission models for the first time, and study the effects of the jet beaming factor in the analyses. Our results show no significant correlation. Even so, the conservative upper limits we set to the contribution of cjSNe to the diffuse astrophysical neutrino flux still allow SNe Ib/c to be the dominant source of astrophysical neutrinos observed by IceCube. We discuss implications to the cjSNe scenario from our results and the power of future neutrino and supernova observations.

I. INTRODUCTION
The discovery of high-energy (HE) astrophysical neutrinos by the IceCube Neutrino Observatory [1,2] has opened a new era of neutrino physics, astrophysics, and multimessenger astronomy. These neutrinos are nearly isotropic on the sky and have energies between 10 TeV to above PeV, suggesting their sources to be the extragalactic populations of extreme cosmic accelerators. Searching for HE neutrino sources is also crucial to identify sources of their parent HE cosmic rays and it offers unique opportunities to understand the acceleration mechanisms of the sources.
Compelling evidence for a few point sources has been reported: blazar TXS 0506+056 [18,19], Seyfert II galaxy NGC 1068 [20], and TDE candidates AT2019dsg, AT2019fdr, and AT2019aalc [21][22][23]. However, stacking analyses show that none of the aforementioned types of sources contribute a major fraction of the all-sky neutrino flux (a possible exception might be the non-jetted AGN if the observed 2.6σ significance is interpreted as a sig-arXiv:2210.03088v1 [astro-ph.HE] 6 Oct 2022 nal [24]). Furthermore, the measurement of the extragalactic gamma-ray background by the Fermi Large Area Telescope (LAT) [25] has set robust constraints on gamma-ray-bright sources as dominant HE neutrino emitters, as a comparable flux of gamma rays are expected to be coproduced with neutrinos following pp and pγ interactions of cosmic rays [26]. Thus, it is likely that there is a class of HE neutrino sources opaque to GeV -TeV gamma rays, in which only neutrinos can freely escape from the cosmic-ray accelerators [27][28][29][30].
Choked-jet supernovae (cjSNe) are promising as such hidden neutrino sources [31][32][33][34][35][36]. Figure 1 shows an illustration: in this scenario, the jets are stalled or "choked" inside the progenitor envelopes or circumstellar materials as they are not powerful enough. Gamma rays produced from cosmic rays accelerated in the choked jet are attenuated through optically thick environments below the stellar photosphere, leaving HE neutrinos as primary signals. Previous theoretical studies have shown that cjSNe could even explain all the HE astrophysical neutrinos observed by IceCube [35,[37][38][39].
In addition, cjSN provides a unified scenario for Type Ib/c supernovae (SNe Ib/c), hypernovae, and the different types of GRBs, in which jet properties and shockbreakout conditions are crucial in making the difference [38,40]. The link to low-power (LP) GRBs such as low-luminosity (LL) GRBs and ultralong (UL) GRBs further strengthens the role of cjSNe as gamma-ray dark factories of HE neutrinos [35,38,39]. From the observational side, many LP GRBs are missed by current GRB surveys, so current stringent limits on HE neutrino emission from classical, high-luminosity GRBs do not apply. From the theoretical side, the intrinsically weak jets associated with LP GRBs are more ideal for neutrino production than classical GRBs, as powerful jets generally lead to inefficient cosmic-ray acceleration in radiationmediated shocks [35]. Therefore, HE neutrinos also provide us with an essential tool to study the cjSN scenario and the observed GRB-SN connection [41][42][43][44].
In this paper, we search for HE neutrinos from cjSNe and study the theoretical implications. Previous studies [45][46][47][48][49] have performed searches using early datasets of IceCube and found no association of neutrinos with supernovae. Here we use 10 years of IceCube neutrino data [50,51], in which high-quality events are recorded and have never been analyzed for cjSNe. We perform an unbinned-maximum-likelihood analysis to search for a statistical correlation between neutrinos and SNe Ib/c, as cjSNe can in principle be observed as SNe Ib/c, where progenitor stars are more massive and typically enclosed by denser extended materials. Our analyses do not find any excess of neutrinos from SNe Ib/c with respect to the background, from which we set upper limits on cjSN models and their contribution to the total astrophysical neutrino fluxes observed by IceCube. Figure 2 illustrates how muon-neutrino signals from all SNe Ib/c in our analyzed sample would look in IceCube, assuming the highest cjSN flux allowed by our analysis.  neutrino spectrum with parameters {ECR, fjet} (detailed in Sec. IV) defined by the upper limit from our analysis. In comparison, we also show the effective number of background events associated with the supernovae [within the effective size of time (20 days) and spatial window (9 deg 2 ); in our analyses, we use more events than shown here (Sec. III B)]. The sensitivity of our analysis is mainly driven by the neutrino data with energy above few tens of TeV, where the background becomes negligible compared to the signal.
For the first time, we take into account realistic models of neutrino emission in LP GRBs, instead of simply assuming HE neutrinos to follow power-law spectra. This is important, as the astrophysical neutrino flux may originate from multiple source populations with various neutrino spectra. Moreover, as the current LP GRB sample is highly incomplete and model uncertainties of LP GRBs are largely unconstrained, our survey in model parameters using SNe Ib/c catalogs is meaningful. Finally, our conservative limits show that, for most cjSN models we consider, SNe Ib/c can still account for 100% of the Ice-Cube diffuse neutrino flux. This paper is organized as follows. In Sec. II, we describe the neutrino dataset and the supernova sample we use and introduce the cjSN models we consider. In Sec. III, we discuss our likelihood formalism and how we obtain the correlation significance from background simulation. In Sec. IV, we detail the procedures of setting upper limits, including how we simulate neutrino signals from cjSNe and how we look for an excess of signals among background fluctuations. In Sec. V, we show our results from single-source and stacking analyses, and present constraints on cjSN models and their contributions to the HE astrophysical neutrino fluxes observed by IceCube. We further discuss the implications of cjSNe as the origin of HE neutrinos. In Sec. VI, we conclude our findings with a future roadmap. The IceCube Neutrino Observatory detects neutrinos through the Cherenkov photons emitted by relativistic charged particles produced from neutrino interactions in (starting events) and outside (throughgoing events) the detector [52,53]. The Cherenkov photons trigger the nearby digital optical modules and can form two kinds of basic event topologies: elongated tracks formed by muons and showers which looks like a round and big blob formed by electrons (electromagnetic shower) or hadrons (hadronic shower). The track events, which are dominated by throughgoing tracks, have a much better angular resolution (as good as < 1 • ) though worse energy resolution (∼ 200% at ∼ 100 TeV) than the shower events (∼ 10 • -15 • and ∼ 15% above 100 TeV) [54]. Thus, track events are suited to search for point sources.
The data released by the IceCube Collaboration span from April 2008 to July 2018 [50,51]. The same data have been used in the 10-year time-integrated neutrino point-source search by IceCube collaboration [55] and searching for high-energy neutrino emission from radio bright AGN [56]. In total, there are 1,134,450 muontrack events. The information for each track is provided, including arrival time, angular direction, angular uncertainty, and reconstructed energy. The arrival time is given in precision of 1 × 10 −8 days (8.6 × 10 −4 s). These ten years of data are grouped into five samples corresponding to different construction phases of IceCube and instrumental response functions, including 1) IC40, 2) IC59, 3) IC79, 4) IC86-I, and 5) IC86-II to IC86-VII. The numbers in the names represent the numbers of strings in the detector on which digital optical modules are deployed. Distributions of these events in the sky can be found in Figs. 1 and 2 in Ref. [56]. We use the events with declination (Dec) between −10 • and 90 • for the following reasons. First, the events from Dec < −10 • (southern sky with respect to IceCube) have much higher backgrounds from atmospheric muons [50]. Second, we find that the given smearing matrices from simulations have too low statistics to obtain good enough energy PDFs for our analysis (Sec. III B).
We also process the 19 × 2 double-counted tracks in the dataset found in Ref. [57] (listed in its table III). These events arise from an internal reconstruction error that identifies some single muons crossing the dust layer as two separate muons arriving at the same time and closely in direction [57]. This would affect neutrinosource searches, especially transients, as finding two associated events versus one would be quite different. Thus, we combine the 19 misreconstructed pairs into 19 single events by averaging the directions and summing up the reconstructed energies.

B. Supernova sample
The supernova sample we use for our analysis is from combining SNe Ib/c from the Open Supernova Catalog [58], the Weizmann Interactive Supernova Data Repository (WISeREP) [59], and the All-Sky Automated Survey for Supernovae (ASAS-SN) [60][61][62][63]. These catalogs have collected more than 36,000, 20,000, and 1,300 supernovae, respectively, from a variety of astronomical surveys and existing archives. We further compare our combined supernova sample with the publicly available catalog of bright supernovae [64,65] and incorporate those that are missed in the above.
Sometimes a supernova is independently discovered by different groups and thus has multiple aliases. This leads to a small fraction of potentially duplicate sources in our sample. To avoid double counting, we first search for the supernovae with an angular distance smaller than 0.1 • . We then merge these supernovae if they are classified as the same type and the difference in their maximal brightness time and redshift are less than 30 days and 10%, respectively. The examples of supernova pairs satisfying our criteria are {SN 2010O, SN 2010P} and {SN 2016coi, ASASSN-16fp}. As these potentially duplicate sources have very similar observational properties, we remove one of them from our sample.
Finally, we keep the supernovae in our sample only if they were observed at Dec ≥ −10 • and within the uptime of IceCube between April 2008 and July 2018 to match our selected data (Sec. II A). In total, our final sample consists of 386 SNe Ib/c, including 30, 36, 36, 36, and 248 for IC40, IC59, IC79, IC86-I, and IC86-II -VII, respectively. We provide the details of our supernova sample at this URL .

C. cjSN models for neutrino emission
We assume choked jets to be nearly calorimetric sources (except the suppression factor) so that neutrinos are produced by all available energy E CR in cosmic rays. The all-flavor neutrino spectrum from a single burst of supernova is thus given by [67,68] where the factor 3/8 is the fraction of energy taken away by neutrinos from charged pions produced in the pγ interactions; the energy-dependent suppression factor f sup (ε p ) accounts for the meson and muon cooling processes that depend on the detailed modeling of choked jets; and the meson production efficiency min[1, f pγ ] is set to 1 for the choked jets as long as the minimal cosmic-ray energy is larger than the pion production threshold. As the average fraction of energy transferred from a parent proton to each neutrino after an interaction is 1/20, A comparison of single-burst neutrino spectra from the cjSN models we consider. Here, the LLGRB-PE model involves prompt neutrino emission from low-luminosity GRBs [66], and the ULGRB-CS model involves neutrinos from cosmic rays accelerated at the collimation shock of ultra-long-duration GRBs [35]. All spectra are normalized to have the same isotropic equivalent cosmic ray energy, ECR = 10 51 erg. Both panels use the same line style for each model while in the right panel each spectrum is weighted by the IceCube effective area (averaged over Dec ≥ −10 • ). The realistic modeling of cjSNe takes into account a variety of cooling mechanisms, leading to the suppression in the neutrino spectra and different energy distributions of neutrino events detected by IceCube.
we have ε ν ≈ (1/20)ε p . R p denotes the bolometric correction factor for the cosmic-ray spectrum. The cosmic rays are expected to follow a power-law spectrum [69] (i.e., dN p /dε p ∝ ε −γ p ) as they are typically accelerated through the first-order Fermi process [70] in the shock. In this case, we have R p (ε p ) = ln (ε max p /ε min p ) for γ = 2 and When producing our results in Sections V B and V C, we take (ε min p , ε max p ) = (2 × 10 3 , 2 × 10 10 ) GeV as this leads to efficient neutrino production within (ε min ν , ε max ν ) ≈ (10 2 , 10 9 ) GeV, the energy range that could be detected by IceCube.
We also consider two well-motivated cjSN models [35,66]. Both models assume power-law parent cosmicray spectra with γ = 2.0 but the neutrino spectra are different due to various processes of meson and muon cooling in different shocked regions of the jet. In reality, f pγ in Eq. (1) is below the unity at low energies and that f sup < 1 is possible depending on cjSN parameters.
Our first class of models assumes the neutrino spectrum to follow a power law with the same spectral index (γ) as the parent proton spectrum. This neglects all complicated mechanisms that could lead to a nontrivial neutrino spectrum. We consider γ = 2.0, γ = 2.5, and γ = 3.0 as three benchmark models, which match the best fits to the 10 years of muon-track events [71], a combination of track and shower events [72,73], and the 7.5 years of the high-energy starting events (HESE) [74], respectively.
Our second class of models takes into account more realistic modeling that link cjSNe with LL GRBs and UL GRBs. We consider two realistic models from different detailed considerations of jet propagation inside the progenitor star and energy losses of particles: 1) the LLGRB-PE model (prompt ν emission from LL GRBs) [66], and 2) the ULGRB-CS model (neutrinos from cosmic rays accelerated at the collimation shock of UL GRBs) [35]. Figure 3 left panel shows that for the same E CR , the all-flavor neutrino spectra of a supernova burst from different models can differ by several orders of magnitude at certain energies. The y-axis is shown as ε 2 ν dN ν /dε ν = 2.3 −1 ε ν dN ν /d log ε ν so that the area under each curve is proportional to the total energy of neutrinos. Figure 3 right panel shows the spectra of detectable muon neutrinos in IceCube, calculated from the curves in the left panel and the average IceCube effective areaĀ eff over Dec ≥ −10 • . For comparison, here we assume that neutrinos are evenly distributed among all flavors. The yaxis is shown as ε νĀeff dN ν /dε ν = 2.3 −1Ā eff dN ν /d log ε ν so that the area under each curve is proportional to the total number of detected neutrinos. For the same E CR , power-law models (in which there are no cooling effects) or realistic models with harder emission spectra between 10 TeV and 1 PeV (where IceCube has the best sensitiv-ity) tend to produce more neutrino events in detectors.

A. Likelihood function and test statistic
To search for a possible correlation between IceCube events and supernovae, we use an unbinned maximumlikelihood method. A similar formalism was commonly used to search for transient neutrino sources by Ice-Cube [75][76][77][78][79][80][81][82][83] and other experiments [84]. The likelihood function is defined as where k labels the five data samples we use in our analysis, N k SN number of supernovae observed during the period of each data sample, and N j is the total number of muon-track events appearing in the temporal and spatial window of the jth supernova. Moreover, s j and b j are the signal and background event rates for the jth supernova, respectively. In our analysis, s j is the parameter to be determined by maximizing the likelihood function, and b j is computed by scaling the number of track events in the supernova's spatial window and outside the temporal window by the ratio between the total time outside and inside its time window. Details of the windows are in Sec. III B.
The L k ij (s j , b j ) in Eq. (2) is the individual likelihood function of a track event i from the kth data sample associated with the jth supernova, with S k ij the signal probability density function (PDF) and B k ij the background PDF, details in Sec. III B. The p j (s j , b j ) in Eq. (2) is a weighting factor, to take into account the Poisson fluctuations of N j . The test statistic (TS) is then defined as the ratio of the likelihood function to its value under the null hypothesis From Eqs.
(2) to (4) we have where TS k j is the test statistic for a single source covered by the kth data sample. Comparing the TS values from the real with simulated data (details in Sec. III C), we can get the significance (or p-value) of the correlation between the supernova sample and the real data.

B. Signal and background PDFs
Both signal and background PDFs for a track event are multiplications of temporal, spatial, and energy PDFs, These PDFs describe the expected temporal, spatial, and energy distributions of the events if they are from the sources or the background. For the background PDFs, we calculate them directly from data as they have enough statistics.
The signal temporal PDF is taken as a Gaussian distribution in T SN − T ν , the difference between the time of a supernova's maximum brightness and the arrival time of an event, where σ T = 4 days and λ T = 13 days given below.
In the standard scenarios of cjSNe, high-energy neutrinos and gamma rays are produced simultaneously via the interaction of relativistic jet inside the extended stellar envelope. As it is found that SNe Ib/c correlated with GRBs typically reach their maximal optical brightness after 13 days from the prompt gamma-ray emission [85], we take λ T = 13 days in Eq. (10).
When searching for temporal correlation, the time window of each supernova is determined by the 99% confidence interval of Eq. (10).
This gives T SNj − T i ∈ [λ T − 2.58σ T , λ T + 2.58σ T ], the width of which is 20 days. The events outside the window are almost impossible to come from the supernova.
For the background temporal PDF, we use a uniform distribution inside the time window and zero outside, because the background event rate from high-energy atmospheric neutrinos and muon interactions are constant. The signal spatial PDF is given by the Fisher-Bingham (Kent) distribution [86][87][88], where µ ij = cos(∆ψ ij ), with ∆ψ ij the angular distance between the track i and supernova j, and σ νi is the reconstructed angular error of a track [50,51]. Equation (12) can be understood as a generalization of the 2D Gaussian distribution with the standard deviation σ νi on a sphere. For small ∆ψ ij , it reduces to a 2D Gaussian distribution. Similar to the time window, we search for spatial correlation through the spatial window defined by the 99% confidence interval of Eq. (12), which gives µ ij ≥ µ 99 (σ νi ) where The background spatial PDF is a function of Dec, δ i , as distribution with respect to right ascension (RA) is nearly isotropic because IceCube locates at the South Pole. We use the same sliding-window method as in Ref. [56] to calculate the background spatial PDF. We normalize the background spatial PDF over a Dec from −10 • to 90 • , matching the data used in our analysis. The results can be found in Fig. 2 of Ref. [56].
The signal energy PDF of an event in sample k can be calculated by where φ j (E ν ) is the expected muon-neutrino fluence from source j (Sec. IV A), A k eff (E ν |δ i ) the effective area, and P k is the reconstructed muon-energy (energy-proxy; E prx i ) distribution for a specific E ν and δ i , given in Refs. [50,51].
The background energy PDF, B k E, ij (E prx i |δ i ), is a normalized distribution of E prx i as a function of δ i and is obtained using a sliding-window method with the window size of sin δ i ± 0.05 and log 10 (E prx i ) ± 0.2.

C. Significance and simulation
To get the p-value (or the statistical significance) for the correlation between the events and supernovae, we need to simulate a large number of synthetic neutrino datasets and calculate the TS. The p-value quantifies the probability that a correlation is observed due to background alone. In our analysis, it is defined as the fraction of the TS of simulated datasets that are larger than the TS of the real data (TS obs ) and it can be converted to significance under the standard normal distribution.
The simulation is detailed as follows. We simulate 10 5 synthetic datasets so that the results converge. For each dataset, it has the same fives phases and for each phase, we simulate the same number of events as the real data. For a simulated event in a specific phase of a specific dataset: • The arrival time and RA are randomly chosen from uniform distributions over the uptime of the phase and [0, 2π), respectively.
• The Dec and the energy proxy are randomly chosen according to the background spatial and energy distributions of the real dataset, given in Sec. III B.
• Given the Dec and energy proxy of a simulated event, the angular error is drawn from the distribution of the angular errors of the events with similar Dec and proxy energy in the real dataset.
• We keep only the events with δ i −10 • , matching our selection for the real data.
The signal and background PDFs of the events in the simulated datasets are calculated in the same way as those in the real dataset.

IV. SETTING UPPER LIMITS
In this section, we detail the procedures of setting upper limits on the parameters of cjSN models. All the models are characterized by the same parameters, the isotropic equivalent cosmic ray energy, E CR , and the fraction of supernovae that may contribute signal neutrinos, f jet . The former represents the total energy budget of which cosmic rays can efficiently produce signal neutrinos in the energy range of interest, and the latter accounts for the uncertainties in the fraction of supernovae developing choked jets towards the Earth. The constraints for each model are set from signal-injected simulations, detailed below.
First, we randomly select f jet N k SN supernovae from our sample and calculate their neutrino spectra for a specific model in Sec. II C. After convolving with the effective area, we get the detectable neutrino spectra, from which we simulate the neutrino events. Second, we simulate each muon-track event from each parent neutrino event using the smearing matrix [50,51], which takes into account the effects of neutrino interactions, muon energy losses, and detector efficiency. We detail the calculation of the number of detectable signal events from individual supernovae in Sec. IV A, and the simulations of neutrino events and track events in Sec. IV B. In every signal-injected simulation, which corresponds to a specific {E CR , f jet } and specific model, the above procedure is done for each supernova in our sample and we inject the simulated signal track events into a randomly selected simulated neutrino dataset (Sec. III C). Next, we calculate the TS value for each signal-injected simulation following the procedure in the previous section and compare them with the TS-values distribution of the simulated datasets (Sec. III C).
The above procedure is repeated for different {E CR , f jet } and models. For each model, to produce our limits, we scan over {E CR , f jet } parameter space, in which the TS-value distribution of each point is obtained by performing 10 4 simulations. We have checked that our results converge well. The sensitivity of our analysis to a model strongly depends on its total signal-neutrino flux. For the scenarios of cjSNe, this is closely related to E CR , as the neutrinos are mainly produced through the charged pion decay following the meson production from the pγ interactions of accelerated protons in the jet [67]. For the kth data sample, we calculate the number of detected signal neutrino events by the sum of contributions from the individual supernovae, which is where, with the flux equally distributed in flavors, the muon-neutrino fluence from a single supernova burst at luminosity distance d j (redshift z j ) is related to Eq. (1) by the observed neutrino energy is E ν = (1 + z j ) −1 ε ν , and the Dec of the signal neutrinos is the same as that of their parent supernova, δ SNj . E min ν and E max ν are determined by the energy range of the accelerated protons at the source.

B. Simulate detected track events
The detected track events from the neutrino events are simulated following the following procedure.
• The arrival time T i is randomly chosen from the signal temporal PDF in Eq. (10).
• The Dec of the neutrino event, δ i , is set to be that of the supernova, δ SNj .
• The energy of a detected signal neutrino arriving at Earth, E det νi , is randomly drawn from the energy spectrum of detected neutrinos; that is, (15). • Given the E det νi and δ i of a neutrino event, the proxy energy E prx i and angular uncertainty of its daughter track event are randomly drawn following the smearing matrix.
• Given the angular uncertainty, the angular distance, RA, and Dec of the track are drawn following the signal spatial PDF in Eq. (12).

V. RESULTS AND DISCUSSION
In this section, we present our results. In Sec. V A, we show our single-source and stacking analysis results. Then we show the upper limits on the choked-jet model parameters in Sec. V B and their contributions to the diffuse high-energy neutrino flux in Sec. V C.

A. Single-source and stacking analysis results
In our single-source analysis, the pre-trial p-value for individual supernovae (p k pre, j ) is determined by the probability of finding TS k bkg, j ≥ TS k obs, j in the test statistic distribution under background-only hypothesis. However, a low pre-trial p-value does not necessarily mean a discovery but could be caused by the statistical fluctuations of the background. This can be taken into account by the post-trial p-value, which can be calculated by that none of the supernovae in our sample have significant neutrino emission. The highest pre-trial significance is 2.9σ (SN 2009hy), but its post-trial significance is only 1.5σ. Here we only report the p-values based on the power-law spectrum with γ = 2.0, as we find that changing the neutrino spectrum only marginally affects the significance of each source. We provide the p-values for each supernova given all cjSN models at this URL . Table I shows our stacking analysis results. The posttrial p-values are calculated by p post = 1 − (1 − p pre ) 5 , where 5 is the total number of spectra/models we test. The results show that for all the models we consider, the data are consistent with the background-only hypothesis.  Also shown in the upper right corner is the existing 90% confidence level limit of γ = 2.0 obtained by previous work [47] using one year of track events. Realistic models are generally allowed to have larger total energy in cosmic rays, as they are less efficient in producing detectable neutrinos.
eter space above the curves is excluded. Among all the scenarios we consider, the power-law model with γ = 2.0 have the strongest constraints, because it produces the most signal events in the energy range to which IceCube detectors are most sensitive. For the same reason, powerlaw models with softer spectra are less constrained, and the limit on the ULGRB-CS model is relatively weak, as it produces the least amount of neutrinos in that energy range.
Our limit on E CR for the power-law model with γ = 2.0 is stronger than the existing one-year limit in Ref. [47] by more than an order of magnitude at f jet = 1. The improvement is expected because our supernovae sample covered by the 10-year IceCube dataset is more than ten times larger than that used in Ref. [47].
Theoretically, E CR and to f jet are completely degenerate. In that case, the contours should all follow ∼ 1/E CR . However, there is a small deviation at small f jet . This is due to limited number of supernovae in our sample, and for small f jet the uncertainty of our limits on E CR is larger. Also, it should be noted that our limits are mainly driven by nearby supernovae. Thus, the degeneracy holds when k f jet (N k SN · f nb ) > 1, where f nb is the fraction of nearby sources in our samples. During our signal-inject procedure, we find that more than 50% of the signal events are contributed by the supernovae with luminosity distance d j 16 Mpc, which account for 2% of all the supernovae in our sample, so we set f nb = 0.02. With k N k SN = 386 we expect that the limits are obtained for f jet 0.1. This implies that with the present supernova sample the limits on cjSNe models are rather meaningful only if the jets have relatively wide opening angles, θ j 0.3 rad. To test cjSN models with smaller opening angles, e.g., θ j ∼ 0.1 rad, better samples with more nearby supernovae are required.

C. Upper limits on diffuse neutrino flux
Now we convert our limits on E CR and f jet to upper limits on the contribution to diffuse neutrino flux from all the cjSN models. The diffuse flux can be calculated by (e.g., [27]) where ε ν ≡ (1 + z)E ν , dρ(z)/dL γ is the luminosity function that takes into account the distribution of on-axis choked-jet supernova rate density ρ(z) = f jet R(z) with respect to the observed GRB luminosity, H 0 is the Hubble constant, and {Ω Λ , Ω m } are the cosmological energy density parameters. Combining the above equation with Eq. (1) gives the all-flavor diffuse neutrino flux by [27,89] where ξ z and R L ≡ R(0) are the redshift evolution factor [67,89] and the true local rate density (i.e., volumetric rate) of the supernovae that harbor choked jets (including off-axis ones), respectively. Eq. (19) shows that an upper limit on f jet E CR directly constrains the maximal neutrino flux we can detect. Figure 7 shows the upper limits (at 90% confidence level) on the all-flavor diffuse neutrino flux from SNe Ib/c for different cjSN models. The bounds are set via the E CR limits for f jet = 1 from Fig. 6. (Note that our limits are sensible for f jet 0.1, as discussed below.) To be conservative and consistent with our supernova sample, we take R L = R Ib/c ∼ 2 × 10 4 Gpc −3 yr −1 , the local rate density of SNe Ib/c [90,91]. We find that most of the models can still explain all the astrophysical neutrino fluxes. The harder the spectrum, the stronger constraints we get. Although we focus on cjSN models, SNe Ib/c may be aided by newborn pulsar winds, which could be emitters of neutrinos in the PeV -EeV range [92][93][94]. Our analysis results can be applied to such models but the limits are still above the model predictions.
Our limits could potentially be improved by selecting a subsample of nearby sources that dominates the total neutrino flux. This would reduce the total background in the analysis and lead to better sensitivity. For our main results, we do not select any subsamples in the analyses, so that our limits are regarded as robust and 'conservative' constraints on cjSN models; also, we caution that the stacking analysis with a small subsample will greatly decrease the sensitivity for f jet < 1. Even if all supernovae harbor jets, only a fraction (f jet ) of them point to us. If none of the on-axis events are included in the nearby sample, sensible constraints cannot be placed as the viewing angle of choked jets is unknown.
It is worth noting that, although we constrain realistic cjSN models motivated by various types of rarer types of supernovae or GRBs, one should not simply downscale our diffuse flux limits by applying the local rate density of rarer source classes to Eq. (19). This is because our SN Ib/c sample does not necessarily include any GRBs or hypernovae. If hypernovae or GRBs (including LL GRBs) only compose a small fraction of the entire supernova sample (i.e., R L < R Ib/c ), our approach described in Sec. (V B) may not give a sensible limit. This is the case even if all SNe Ib/c harbor jets (i.e., R L = R Ib/c ) because only a fraction of them have jets pointing to us. In the current analysis with ∼ 400 supernovae, we can apply the limits down to f jet ∼ 0.1 and more supernova samples are necessary to probe cjSN models with smaller  [72,73]. The data points with error bars and the red band represent the best-fit results from the 2010 -2015 shower data of IceCube [73]. Because of the high local rate density of supernovae, most cjSN models we consider remain viable to explain the diffuse flux of astrophysical neutrinos measured by IceCube.
jet opening angles. This point is important not to misinterpret the results of the stacking analyses. The choice of the source rate density in the stacking analysis presented here must always be consistent with the sample of sources one takes into account.
Previous studies indicated that most astrophysical neutrino sources should be optically thick to gamma rays [27]. Among promising gamma-ray dark candidates, core-collapse supernovae with jets choked by surrounding materials are especially interesting [35,37,38], as they 1) naturally explain the potential connections between SNe Ib/c, hypernovae, and classical, LL, and UL GRBs, and 2) are efficient in producing HE astrophysical neutrinos.
In this paper, we studied whether cjSNe can be the dominant source contributing to the diffuse astrophysical neutrinos detected by IceCube. We used the unbinnedmaximum-likelihood method to search for the associations between IceCube neutrinos and SNe Ib/c. Compared to existing searches, our analysis takes advantage of 10 years of IceCube neutrino data. We also collected SNe Ib/c from publicly available supernova catalogs covering the same period. Importantly, for the first time, we looked for the neutrino signals of realistic cjSN models in IceCube data. These models take into account a variety of time-dependent cooling processes for cosmic rays and mesons in LP GRBs, which are well-motivated but have never been constrained by data.
For all the cjSN models we considered, our singlesource analysis (Fig. 5) and stacking analysis (Table I) found no significant correlation between the muon-track events and any SNe Ib/c in our source sample with respect to the backgrounds. Thus, we put upper limits on, E CR , the total energy budget of cosmic rays and, f jet , the fraction of SNe Ib/c that has choked jets pointing to us for different cjSN models. Our 10-year limit improves the prior 1-year limit in Ref. [47] by more than an order of magnitude (Fig. 6). Moreover, the limits of both parameters are generally sensitive to the efficiency of a cjSN model in producing detectable neutrinos by IceCube (see the right panel of Fig. 3), which can vary by a factor up to ∼ 5 among the models we consider.
Future prospects of testing cjSNe as neutrino sources are very promising. From Fig. 7, we clearly see that most of our conservative limits obtained from 10 years of data are close to the measured astrophysical neutrino flux in IceCube. Given that IceCube-Gen2 [129] will be at least 5 times more sensitive than IceCube, IceCube-Gen2 will effectively strengthen our stacking limits and critically constrain most cjSN models with just 2 years of operation. Additionally, there will be much more high-quality data from other observatories such as KM3NeT [130], Baikal-GVD [131], P-ONE [132], and TRIDENT [133], which are sensitive to different parts of the sky. Importantly, the Vera C. Rubin Observatory's upcoming Legacy Survey of Space and Time (LSST) [134,135] will increase the detection rate of supernovae by more than an order of magnitude (> 10 5 core-collapse supernovae per year [136,137]). The improvement of the redshift completeness of nearby supernova samples at low redshift (z 10 −2 ) will dramatically boost the sensitivities of future stacking analyses; furthermore, an increase of the observed supernovae at higher redshift will be beneficial to constraining the cjSN models with f jet 0.1. This is crucial because the typical beaming factor of GRB jets is believed to be f jet ∼ 0.01. Together, core-collapse supernovae with choked jets have reached the fork of being discovered or ruled out as dominant neutrino sources. Our insights into the interrelationship between the most elusive high-energy particles and the stellar explosions within dense environments will drastically grow in the near future.