Could Nearby Star-forming Galaxies Light Up the Pointlike Neutrino Sky?

IceCube has been measuring high-energy neutrinos from tens of TeV up to PeVs (Aartsen et al. 2020; Albert et al. 2020) during the last decade and yet their origin is still unclear. These measurements certainly reveal the presence of neutrino emitters in our universe, even though the only clear multimessenger indication of a neutrino source has been related to the blazar TXS 0506+056 (Aartsen et al. 2018). In this context, star-forming and starburst galaxies (SFGs and SBGs) are abundant sources that could be responsible for at least a part of such observations (Ambrosone et al. 2021). Indeed, SFGs and SBGs are galaxies endowed with an intense star-forming activity, and therefore they have a particularly high gas density (n_ISM ∝ 10–100 cm⁻³), strong magnetic fields (B ∼ 10² μG), and supernova explosion rates of the order of 0.01–1 yr⁻¹ (Thompson et al. 2006). For these reasons, they are supposed to exhibit strong turbulence which should be able to trap high-energy cosmic-rays inside their environment, giving them an enhanced probability to produce high-energy neutrinos and gamma-rays through hadronic collisions (Peretti et al. 2019). In this regard, many authors have studied the properties of the diffuse gamma-ray and neutrino fluxes of these sources (Loeb & Waxman 2006; Murase et al. 2013; Tamborra et al. 2014; Senno et al. 2015; Murase & Waxman 2016; Bechtol et al. 2017; Sudoh et al. 2018; Ajello et al. 2020; Peretti et al. 2020, 2021; Ambrosone et al. 2021; Blanco & Linden 2021). Different cosmic-ray transport scenarios have been used to describe the observed gamma-ray emission from SFGs and SBGs, disentangling a calorimetric behavior from a diffusion-dominated regime (Peretti et al. 2019; Krumholz et al. 2020; Müller et al. 2020; Ha et al. 2021; Owen et al. 2021; Shimono et al. 2021; Werhahn et al. 2021a, 2021b; Xiang et al. 2021b). However, in order to explain a sizable portion of measured high-energy neutrino flux with SFGs and SBGs, the deep universe must be considered up to redshift ∼4–5 because of their dimness (Ajello et al. 2020; Peretti et al. 2020; Ambrosone et al. 2021). The low gamma-ray luminosity of SFGs and SBGs typically represents a bound for their contribution to the observed astrophysical neutrinos as a pointlike component. Currently, only a dozen of these sources have been cataloged as gamma-ray pointlike sources using the Fermi-LAT data and only a few of them have been observed through imaging Cherenkov telescopes.

In this Letter, for the first time, we employ a multimessenger and multiwavelength approach to assess the ability of current and upcoming neutrino telescopes to observe such galaxies as neutrino pointlike sources. This is a question of utmost importance after the IceCube collaboration reported a 2.9σ excess of neutrino events coming from the direction of NGC 1068 analyzing 10 yr data (Aartsen et al. 2020). Even though additional active-galactic-nuclei (AGN) activity could be present, observations like this one suggest that star-forming activity could trace the hadronic emission and produce pointlike excesses in the TeV sky observed by IceCube as well as in the future skymap of the upcoming KM3NeT/ARCA and IceCube-Gen2 telescopes.

In this work, following the scenario proposed by Peretti et al. (2019) which considers the transport of high-energy protons and electrons, we analyze the 10 yr Fermi-LAT data provided by Ajello et al. (2020), taking into account the spectral energy distributions (SEDs) of 13 SFGs and SBGs. We describe these gamma-ray observations through the hadronic and leptonic processes at work in the star-forming regions and determine the corresponding neutrino fluxes. For each SBG, we require the star formation rate (SFR) to be consistent with the one derived from infrared (IR) and ultra-violet (UV) observations (Kornecki et al. 2020) within a maximal level of discrepancy, thus making our predictions more robust. Hence, we compare the most likely neutrino-flux normalizations at 1 TeV with IceCube, IceCube-Gen2, and KM3NeT sensitivities (see Figure 1). Besides this, for the brightest sources we also show in Figure 2 the prospects of the Cherenkov Telescope Array (CTA).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In this section, we summarize the emission model adopted for the neutrino and gamma-ray spectra of SFGs and SBGs. We adapt the model provided by Peretti et al. (2019; see also Peretti et al. 2020 and Ambrosone et al. 2021 for further details). We assume that the star-forming activity is limited to a small spherical region of the galaxy called the starburst nucleus (SBN). Hence, we solve the leaky-box equations to determine the energy distribution of high-energy cosmic-rays (protons and electrons) injected by supernovae explosions. The solution of the leaky-box equations can be written as

