Constraints on Populations of Neutrino Sources from Searches in the Directions of IceCube Neutrino Alerts

Beginning in 2016, the IceCube Neutrino Observatory has sent out alerts in real time containing the information of high-energy ( E  100 TeV ) neutrino candidate events with moderate to high (  30% ) probability of astrophysical origin. In this work, we use a recent catalog of such alert events, which, in addition to events announced in real time, includes events that were identi ﬁ ed retroactively and covers the time period of 2011 – 2020. We also search for additional, lower-energy neutrinos from the arrival directions of these IceCube alerts. We show how performing such an analysis can constrain the contribution of rare populations of cosmic neutrino sources to the diffuse astrophysical neutrino ﬂ ux. After searching for neutrino emission coincident with these alert events on various timescales, we ﬁ nd no signi ﬁ cant evidence of either minute-scale or day-scale transient neutrino emission or of steady neutrino emission in the direction of these alert events. This study also shows how numerous a population of neutrino sources has to be to account for the complete astrophysical neutrino ﬂ ux. Assuming that sources have the same luminosity, an E − 2.5 neutrino spectrum, and number densities that follow star formation rates, the population of sources has to be more numerous than 7 × 10 − 9 Mpc − 3 . This number changes to 3 × 10 − 7 Mpc − 3 if number densities instead have no cosmic evolution.


