Search for Continuous and Transient Neutrino Emission Associated with IceCube's Highest-energy Tracks An 11 yr Analysis

IceCube alert events are neutrinos with a moderate-to-high probability of having astrophysical origin. In this study, we analyze 11 yr of IceCube data and investigate 122 alert events and a selection of high-energy tracks detected between 2009 and the end of 2021. This high-energy event selection ( alert events + high-energy tracks ) has an average probability of 0.5 of being of astrophysical origin. We search for additional continuous and transient neutrino emission within the high-energy events ’ error regions. We ﬁ nd no evidence for signi ﬁ cant continuous neutrino emission from any of the alert event directions. The only locally signi ﬁ cant neutrino emission is the transient emission associated with the blazar TXS 0506 + 056, with a local signi ﬁ cance of 3 σ , which con ﬁ rms previous IceCube studies. When correcting for 122 test positions, the global p -value is 0.156 and compatible with the background hypothesis. We constrain the total continuous ﬂ ux emitted from all 122 test positions at 100 TeV to be below 1.2 × 10 − 15 ( TeV cm 2 s ) − 1 at 90% con ﬁ dence assuming an E − 2 spectrum. This corresponds to 4.5% of IceCube ’ s astrophysical diffuse ﬂ ux. Overall, we ﬁ nd no indication that alert events in general are linked to lower-energetic continuous or transient neutrino emission.


INTRODUCTION
The IceCube Neutrino Observatory (Aartsen et al. 2017a) is a Cherenkov detector using a cubic kilometer of Antarctic ice at the geographic South Pole to primarily (but not exclusively) study high-energy astrophysical neutrinos.Its duty cycle is greater than 99% (Aartsen et al. 2017a), and its field of view covers the full sky while being most sensitive to high-energy neutrino events near the celestial equator.This makes IceCube ideal for surveying the sky (Aartsen et al. 2017b).As part of the realtime program, IceCube alerts other telescopes upon detection of a neutrino event with a high probability of being of astrophysical origin, which can then trigger follow-up observations (Aartsen et al. 2017b;Blaufuss et al. 2019;Kintscher 2016).
On the 22nd of September 2017, IceCube detected a neutrino of likely astrophysical origin (IceCube-170922A 1 ).This triggered multi-wavelength follow-up observations, which detected a flaring blazar (TXS 0506+056) at the reconstructed origin direction of IceCube-170922A (Aartsen et al. 2018a).This correlation is significant at a 3σ level (Aartsen et al. 2018a).Additionally, a neutrino flare was identified originating from the same direction between September 2014 and March 2015 with a significance of 3.5σ (Aartsen et al. 2018b).
This detection demonstrates that IceCube alerts can point to neutrino source candidates due to their high probability of being of astrophysical origin, and we aim to investigate the origin directions of other IceCube alerts.A preliminary search showed no indication of continuous neutrino emission (Karl 2019).However, the IceCube alert criteria have since been updated (Blaufuss et al. 2019;Abbasi et al. 2023a).The IceCube data have also been reprocessed with improved calibration of the optical sensors (Abbasi et al. 2021a;Aartsen et al. 2020).This leads to improved energy and direction reconstruction compared to previous results in Aartsen et al. (2018b); Abbasi et al. (2021b); Aartsen et al. (2020Aartsen et al. ( , 2015)); Abbasi et al. (2021d,c).A first analysis benefiting from this new processing (Abbasi et al. 2022a) detected the neutrino signal from the Seyfert II galaxy NGC 1068 with a significance of 4.2σ (compared to 2.9σ in Aartsen et al. (2020)).A large part of the increase (0.9σ) is due to improved data processing and calibration.More details about effects on data are discussed in appendix B and in the supplementary material of Abbasi et al. (2022a).
In this work, we analyze 11 years of reprocessed IceCube data (through-going muon tracks, see Table 1) and search for an excess of neutrino-induced muons.We apply a conservative lower limit on the angular uncertainty of 0.2 degrees, whereas the median angular resolution is 0.57 degrees (compared to a median angular resolution of 0.59 degrees before the reprocessing).We identify possible neutrino production sites by looking at the origin of high-energy neutrinos that have a high probability of being of astrophysical origin.IceCube's highest energy neutrinos with the largest astrophysical purity are events from the new selection of IceCube alerts published in the so-called "gold" alert channel (Blaufuss et al. 2019;Abbasi et al. 2023a).Additionally, we extend the list by including 18 high-energy events from Abbasi et al. (2022b) that were confirmed to be likely astrophysical events by a new event classifier (Kronmueller & Glauch 2019).Since we use a combination of IceCube alert events and high-energy tracks identified retrospectively, we will refer to our event selection as "alert+ events" for brevity.All IceCube data used in this work (lower and high-energy events) have been reprocessed.
In this work, we excluded alert+ events within 30 degrees of the geographic poles (affecting three events), for which we have smaller statistics for the background.Other IceCube analyses have applied different declination cuts (for example including all events up to 81 • declination (Abbasi et al. 2022a) or up to 82 • (Aartsen et al. 2020)).Additionally, we removed alert+ events with large uncertainties (≥ 100 square degrees, affecting two events).As a result, our final sample consists of 122 high-energy events (104 IceCube alert events and 18 high-energy tracks, listed in Table 3), detected between 2009 and the end of 2021.On average, our selected alert+ events have a probability to be astrophysical of ≳ 0.5.The probability to be astrophysical is spectrum dependent and based on the muon neutrino spectrum measured by IceCube (Haack & Wiebusch 2017;Abbasi et al. 2022b).The median angular resolution (90% Figure 1.Sky map in right ascension and declination (epoch=J2000) with the arrival directions of events fulfilling the IceCube alert criteria (with the highest probability to be of astrophysical origin) we investigate in this work.The events were detected between August 2009 and the end of 2021.The shaded regions represent the 90% uncertainty region of the reconstruction.
uncertainty regions) of alert+ events is 2.1 degrees.In Figure 1, we show a map of all arrival directions and their 90% uncertainty regions of IceCube alert+ events investigated in this work.These events provide positions of interest analogous to a catalog of possible neutrino sources.Since IceCube alert+ events trigger this analysis, we remove the respective alert+ event from the 11 years of IceCube data when running the analysis.We present the analysis method in Section 2 and the results in Section 3.
Table 1.Overview of the improved and reprocessed data samples in this analysis.The columns list the configuration of the detector ("IC" and the number of deployed strings), the uptime (livetime) of the detector in days, the number of events in each sample, and the start and end dates of the data subset.