$$\begin{eqnarray}&&{F}_{p,e}={Q}_{p,e}{\left(\displaystyle \frac{1}{{T}_{\mathrm{adv}}}+\displaystyle \frac{1}{{T}_{\mathrm{loss}}}+\displaystyle \frac{1}{{T}_{\mathrm{diff}}}\right)}^{-1},\end{eqnarray} \tag{ 1 }$$

where F_p,e and Q_p,e are, respectively, the distribution function and the injection rate for protons and electrons. The cosmic-ray distribution stems from a balance of different processes with characteristic timescales: the advection timescale T_adv = R_SBN/v_wind, with R_SBN being the SBN radius and v_wind the dust-wind velocity, the energy loss timescale T_loss, and the diffusion timescale T_diff. For protons, energy losses are dominated by ionization, proton–proton, and Coulomb interactions, whereas for electrons they are dominated by ionization, inverse Compton scatterings, and synchrotron interactions (Peretti et al. 2019). For the cosmic-ray diffusion, we consider a Kolmogorov-like model with diffusion coefficient D(E) ∝ E^1/3, so that T_diff(E) ∝ E^−1/3. The diffusion timescale is typically larger than the other ones due to the large magnetic fields and the high levels of turbulence of SBGs and SFGs. This implies that the gamma-ray and neutrino emission is only marginally affected by the specific choice of the diffusion model.

We assume that protons are injected with a power-law spectrum in momentum space with spectral index Γ+2 and an exponential cutoff at 10 PeV, in agreement with the combined fit of IceCube and Fermi-LAT diffuse data (Ambrosone et al. 2021). The precise value of the cutoff energy is irrelevant for our purposes, since the gamma-ray flux is suppressed already at much smaller energies by internal and external absorption. The proton spectrum is directly proportional to the star formation rate ${\dot{M}}_{* }$ and normalized by requiring that each supernova releases into protons 10% of its total explosion kinetic energy (∼10⁵¹ erg). On the other hand, electrons have the same spectral index, a Gaussian cutoff at 10 TeV, and a normalization fixed to one-fiftieth of that of protons, similar to what is inferred for our Galaxy (Torres 2004; Peretti et al. 2019).

Neutrinos are emitted through the decay of charged pions (π → μ ν_μ , μ → e ν_e ν_μ ) that are produced in hadronic interactions of the injected protons with the interstellar gas. We determine the neutrino production rate by assuming that pions always carry 17% of their parent proton energy (Kelner et al. 2006). Gamma-rays are emitted in hadronic processes through neutral pion decays (π⁰ → γγ) and in leptonic processes through bremsstrahlung and inverse Compton scatterings of electrons. We also take into account the contribution of secondary electrons produced by the decay of charged pions. We emphasize that both the hadronic and leptonic components examined in this study are dictated by the star formation rate. Investigating the possible presence of additional components with a different origin (e.g., related to AGN activities) is beyond the scope of this work. The gamma-ray emission from inverse Compton scattering depends on the density of background photons acting as a target. In Peretti et al. (2019), such radiation is modeled as a superposition of blackbody spectra. However, as also emphasized there, we can approximate it as a monochromatic spectrum with an energy equal to the peak energy _peak = 0.1 eV. Therefore, the inverse Compton spectrum can be simply characterized by the energy density U_rad of the background photons. Analogously to Ambrosone et al. (2021), we consider all sources to have a similar spectral shape for the background photons equal to the M82 best-fit background spectrum reported by Peretti et al. (2019). The normalization for the different sources is self-consistently determined by the radiation energy density of each source U_rad (see below). Finally, we account for internal and external gamma-ray absorption due to pair-production processes with background photons. For the latter, we consider as a target CMB photons as well as the extragalactic background light model reported in Franceschini & Rodighiero (2017). We point out that gamma-ray absorption is relevant at energies larger than a few TeV, and therefore it affects only the forecast predictions for the CTA telescope.

