Search for High-Energy Neutrinos from Ultra-Luminous Infrared Galaxies with IceCube

Ultra-luminous infrared galaxies (ULIRGs) have infrared luminosities $L_{\mathrm{IR}} \geq 10^{12} L_{\odot}$, making them the most luminous objects in the infrared sky. These dusty objects are generally powered by starbursts with star-formation rates that exceed $100~ M_{\odot}~ \mathrm{yr}^{-1}$, possibly combined with a contribution from an active galactic nucleus. Such environments make ULIRGs plausible sources of astrophysical high-energy neutrinos, which can be observed by the IceCube Neutrino Observatory at the South Pole. We present a stacking search for high-energy neutrinos from a representative sample of 75 ULIRGs with redshift $z \leq 0.13$ using 7.5 years of IceCube data. The results are consistent with a background-only observation, yielding upper limits on the neutrino flux from these 75 ULIRGs. For an unbroken $E^{-2.5}$ power-law spectrum, we report an upper limit on the stacked flux $\Phi_{\nu_\mu + \bar{\nu}_\mu}^{90\%} = 3.24 \times 10^{-14}~ \mathrm{TeV^{-1}~ cm^{-2}~ s^{-1}}~ (E/10~ \mathrm{TeV})^{-2.5}$ at 90% confidence level. In addition, we constrain the contribution of the ULIRG source population to the observed diffuse astrophysical neutrino flux as well as model predictions.