Year
Livetime We use an unbinned likelihood approach as presented by Braun et al. (2008).In this work, we investigate two source types: continuous sources and transient sources.We compare two hypotheses (each with a set of parameters ⃗ θ) • Background hypothesis H 0 ( ⃗ θ 0 ): The background comprises atmospheric neutrinos, atmospheric muons (remaining after event selection cuts), and diffuse astrophysical neutrinos.The flux is uniform in time and right ascension.
• Signal hypothesis H 1 ( ⃗ θ 1 ): There is a signal component additional to the atmospheric background and the average diffuse astrophysical neutrino emission.The signal neutrinos cluster around their source (subscript S) at right ascension, declination ⃗ x S = (α S , δ S ).The energy spectrum of the emitted flux is an unbroken power law: dϕ dEν ∝ E −γ ν .In the specific case of a transient source hypothesis (see Section 2.2), the neutrino emission has a Gaussian time profile with mean µ T and width σ T .
We remove the high-energy alert+ events that triggered this analysis from the data set.Hence, we look for additional neutrino emission from the direction of the high-energy alert+ events.We then maximize the likelihood, L, and compute the likelihood ratio . (1) The likelihood maximization varies the expectation value of the number of detected signal neutrinos, n S , and the emitted energy spectral index, γ.We allow values for γ between 1.5 and 4. For the background hypothesis, n S is fixed to 0. The likelihood is the probability density of observing the data given a specific hypothesis.The probability density of observing an event, i, is a sum of its probability to be signal, S i , or background, B i : where N is the total number of detected events (signal and background combined).
We define the test statistic, TS, as for a signal hypothesis with the best-fit value of nS neutrinos (the "ˆ" denoting the best-fit of a parameter) as the mean number of neutrinos we expect to detect from the neutrino source.
The investigated source candidates have directional uncertainties (see Figure 1).However, we assume potential sources are smaller than the best resolution of 0.2 degrees in our data (TXS 0506+056 has an angular size of ∼ 2.6 arc-seconds).Hence, we fit the best point-source position within a reconstructed 90% uncertainty region by dividing the region in a grid with steps of 0.2 • in right ascension and declination, the best angular uncertainty for events used in this study.The likelihood is optimized at each grid point.The grid point yielding the best result (i.e., the highest TS value) is subsequently considered the point-source position.
This procedure is run on different realizations of background data ∼ 10 4 times.The background data are all 11 years of muon tracks with randomly assigned right ascensions.In the final step, we calculate the test statistic value, TS data , for the true data and compare this with the simulated background test statistic distribution.The local p-value is the probability of getting this TS data (or a larger value) from a random background realization.This procedure is repeated for all remaining regions in the sky, yielding 122 local p-values.From these 122 values, we take the most significant local p-value, p 0 , to identify the most significant source.
As a next step, we correct the significance for having tested 122 regions in the sky.Considering only background realizations, we take the most significant p-value out of 122 positions for each realization and generate a distribution of best local p-values, p 0,BG .The final global p-value of our analysis is the probability for p 0,BG to be at least as significant as the p-value we got from our real data, p 0 .Since we are investigating only a limited number of points (122), weaker neutrino emissions have a higher significance in this analysis than in an all-sky scan, for example, in Abbasi et al. (2022a).
When testing the method with Monte Carlo simulations, the best-fit number of signal neutrinos, nS , and source spectral index, γ, show a bias compared to the true simulated source properties.For sources with simulated hard spectral indices (i.e., γ = 2), there is a tendency to fit slightly softer spectra and a slightly larger number of signal neutrinos.For simulated sources following softer spectral indices (i.e., γ = 3), the tendency is reversed to fitting slightly harder spectral indices and smaller numbers of signal neutrinos.Appendix A presents a more in-depth discussion of this bias.Correcting the bias is not straightforward, and we have decided not to include an, at best, incomplete correction.Hence, the best-fit fluxes are only indicative.This bias does not affect the flux limits since they are based on simulated fluxes where the true source strength is known.