The neutrino and gamma-ray emission mainly depends on the spectral index Γ and on the star formation rate. For the density n_ISM of the interstellar gas (target of the proton interactions), we rely on the Kennicutt relation (Kennicutt 1998; Kennicutt & Evans 2012; Kennicutt & De Los Reyes 2021). This relation connects the surface density of SFR, ${{\rm{\Sigma }}}_{\mathrm{SFR}}={\dot{M}}_{* }/(\pi {R}_{\mathrm{SBN}}^{2})$, and the gas surface density, Σ_gas = m_p n_ISM R_SBN. For simplicity, we consider for all the galaxies R_SBN = 200 pc, which is an average size of the circumnuclear regions (Peretti et al. 2019). By using the benchmark value n_ISM = 175 cm⁻³ for ${\dot{M}}_{* }=5\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ (Peretti et al. 2019), we obtain the scaling relation

$$\begin{eqnarray}&&{n}_{\mathrm{ISM}}=175{\left(\displaystyle \frac{{\dot{M}}_{* }}{5\,{M}_{\odot }\,{\mathrm{yr}}^{-1}}\right)}^{2/3}\,{\mathrm{cm}}^{-3}.\end{eqnarray} \tag{ 2 }$$

For the background photon density U_rad, we assume a direct proportionality to the SFR, which is expected to be tightly related to the infrared (IR) luminosity (Kennicutt 1998; Inoue et al. 2000; Hirashita et al. 2003; Yuan et al. 2011; Kennicutt & Evans 2012; Kennicutt & De Los Reyes 2021). In particular, we consider the scaling relation

$$\begin{eqnarray}&&{U}_{\mathrm{rad}}=2500\left(\displaystyle \frac{{\dot{M}}_{* }}{5\,{M}_{\odot }\,{\mathrm{yr}}^{-1}}\right)\,\mathrm{eV}\,{\mathrm{cm}}^{-3},\end{eqnarray} \tag{ 3 }$$

where the reference values, U_rad = 2500 eV cm⁻³ and ${\dot{M}}_{* }=5\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$, are obtained from what has been inferred for M82 by Peretti et al. (2019, 2020). Such a relation is generally satisfied by high-SFR galaxies which behave as good calorimeters. On the other hand, low-SFR galaxies have been shown to deviate from the assumed linear SFR–IR relation, as also recently pointed out by Kornecki et al. (2020), Werhahn et al. (2021b). Nevertheless, for these galaxies we have verified that Equation (3) is a good and conservative approximation: indeed the gamma-ray emission depends only slightly on U_rad since the leptonic component is in general subdominant in our model. Moreover, a smaller value for U_rad, as observed for low-SFR galaxies (Kornecki et al. 2020), would further reduce the leptonic production in favor of the hadronic one. Finally, for all the galaxies, we adopt the average values B = 200 μG (Thompson et al. 2006) and v_wind = 500 km s⁻¹ (Peretti et al. 2020) for the magnetic field and the dust-wind velocity, respectively.