INTRODUCTION
Nearly a decade after the detection of a diffuse flux of astrophysical neutrinos (Aartsen et al. 2013), the origins of this flux largely remains a mystery.However, to that end, IceCube -a cubic kilometer neutrino telescope operating at the geographic South Pole -found compelling evidence that a blazar, TXS 0506+056, is a source of high-energy neutrinos (Aartsen et al. 2018a,b), though this object alone can only explain a small portion of the diffuse flux.Pinpointing more sources of cosmic neutrinos, or understanding the populations of sources which contribute to the overall measured diffuse flux (Aartsen et al. 2015(Aartsen et al. , 2020a;;Abbasi et al. 2021a,b) could prove pivotal in understanding the processes behind the acceleration and propagation of high-energy cosmic rays.
The identification of the blazar TXS 0506+056, as well as some more recent claims of possible neutrino sources, e.g.Stein et al. (2021); Franckowiak et al. (2020), were enabled in part because of the correlations of these sources with public IceCube neutrino alerts (Aartsen et al. 2017a).In addition to searching for correlations with individual alerts, some have looked for correlations between catalogs of sources and neutrino alerts, e.g.Plavin et al. (2020).These high-energy neutrino alerts are often used to trigger follow-up observations because of their significant probabilities of astrophysical origin.By only looking at neutrino candidate events with high estimated initial energies (E ν 100 TeV), the typically overwhelming backgrounds from atmospheric cosmic-ray interactions can be suppressed.
While it is clear that using high-energy neutrino alerts to trigger multi-wavelength (MWL) observations is a promising and fruitful way to identify cosmic neutrino sources, it is not without its limitations.First, it is not clear which astrophysical objects are sources of highenergy neutrinos.This, in combination with the nonnegligible localization uncertainties for neutrino events, can lead to source confusion when using pointed MWL observations, especially as recent limits suggest that any neutrino source population responsible for a significant fraction of the diffuse flux must be fairly numerous (Aartsen et al. 2019a).Second, there is not a consensus on the types of MWL emission that one expects to see with high-energy neutrinos.For example, works motivated by the IC170922A event in 2017 found to be coincident with a gamma-ray flare of TXS 0506+056 providing evidence for a connection (Aartsen et al. 2018b), have studied for potential correlation between neutrino emission and gamma-ray emission in blazars.These works have yielded constraining upper-limits, e.g.Murase & Waxman (2016); Aartsen et al. (2017b); Murase et al. (2018); Yuan et al. (2020); Oikonomou et al. (2019), with more recent studies suggesting that correlation with time-dependent behavior at other frequencies might be more telling for correlating with neutrino emission (Kun et al. 2021).Third, try-ing to consistently model neutrino emission from sources based on alert and MWL observations are subject to the Eddington bias (Strotjohann et al. 2019), because alert events are likely from sources where, although the joint contribution from a population of sources to the diffuse flux might be large, alert events from individual sources are likely Poisson overfluctuations.Lastly, pointed MWL follow-up observations for alert events can be expensive, and not all alerts are followed up by instruments either because of low event-by-event astrophysical probabilities or because of observational constraints.
In this work, we rely on a different method of using neutrino alerts to search for astrophysical neutrino sources -by following up neutrino candidate alert events using lower-threshold and higher statistics neutrino data.We use the word "followup" to denote performing analyses that use other neutrino data to search for astrophysical signals in the vicinity of the neutrino candidate alert.In addition to circumventing many of the difficulties that surround using MWL observations to followup neutrino alerts, this search strategy complements the all-sky searches and catalog based searches by reducing the number of unique locations on the sky which need to be investigated, without requiring a fixed hypothesis on the astrophysical source class.This large trials-factor from the "look-elsewhere effect" typically degrades the sensitivity of searches that look for possible neutrino sources at every location on the sky.This could be one of the reasons why no source has been detected at above the 3σ level from these types of searches (Aartsen et al. 2020b).
In this paper, we report the results of searches for neutrino sources in the directions of IceCube neutrino candidate alerts.In Section 2, we describe the data samples used and then we outline the analysis techniques to analyze these data in Section 3. In Section 4, we show how we combine the results from the individual analyses to search for an overall excess of lower-energy neutrinos in the direction of high-energy neutrino alerts.After discussing the results of the analysis in Section 5, we show how these results can be used to constrain populations of neutrino sources in Section 6.

DATA SAMPLES
The IceCube Neutrino Observatory is a gigaton-scale Cherenkov detector embedded in the ice at the geographic South Pole (Aartsen et al. 2017c).The detector consists of 5,160 Digital Optical Modules (DOMs) dispersed on 86 "strings" arrayed in a hexagonal grid and deployed at depths of 1450-2450 m beneath the ice surface.Each DOM contains a 10 inch photomultiplier tube suited to detect optical Cherenkov photons (Ab- 90% cont.Although sensitive to all flavors of neutrino interactions, this study relies only on muon track events from muon-neutrino charged current interactions as well as a 10% contribution from muonic tau decays from charged current tau-neutrino interactions.These "track" events enable a better angular resolution than the other event type, "cascades" (from charged current electron-and tau-neutrino interactions or neutral current interactions of all flavors), at the cost of a poorer energy resolution.The angular resolution of track events is preferable when searching for point sources in the region of the sky where most of the alert events are detected.
In addition to neutrinos from astrophysical sources, IceCube detects many neutrinos and muons from cosmic-ray interactions in the atmosphere.In the southern celestial hemisphere, the events detected by IceCube are dominated by atmospheric muons, with events from atmospheric neutrinos still occurring at rates a few orders of magnitude larger than astrophysical neutrinos.In the northern celestial hemisphere, the atmospheric muons are attenuated by the Earth, and the rate is dominated by atmospheric neutrinos.
The analysis presented here leverages the strengths of two different IceCube data streams: the alert-event stream and the "gamma-ray follow-up" (GFU) stream (Aartsen et al. 2016).Both of these event selections try to isolate neutrino candidate events with low latency, enabled by a realtime alert infrastructure that began sending alerts publicly in April 2016, and are described in full in Aartsen et al. (2017a).However, the selection criteria used to identify alert events was revisited in 2019 (Blaufuss et al. 2020) to expand the alert program, and the alert stream now consists of two unique channels: "Gold" events which have an average astrophysical signal purity above 50%, and "Bronze" which have an average astrophysical signal purity above 30%.These event-by-event astrophysical purities are calculated by finding the event "signalness", S, which is the ratio of the expected number of events from signal to the expected total number of events (signal plus background) at a given declination with energies greater than the reconstructed energy of the event (Blaufuss et al. 2020).The final rate is approximately 10 events per-year in the "Gold" selection, and 30 events overall in the "Gold" and "Bronze" selection.Signalness is dependent on the assumed spectral index.To avoid this dependence, the alert stream effective area, which is a function of the energy, is used for the analysis while treating all events to be on the same footing.
Whereas the alert stream is optimized for finding individual events with moderate-to-high probability of as-trophysical origin, the GFU sample is focused on optimizing sensitivity to short timescale transients.When searching for transient neutrino emission, the effective background of the analysis is reduced because each analysis only looks at a narrow swath of sky and for a limited period of time.Because of this, the cuts for the event selection can be looser than for the alert stream.The final all-sky rate is ∼ 6.7 mHz (approximately 2 × 10 5 events per year).While the vast majority of this is from atmospheric backgrounds, the effective area for astrophysical neutrinos is significantly larger with the GFU sample than with the alert sample.For example, consider the quantity which is the expected number of detected events for a source at a declination δ with an E −γ spectrum with normalization φ 0 in a given event stream (either the alert stream or GFU).The term A eff stream describes the effective area of the event selection.Then, for a source with γ = 2.5 (γ = 2.0), the ratio of the expected number of events in the GFU stream to the alert stream is 97 (13) for a source at δ = 0 • and 270 (40) for a source in the northern sky.This means that if a source, in the same direction of the alert, is bright enough for there to be an expectation of observing one alert event, then a search with the GFU event selection could result in tens to hundreds of lower-energy events coincident with the source.This also helps to remove some of the aforementioned Eddington bias, as using a larger statistics sample can provide a better estimate of the true source flux.It is worth noting that most (88%) of the alert events are included in the GFU sample.When we perform a followup for each alert, we remove the alert that prompted the analysis from the followup.
In addition to events detected in real time after the creation of the alert selection criteria, archival IceCube data dating back to 2011 have been processed with the same sets of cuts used for the alert and GFU streams.Overall, our datasets span May 13, 2011 to December 31, 2020, and the final alert stream has 275 events and the GFU selection has 1.8 × 10 6 events from this period.
Complementing the reconstructions that were applied to each event in the GFU stream, described in Aartsen et al. (2017a), alert events have an additional reconstruction applied to them.This includes effects from systematic uncertainties after they have been initially circulated to the public.Because of the computational constraints of this reconstruction, it is reserved for only the events passing alert cuts.In order to calculate uncer-tainty contours for the directional reconstruction of each alert event, they are compared against re-simulations of other high-energy track events which varied the allowed models of the optical properties of the deep glacial ice.This process is described in full in Abbasi et al. (2021c), and it allows one to quote 50% and 90% containment contours for the directional reconstruction of the event.This process was first applied for the event IceCube-160427A, for which a likely unrelated coincident supernova was found by Pan-STARRS (Kankare et al. 2019).For events detected in real time, bounds to these these contours are typically circulated to the community within a few hours of the initial detection of the event.
Figure 1 shows the sample of alert-quality events that we analyze with our followup analyses by showing the 50% and 90% contours calculated using the procedure outlined above.A few events from the full catalog of alert events are excluded from our analysis because the computation time for each followup scales with the alert uncertainty size, and some alert events had exceptionally large angular uncertainties which prevented us from performing enough pseudo-experiments to characterize our background expectations.Although those events are excluded from the Figure 1, we tabulate them in the full table of results in Appendix A. Additionally, when searching for emission on short timescales, some alerts are excluded just from the transient analyses we perform if there was a significant period of detector downtime near the time of the alert.These are also mentioned in the Table in Appendix A.

ANALYSIS METHOD
The analyses we perform rely on an unbinned maximum likelihood technique that is well-established in many searches (Braun et al. 2008).For recent examples of such analyses, see e.g.Aartsen et al. (2020bAartsen et al. ( , 2019a)).For each of these analyses, we are searching for GFU events that are spatially clustered and coincident with alert events.This is done by checking if the alert event which prompted the followup is also in the GFU sample, and if it is, we remove it from the search.As we do not know if populations of neutrino sources that are responsible for the diffuse flux are transient sources or sources that emit over long timescales, we perform several analyses that test different temporal hypotheses.Specifically, we perform likelihood analyses on three different timescales: (1) centered on the alert time and searching for spatially coincident events in a time-range of ±500 s, (2) centered on the alert time and searching for spatially coincident events in a time-range of ±1 day, and (3) searching for spatially coincident events during the entire livetime of the GFU sample.The durations of the transient analyses were chosen to strike a balance between theoretical and empirical constraints.The shorter timescale (±500 s) has been suggested as a conservative approximate timescale for neutrino emission from compact binary mergers (Baret et al. 2011;Aartsen et al. 2020c), and has been used in many IceCube searches for transients.The 2 day timescale reflects the longest timescales of neutrino emission proposed for neutrino emission from some short timescale transients, e.g.Fast Radio Bursts (Metzger et al. 2020), and it is also the maximal time window in which our analysis remains sensitive to individual coincident events.It is worth noting that searches that look for GFU events that are clustered in time but are not necessarily coincident in time with the alert event (so-called "flare analyses") are not performed in this work, although a dedicated analysis repeating the flare analysis that identified the 2014-2015 flare from TXS 0506+056 are in development (Abbasi et al. 2021d).
At the core of each analysis presented here is the same likelihood framework, with some differences which we outline below.Given a location on the sky with Equatorial coordinates, Ω = (α, δ), the likelihood consists of the weighted sum of a signal probability distribution function (PDF), S, and a background PDF, B, and is given by The index i iterates over all N neutrino candidate events in the GFU sample and n s is the number of signal neutrino events.The signal and background PDFs, S and B, are functions of four observables associated with each event: the reconstructed right ascension and declination, Ω i = (α i , δ i ), the reconstructed energy, E i , and the angular uncertainty σ i .Both S and B are themselves products of two terms: energy and spatial PDFs.
The signal energy PDF is a declination-dependent reconstructed energy distribution, where the underlying neutrino flux is modeled as a power-law, where φ 0 is the flux normalization and γ is the spectral index.The spatial term of the signal PDF is modeled as a Gaussian with width σ i , given by the angular uncertainty estimator of each neutrino candidate event in the GFU sample.The energy and spatial components of the background PDF are determined as functions of the reconstructed energy, E i and declination, δ i , for each event.A more thorough description of how the signal and background PDFs are calculated is provided in Aartsen et al. (2017d).
The difference in the transient analyses and the timeintegrated analysis is encapsulated in the λ term in Eq. 2. When searching for short timescale transient emission, the data is divided into the time period of interest, "on-time" data, and the remaining, "off-time", data.We then set λ to a Poisson probability, λ = where n b is the expected number of background events on the entire sky in the time-window estimated using the surrounding off-time data.This "extended likelihood" methodology (Barlow 1990) has been a feature of many analyses searching for short timescale neutrino emission (Abbasi et al. 2021e;Aartsen et al. 2017e, 2020d.For the time-integrated analysis looking for steady emission over the entire livetime of the GFU sample, we set λ = 1.
The best-fit signal parameters at a given source position are then obtained using the maximum-likelihood method.For the transient analyses, we only maximize the likelihood with respect to n s , and we keep γ fixed to 2.5.This is because in transient analyses we are looking for few individual coincident events, and it is not feasible to fit a spectral index from the observation of a single event.For the time-integrated analysis, we maximize with respect to both n s and γ.Using a likelihood ratio test, we calculate the point-source test statistics T S ps , as where ns is the best-fit n s and γ is the best-fit γ (or fixed to 2.5 for the transient analyses).Here, the null hypothesis is defined as n s = 0, representing the case of no neutrino source in the direction Ω.
Thus far, we have described how to search for a source at a location Ω.However, we would like to search for sources in the direction of IceCube alert events, and these alert events are not perfectly localized.In order to do this, we first divide the sky into a grid using the HEALpix scheme (Górski et al. 2005) and calculate T S ps on each point of this grid, using a resolution of approximately 0.2 • .
We then use the likelihood scan of each alert event to create a skymap.Abbasi et al. (2021c) details how 50% and 90% containment contours are calculated based on finding critical likelihood values.However, directly normalizing the likelihood space as a function of location on the sky would neglect these critical values.In order to obtain a skymap that reflects the uncertainties dictated . Schematic of the likelihood analysis described in Section 3.For each neutrino alert candidate event, the local point source T S-map is calculated after excluding the alert event from the dataset (left).From this map we then subtract a term that is constructed from the likelihood map derived from a more sophisticated and dedicated reconstruction of the alert event (middle) to produce an overall map (right) from which the maximum is used as the analysis test statistic, T S, corresponding to the location denoted here with a white cross hair.The alert shown here, Run120027:Event12133428, was chosen just for visualization purposes, and black contours show 50% (solid) and 90% (dashed) containment derived from the map displayed in the middle.This shows a time-integrated followup, but we repeat the procedure for the shorter timescale followups as well.
by systematic resimulation studies, we apply an orderpreserving transformation to re-calibrate the likelihood space such that when we normalize the likelihood as a function of location on the sky, approximately 90% of the skymap lies within the quoted 90% systematic contour (Abbasi et al. 2021a, Appendix I).We then use an algorithm which effectively fits for the position of a point-like source in the environment of the alert direction.This is done, once we have a skymap normalized as a PDF as a function of location on the sky for an alert event, P j ( Ω).We include the skymap as a spatial constraint by multiplying it to the neutrino likelihood function via L → L • P j ( Ω).Because P j ( Ω) is independent of the variables which are floated when maximizing the likelihood, n s , γ, the constraint manifests as adding a logarithmic term to the pointsource test statistic defined for each grid point.Finally, the test statistic from each alert followup is given by where w( Ω) = P j ( Ω)/ max Ω (P j ( Ω)).Normalizing the skymap term, w, based on the maximum value of the skymap is just a choice of convention, and different convention choices add a constant term to the test statistic which can be omitted when calculating significances and p-values based on pseudo-experiments which are unique to each individual alert followup (see below).The analysis process is shown schematically in Figure 2, which shows the method to calculate the maximum test statistic for each alert followup.This is done by calculating the local point-source test statistic map and adding to it a logarithmic penalty term from the skymap.This procedure was initially developed in the context of searching for joint sources of ultra-high-energy cosmic rays and neutrinos (Albert et al. 2022) and has been used to search for neutrinos coincident with ANITA events (Aartsen et al. 2020e), gravitational waves (Aartsen et al. 2020c), and gamma-ray bursts (Abbasi et al. 2022a).We then repeat this followup procedure for every alert in our catalog and for each of the 3 timescales described above.
In order to calculate significances for each of the alert followups, we compare each observed test statistic to pseudo-experiments with scrambled data.For all pseudo-experiments, the PDFs in the likelihood and the alert skymaps are kept fixed.We then randomize the GFU data in right ascension.We also perform tests to see how well the analysis is able to recover simulated signal.We randomize the arrival directions in the same method as we do for background-only pseudoexperiments, but in these cases, we also inject signal events assuming a true source position that is sampled from the skymap of the alert event.We can summarize the analysis performance by calculating a "sensitivity" for each followup, defined as the expected median onesided Neyman upper limit at 90% confidence under the assumption of the null hypothesis (no source in the direction of the alert event).As some of the alert events have extremely large localization uncertainties, this increases the effective background and degrades the sensitivity of the analysis.For the short timescale analyses, this results in a sensitivity that can be up to 10% worse than in the case of no localization uncertainty, when assuming the source has an E −2.5 spectrum.For the timeintegrated case, the sensitivity (and correspondingly, the upper limits in the cases of non-detections) when searching for sources coincident with alert skymaps is any-where from 10% to a factor of 2 worse than the localized point source case.In the case of potential detections, the greater effective background as compared to the localized point source scenario can have a larger effect on the signal strength required to confidently detect a source.For example, for most of the alert skymaps, the signal strength required to detect a source at the 3σ level, before trials correction, using this skymap construction is at least a factor of 2 higher than the signal required to detect a source in the perfectly localized point source scenario.

POPULATION ANALYSIS
Thus far, we have described how to search for significant correlations of GFU events with individual alert events.We perform the analysis described above for each of the alerts in the catalog and for each of our three analysis timescales.However, we would like to search for signals from populations of sources which might manifest as multiple alert followups that are individually not significant, but when combined, the total flux is detectable.
In order to accomplish this, we begin by calculating pvalues for each of the individual followups using the procedure outlined in Section 3.Then, for each time window, we sort the p-values into a list, p t,1 , p t,2 , . . ., p t,N , where the index t identifies the time-window1 .Under the background hypothesis (n s = 0), these N p-values are expected to be uniformly distributed between 0 and 1.The probability that k or more p-values are smaller than or equal to p t,k is thus given by the binomial probability: We then evaluate this probability for all possible number of sources, k, to calculate the overall binomial pvalue, α = min k α k .In order to account for the fact that we performed multiple tests by finding the most significant number of sources, k, we repeat this entire procedure on ensembles of background pseudoexperiments.We can then calculate an overall analysis p-value for each time window by comparing α to the distribution of this value that we obtain from these pseudo-experiments.When performing these ensembles of background pseudo-experiments, we use the same sets of scrambled data for all of the alert followup analyses.This ensures that any correlations between overlapping alert contours is properly accounted for when calculating an overall analysis p-value.
In Section 2, we discuss how each alert has a corresponding "signalness", S. Preliminary versions of the analysis attempted to incorporate this signalness parameter into the overall population test statistic, so as to give more weight to alert events that have a higher probability of astrophysical origin.However, the signalness of each alert event is calculated under the assumption that the diffuse astrophysical neutrino flux is described by the single power law fit reported in Haack & Wiebusch (2018), and thus the signalness is sensitive to changes in the spectral shape of the diffuse astrophysical neutrino flux.For this reason, we instead treat all alert events on equal footing, regardless of signalness.

ANALYSIS RESULTS
The pre-trials p-values for all of the analyses are displayed in Figure 3, and we also tabulate these values as well as the best-fit information from the time-integrated analysis in Table 3 in Appendix A. We do not include the best-fit information for the short-timescale analyses in Table 3, as most of the best-fit number of events are ns = 0 and because we do not fit for the spectral index in those analyses.
No individual alert followup is significant, especially after accounting for trials correction.The most significant followup comes from the ±500 s followup of the alert event 118342:24578488, which was detected on 2011-06-16.A single spatially coincident event with a reconstructed energy of ∼750 GeV that arrived 346 s after  .Binomial scan (left) and distribution of binomial p-values (right) calculated from pseudo-experiments using data scrambled in right ascension for the time-integrated analysis, as described in Section 4. In the binomial scan, the results for all alert followups is ordered by decreasing significance, and for each possible number of sources, k, we calculate the probability of obtaining k results each with p-values less than or equal to p k .The global minimum corresponds to k = 23 and a binomial p-value of α = 6.1 × 10 −4 (black line).The observed binomial p-value from data is shown with a black line, compared to the expected background distribution in blue, and results in an overall p-value of 1.8% for the time-integrated analysis.
the alert event drives the significance, yielding ns = 1 and a pre-trials p-value of 6 × 10 −4 .After accounting for the fact that we performed nearly 275 analyses for each of the three time windows, this is consistent with background expectations, with a trials corrected p-value of approximately 40%.
It is worth noting that the alert event that was identified as coincident with TXS 0506+056, IceCube-170922A, does not yield a significant result in this analysis (130033:50579430 in Table 3 with pre-trial p value of 0.36 for the time-integrated case).However, this is completely consistent with what was reported in Aartsen et al. (2018a), mainly because of the temporal hypothesis that is being tested.The analysis in Aartsen et al. (2018a) was searching for events that were clustered in time, and the clustering of events in 2014-2015 led to a more significant result than if the temporal hypothesis had been searching for time-integrated emission, as we do here.An analysis that performs the flare search that was done in Aartsen et al. (2018a) on the entire catalog of alert events is beyond the scope of this work, because of the computational constraints that come from the added dimensionality of maximizing the likelihood with respect to parameters that describe the temporal signal hypothesis.Additionally, that analysis was performed at the location of the source TXS 0506+056, whereas this analysis has increased background because we consider all locations within the uncertainty contour of the alert event.We do still fit ns = 17.5 events, and the best-fit location is less than 0.2 degrees from the object TXS 0506+056.When comparing against other time-integrated analyses (Aartsen et al. 2020b), we fit a slightly softer spectrum (γ = 2.75), mainly because we exclude the alert event from our sample when performing the followup, and if it were included it would shift the spectral fit to be harder.
In addition to performing the individual followups, we can also perform the population tests outlined in Section 4 for each timescale.The results for each of these tests is shown in Table 1.The most significant result comes when searching for steady neutrino sources over the entire GFU livetime, which yields an analysis p-value of 0.018 when comparing the observed α (6.1 × 10 −4 ) to a distribution generated from ensembles of pseudoexperiments.We find this to be consistent with expectations from background, especially after considering that we perform three tests, one for each time window, which is not reflected in the quoted p-value above (treating the time-windows as independent would result in a trials corrected p-value of approximately 5%).The comparison of the observed value and the background expectation is shown in Figure 4, where we also show how we calculate α, namely, by calculating α k for all possible number of sources, k, and returning the minimum.The analysis finds a best-fit number of sources of k = 23 (α 23 = 6.1 × 10 −4 ) for the time-integrated analysis.

CONSTRAINTS ON POPULATIONS OF SOURCES
In Section 1, we have described how this analysis is model-independent in that it does not rely on searching for neutrino emission from a specific astrophysical class of objects.Because of this, we can use the results of this analysis to constrain a wide variety of populations of potential astrophysical neutrino sources.
In order to constrain these populations of possible neutrino sources, we must first simulate how they would appear in the analysis.To do this, we begin by using the publicly available FIRESONG code (Tung et al. 2021).FIRESONG takes a few variables describing the population as inputs: a neutrino luminosity function characterizing the distribution of intrinsic luminosities of sources, a local number or rate density, a cosmic evolution model as a function of redshift, and an assumed spectral shape for the neutrino emission.Given these inputs, FIRESONG simulates a population of neutrino sources and calculates the neutrino flux at Earth from each of these objects.FIRESONG can simulate either transient or steady neutrino sources.We simulate transient sources when considering the populations we can constrain with the short timescale analyses, and we simulate steady sources when considering populations we can constrain with the time-integrated analysis.Unless stated otherwise, we simulate all populations assuming a spectral index of γ = 2.5, which is consistent with recent measurements of the diffuse astrophysical neutrino flux (Aartsen et al. 2020a).
Once we have a list of simulated sources, each with a flux and declination, φ l , δ l , we then determine which sources would yield an alert event.In order to do this, we calculate N alert l (δ l , γ = 2.5, φ l ) according to Equation 1.As this number is an expectation, we then randomly generate the actual number of observed alerts by sampling from a Poisson distribution with mean N alert l (δ l , γ = 2.5, φ l ) for each simulated source.We are most interested in the regime where N alert l (δ l , γ = 2.5, φ l ) is much less than one, because previous limits on neutrino source populations suggest there are more astrophysical neutrino sources than there are observed alert events from astrophysical sources (Aartsen et al. 2019a), and therefore we do not inject more than one alert event when an individual source could yield more than one alert event.However, in these cases, there are many additional events that are still injected into the GFU selection for the first alert, and thus they are still distinguishable from background-like populations in the analysis.
In Figure 5, we show one particular realization when simulating source populations whose number densities are assumed to track star-formation rates (Madau & Dickinson 2014).We simulate two different populations, where each of these populations would completely saturate the diffuse neutrino flux.The difference between the two simulated populations is in their number densities and per-source luminosities (explained below).For one population, we simulate a relatively rare population of sources (local number density of 10 −8 Mpc −3 ), and for the other, a relatively high density population of sources (local number density of 10 −6 Mpc −3 ).
In order for both of these populations to saturate the diffuse flux, which is equivalent to ensuring that the expectation for the best fit flux normalization matches the observed rate of alerts with astrophysical origin, the population with the smaller number density must have, on average, brighter neutrino sources.This is highlighted in flux distribution shown in Figure 5, as the population with fewer sources extends to higher fluxes.In addition to showing the flux distribution of every source in the population, we also show the subset of sources which resulted in an alert event, for a given realization.
After we have calculated those sources which will yield alert events in a particular realization, we then find how many events in the GFU sample we expect to observe from those sources, i.e. we calculate N GFU l (δ l , γ = 2.5, φ l ) according to Equation 1, and then fluctuate each of these numbers in the same way that we do for the alert event observations.Once we have a list of sources which caused alert events and we also have the number of additional events in the GFU sample that we should observe in a given pseudo-experiment, we inject these additional events on top of our scrambled data to perform a single pseudo-experiment.In addition to alerts from sources, there are also alerts that are the result of atmospheric backgrounds.For these alert events, we Poisson fluctuate the rates cited in Blaufuss et al. (2020) to find the number of expected alerts from atmospheric backgrounds, and we include these in our list of alert events for a psuedo-experiment, but for these alert events, we do not inject any additional events into the GFU sample, and we perform followups for these alerts as well.We calculate all analysis parameters (test statistics as well as stacked binomial p-values) to find if the The population with a smaller density (blue) has alert events that have, on average, higher fluxes.The dashed and dotted black lines show the E −2.5 fluxes needed from a source at δ = 0 • to have an expectation of detecting one alert event or one event in the large-statistics GFU sample, respectively.The higher fluxes from the blue population lead to sources in the upper-left corner of the figure, which produce events in both the alert and the GFU sample.The cumulative log(N )−log(S) distribution is given for the low density population (solid blue) and high density population (dashed green) on the right.
simulated population of sources is distinguishable from our observed data.
The motivation for this analysis lies in the differences between the flux distributions for the sources which cause alert events.Those sources which cause alert events and come from low density source populations have, on average, higher fluxes than sources which cause alert events and come from high density source populations.This means that, although we may be in the Eddington bias regime for alert events (i.e.N alert l (δ l , γ, φ l ) < 1), we may not be in the Eddington bias regime for events from the selection with a larger effective area (i.e.N GFU l (δ l , γ, φ l ) > 1).If, when we search for sources in the direction of alert events, there are additional GFU events coming from a source, then this allows us to more accurately calibrate the flux of that source.This comparison is made more explicit in Figure 6, where we show the same simulated populations of sources as we did in Figure 5.However, here we also draw attention to the flux normalizations where there is an expectation of observing exactly one event in the alert sample (GFU sample) for a source at δ = 0 • and with an E −2.5 spectrum, i.e.N alert (δ = 0 For the populations which we simulate, there is a qualitative difference between the sources which cause alerts in the low-density and high-density populations.Namely, the majority of alert events from the rare population will be accompanied by lower-energy events in the GFU sample, whereas only a small fraction of alert events will have additional detected events when the population is more numerous.From Figure 6, it becomes clear that the some low density populations can be distinguished from high density populations using the population analysis described in Section 4. This is be-cause the additional lower-energy GFU events which accompany the alerts in the low density population scenario could be identified in each of the alert followups.This, when stacked together using the population analysis, would be more signal-like than the case of the high density population.
In order to quantify how distinguishable certain populations are from our observed data, we repeat this process of injecting source populations into the analysis for a variety of different source input parameters.For each of the spots in the input source parameter space, we calculate expected distributions of binomial p-values, α, and we compare this to our observed values for both steady source hypotheses as well as transient source hypotheses.
Our resulting per-flavor limits on steady neutrino source populations are shown in Figure 7.Our per-flavor limits on transient source populations, which we calculate based off of our results from the short-timescale analyses, are displayed in Figures 8 and 9.For the transient analysis, we inject populations of transients where the emission timescale matches the analysis duration in the observer frame, and thus the ±500 s results can be used to constrain any transient sources where the true emission in the observer frame is less than or equal to 500 s, and the ±1 day timescale can be used to constrain any transient source population with true durations in the observer frame less than or equal to 1 day (we will call these minute-scale and day-scale transients, respectively).These limits are compared to the median upper limits that one could expect under the assumption of the null-hypothesis.In Figures 7, 8 and 9, the faded band represents how much the upper limit (90% CL) can typically fluctuate under the assumption of the null hypothesis.In case of an underfluctuation, previous analyses have quoted the sensitivity as the upper limit in these cases (see for e.g.Abbasi et al. 2022b).However, the binomial analyses for all three time windows analyzed here results in pre-trial p-values that are less than p=0.5, and thus we quote upper limits that are calculated using the classical Neyman approach.
The input parameter space of potential source populations is highly-dimensional, so we choose some benchmark models for these limits.Namely, we assume that source densities track star-formation rates as measured in Madau & Dickinson (2014) and we assume that all sources have a "standard candle" (SC) luminosity function, i.e. one in which all of the sources have the same luminosity and the observed flux from an individual source is solely determined by its distance.Under these assumptions, the time-integrated analysis shows that if a population of standard candle sources whose densities track star-formation rates were to be solely responsible for the diffuse astrophysical neutrino flux, then their local density must be greater than 7 × 10 −9 Mpc −3 at the 90% confidence level.For this time-integrated analysis, the fact that our observed result was more signal-like than our median expectation from the null hypothesis is reflected in the weakening of our limits with respect to the median expectation, as shown in Figure 7.For transient source populations, we find that rare populations of transient sources (rate density less than 10 −9 Mpc −3 yr −1 ), can be responsible for no more than 6% (14%) of the diffuse flux for minute-scale (dayscale) transients.We also extrapolate these limits to higher densities in the same manner as Aartsen et al. (2019a).This tells us that populations of transient sources with these same population parameters cannot be solely responsible for the diffuse flux unless their local rate densities are greater than 8 × 10 −5 Mpc −3 yr −1 (1×10 −5 Mpc −3 yr −1 ) for minute-scale (day-scale) transient sources.
We explore how our limits change by altering these assumptions, and the results of that analysis are tabulated in Table 2. Overall, we find that our limits are not extremely sensitive to changes in the assumed luminosity function.However, changes in the cosmic evolution can have a large effect on the limits.In general, the more rapidly the number of sources grows as a function of redshift, the less constraining this analysis will be on the local density, as rapidly growing source populations have a larger total number of sources at higher redshifts that could result in alert events.
Although we choose models as generic as possible for source evolution (to reflect the fact that we do not know the true properties of the source classes responsible for the diffuse flux), we also inject a more rapidly evolving population into the time-integrated analysis, to see how these limits might change for even more strongly evolving populations.For this example, we inject a toy-model that was used to describe a generic AGN-like evolution, in which the density of sources ρ(z) is described as a piecewise function of z, where z a = 1.7 and z b = 2.7.The observational study underlying this parameterization is originally from Hasinger et al. (2005) and was later reduced to this one-dimensional simplification in Stanev (2008).For this case, treating the luminosity function again as a standard candle population, we find we can only exclude Figure 7. Constraints on the per-source per-flavor neutrino luminosity between 10 TeV and 10 PeV from populations of steady neutrino sources with E −2.5 spectra, and whose densities track star-formation rates.The faded line and band show the analysis sensitivity and the 1σ (68%) expected fluctuation of a one-sided Neyman upper-limit under the null hypothesis.The data are inconsistent with a sole population of sources with the same luminosity being responsible for the diffuse flux (shown with uncertainties as the blue shaded regions) unless it has a density greater than 7 × 10 −9 Mpc −3 .The left hand side shows these constraints in the density/luminosity plane, as in Kowalski (2015), whereas the figure on the right scales the luminosity by density, which is proportional to the energy density, to focus on the most-relevant section of the parameter space.This analysis (±500 s, 9.6 yr.) Figure 8. Constraints on the per-source emitted energy per-flavor between 10 TeV and 10 PeV from populations of transient neutrino sources with E −2.5 spectra whose rate densities track star-formation rates.The faded line and band show the analysis sensitivity and the 1σ (68%) expected fluctuation of a one-sided Neyman upper-limit under the null hypothesis.Upper limits (90% CL) are inconsistent with rare populations (rate density less than 10 −9 Mpc −3 yr −1 ), of standard candle transients producing more than 6% of the diffuse flux (shown with uncertainties as the blue shaded regions) for minute-scale transients.
populations rarer than 6×10 −11 Mpc −3 , which is nearly two orders of magnitude less numerous than the local density of sources we can exclude if sources are assumed to follow star-formation rates.This is due to the greater number of sources at higher redshifts z 1 for the AGNlike evolution compared to the star-formation-like evolution.Although the luminosity function of AGN is known to not be a standard candle function, we treated the luminosity function as such for this example because we know that the effect from varying luminosity functions is second-order when compared to the density evolution.
These limits are the first limits reported by IceCube on populations of sources from searches for neutrino sources in the direction of alert events.Other popu-lation constraints, such as the ones reported in Aartsen et al. (2019a), are at a similar level for standard candle sources, although that analysis constrained sources with harder spectra.Although other neutrino source analyses tend to suffer when searching for neutrino sources with soft spectra, this analysis does not suffer as much because the ratio N GFU l (δ, γ, φ) / N alert l (δ, γ, φ) gets larger for soft spectra, as the alert stream is most sensitive at higher energies, which effectively boosts our signal when looking for correlations with alerts by using lower-energy events.When we inject harder spectra into our analysis, the constraints get weaker than those reported in Aartsen et al. (2019a).As for limits on transient source populations, the limits on the total emit- This analysis (±1 day, 9.6 yr.) Figure 9. Constraints on the per-source emitted energy per-flavor between 10 TeV and 10 PeV from populations of transient neutrino sources with E −2.5 spectra whose rate densities track star-formation rates.The faded line and band show the analysis sensitivity and the 1σ (68%) expected fluctuation of a one-sided Neyman upper-limit under the null hypothesis.Upper limits (90% CL) are inconsistent with rare populations (rate density less than 10 −9 Mpc −3 yr −1 ), of standard candle transients producing more than 14% of the diffuse flux (shown with uncertainties as the blue shaded regions) for day-scale transients.
ted energy agree within 10% with the limits reported in Aartsen et al. (2019b).However, those limits were only on sources with emission timescales less than 100 s in duration, meaning that it was effectively insensitive to transients that emit on longer timescales.Thus, the limits on transient source populations presented here are the first limits reported by IceCube on transient source populations with emission timescales greater than 100 s and up to one day in duration.

DISCUSSION & CONCLUSION
We present a search for neutrino events in the direction of IceCube neutrino candidate alert events.This search strategy is model-independent in that it does not rely on any assumptions of a specific source class being responsible for the diffuse neutrino flux.However, it is complementary to other model-independent searches that are typically accompanied by large trials-factor, be-cause we only need to search for emission coincident with the smaller statistics sample of alert quality events.
We search for neutrino emission on short timescales coincident with the alert events as well as for emission over the entire livetime of our datasets.No individual followup yields a significant result.We use these results to constrain contributions to the diffuse astrophysical neutrino flux from generic populations of sources.
Our most significant result comes from the timeintegrated population analysis, where we look for joint contributions from multiple subthreshold sources in the direction of alert events.Although consistent with background (p-value of 1.8% before accounting for the three time windows investigated), it cannot be said that there are no sources in the direction of these alert events.This is merely an indication that any sources that are in the directions of these alert events are currently at a level that is too dim to be significantly detected with the cur-rent analysis sensitivity.As of July 2021, we began using the procedure described in Section 3 to search for coincident transient neutrino emission for alerts that were detected in real time, with results frequently circulated via the Gamma-ray Coordinates Network in the hopes of identifying future sources that are producing alerts.
Future improvements might begin to reveal sources that were subthreshold to this analysis.One major improvement could come from a better localization of Ice-Cube alert events.Reducing the localization uncertainty on these events would reduce the effective background of this analysis and limit the number of locations on the sky that we need to investigate.Studies are underway to find the optimal way to construct these localization uncertainty spaces Abbasi et al. (2021c).Additionally, if sources are active in the time-domain and may emit lower-energy neutrinos at different times than we detect high-energy alert events.This would mean, an analysis that can search for this, similar as to what was done for TXS 0506+056, would be more sensitive to these fluxes.Such an analysis is in progress, with preliminary sensitivities reported in Abbasi et al. (2021d).Finally, a better treatment of the point-source likelihood could improve the per-source sensitivity to neutrino sources.Improvements to this analysis are underway, and a population analysis that utilizes these changes could prove fruitful in identifying populations of individually subthreshold neutrino sources (Bellenghi et al. 2021).

ACKNOWLEDGEMENTS
The IceCube collaboration acknowledges the significant contributions to this manuscript from Alex Pizzuto, Abhishek Desai and Justin Vandenbroucke.The authors gratefully acknowledge the support from the following agencies and institutions: USA -U.S.National Science Foundation-Office of Polar Programs, U.S. National Science Foundation-Physics Division, U.S. National Science Foundation-EPSCoR, Wisconsin Alumni Research  Table 3. Results from the individual skymap analyses for each of the time windows analyzed.In addition to providing the best-fit position from the alert event skymaps in Equatorial coordinates (J2000, all coordinates quoted in degrees), we also provide the location which yielded the largest test-statistic for the time-integrated analysis in the (RA, dec.) steady column.For the short time windows, we only provide the p-values from each analysis, and for the time-integrated analysis, we also provide the test-statistic and the best-fit parameters from the likelihood analysis.Those alerts which were excluded from the analyses, for reasons described in the text, are labeled with "Excl." in the relevant columns.Times are quoted in terms of Modified Julian Day (MJD).

Figure 1 .
Figure 1.Skymap in Equatorial coordinates (J2000) of all neutrino candidate alerts used in the analysis described in Section 3. Contours denote the 50% (solid) and 90% (dashed) containment based on rescalings of the likelihood space according to resimulations of the event IceCube-160427A (Kankare et al. 2019, see Section 2 for more details).The color indicates the signalness of each alert event, described in the text.The Galactic Plane and Galactic Center are shown as a black solid line and dot, respectively.basi et al. 2010) as well as read-out and digitization electronics (Abbasi et al. 2009).Neutrinos are detected indirectly via the Cherenkov radiation produced from relativistic charged particles created by deep inelastic neutrino nucleon interactions in the surrounding ice or nearby bedrock beneath IceCube.Although sensitive to all flavors of neutrino interactions, this study relies only on muon track events from muon-neutrino charged current interactions as well as a 10% contribution from muonic tau decays from charged current tau-neutrino interactions.These "track" events enable a better angular resolution than the other event type, "cascades" (from charged current electron-and tau-neutrino interactions or neutral current interactions of all flavors), at the cost of a poorer energy resolution.The angular resolution of track events is preferable when searching for point sources in the region of the sky where most of the alert events are detected.In addition to neutrinos from astrophysical sources, IceCube detects many neutrinos and muons from cosmic-ray interactions in the atmosphere.In the southern celestial hemisphere, the events detected by IceCube are dominated by atmospheric muons, with events from atmospheric neutrinos still occurring at rates a few orders of magnitude larger than astrophysical neutrinos.In the northern celestial hemisphere, the atmospheric muons are attenuated by the Earth, and the rate is dominated by atmospheric neutrinos.

Figure 3 .
Figure 3. Results from all of the individual alert followup analyses.Observed p-values (markers) are compared to expectations from background-only pseudo-experiments (solid lines) for the time-integrated (blue), ±500 s (green), and ±1 day (red) analyses.p-values are shown before accounting for the trials factor accrued from the look elsewhere effect.

Figure 4
Figure4.Binomial scan (left) and distribution of binomial p-values (right) calculated from pseudo-experiments using data scrambled in right ascension for the time-integrated analysis, as described in Section 4. In the binomial scan, the results for all alert followups is ordered by decreasing significance, and for each possible number of sources, k, we calculate the probability of obtaining k results each with p-values less than or equal to p k .The global minimum corresponds to k = 23 and a binomial p-value of α = 6.1 × 10 −4 (black line).The observed binomial p-value from data is shown with a black line, compared to the expected background distribution in blue, and results in an overall p-value of 1.8% for the time-integrated analysis.

Figure 5 .Figure 6 .
Figure 5. Redshift distribution (left) and flux distribution (right) for two different simulated populations of neutrino sources, both of which saturate the total diffuse astrophysical neutrino flux.The population in green (Population A) has a higher local density (10 −6 Mpc −3 ) but a smaller per-source luminosity while the blue population (Population B) has a lower density (10 −8 Mpc −3 ) but higher per-source luminosity.Both populations track the star-formation model from Madau & Dickinson (2014) and assume sources have E −2.5 spectra.For each population, those sources which result in alert events are shown in the corresponding darker colors.Both populations have the same redshift distribution in alert events, as the cumulative flux distribution as a function of redshift is the same for both populations.Alerts from the rarer population (dark blue) correspond to sources with higher fluxes than the higher density population (dark green).

Foundation,
Center for High Throughput Computing (CHTC) at the University of Wisconsin-Madison, Open Science Grid (OSG), Advanced Cyberinfrastructure Coordination Ecosystem: Services & Support (AC-CESS), Frontera computing project at the Texas Advanced Computing Center, U.S. Department of Energy-National Energy Research Scientific Computing Center, Particle astrophysics research computing center at the University of Maryland, Institute for Cyber-Enabled Research at Michigan State University, and Astroparticle physics computational facility at Marquette University; Belgium -Funds for Scientific Research (FRS-FNRS and FWO), FWO Odysseus and Big Science programmes, and Belgian Federal Science Policy Office (Belspo); Germany -Bundesministerium für Bildung und Forschung (BMBF), Deutsche Forschungsgemeinschaft (DFG), Helmholtz Alliance for Astroparticle Physics (HAP), Initiative and Networking Fund of the Helmholtz Association, Deutsches Elektronen Synchrotron (DESY), and High Performance Computing cluster of the RWTH Aachen; Sweden -Swedish Research Council, Swedish Polar Research Secretariat, Swedish National Infrastructure for Computing (SNIC), and Knut and Alice Wallenberg Foundation; European Union -EGI Advanced Computing for research; Australia -Australian Research Council; Canada -Natural Sciences and Engineering Research Council of Canada, Calcul Québec, Compute Ontario, Canada Foundation for Innovation, WestGrid, and Compute Canada; Denmark -Villum Fonden, Carlsberg Foundation, and European Commission; New Zealand -Marsden Fund; Japan -Japan Society for Promotion of Science (JSPS) and Institute for Global Prominent Research (IGPR) of Chiba University; Korea -National Research Foundation of Korea (NRF); Switzerland -Swiss National Science Foundation (SNSF); United Kingdom -Department of Physics, University of Oxford.APPENDIX A. INDIVIDUAL FOLLOWUP RESULTS TABLE

Table 1 .
Results from the binomial tests for each of the time-windows analyzed.We find all results to be consistent with background expectations.

Table 2 .
Constraints on the fractional contribution to the diffuse astrophysical neutrino flux from populations of steady (transient) neutrino sources for various densities (rate densities).Populations either track star-formation rates (SFR) or have no cosmic evolution (No Evolution).