Time-integrated search for continuous sources
We define the signal and background probability density functions (pdfs) S i and B i in a spatial and an energy part (see, e.g., Braun et al. (2008); Abbasi et al. (2011)).The spatial part depends on the source position ⃗ x S and the reconstructed event properties: reconstructed origin ⃗ x i and the angular uncertainty of the reconstructed origin σ i .The energy part depends on the reconstructed muon energy, E i , the reconstructed origin declination, δ i , and the source energy spectral index, γ.The signal pdf for a steady source is hence The energy pdf, S energy , is the probability of detecting a neutrino with reconstructed energy, E i , at declination, δ i , assuming the source emits neutrinos with a spectrum of E −γ .The background pdfs, B i , are defined similarly The spatial term depends only on the event declination, δ i .We assume uniformity in right ascension for the background data due to IceCube's unique position at the South Pole.B energy is derived directly from experimental data.
Searching for neutrino counterparts of the alert+ events, we want to be sensitive to a single strong emission from one (or a few) sources and, additionally, to faint emissions from a larger number of sources.Hence, our search for continuous sources consists of two parts.The first part searches for single strong neutrino emitters.The second part investigates the overall neutrino emission from all 122 positions of interest.In the latter case, we combine the neutrino emission and define a new test statistic value, TS stacked , by summing the test statistic values of all alert+ positions, k, (5) We take the TS k from the individual search, hence we do not correct for overlapping uncertainty regions of alert+ events.