Our main purpose is to determine the high-energy neutrino emission from SFGs and SBGs based on a likelihood analysis of gamma-ray data. A crucial observation must be made: for some of these galaxies, a possible additional source of gamma-rays may be related to AGN activity (Alonso-Herrero et al. 2013; Yoast-Hull et al. 2017; Guo et al. 2019; Ajello et al. 2020; Inoue et al. 2020; Murase et al. 2020; Kheirandish et al. 2021; Xiang et al. 2021a), which is not included in the following analysis. Thus we test the hypothesis that the observed gamma-ray SEDs are saturated by star-forming activity only, and determine the most likely neutrino emission under this assumption. We analyze the gamma-ray SEDs of 13 galaxies observed by Fermi-LAT with 10 yr of observation (Ajello et al. 2020). ⁴ For M82 and NGC 253 we also make use of the data provided by VERITAS (VERITAS Collaboration et al. 2009) and H.E.S.S. (H. E. S. S. Collaboration et al. 2018), respectively. For each galaxy, we pursue a Bayesian approach to assess the most likely values for the two free parameters of the model: the star formation rate ${\dot{M}}_{* }$ and the spectral index Γ of injected protons and electrons. We determine the posterior distribution as

$$\begin{eqnarray}&&p({\dot{M}}_{* },{\rm{\Gamma }}| \mathrm{SED})\propto p(\mathrm{SED}| {\dot{M}}_{* },{\rm{\Gamma }})p({\dot{M}}_{* })\,p({\rm{\Gamma }}),\end{eqnarray} \tag{ 4 }$$

with a Gaussian likelihood function

$$\begin{eqnarray}&&p(\mathrm{SED}| {\dot{M}}_{* },{\rm{\Gamma }})={e}^{-\tfrac{1}{2}\sum _{i}{\left(\tfrac{{\mathrm{SED}}_{i}-{E}_{i}^{2}{{\rm{\Phi }}}_{\gamma }({E}_{i}| {\dot{M}}_{* },{\rm{\Gamma }})}{{\sigma }_{i}}\right)}^{2}}.\end{eqnarray} \tag{ 5 }$$

Here, SED_i is the measured data, where i runs over the energy bins centered around the energy E_i , and σ_i is the observational uncertainty. The data are compared to the gamma-ray emission ${{\rm{\Phi }}}_{\gamma }({E}_{i}| {\dot{M}}_{* },{\rm{\Gamma }})$ predicted by our model. We adopt the source distances provided by Kornecki et al. (2020) and neglect the uncertainty on such values. For all the galaxies we consider the same uniform prior on the spectral index in the range 1.0–3.0. For the SFR, we instead adopt a different uniform prior distribution for each source. We require ${\dot{M}}_{* }$ to be consistent within a factor of three with the corresponding values reported in Kornecki et al. (2020). This choice is made to cover the wide variety of SFR estimates presented in the literature (Groves et al. 2008; Bolatto et al. 2011; For et al. 2012; Rahmani et al. 2016; Yoast-Hull et al. 2017; Ajello et al. 2020; Peretti et al. 2020).

The results of the analysis are summarized in Table 1. For M82, NGC 253, NGC 2146, NGC 2403, the Circinus galaxy, M31, M33, and the SMC, we find a reasonable agreement both with the data and the measured SFR values. We emphasize that in our model all the galaxies are approximated as pointlike sources with a dominant emission limited to the SBN. Hence, the SBN extension intrinsically defines the extension of the gamma-ray source emission. Most of the sources have not been resolved and therefore they are consistent with our starburst nucleus approximation. By contrast, the SMC and M31 have been observed by Ajello et al. (2020) as extended sources with an angular extension of 15 and 045, respectively. Thus, for these sources, we have also performed our analysis by assuming a different core radius in agreement with these observations (i.e., 800 pc and 3000 pc, respectively) and have obtained results similar to those reported in Figure 1. Indeed, the most likely values of the neutrino-flux normalizations change by less than a factor of two, thereby making our conclusions unchanged. We stress that our main goal is to assess the possibility of a future neutrino telescope singularly observing these sources, only using their gamma-ray data and the star formation rate as constraints for their neutrino emissions. In order to robustly constrain the parameters of these extended sources, it would be necessary to conduct a numerical simulation of their astrophysical environment, which we leave for future work.

Table 1. Results of the Likelihood Analysis of Current Gamma-Ray Data

Source	Uniform Prior	Most Likely Values	68% Credible Intervals		χ2/dof
	${\dot{M}}_{* }$	$({\dot{M}}_{* },{\rm{\Gamma }})$	${\dot{M}}_{* }$	Γ