INTRODUCTION
The observation of high-energy astrophysical neutrinos with IceCube (Aartsen et al. 2013(Aartsen et al. , 2014 marked the birth of neutrino astronomy. Numerous studies have been performed searching for the sources of these astrophysical neutrinos, which are cosmic messengers that represent a smoking-gun signature of hadronic acceleration. So far, the blazar TXS 0506+056 (Aartsen et al. 2018a,b) is the sole neutrino source candidate that has been identified with a significance at the 3σ level, and the first indications of neutrino emission have been found from the starburst galaxy NGC 1068 (Aartsen et al. 2020a). Studies from the community have also shown indications of astrophysical neutrinos corre- * also at Università di Padova, I-35131 Padova, Italy † also at National Research Nuclear University, Moscow Engineering Physics Institute (MEPhI), Moscow 115409, Russia ‡ also at Earthquake Research Institute, University of Tokyo, Bunkyo, Tokyo 113-0032, Japan lated with the tidal disruption event AT2019dsg (Stein et al. 2021) and radio-bright active galactic nuclei (AGN; Plavin et al. 2020). However, these results are insufficient to explain the origin of the observed diffuse astrophysical neutrino flux (Aartsen et al. 2015(Aartsen et al. , 2016(Aartsen et al. , 2020bAbbasi et al. 2020a). Further dedicated IceCube studies (see Abbasi et al. 2021 for an overview of the latest IceCube searches) have investigated a wide range of candidate neutrino-source classes and found no evidence for neutrinos originating from such sources. Constraints from these IceCube studies imply that if the neutrino sky is dominated by a single source population, this population likely consists of relatively numerous, low-luminosity sources (Lipari 2008 Fermi -LAT observations provide an upper bound on the contribution of non-blazar sources to the extragalactic gamma-ray background (EGB; Ackermann et al. 2015Ackermann et al. , 2016. This bound implies constraints on the non-blazar source populations that could be responsible for the diffuse neutrino flux (e.g. Bechtol et al. 2017). In particular, the non-blazar EGB bound hints towards neutrino sources that are gamma-ray opaque Vereecken & de Vries 2020).
ULIRGs are predominantly powered by starbursts with star-formation rates 100 M yr −1 (Lonsdale et al. 2006;Rieke et al. 2009;da Cunha et al. 2010;Piqueras López et al. 2016). The strong thermal IR emission of ULIRGs is the result of abundant dust and gas reprocessing higher-frequency radiation. In addition, can show signs of AGN activity (Nagar et al. 2003;Lonsdale et al. 2006;Clements et al. 2010;Fadda et al. 2010). The fraction of ULIRGs containing an AGN is 40%, with higher values reported as a function of increasing infrared luminosity (Veilleux et al. 1995;Veilleux et al. 1999;Goto 2005;Imanishi et al. 2008;Hou et al. 2009). The AGN contribution to the infrared luminosity is typically of the order of ∼10%, although AGN can become the dominant contributors for ULIRGs with a stronger infrared output (Farrah et al. 2003;Veilleux et al. 2009;Imanishi et al. 2010;Nardini et al. 2010;Yuan et al. 2010). Lonsdale et al. (2006) also note that the AGN power could be underestimated due to strong obscuration of the AGN. This result is consistent with the work of Nardini & Risaliti (2011), who find that AGN in a sample of local ULIRGs are Compton-thick, i.e. the AGN are obscured by dust clouds with column densities N H 1.5 × 10 24 cm −2 (Comastri 2004). Observations of Arp 220, i.e. the closest ULIRG, indicate that a possible AGN hosted by this object should be Compton-thick with N H 10 25 cm −2 (Wilson et al. 2014;Scoville et al. 2017).
Hadronic acceleration, and hence neutrino production, could occur both in starburst reservoirs (e.g. Tamborra et al. 2014;Peretti et al. 2019;Ambrosone et al. 2021) and in AGN (e.g. Murase 2017;Inoue et al. 2019;Kheirandish et al. 2021). Since ULIRGs host such environments, this implies that ULIRGs are candidate sources of astrophysical neutrinos. He et al. (2013) have constructed a reservoir model of the starburst regions of ULIRGs, where the enhanced hypernova rates are responsible for hadronic acceleration. They predict a PeV diffuse neutrino flux from ULIRGs that can explain a significant fraction of the diffuse IceCube observations. The generic reservoir model of Palladino et al. (2019) considers a neutrino flux from hadronicallypowered gamma-ray galaxies (HAGS), which include ULIRGs and starburst galaxies with L IR < 10 12 L . For power-law spectra (Φ ν+ν ∝ E −γ ) with spectral indices γ ≤ 2.12, this model can partially explain the diffuse neutrino observations while remaining consistent with the non-blazar EGB bound above 50 GeV (Ackermann et al. 2016;Lisanti et al. 2016;Zechlin et al. 2016). In particular, a population of HAGS could be the dominant contributor to the diffuse neutrino flux above 100 TeV, and it could be responsible for roughly half of the diffuse neutrino flux between 10-100 TeV. In contrast with such reservoir scenarios, Vereecken & de Vries (2020) model the neutrino production in ULIRGs through a beam dump of hadrons accelerated in a Compton-thick AGN. They find that if the AGN is obscured by column densities N H 5 × 10 25 cm −2 , ULIRGs could fit the diffuse neutrino observations without violating the Fermi -LAT non-blazar EGB bound above 50 GeV (Ackermann et al. 2016).
This paper presents an IceCube search for high-energy astrophysical neutrinos originating from ULIRGs. In Section 2 we describe the source selection, the Ice-Cube dataset, and the analysis method used in this search. The results of the analysis are presented in Section 3. The interpretation of these results is given in Section 4, where we discuss the implications for the neutrino flux from the ULIRG source population, and where we compare our results to the models of He et al. (2013), Palladino et al. (2019), and Vereecken & de Vries (2020). Finally, we present our conclusions in Section 5. In this work, we assume a ΛCDM cosmology using the Planck 2015 results (Ade et al. 2016), H 0 = 67.8 km s −1 Mpc −1 , Ω m = 0.31, and Ω Λ = 0.69.

Selection of ULIRGs
The ULIRGs for this analysis are obtained from three different catalogs, primarily based on data from the Infrared Astronomical Satellite (IRAS ; Neugebauer et al. 1984), listed below. The infrared luminosity L IR listed in the catalogs is determined using the IRAS flux measurements at 12 µm, 25 µm, 60 µm, and 100 µm (see e.g. Sanders & Mirabel 1996). The three catalogs are: 1. The IRAS Revised Bright Galaxy Sample (RBGS; Sanders et al. 2003). This catalog contains the brightest extragalactic sources observed by IRAS. These RBGS objects are sources with a 60 µm infrared flux f 60 > 5.24 Jy and a Galactic latitude |b| > 5 • to exclude the Galactic Plane. The RBGS provides the total infrared luminosity between 8-1000 µm for all objects in the sample, containing 21 ULIRGs.
2. The IRAS 1 Jy Survey of ULIRGs ) selected from the IRAS Faint Source Catalog (FSC; Moshir et al. 1992) contains sources with f 60 > 1 Jy. The survey required the ULIRGs to have a Galactic latitude |b| > 30 • to avoid strong contamination from the Galactic Plane. Furthermore, in order to have accessible redshift information from observatories located at Mauna Kea, Hawaii, this survey is restricted to declinations δ > −40 • . The resulting selection is a set of 118 ULIRGs.
3. The ULIRG sample used by Nardini et al. (2010). Their selection is primarily based on the redshift survey (Saunders et al. 2000) of the IRAS Point Source Catalog (PSC; Beichman et al. 1988). This catalog contains objects with f 60 0.6 Jy and covers 84% of the sky. In addition, the authors require the ULIRGs to be observed by the Infrared Spectograph (IRS; Houck et al. 2004) on the Spitzer Space Telescope (Werner et al. 2004). As such, they obtain a sample of 164 ULIRGs.
There exists some overlap between the ULIRGs of the above three catalogs. Therefore, the NASA/IPAC Extragalactic Database 2 (NED) is used to cross-identify 2 Available at ned.ipac.caltech.edu. these sources. This cross-identification was done by obtaining the NED source list of each catalog, and then cross-identifying sources with the same NED information in each list. The result is a selection of 189 unique ULIRGs. For uniformity, NED is also used to obtain the equatorial coordinates and redshift for each object. Since L IR values depend on IRAS flux measurements, which are optimized separately for the different IRAS surveys, they are taken in the following order: • From Catalog 1 if available; • From Catalog 2 if not available in Catalog 1; • From Catalog 3 if not available in Catalog 1 or 2.
Note that Catalog 3 is the only catalog that provides uncertainties on these values. We therefore include all objects from this catalog that are consistent with L IR = 10 12 L within one standard deviation 3 .
To obtain a representative sample of the local ULIRG population, which we define as a sample of ULIRGs that is complete up to a given redshift, we place a cut on the redshifts in our initial ULIRG selection. We find that keeping only those objects with z ≤ 0.13 allows the least luminous ULIRGs (i.e. L IR = 10 12 L ) to be observed, given a conservative IRAS sensitivity f 60 = 1 Jy (see Appendix A, where we also take into account the effect of the limited sky coverage of Catalog 2). The result is a representative sample of 75 ULIRGs with z ≤ 0.13, shown in Fig. 1. Table 3 lists their NED identification, equatorial coordinates, redshifts, fluxes at 60 µm, and total IR luminosities.

Detector and Data Set
The IceCube Neutrino Observatory is a 1 km 3 optical Cherenkov detector located at the geographic South Pole (Aartsen et al. 2017b). The detector is buried in the ice at depths between 1450 m and 2450 m below the surface. It consists of 5160 digital optical modules (DOMs), which are distributed over 86 vertical strings. Each DOM contains a photomultiplier tube and onboard read-out electronics (Abbasi et al. 2009(Abbasi et al. , 2010. When a high-energy neutrino interacts with the ice or the bedrock in the vicinity of the detector, secondary relativistic particles are produced that emit Cherenkov radiation, which can be registered by the DOMs. The number of photons collected by the DOMs gives a measure of the energy deposited by these secondary charged particles. This feature combined with the arrival time of the light at the DOMs is used to reconstruct the particle trajectories, which are subsequently used to determine the arrival direction of the original neutrino.
IceCube is sensitive to all neutrino flavors, although neutrinos and antineutrinos can generally not be distinguished. Neutral-current interactions of the three flavors, and charged-current interactions of electron neutrinos (ν e +ν e ) and tau neutrinos (ν τ +ν τ ) are observed by means of their particle cascades. These cascades lead to a rather spherical pattern of recorded photons inside the detector, with a relatively poor angular resolution ( 8 • ; see Aartsen et al. 2019c). However, charged-current interactions of muon neutrinos (ν µ +ν µ ) produce muons, which can leave extended track-like signatures in the detector 4 . The typical angular resolution of these tracks is 1 • for muons with energies 1 TeV (Aartsen et al. 2017c(Aartsen et al. , 2020a. For this analysis we use the IceCube gamma-ray follow-up (GFU) data sample (Aartsen et al. 2017c), which contains well-reconstructed muon tracks with 4π sky coverage. The track selection is performed separately for both hemispheres due to their differing background characteristics. The background for astrophysical neutrinos originates from cosmic-ray induced particle cascades in Earth's atmosphere. Atmospheric neutrinos produced in such cascades are able to reach the detector from both hemispheres, leading to a background event rate at the mHz level. Compared to astrophysical neutrinos, atmospheric neutrinos have a relatively soft spectrum that dominates at energies 100 TeV (Abbasi et al. 2011;Aartsen et al. 2015Aartsen et al. , 2016Aartsen et al. , 2020b4 Tracks can also be the signature of taus produced in chargedcurrent interactions of tau neutrinos. However, such tracks only become detectable at the highest tau energies Eτ , since the decay length of the tau is roughly 50 m × (Eτ /PeV) (see also Abbasi et al. 2020b). Abbasi et al. 2020a). Atmospheric muons produced in the Northern hemisphere are attenuated by Earth before reaching IceCube. In the Southern hemisphere, however, they are able to penetrate the Antarctic ice and leave track signatures in the detector. These atmospheric muons have a median event rate of 2.7 kHz, requiring a more stringent event selection in the Southern sky. The event selection criteria result in a GFU track sample at the atmospheric neutrino level, with an all-sky event rate of 6.6 mHz. Note that this is still various orders of magnitude above the expected rate of astrophysical neutrinos at the µHz level. The selection efficiency of the GFU sample is determined by simulating signal neutrinos according to an unbroken power-law spectrum Φ ν+ν ∝ E −2 ν between 100 GeV < E ν < 1 EeV. In the Northern hemisphere, the selection efficiency is ∼50% at E ν ∼ 100 GeV and reaches ∼95% for E ν 100 TeV. In the Southern hemisphere, the selection efficiency is ∼5% at E ν ∼ 10 TeV and exceeds 70% for E ν 1 PeV. The GFU data used in this search contains over 1.5×10 6 events, spanning a livetime of 2615.97 days ∼ 7.5 years of the full 86-string IceCube configuration between 2011-2018.

Analysis Method
In this work, we search for an astrophysical component in the IceCube data which is spatially correlated with our representative sample of 75 ULIRGs. This astrophysical signal would be characterized by an excess of neutrino events above the background of atmospheric muons and atmospheric neutrinos. To search for this excess, a time-integrated unbinned maximum likelihood analysis is performed (Braun et al. 2008). In addition, the ULIRGs are stacked in order to enhance the sensitivity of the analysis (see Achterberg et al. 2006).
The unbinned likelihood for a search stacking M candidate point sources using N data events is given by Here, n s denotes the number of signal events, and γ denotes the spectral index of the signal energy distribution, which is assumed to follow an unbroken power-law spectrum, Φ ν+ν ∝ E −γ . These fit parameters are restricted to n s ≥ 0 and γ ∈ [1, 4] but allowed to float otherwise. For each event i, the probability density functions (PDFs) of the background, B i , and the signal for source k, S k i , are evaluated 5 . The parameter w k represents the stacking weight of source k.
Both the signal PDF and background PDF are separated into spatial and energy components. The background PDF is evaluated as with δ i and E i the reconstructed declination and energy of event i, respectively. The spatial component is uniform in right ascension due to the rotation of the IceCube detector with Earth's axis, resulting in a factor 1/2π. The spatial PDF, B, and the energy PDF, E B , are constructed from experimentally obtained background data in reconstructed declination and energy. Analogously, the signal PDF of candidate point source k is evaluated as For the spatial PDF of the signal, S, we take a two-dimensional Gaussian centered around the location of the source, x k . The Gaussian is evaluated using the reconstructed location of event i, x i , and its angular uncertainty, σ i , as the standard deviation. The estimation of the angular uncertainty does not take into account systematic effects, such that we follow Aartsen et al. (2020a) and impose a lower bound σ i ≥ 0.2 • . Since σ i is much larger than the uncertainty on the source location, the latter may safely be neglected. The energy PDF of the signal, E S , is modeled using a power-law spectrum with spectral index γ ∈ [1, 4], where we fit a single value of γ for the stacked flux of all ULIRGs. This modeling is motivated by the diffuse neutrino observations, which are currently best described by a power law (e.g. Abbasi et al. 2020a;Stettner 2020).
The stacking weight of source k ∈ {1, 2, ..., M } is constructed as w k = r k t k / M j=1 r j t j . Firstly, it depends on the detector response for an unbroken E −γ power-law spectrum at the source declination δ k . This dependence is reflected by the weight term r k ∝ A eff (E, δ k )E −γ dE, where A eff is the effective area of the detector (Aartsen et al. 2017c). Secondly, the stacking weight depends on a theoretical weight t k , which is modeled according to the hypothesis tested by the analysis. In this work, we make the generic assumption that the total IR luminosity between 8-1000 µm, L IR , is representative for the power of a possible hadronic accelerator in ULIRGs. This assumption is motivated by the fact that the total IR luminosity is directly related to the star-formation rate (see e.g. Rieke et al. 2009), and that a possible AGN contribution likely increases with 5 In the unbinned likelihood formulation of Eq. (1), the PDFs are assumed to be continuous (Braun et al. 2008). This continuity is approximated by choosing a PDF bin size that is well below the angular resolution ( 1 • ) and energy resolution (log(E/GeV) ∼ 0.3) for track-like events.
IR luminosity (see e.g. Nardini et al. 2010). Additionally, we assume that the neutrino production is identical in all candidate sources. Hence, we set t k = F IR,k , which is the total IR flux (8-1000 µm) of ULIRG k. The total infrared flux is computed as F IR = L IR /(4πd 2 L ), where L IR is taken from the respective ULIRG catalog (see Section 2.1), and where d L is the luminosity distance determined from the redshift measurements. We note that the theoretical weight is predominantly driven by the luminosity distance, since the IR luminosity values span less than an order of magnitude.
Subsequently, a test statistic is constructed as TS = 2 log L(n s =n s , γ =γ)/ log L(n s = 0). Here,n s andγ are those values that maximize the likelihood given in Eq. (1), where a singleγ is fitted for the full ULIRG sample. This TS is used to differentiate data containing a signal component from data being compatible with background. The background-only TS distribution is determined by randomizing the GFU data in right ascension 10 5 times, and calculating the TS in each iteration, which is called a trial. However, ∼10 8 trials are required to obtain a background-only TS distribution that is accurate up to the 5σ significance level. As this is computationally unfeasible, we fit a χ 2 PDF to the 10 5 background-only trials in order to approximate the background-only TS distribution (Wilks 1938). This PDF is then used to compute the one-sided p-value corresponding with a certain observed TS obs , which is the unique TS value associated with the unscrambled GFU data. The p-value reflects the compatibility of TS obs with the background-only scenario. This p-value is determined by evaluating the background-only survival function at TS obs , i.e. the integral of the backgroundonly PDF over all TS ≥ TS obs .
To test the performance of the analysis, we simulate muon neutrinos with a true direction exactly at the location of our 75 selected ULIRGs. The relative number of neutrinos simulated from an ULIRG location is proportional to the stacking weight of that source. These neutrinos, with energies 100 GeV < E ν < 100 PeV, are generated from an unbroken power-law spectrum, which is normalized at an energy 6 E 0 = 10 TeV. Here we assume that all sources have identical accelerator properties, such that we use a single spectral index γ to represent the flux of all ULIRGs. In such a simulation, the neutrino produces a track signature in the detector, i.e. a pseudosignal event. The angle between the reconstructed direc- E 0 = 10 TeV 90% sensitivity 3σ discovery potential 5σ discovery potential Total number of ULIRG neutrinos 90% sensitivity 3σ discovery potential 5σ discovery potential (b) Figure 2. Sensitivity, 3σ, and 5σ discovery potentials as a function of the spectral index γ, for injected pseudo-signal following an unbroken E −γ power-law spectrum. Panel (a) shows these quantities in terms of the flux at the normalization energy E0 = 10 TeV. Panel (b) shows them in terms of the mean number of injected pseudo-signal events. This number is found by convolving the injection spectrum with the effective area of the detector, and integrating it over energy and detector livetime.
tion of the track and the true direction of the neutrino creates a spread of these pseudo-signal events centered around the ULIRG locations. The number of pseudosignal events are drawn from a Poisson distribution with a given mean. For various values of this mean pseudosignal, we determine the corresponding TS distribution by performing 10 3 trials. We define the sensitivity of the analysis at 90% confidence level as the mean number of pseudo-signal events required to obtain a p-value ≤ 0.5 in 90% of the trials. In addition, we define the 3σ (resp. 5σ) discovery potential as the mean number of pseudo-signal events required to obtain a p-value ≤ 2.70 × 10 −3 (resp. ≤ 5.73 × 10 −7 ) in 50% of the trials. Panel (a) of Fig. 2 shows these quantities in terms of the stacked muon-neutrino flux evaluated at 7 E 0 = 10 TeV as a function of the spectral index γ. Panel (b) shows the same quantities in terms of the total mean number of pseudo-signal events. We determine this number by integrating the injection spectrum, after convolving it with the detector effective area, over energy and detector lifetime. We find that the sensitivity and discovery potentials are more competitive for harder spectra. This dependence is expected since harder spectra are more easily distinguished from the atmospheric background, which is well-described by an E −3.7 power-law spectrum (Abbasi et al. 2011).
In a more realistic scenario, explored by Ambrosone et al. (2021), the spectral index of the power-law spec-7 The shape of the sensitivity and discovery potentials plotted in panel (a) of Fig. 2 is a direct consequence of the evaluation of the flux at E 0 = 10 TeV.
trum may vary from object to object. For ULIRGs, such a variation could e.g. be the result of a different relative starburst-versus-AGN contribution to the neutrino flux, since the AGN contribution seems to increase with IR luminosity (Nardini et al. 2010). The stacked neutrino flux of our selected ULIRGs would therefore be a superposition of these spectra. This effect is also referred to as spectral-index blending (Ambrosone et al. 2021). In this scenario, the single spectral indexγ fitted in our analysis is therefore slightly biased towards the sources providing most of the neutrino flux at Earth. Our sensitivity in the spectral-index blending scenario will worsen slightly compared to the values given in Fig. 2. However, this reduction of the sensitivity is mitigated by the fact that the the fitted spectral indexγ as well as the sensitivity are mostly driven by the strongest neutrino sources.

RESULTS
We report the results of the stacking search for neutrino emission from our representative sample of 75 ULIRGs using 7.5 years of IceCube data, with the total IR flux of each individual source as a stacking weight. The analysis yields a best fit for the number of signal eventsn s = 0, such that the best fit for the spectral index,γ, remains undetermined. The fitted results correspond with a TS = 0, and p-value = 1. These observations are consistent with the hypothesis that the data is compatible with background. We therefore set upper limits on the stacked muon-neutrino flux of the 75 ULIRGs considered. Since the analysis yields a TS = 0, the upper limits at 90% confidence level (CL) are set equal to the 90% sensitivity, shown in Fig. 2 for unbroken E −γ power-law spectra. Table 1 lists the upper limits, Φ 90% νµ+νµ , for some specific values of the spectral index γ at the normalization energy E 0 = 10 TeV.

Limits on the ULIRG Source Population
The upper limits Φ 90% νµ+νµ of the ULIRG stacking analysis can be translated into an upper limit on the diffuse muon-neutrino flux originating from all ULIRGs with redshift z ≤ 0.13, computed as Φ z≤0.13 νµ+νµ = c Φ 90% νµ+νµ /(4π sr). Here, c = 1.1 which corrects for the completeness of the ULIRG sample (see Appendix A for more details). This result can be extrapolated to an upper limit on the contribution of the total ULIRG source population to the diffuse neutrino flux. For this extrapolation we assume that all ULIRGs have identical properties of hadronic acceleration and neutrino production over cosmic history. We follow the method presented in Ahlers & Halzen (2014) to estimate the neutrino flux of all ULIRGs up to a redshift z max as The redshift evolution factor ξ z effectively integrates the source luminosity function over cosmic history up to a redshift z. For an unbroken E −γ power-law spectrum, ξ z becomes energy-independent 8 . It can then be computed directly given a parameterization H(z) of the source evolution with redshift (Vereecken & de Vries 2020). More details and numerical values of ξ z are provided in Appendix B. The ULIRG luminosity function has been a topic of several studies Blain et al. 1999;Genzel & Cesarsky 2000;Cowie et al. 2004;Chapman et al. 2005;Le Floc'h et al. 2005;Lonsdale et al. 2006;Caputi et al. 2007;Reddy et al. 2008;Magnelli et al. 2009;Clements et al. 2010;Goto et al. 2011;Magnelli et al. 2011;Casey et al. 2012;Gruppioni et al. 2013). Qualitatively, it is observed that ULIRGs have a rapidly increasing source evolution up to a redshift z ∼ 1, followed by a relative flattening of the source evolution at higher redshifts. Unless stated otherwise, in this work we follow Vereecken & de Vries (2020) by parameterizing the ULIRG redshift evolution as H(z) ∝ (1 + z) m , with m = 4 for z ≤ 1 and m = 0 (i.e. a flat source evolution) for 1 < z ≤ 4.
Using Eq.
(2), we are now able to determine the upper limits at 90% CL on the diffuse muon-neutrino flux of the ULIRG population up to a redshift z max = 4.0, which is chosen following Palladino et al. (2019). The integral limits for unbroken E −2.0 , E −2.5 , and E −3.0 power-law spectra are shown in panel (a) of Fig. 3. Each of these limits is plotted in their respective 90% central energy region, defined as the energy range in which events contribute to 90% of the total sensitivity 9 . We find that the limits under the E −2.0 and E −2.5 assumptions exclude ULIRGs as the sole contributors to the diffuse neutrino observations up to energies of ∼3 PeV and ∼600 TeV, respectively. The E −3.0 limit mostly constrains the diffuse ULIRG flux at energies below the diffuse observations. Panel (b) of Fig. 3 shows the quasi-differential limits on the diffuse neutrino flux from ULIRGs. These quasi-differential limits are determined by computing the E −2.0 limit in each bin of energy decade between 100 GeV and 10 PeV. We find that the quasi-differential limits exclude ULIRGs as the sole contributors to the diffuse neutrino observations in the 10-100 TeV and 0.1-1 PeV bins, respectively.
The above results do not constrain the possible contribution of LIRGs (10 11 L ≤ L IR < 10 12 L ) to the diffuse neutrino flux. Within z < 1, the IR luminosity density of LIRGs evolves roughly the same with redshift as the IR luminosity density of ULIRGs, although the former is ∼10-50 times larger (Le Floc'h et al. 2005;Ca-9 The 90% central energy range of the sensitivity is determined by reducing (resp. increasing) the maximum (resp. minimum) allowed energy of injected pseudo-signal events in steps of log 10 (E/GeV) = 0.2, and computing the sensitivity in each step. The maximum energy bound (resp. minimum energy bound) is the value for which the sensitivity flux increases with 5% compared to the sensitivity flux integrated over all energies.

E ν [TeV]
10 −12 10 −11 10 −10 ULIRG evolution, z max = 4.0 HESE (7.5 yr) ν µ +ν µ (9.5 yr) E −2.0 limit ULIRGs E −2.5 limit ULIRGs E −3.0 limit ULIRGs   Casey et al. 2012). Assuming that the correlation between the total IR and neutrino luminosities holds down to L IR ≥ 10 11 L , the contribution of LIRGs to the neutrino flux is expected to dominate with respect to ULIRGs by a factor ∼10-50. Under this assumption, the diffuse neutrino observations would likely provide more stringent constraints than our stacking analysis. Nevertheless, we note that the AGN contribution to the total IR luminosity seems to be smaller for LIRGs than for ULIRGs (Nardini et al. 2010). This difference of the AGN contribution suggests that the hadronic acceleration properties of LIRGs may differ from those of ULIRGs. A dedicated study searching for high-energy neutrinos from LIRGs would provide further insights on the possible acceleration mechanisms within these objects.

Comparison with Model Predictions
First, we compare our results with the cosmic-ray reservoir model developed by He et al. (2013). They consider hypernovae as hadronic accelerators which are able to produce O(100 PeV) protons. He et al. (2013) suggest that the enhanced star-formation rate in ULIRGs leads to an enhanced hypernova rate. This enhanced hypernova rate results in ULIRG neutrino emission consistent with an E −2.0 power-law spectrum, with a cutoff at ∼2 PeV. They predict a diffuse neutrino flux originating from ULIRGs up to a redshift z max = 2.3. Hence, we compare the prediction by He et al. (2013) to our E −2.0 upper limit (90% CL) on the diffuse ULIRG neutrino flux up to a redshift z max = 2.3. This comparison is presented in panel (a) of Fig. 4. We find that our upper limit at 90% CL using 7.5 years of IceCube data is at the level of the predicted flux of He et al. (2013). A followup study using additional years of data is required to further investigate the validity of this model. Second, we consider the work of Palladino et al. (2019). Their study applies a multimessenger approach to construct a generic cosmic-ray reservoir model of hadronically powered gamma-ray galaxies (HAGS). Candidate HAGS include ULIRGs and starburst galaxies with L IR < 10 12 L . Palladino et al. (2019) model the diffuse neutrino and gamma-ray emission of a population of HAGS up to a redshift z max = 4.0 according to a power-law spectrum with an exponential cutoff. This spectrum is fitted to the diffuse IceCube observations given in Haack & Wiebusch (2018). They find that spectra with indices γ ≤ 2.12 are able to explain a significant fraction of the diffuse neutrino observations without violating the non-blazar EGB bound above 50 GeV (Ackermann et al. 2016;Lisanti et al. 2016;Zechlin et al. 2016). In panel (b) of Fig. 4 we compare the baseline HAGS model with spectral index γ = 2.12 to our E −2.12 upper limit (90% CL) on the diffuse ULIRG neutrino flux up to a redshift z max = 4.0. To deter-   Yüksel et al. 2008). We find that our limits exclude ULIRGs at 90% CL as the sole HAGS that can be responsible for the diffuse neutrino observations. However, it should be noted that this result does not have any implications for other candidate HAGS, such as starburst galaxies with L IR < 10 12 L . Last, we make a comparison with the beam-dump model of Vereecken & de Vries (2020). This model considers Compton-thick AGN, i.e. AGN obscured by clouds of matter with hydrogen column densities N H 1.5 × 10 24 cm −2 (Comastri 2004), as hadronic accelerators. The neutrino production is modeled through the pp-interactions of accelerated hadrons with cloud nuclei. Vereecken & de Vries (2020) predict a diffuse neutrino flux from ULIRGs, for which they consider two methods to normalize the proton luminosity. In one approach, it is normalized to the IceCube observations of Kopper (2018). The authors find that ULIRGs can fit the Ice-Cube data without exceeding the Fermi-LAT non-blazar EGB bound above 50 GeV (Ackermann et al. 2016) for column densities N H 5 × 10 25 cm −2 . A second approach is found by normalizing the ULIRG proton luminosity to the radio luminosity using the observed ULIRG radio-infrared relation (Sargent et al. 2010). In this case, the modeling mainly depends on the electron-to-proton luminosity ratio, f e (see e.g. Merten et al. 2017).
Here, we fit the results of Vereecken & de Vries (2020) for an E −2.0 spectrum to our E −2.0 upper limit (90% CL) on the diffuse neutrino flux from ULIRGs up to a redshift z max = 4.0. Subsequently, we use the method presented in Vereecken & de Vries (2020) to estimate the ULIRG radio luminosity using the radio-infrared relation of Sargent et al. (2010). For this estimation, we take a conservative value L IR = 10 12 L , and we assume that AGN contribute to 10% of this infrared luminosity based on the results of Nardini et al. (2010). This allows us to set a lower limit on the electron-to-proton luminosity ratio, f e 10 −3 , by fixing all other parameters in the model of Vereecken & de Vries (2020). It should be noted that among these parameters are quantities with large uncertainties, such as the electron-to-radio luminosity ratio (Becker Tjus et al. 2014). Hence, the lower limit f e 10 −3 should be regarded as an order-ofmagnitude estimation. Nevertheless, we point out that this limit on the electron-to-proton luminosity ratio lies within the same order of magnitude as the lower limits provided by Vereecken & de Vries (2020) for several obscured flat-spectrum radio AGN that were also studied with IceCube (Maggi et al. 2016(Maggi et al. , 2018.

CONCLUSIONS
ULIRGs have IR luminosities L IR ≥ 10 12 L , making them the most luminous objects in the IR sky. They are mainly powered by starbursts, possibly combined with an AGN contribution that likely increases with IR luminosity. These starburst and AGN environments are plausible hosts of hadronic accelerators, suggesting that neutrino production occurs in such environments. ULIRGs are therefore candidate sources of high-energy astrophysical neutrinos. Studies by He et al. (2013), Palladino et al. (2019), and Vereecken & de Vries (2020) suggest that the ULIRG population could be responsible for a significant fraction of the diffuse astrophysical neutrino flux observed with the IceCube Neutrino Observatory at the South Pole.
In this work, we presented a stacking search for astrophysical neutrinos from ULIRGs using 7.5 years of IceCube data. A representative sample of 75 ULIRGs with redshifts z ≤ 0.13 was obtained from three catalogs based on IRAS data Sanders et al. 2003;Nardini et al. 2010). An unbinned maximum likelihood analysis was performed with the IceCube data to search for an excess of an astrophysical signal from ULIRGs above the atmospheric background. No such excess has been found. The analysis yields a p-value = 1, which is consistent with the hypothesis that the data is compatible with background. Hence, we report upper limits (90% CL) on the neutrino flux originating from these 75 ULIRGs. The stacked flux upper limit for an unbroken E −2.5 power-law spectrum is Φ 90% νµ+νµ = 3.24 × 10 −14 TeV −1 cm −2 s −1 (E/10 TeV) −2.5 .
We studied the implication of our null result on the contribution of the ULIRG source population to the Ice-Cube diffuse neutrino observations. The integral limits for unbroken E −2.0 and E −2.5 power-law spectra exclude ULIRGs as the sole contributors to the diffuse neutrino flux up to energies of ∼3 PeV and ∼600 TeV, respectively. The quasi-differential limits exclude ULIRGs as the sole contributors to the diffuse neutrino flux in the 10-100 TeV and 0.1-1 PeV energy bins. We remark that these results do not constrain the possible diffuse neutrino contribution of the less luminous but more numerous LIRGs (10 11 L ≤ L IR < 10 12 L ), which have physical properties that are similar to those of ULIRGs (see e.g. Lonsdale et al. 2006).
Finally, we compared our upper limits at 90% CL with three models that predict a diffuse neutrino flux from ULIRGs. First, the prediction of He et al. (2013) is in tension with our E −2.0 upper limit, although a follow-up study with additional years of data is required to test the validity of this model. Second, our E −2.12 upper limit excludes ULIRGs as the sole population of HAGS, pro-posed by Palladino et al. (2019), that are responsible for the diffuse neutrino observations. We note that this result does not constrain the possible diffuse neutrino contribution from other candidate HAGS, such as starburst galaxies with L IR < 10 12 L . Third, we report a lower limit on the electron-to-proton luminosity ratio (Merten et al. 2017), f e 10 −3 , in the beam-dump model of Vereecken & de Vries (2020). This lower limit was determined by fitting the model to our E −2.0 upper limit, while fixing all other model parameters. These parameters include quantities with large uncertainties, such as the electron-to-radio luminosity ratio (Becker Tjus et al. 2014). Therefore, the lower limit on the electronto-proton luminosity ratio presented here should be regarded as an order-of-magnitude estimation. A dedicated search for high-energy neutrinos from Comptonthick AGN could provide more insights on the possible neutrino production in such beam-dump scenarios. A representative sample of the local ULIRG population is found by determining the redshift up to which the least luminous ULIRGs (i.e. L IR = 10 12 L ) can be observed, given a conservative IRAS sensitivity f 60 = 1 Jy. For this redshift determination we use the observed correlation between f 60 and the total IR flux F IR = L IR /(4πd 2 L ) of our initial ULIRG selection (see Section 2.1), which is shown in Fig. 5. Here, d L is the luminosity distance, which is determined from the redshift measurements. Since the f 60 measurements are optimized separately for the different IRAS surveys, these are taken in the following order: • From the RBGS (Sanders et al. 2003) if available; • From the FSC (Moshir et al. 1992) if not available in the RBGS; • From the PSC (Beichman et al. 1988) if not available in the RBGS or FSC.
Note that the data from the FSC and PSC are obtained from NED.
To determine the F IR value corresponding with f 60 = 1 Jy, we perform a log-linear fit on the f 60 and F IR catalog data of the form with best-fit parameters a = (95.85 ± 0.81) × 10 −2 and b = 12.645 ± 0.094. This fit is shown in Fig. 5. Subsequently, using Eq. (A1), the luminosity distance d L can be determined for which L IR = 10 12 L given f 60 = 1 Jy. We find that d L ∼ 700 Mpc, which corresponds to a redshift z ∼ 0.143. However, the uncertainties on L IR were not taken into account, since they are not provided by Catalogs 1 & 2. Therefore, a conservative redshift z = 0.13 is adopted up to which we are confident that the ULIRG selection is representative for the local ULIRG population. This conservative redshift cut results in a final selection of 75 ULIRGs with z ≤ 0.13. The redshift cut at z = 0.13 effectively corresponds with a flux constraint f 60 1 Jy. However, the final sample of 75 ULIRGs likely misses sources with Fit ULIRGs Figure 5. Correlation between the flux density at 60 µm, f60, and the total IR flux between 8-1000 µm, FIR, for our initial selection of 189 ULIRGs. The log-linear fit through these data points is also shown.
1 Jy ≤ f 60 < 5.24 Jy, called "1-Jy sources" hereafter. This lack of sources is due the fact that Catalog 2, although complete, only covers 40% of the sky. In addition, the required Spitzer observations limit the coverage of Catalog 3. The final ULIRG sample contains 37 1-Jy sources from Catalog 2, and 15 1-Jy sources from Catalog 3 which are located in the complementary 60% of the sky. Hence, we likely miss ∼40 1-Jy sources in our final ULIRG selection. The effect of these missing 1-Jy sources is estimated by computing their contribution to the cumulative stacking weight of the analysis (see Section 2.3). For this estimation, 40 1-Jy sources are simulated evenly over the part of the sky not covered by Catalog 2, where each declination will directly determine the corresponding detector response. Each of these sources is given a theoretical weight equal to the median value of the total IR flux of the 37 ULIRGs taken from Catalog 2. After repeating this simulation 10 4 times, we find that the median contribution of 40 missing 1-Jy sources to the cumulative stacking weight is roughly 10% for all spectra considered in this work. Our sample of 75 ULIRGs can therefore still be regarded as a representative sample of the local ULIRG population. We compute the upper limit on the diffuse neutrino flux from all ULIRGs within z ≤ 0.13 as Φ z≤0.13 νµ+νµ = c Φ 90% νµ+νµ /(4π sr). Here, Φ 90% νµ+νµ is the upper limit on the stacked flux from the 75 analyzed ULIRGs (see e.g. Table 1), while the completeness correction factor c = 1.1 takes into account the missing contribution of ∼40 missing 1-Jy sources.

B. REDSHIFT EVOLUTION FACTOR
A redshift evolution factor ξ z (Ahlers & Halzen 2014) was introduced to estimate the muon-neutrino flux corresponding with a certain fraction of the ULIRG source population (see Section 4.1). For an unbroken E −γ power-law spectrum, it takes the energy-independent form (e.g. Vereecken & de Vries 2020) Table 2 lists the values of ξ z for different combinations of the redshift z up to which Eq. (B2) is integrated and the spectral index γ. These values are given for three parameterizations of the source evolution H(z) ∝ (1 + z) m : • ULIRG evolution, where m = 4 for z ≤ 1 and m = 0 for 1 < z ≤ 4 (Vereecken & de Vries 2020).
• Flat evolution, where m = 0 for all z.      Note-Column (1): Primary name identification in the NASA/IPAC Extragalactic Database (NED). Columns (2)-(4): Equatorial J2000 coordinates and redshift taken from NED. Column (5): Luminosity distance determined from the redshift measurements using the cosmological parameters in Ade et al. (2016). Columns (6)- (7): Infrared flux at 60 µm and its uncertainty. Taken from the IRAS Revised Bright Galaxy Sample (RBGS; Sanders et al. 2003) if available; taken from the IRAS Faint Source Catalog (FSC; Moshir et al. 1992) if not available in the RBGS; taken from the IRAS Point Source Catalog (PSC; Beichman et al. 1988) if not available in the RBGS or FSC. Data from the FSC and PSC is obtained from NED. Columns (8)-(9): Total infrared luminosity between 8-1000 µm and the ULIRG catalog from which this value is taken.