Transient sources
For transient sources, we multiply a temporal pdf with the previously defined spatial and energy pdfs in equations ( 3) and (4) (Braun et al. 2010).We assume a Gaussian-shaped time profile centered around µ T with width σ T for the signal part.The temporal signal pdf becomes with t i as the time the event was detected.The background expectation is a constant rate over the whole data taking time, t data : The search for time-dependent sources adds another optimization step for the best flaring time.This introduces a bias towards shorter flares since the number of possible shorter flares is larger than the number of possible longer flares.We correct for this effect by multiplying the test statistic by a marginalization factor, √ 2πσ T 300 days (Braun et al. 2010).Here, 300 days is the maximal flaring time.Longer time scales would result in worse sensitivity than the time-integrated search.We assume a minimal σ T of 5 days to ensure the background uniformity in right ascension.
Conventional methods to find neutrino flares as in Aartsen et al. (2015); Abbasi et al. (2021d); Aartsen et al. (2018b); Abbasi et al. (2021c) apply a brute-force scan of all possible time-intervals between events where the ratio of Equation (3) over Equation (4) exceeds a certain threshold.This is computationally expensive.The computational cost can be reduced by increasing the required threshold and hence reducing the possible number of intervals scanned.We want to include as few biases as possible, and if following conventional approaches, we would apply the same threshold as in Aartsen et al. (2018b), where the ratio had to be ≥ 1.However, Aartsen et al. (2018b) performs this search only on one position in the sky.In our case, this would mean scanning the uncertainty region of 122 alert+ events in steps of 0.2 degrees and, at each step, evaluating every possible time window between 5 and 300 days in 11 years for neutrino emission.This proved to be computationally unfeasible.To overcome this problem, we investigated new approaches (Karl et al. 2021;Karl 2022;Karl et al. 2023), which do not rely on thresholds, such as a different test statistic to evaluate if an emission is time-dependent (Eller & Shtembari 2021) or finding an analytical description of the test statistic such that we would not need to simulate a large number of background and signal models.
Here, we have applied an unsupervised learning algorithm looking for clustering in data: expectation maximization (Dempster et al. 1977).This is the first time we apply expectation maximization on IceCube data and use it to fit the best time of transient neutrino emission.
The procedure is as follows (Karl & Eller 2023).For a source position to be tested (grid point), we assume a two-component mixture model for the temporal distribution of our data (a neutrino flare in the form of a Gaussian signal and uniform background).As a starting flare, we choose a single very broad flare, extending beyond the whole data-taking period.For each event, we compute the probability of it belonging to the neutrino flare (the membership probability).These probabilities are then used to improve the flare parameters iteratively.In the calculation of the membership probability for event i, we include the pdf values for the spatial and energy signal and background pdfs (as in Equation (3) and Equation ( 4)) as event weights.The membership probability is: and at each iteration the mean time, µ T , and the width, σ T , are recalculated using and The quantity n flare scales the Gaussian temporal pdf according to the expected number of signal events.However, n flare is only used to determine µ T and σ T ; n S is fitted independently once we determine the time pdf of the neutrino flare.We stop the iterations when there is no change in the likelihood in the past 20 iterations or once 500 iterations have been performed.
The signal weight, S i /B i , depends on the assumed source spectral index, γ.We want to avoid favoring a specific index; hence we run expectation maximization for different fixed spectral indices, γ EM , between 1.5 and 4 in steps of 0.2 (Karl & Eller 2023).We get an optimized time pdf for each γ EM .We then optimize the test statistic as in Equation ( 2) with the signal and background pdfs, including the temporal pdfs for each γ EM .In this step, we fit n S and γ while keeping the temporal pdf with μT (γ EM ) and σT (γ EM ) fixed.The flare yielding the highest TS value is then the best-fit flare for this grid point.For each alert+, we repeat this procedure at every grid point in the uncertainty region.The point with the most significant result is then the preferred source location.For the background TS distribution, we shuffled the event times and calculated the new right ascension values based on the event azimuths and the shuffled times.