M82	3.0–30	$\left(4.5,2.30\right)$	$\left[4.3,4.6\right]$	$\left[2.27,2.33\right]$	1.24
NGC 253	1.4–17	$\left(3.3,2.30\right)$	$\left[3.14,3.40\right]$	$\left[2.28,2.32\right]$	1.32
ARP 220	60–740	$\left(740,2.66\right)$	$\left[492,740\right]$	$\left[2.51,2.68\right]$	1.52
NGC 4945	0.35–4.15	$\left(4.15,2.30\right)$	$\left[4.05,4.15\right]$	$\left[2.23,2.32\right]$	1.52
NGC 1068	5–93	$\left(16,2.52\right)$	$\left[13,20\right]$	$\left[2.45,2.65\right]$	0.65
NGC 2146	3–57	$\left(15,2.50\right)$	$\left[9,27\right]$	$\left[2.44,2.88\right]$	0.50
ARP 299	28–333	$\left(28,2.15\right)$	$\left[28,200\right]$	$\left[1.40,1.90\right]\cup \left[2.77,3.00\right]$	0.18
M31	0.09–0.90	$\left(0.34,2.40\right)$	$\left[0.31,0.40\right]$	$\left[2.29,2.61\right]$	0.52
M33	0.09–0.90	$\left(0.44,2.76\right)$	$\left[0.19,0.56\right]$	$\left[2.57,2.96\right]$	0.44
NGC 3424	0.4–5.4	$\left(5.4,2.22\right)$	$\left[2.5,5.4\right]$	$\left[1.92,2.67\right]$	1.63
NGC 2403	0.1–1.2	$\left(0.75,2.12\right)$	$\left[0.58,0.96\right]$	$\left[1.92,2.36\right]$	0.38
SMC	0.008–0.090	$\left(0.038,2.14\right)$	$\left[0.037,0.039\right]$	$\left[2.13,2.16\right]$	1.90
Circinus Galaxy	0.1–8.1	$\left(6.6,2.32\right)$	$\left[6.2,7.8\right]$	$\left[2.15,2.45\right]$	0.92

Note. The columns report the source name, the SFR prior, the most likely values of the two parameters, the 68% maximum posterior density credible intervals of the marginal distributions, and the reduced chi-squared values considered as an estimate of the goodness of the fit. The star formation rate ${\dot{M}}_{* }$ is in units of M_⊙ yr⁻¹.

Download table as:  ASCIITypeset image

For NGC 4945, NGC 3424, and ARP 220, we find a lower limit for the SFR, implying the need for additional contributions to the gamma-ray production. Indeed these sources host an AGN (Yoast-Hull et al. 2017; Ajello et al. 2020; Xiang et al. 2021a). A special case is ARP 299, for which we find an upper limit for the SFR. Most of the gamma-ray measurements of this source are only flux upper limits, thus jeopardizing the likelihood procedure. Finally, for NGC 1068, we find a good agreement both with the data and the SFR value given by Kornecki et al. (2020). However, as will be discussed later, the corresponding neutrino emission predicted by our model is much lower than the one released by IceCube as a 2.9σ excess (Aartsen et al. 2020). Therefore, while the gamma-ray flux could well be explained by the star-forming activity of the source, the explanation of its potential neutrino emission might require an intense AGN hot corona activity (Ajello et al. 2020; Inoue et al. 2020; Anchordoqui et al. 2021; Kheirandish et al. 2021).