Continuous sources
The search for the strongest single continuous source yields a global p-value of 0.98 and is compatible with the background hypothesis.We determine an upper flux limit by simulating neutrino emission with an E −2 spectrum.The upper flux limit is the flux for which 90% of the corresponding test statistic distribution lies above the test statistic value of the strongest single continuous source.We get an upper flux limit (for muon neutrinos and antineutrinos) at 90% confidence level for the most significant position of Φ νµ+νµ,single 90%,100TeV = 6.9 × 10 −17 (TeV cm 2 s) −1 .In general, the energy-dependent flux, Φ(E), of this flux limit is Φ The acceptance for the simulated flux has a limited range in energy.We define the energy range for the flux limit as the central 90% quantile of detected simulated events.In this case, we limit the flux from 0.9 TeV to 483 TeV.Table 4 lists the results for all 122 investigated regions.
For the combined emission of all sources, we get a p-value of 8%, which is also compatible with the background hypothesis.We determine the 90% confidence level upper flux limit, Φ νµ+νµ,stacked 90%,100TeV , by simulating an increasing number of sources emitting a weak flux, ϕ 1 , corresponding to one neutrino coming from a source at the celestial equator -IceCube's most sensitive region for detecting neutrinos at the highest energies -in 11 years (ϕ 1 = 4.502 × 10 −18 (TeV cm 2 s) −1 ).We repeat the simulation ∼ 10 4 times for each combined flux and create a TS stacked distribution.Based on this distribution, we determine the combined flux strong enough to yield a higher test statistic value than our result with 90% probability.
The upper limit of emission additional to the alert+ events is Φ νµ+νµ,stacked 90%,100TeV = 4.2×10 −16 (TeV cm 2 s) −1 for a spectral index of γ = 2 and within the energy range from 4.2 TeV to 3.6 PeV.For comparison, the diffuse astrophysical neutrino flux is Φ diffuse,100TeV = 1.44 × 10 −15 (TeV cm 2 s sr) −1 in the range of 15 TeV to 5 PeV (Abbasi et al. 2022b) with a  (Abbasi et al. 2022b)).The dotted purple line shows the 90% confidence level upper flux limit for the spectral index of the diffuse flux (γ = 2.37) between 0.6 TeV and 1 PeV.Left: The upper flux limit, excluding the alert+ events in the analyzed data, is 1.6% (γ = 2) of the astrophysical diffuse flux in the overlapping energy range, and 1.5% when assuming the same spectral index (γ = 2.37) as for the astrophysical diffuse flux.Right: The upper flux limit, including the alert+ events in the data, is 4.5% of the astrophysical diffuse flux in the overlapping energy range for γ = 2, and 8% of the diffuse flux when assuming the same spectral index (γ = 2.37) as for the astrophysical diffuse flux.
spectral index of γ = 2.37.Integrating over the energy range where both the diffuse flux and Φ νµ+νµ,stacked 90%,100TeV overlap, Φ νµ+νµ,stacked 90%,100TeV corresponds to 1.6% of the astrophysical diffuse flux.To constrain the maximal possible emission from the alert+ regions, including the highest-energy events, we include the alert+ events just for the following limit.Thus, considering the total emission of all 122 regions, including alert+ events, we get an upper flux limit of Φ νµ+νµ,with alerts 90%,100TeV = 1.2 × 10 −15 (TeV cm 2 s) −1 for the energy range from 4.2 TeV to 3.6 PeV, which corresponds to 4.5% of the diffuse astrophysical neutrino flux where both fluxes overlap in their energy range (see Figure 2).We repeat the upper flux limit calculation with the same spectral index as for the astrophysical diffuse flux and get a limit of 2.1×10 −16 (TeV cm 2 s) −1 (1.5% of the astrophysical diffuse flux) excluding alert+ events, and 1.1×10 −15 (TeV cm 2 s) −1 including alert+ events (8% of the astrophysical diffuse flux) at 100 TeV.For γ = 2.37, the energies of the simulated detected events range from 0.6 TeV to 1 PeV.This energy range differs from the previous range for γ = 2.The energy distribution of the signal events depends on the simulated energy spectral index.There are more neutrinos in lower energies if the simulated energy spectrum is softer compared to a harder emission.
The lack of lower energy neutrino emission (compared to IceCube alert+ events) could be caused by various scenarios.It is, for example, possible that some sources flare in neutrinos, emitting mainly high-energy neutrinos.Another possibility might be a hard neutrino emission, i.e., γ ≤ 1 (for example, models proposed in Waxman & Bahcall (1999); Padovani et al. (2022)).The atmospheric background would dominate the lower-energy neutrino emission.The higherenergy neutrino emission would be detected as single high-energy events, given IceCube's effective area (Aartsen et al. 2020;Abbasi et al. 2021b).This matches our observation.However, there are many different scenarios that agree with this work.In these cases, different source populations or states produce different neutrino spectra compared to one continuous power law.Another possible scenario including a source population emitting single power-laws is described in Abbasi et al. (2023b).Our result agrees with the high-density scenario presented in Section 6. There, a high-density source population with low individual fluxes (with an E −2.5 energy spectrum) is the origin of alert events.Due to the sheer number of sources, we would be able to detect flux fluctuations in high-energies as alert events without a detectable lower-energy component.In lower energies, the flux would be too low to be detected and it would require a simultaneous fluctuation in both lower and higher energies such that both components could be detected from the same object.

Transient sources
In our search for transient sources, we look for the most significant transient neutrino emission.Out of all the investigated 122 alert+ origins, the most significant transient emission is the neutrino flare with the seed alert IceCube-170922A, which is associated with the blazar TXS 0506+056.Our search yields a local p-value of 0.14% (or a significance of 3σ).The main differences between the search in Aartsen et al. (2018b) and this work are: • We have no external trigger in this work, whereas Aartsen et al. (2018b) was triggered by the observation of a flaring blazar.
• We use 11 years of recalibrated IceCube muon data, improving directional and energy reconstruction.For a discussion of how the contributing events are affected, see Appendix B.
• We include a fit for the best source position and use expectation maximization to identify the time of the neutrino flare.
The corresponding flare is centered around a mean flare time μT = 57001 +38 −26 MJD and has a width of σT = 64 +35 −10 days.These properties agree with Aartsen et al. (2018b), as shown in Figure 3.When correcting for the look-elsewhere effect, the global p-value is p global = 0.156, which is not significant.Such a trial correction does not apply for the search reported in Aartsen et al. (2018b).Table 5 lists all results for the investigated regions.
The best-fit parameter can yield insight into the source emission.However, as mentioned in Section 2 and Appendix A, the best-fit results and the resulting flux estimations are biased.The best-fit result of the number of neutrinos in the neutrino flare is nS = 12 +9 −7 with a spectral index of γ = 2.3 ± 0.4.This corresponds to an average flux of Φ νµ+νµ 100TeV = 1.1 +0.9 −0.8 × 10 −15 (TeV cm 2 s) −1 in the energy range of 3.5 TeV to 213 TeV during the period of the neutrino flare.The corresponding single flavor neutrino and anti-neutrino fluence, the flux integrated over the flaring period (μ T − 2σ T to μT + 2σ T ), is J νµ+νµ 100TeV = 1.2 +1.0 −0.8 × 10 −8 (TeV cm 2 ) −1 .This flux estimation also agrees with Aartsen et al. (2018b), as shows the all-flavor neutrino flux (three times Φ νµ+νµ 100TeV ) in Figure 4.In Appendix B, we compare the events contributing to the neutrino flare in this analysis and previous works and explain why the errors differ.
For transient emission, the lack of additional lower-energy neutrino emission (besides the reported local evidence associated with TXS 0506+056) can imply various scenarios.One is that neutrino flares occur rarely or might not necessarily be connected to the production sites of high-energy neutrinos.Similarly to Section 3.1, it could also indicate that these neutrino sources emit a very hard energy spectrum, for example, with γ ≤ 1.

CONCLUSION
Our study focused on the origin of IceCube's highest energy events, or alert+ events, to identify potential sources of additional neutrino emission.To achieve this, we systematically scanned the 90% uncertainty contours of reconstructed alert+ events, with a resolution of 0.2 degrees, to determine the most significant source position.We assumed that the emission followed a power-law distribution, ∝ E −γ , with γ ranging from 1.5 to 4.
Our analysis found no evidence for continuous emission from a single source, as the data were consistent with the background assumption.Therefore, we placed a constraint on the overall combined flux from all positions, which was found to be 1.6% of the diffuse astrophysical neutrino flux observed by IceCube (for γ = 2).If we included the alert+ events in the analysis, we could constrain all expected emissions from their respective directions to no more than 4.5% of the diffuse astrophysical neutrino flux (for γ = 2).For a source spectral index similar to the diffuse astrophysical neutrino flux (γ = 2.37), we constrain the overall combined flux to be less than 1.5% (excluding the alert+ events) and less than 8% (including the alert+ events) of the diffuse astrophysical neutrino flux.This indicates that different source populations or states produce different neutrino spectra compared to one continuous power law.
Our investigation confirmed the neutrino flare associated with the blazar TXS 0506+056 as the most significant transient emission from all investigated positions, with a local significance of about 3σ.When we corrected for the look-elsewhere effect in this analysis, the global significance was 15.6%, consistent with the background expectation.The parameters of the neutrino flare in this study using recalibrated data agreed with previously published results.We identified a Gaussian time window with a center at 57001 +38 −26 MJD, and a width of 64 +35 −10 days as the best fit and estimated that 12 +9 −6 neutrinos were detected during the flare with a best-fit spectral index of γ = 2.3 ± 0.4.This corresponds to a single flavor neutrino fluence of J νµ+νµ 100TeV = 1.2 +1.0 −0.8 × 10 −8 (TeV cm 2 ) −1 and an average flux of Φ νµ+νµ 100TeV = 1.1 +0.9 −0.8 × 10 −15 (TeV cm 2 s) −1 during the 2σ T time window.However, we find no other alert+ event with a similar local significance.TXS 0506+056 remains the only source candidate where we find a connection of a high-energy alert and a lower energetic neutrino emission.
For neither continuous nor transient emission did we find evidence of a lower energy neutrino component.This can be explained in various scenarios.One is a hard neutrino spectrum with γ ≤ 1.In such a scenario, atmospheric background noise would dominate the lower energy range, while the higher energy range would yield single high-energy events.It could also be caused by a high-density source population as investigated in Abbasi et al. (2023b), where high-energy events are the result of fluctuations from a large population of sources with individually weak fluxes.In this case, the lower energy flux would still be too low to be detected.Our finding also suggests neutrino flares may be rare or produced at different sites than IceCube alert+ events or that there are sources mainly emitting high-energy neutrinos.
actual simulated source is smaller than 0.3 degrees for a flux resulting in 5 signal neutrinos on average.This also means that the best-fit n S will be larger than 0 in many cases with no neutrino source since the algorithm will find the position with the largest background fluctuation.Hence, correcting this bias is not straightforward, and this analysis is mainly sensitive to strong neutrino sources.
For transient sources, the bias is smaller.In the same example as above, 10 neutrinos with γ = 2 emitted over a period of σ T ≈ 55 days are a much stronger signal compared to 10 neutrinos over 11 years.Hence in this specific case, the mean best-fit nS = 12 and the best-fit γ = 2.1.However, we still face the case that background fluctuations can dominate weak neutrino emission (in the case of σ T ≈ 55 days, anything below 5 neutrinos is difficult), which makes correcting this bias challenging.We have decided not to include an, at best incomplete, correction in this work.For now, measurements of point-source fluxes are only possible with the KDE approach.