The results of the likelihood analysis of current gamma-ray data are hence employed to estimate the high-energy muon neutrino flux from each source. The normalization of such fluxes at 1 TeV is shown in Figure 1, where the diamonds represent the predictions for the most likely source parameters and the bands cover the 68% interval of the marginal flux distributions. In the plot, our predictions are compared with the IceCube pointlike sensitivity (Aartsen et al. 2020) as well as with those of the upcoming neutrino telescopes KM3NeT/ARCA (Aiello et al. 2019) and IceCube-Gen2 with different observation times. The latter is estimated to be at least five times better than the current IceCube one (Aartsen et al. 2021). Remarkably, we find that the cores of the SMC and the Circinus galaxy could be potentially detected by KM3NeT thanks to its higher sensitivity to the southern sky. In the northern sky, a promising source is ARP 299 which might be within reach of IceCube-Gen2. Other galaxies such as NGC 4945, M31, NGC 2403 and M82, are instead just a factor of ∼3 below the KM3NeT and IceCube-Gen2 pointlike sensitivities within their 68% credible intervals. These considerations are likely to be representative of the real detection potential for these sources. It is worth noting that the experimental sensitivities are determined for an E⁻² neutrino spectrum and therefore the conclusions might change for softer spectra. However, as reported in Table 1, the most likely spectral index for most galaxies is close to 2. Moreover, the largest values for the neutrino-flux normalization reported in Figure 1 are obtained in case of hard spectra, thus supporting the comparison with the reported pointlike sensitivities.

For the brightest neutrino galaxies, we also show in Figure 2 the most likely very-high-energy gamma-ray flux to be compared with the CTA differential sensitivity, for which we consider an observation time of 50 hr toward the direction of the galaxy (Acharya et al. 2018). The cutoff of the gamma-ray flux is due both to internal and external absorption. A few sources are expected to be above the CTA sensitivity in the range between 100 GeV and 10 TeV; for the SMC, this range extends to 100 TeV. This is particularly relevant, since one of the main limitations in our study is the lack of gamma-ray measurements in the range above 1 TeV (apart from the VERITAS and H.E.S.S. data on M82 and NGC 253). Therefore, our results suggest that CTA will allow us to draw more robust conclusions on the gamma-ray production in SFGs and SBGs. It is worth emphasizing that, under the assumption of a dominant hadronic production linked to the star-forming activity, at high energies the gamma-ray and neutrino emission are directly related to each other through the SFR which is derived from the UV and IR observations.

In the near future, upcoming neutrino telescopes will potentially observe nearby star-forming and starburst galaxies as pointlike sources. In particular, they could pose a solid link between the hadronic emission of these sources, supposed to dominate the GeV–TeV gamma-rays, and the expected intense star-forming activity as obtained from IR and UV observations.

As the brightest sources predicted by our emission model are in the southern sky, a leading role will be played by neutrino telescopes located in the northern hemisphere. They include KM3NeT/ARCA as well as the planned Baikal-GVD (Avrorin et al. 2020) and P-ONE (Agostini et al. 2020) telescopes. However, IceCube-Gen2 will have a better pointlike sensitivity in the northern sky, thus remaining a crucial observer for the sources positioned in this portion of the sky. In any case, the advent of a global neutrino network would be necessary to observe most of the pointlike emission from SFGs and SBGs predicted in this work by increasing the available worldwide effective volume of neutrino telescopes.

While our results are mainly derived from the star-forming activity, we cannot exclude a possible AGN emission counterpart for some of the galaxies selected. However, additional AGN activity would typically exhibit a flaring behavior unless the related duty cycle is very high. Therefore, the neutrino measurements over large observational periods would be dominated by the star-forming steady component, thereby making our detection prospects more robust. In this regard, a crucial role will be also played by CTA that will potentially measure the gamma-ray emission above tens of TeV energies for some of the galaxies investigated. This will allow us to place better constraints on the sources' production mechanism as well as on the parameters of the emitting cores.

This work indicates that the next decade will be decisive to assess whether the local star-forming activity can be a tracer of pointlike neutrino production and it highlights the importance of pivoting the presented scenarios through the low-energy thermal emissions.

This work was partially supported by research grant number 2017W4HA7S "NAT-NET: Neutrino and Astroparticle Theory Network" under the program PRIN 2017 funded by the Italian Ministero dell'Università e della Ricerca (MUR). The authors also acknowledge the support of research project TAsP (Theoretical Astroparticle Physics) funded by the Istituto Nazionale di Fisica Nucleare (INFN).