B. TRANSIENT SOURCES ANALYSIS
Figure 5 shows the p-value map of the scanned region around IceCube-170922A on the left.The most significant position is within 0.5 • from TXS 0506+056.The right panel of Figure 5 shows a histogram of the angular distance of events from TXS 0506+056.There is a clustering of events around the source position.The signal events for this plot are simulated according to the best-fit result of the likelihood ratio test (n S = 12, γ = 2.3).The background distribution is scrambled data in right ascension.The signal on top of the background flux matches the observed data.To determine the uncertainties of the best-fit values, we run a likelihood scan over the parameter space and use Wilk's theorem (Wilks 1938) to determine the 68% and 90% contours (see Figure 6).These contours are relevant for the two-dimensional uncertainties of the flux as in   Table 2 lists the top 14 contributing events to the neutrino flare, sorted by their S i /B i value multiplied by S T .We compare this with a previous data sample (Abbasi et al. 2021b) (for events also included in that sample) to emphasize how the updated photomultiplier calibration affects the reconstructed direction, angular error, and energy.
The improved directional and energy reconstruction has changed the contributing events compared to previous analyses (Aartsen et al. 2018b;Abbasi et al. 2021b).Most of the significance is caused by the two most contributing events, which remain the same (see also Karl (2022)).However, their position is shifted, and their energy is changed.For the remaining events, the contributing order has changed, or the events themselves differ.Figure 7 shows the position and energy of the 14 most contributing events to the neutrino flare from the previous data set (left) and the improved data used in this work (right).The event with the largest error region (σ = 1.9 • ) on the right panel is also included in the left panel.However, the uncertainty was underestimated in the previous data sample (σ = 0.36 • ), and its position has shifted.2).Right: The 14 most contributing events from the old data sample (Abbasi et al. 2021b).

C. ICECUBE ALERT+ EVENTS
Table 3.All alert events (track name starting with "IC") and highenergy tracks (track name starting with "DIF" (selected from Abbasi et al. ( 2022b)), "EHE" (extremely-high-energy), or "HESE" (highenergy-starting event)) investigated in this work.The track name includes the time of detection in the format yymmdd.In the case of alert events, the letter "A" or "B" is used to distinguish events detected on the same day.The time is the detection time in MJD, and R.A. and Dec list the best reconstruction coordinates with 90% confidence level uncertainties.

Figure 2 .
Figure2.90% confidence level upper flux limits assuming for all source candidates combined (dashed orange line) valid in the energy range of 4.2 TeV to 3.6 PeV and a neutrino emission following E −2 .The green line is the diffuse astrophysical neutrino flux (Φ diffuse,100TeV = 1.44 × 10 −15 • 4π (TeV cm 2 s) −1 ) in the range of 15 TeV to 5 PeV(Abbasi et al. 2022b)).The dotted purple line shows the 90% confidence level upper flux limit for the spectral index of the diffuse flux (γ = 2.37) between 0.6 TeV and 1 PeV.Left: The upper flux limit, excluding the alert+ events in the analyzed data, is 1.6% (γ = 2) of the astrophysical diffuse flux in the overlapping energy range, and 1.5% when assuming the same spectral index (γ = 2.37) as for the astrophysical diffuse flux.Right: The upper flux limit, including the alert+ events in the data, is 4.5% of the astrophysical diffuse flux in the overlapping energy range for γ = 2, and 8% of the diffuse flux when assuming the same spectral index (γ = 2.37) as for the astrophysical diffuse flux.

Figure 3 .
Figure 3. Logarithm of the signal-over-background ratio, log 10 Si/Bi, distribution of individual events, i, versus their detection time, ti, between 2012 and 2016.The log 10 Si/Bi values are for the best-fit position (close to TXS 0506+056) and the best-fit spectral index.The color indicates the reconstructed muon energy, Erec,µ, increasing from light to dark.The black-dashed line shows this work's best-fit time pdf ST (with the y-axis on the right).It agrees with the grey-dashed pdf of Aartsen et al. (2018b).

Figure 5 .
Figure 5. Left: P-value map of the alert region of IceCube-170922A.The grey dot indicates the reconstructed direction of IceCube-170922A, and the grey contour shows the 90% uncertainties of the reconstruction.The red cross marks the best-fit position of the position scan (0.6 • from the reconstructed alert position).The star shows the location of TXS 0506+056.All black bins have p-values close to 1. Right: Number of events at binned squared angular distances, Ψ2 , between TXS 0506+056 and the reconstructed event directions during the neutrino flare (57001 MJD ± 2 × 64 days).Scrambled data in right ascension provides the background (blue), and Monte Carlo simulations for the best-fitted flux (nS = 12 and γ = 2.31) yield the signal (orange).The grey line combines the background with the signal and matches the data points (black).The data are shown with 68% uncertainties.

Figure 4 .
For the time, we determine the profiled change of the test statistic for different µ T and σ T .The best n S and γ are fitted for each value.The 68% uncertainties determined by a profiled change of the test statistic are μT = 57001 +38 −26 MJD and σT = 64 +35 −10 days.The one-dimensional errors for fluence, number of signal neutrinos, n S , and spectral index γ, are determined with the profiled change of the test statistic where the mean flaring time and the flare width are kept fixed to the best-fit values.For the signal fluence the 68% uncertainties are are J νµ+νµ 100TeV = 1.2 +1.0 −0.8 × 10 −8 (TeV cm 2 ) −1 and for n S and γ we get nS = 12 +9 −7 and γ = 2.3 ± 0.4.

Figure 7 .
Figure 7. Position and energy (color) of the 14 most contributing events to the TXS 0506+056 neutrino flare.The circles show the uncertainty of the directional reconstruction, σi.Left: The 14 most contributing events from the data sample used in this work (see Table2).Right: The 14 most contributing events from the old data sample(Abbasi et al. 2021b).

Table 2 .
Left: Top 14 events with the strongest contribution to the neutrino flare of TXS 0506+056, sorted by significance.
Right: The same events in the data sample published inAbbasi et al. (2021b).The last row states the ranking of the contribution in previous analyses.The data set used in this work has improved directional and energy reconstruction.Some events have shifted in position and have slightly different energies.

Table 4 .
Results of the individual time-integrated analysis sorted by significance.The first column contains the index of the alert+ event as in Table3.The two following columns list the best-fit position of this work.The fourth and fifth columns contain the best-fit parameter of the likelihood optimization ns and γ.The sixth column shows the local p-values, and the seventh column the 90% confidence level upper flux limits.The central 90% quantiles of neutrino energies of the detected simulated events for computing the flux limit are listed in columns eight and nine.They define the range in which the flux limit is valid.The global p-value for the time-integrated single-source search is 0.98.

Table 5 .
Results of the time-dependent analysis sorted by significance.The first column contains the alert+ index as in Table3.The next two columns list the best-fit position.The fourth and fifth columns contain the best-fit parameter of the likelihood optimization ns and γ.The sixth and seventh column list the best-fit results for the Gaussian time window with mean µT and width σt.The last column shows the local p-values.The global p-value for the time-dependent analysis is 0.156.