Neutrinos from LHAASO Sources

Abstract

We briefly review the main results of the IceCube Neutrino Observatory one decade after the discovery of cosmic neutrinos. We emphasize the importance of multimessenger observations, most prominently for the discovery of neutrinos from our own Galaxy. We model the flux from the Galactic plane produced by Galactic cosmic rays interacting with the interstellar medium and discuss the perspectives of understanding the TeV-PeV emission of the Galactic plane by combining multimessenger observations. We draw attention to the interesting fact that the neutrino flux from the Galaxy is not a dominant feature of the neutrino sky, unlike the case in any other wavelength of light. Finally, we review the attempts to identify PeVatrons by confronting the neutrino and gamma-ray emission of Galactic sources, including those observed by LHAASO. We end with a discussion of searches for neutrinos from LHAASO’s extragalactic transient source gamma-ray burst 221009A.

1 IceCube: The First Decade of High-Energy Neutrino Astronomy

Conceptualized more than half a century ago, high-energy neutrino astronomy became a reality when the IceCube project transformed a cubic kilometer of transparent natural Antarctic ice one mile below the geographic South Pole into the largest particle detector ever built [1]. IceCube discovered neutrinos in the $\rm TeV\sim PeV$ energy range originating beyond our galaxy, opening a new window for astronomy [2, 3]. The observed energy density of neutrinos in the extreme universe exceeds the energy in gamma rays [4]; it even outshines the neutrino flux from the nearby sources in our own galaxy [5]. The Galactic plane only appears as a faint glow at a level of ten percent of the extragalactic flux, in sharp contrast with what is found for any other wavelength of light, where the Milky Way is the dominant feature in the sky. This observation implies the existence of sources of high-energy neutrinos in other galaxies that are not present in our own [6].

After accumulating a decade of data with a detector with gradually improved sensitivity, the first high-energy neutrino sources emerged in the neutrino sky: the active galaxies NGC 1068, NGC 4151, PKS 1424+240 [7, 8], and TXS 0506+056 [9, 10]. The observations point to the acceleration of protons and the production of neutrinos in the obscured dense cores surrounding the supermassive black holes at their center, typically within a distance of only $10\sim 100$ Schwarzschild radii [7]. TXS 0506+056 had previously been identified as a neutrino source in a multimessenger campaign triggered by a neutrino of 290 TeV energy, IC170922, and by the independent observation of a neutrino burst from the same source in archival IceCube data in 2014 [9, 10].

The rationale for searching for cosmic-ray sources by observing neutrinos is straightforward: in relativistic particle flows—for instance, onto black holes—some of the gravitational energy released in the accretion of matter is transformed into the acceleration of protons or heavier nuclei. These subsequently interact with radiation and/or ambient matter in the vicinity of the black hole to produce pions and other secondary particles that decay into neutrinos. Both neutrinos and gamma rays are produced with roughly equal rates. While neutral pions decay into two gamma rays, $\pi^{0}\to\gamma+\gamma$, the charged pions decay into three high-energy neutrinos ($\nu$) and antineutrinos ($\bar{\nu}$) via the decay chain $\pi^{+}\to\mu^{+}+\nu_{\mu}$ followed by $\mu^{+}\to e^{+}+\bar{\nu}_{\mu}+\nu_{e}$; see Fig. 1. Based on this simplified flow diagram, we expect equal fluxes of gamma rays and muon neutrinos. The flow diagram of Fig. 1 implies a multimessenger interface between the pionic gamma-ray flux and the three-flavor neutrino flux:

$$E_{\gamma}^{2}\frac{dN_{\gamma}}{dE_{\gamma}}\approx\frac{4}{K_{\pi}}\,\frac{1}{3}E_{\nu}^{2}\frac{dN}{dE_{\nu}}\Big{|}_{E_{\nu}={E_{\gamma}}/2},$$ (1)

where $K_{\pi}$ is the ratio of charged and neutral pions produced, with $K_{\pi}\approx 2(1)$ for $pp(p\gamma$) interactions.

Figure 2 verifies this relation for $pp$ and He-$p$ interactions. We have shown that 1) the conversion between secondary neutrinos and gamma rays expressed in Eq. 1 is accurate up to $\sim 10\%$ level, and 2) in a hadronuclear interaction, the chemical composition of the primary cosmic rays barely matters as the $\gamma-\nu$ relation only relies on the production of pions [11].

This powerful relation connects neutrinos and pionic gamma rays with no reference to the cosmic-ray beam producing the neutrinos; it simply reflects the fact that a $\pi^{0}$ decays into two gamma rays for every charged pion producing a $\nu_{\mu}+\bar{\nu}_{\mu}$ pair. Also, from the fact that in the photoproduction process 20% of the initial proton energy is transferred to the pion, we anticipate that the gamma ray carries one-tenth of the proton energy and the neutrino approximately half of that. Because cosmic neutrinos are inevitably accompanied by high-energy photons, neutrino astronomy is a multimessenger astronomy.

Figure 1: Flow diagram showing the production of charged and neutral pions in $p\gamma$ interactions. The circles indicate equal energy going into gamma rays and into pairs of muon neutrinos and antineutrinos, which IceCube cannot distinguish between. Because the charged pion energy is shared roughly equally among four particles and the neutral pion energy among two photons, the photons have twice the energy of the neutrinos. Unlike the weakly interacting neutrinos, the gamma rays and electrons may lose energy in the target and will lose additional energy in the background light (EBL) before reaching Earth.

Figure 2: Energy spectra of secondary gamma rays (red) and neutrinos (blue) from the interaction between proton (dark colors)/helium (light colors) primaries and rest-mass proton targets. The primary spectrum in use is $dN/dE=E^{-2}\,\exp(-E/10\,{\rm PeV})$, with $E=E_{p}=E_{\rm He}/4$ being the energy per nucleon. We have set $E_{\rm min}=2.5\,\rm GeV$. The secondary spectra are computed using the aafragpy package [12] with cross sections from [13, 14, 15]. For comparison, the dashed curves show the neutrino flux converted from the gamma-ray flux using Eq. 1.

After accumulating a decade of data, the observed flux of neutrinos reaching us from the cosmos is shown in Fig. 3. It has been measured using two different methods to separate the high-energy cosmic neutrinos from the large backgrounds of cosmic-ray muons, 3,000 per second, and neutrinos, one every few minutes, produced by cosmic rays in the atmosphere. The first method identifies muon neutrinos of cosmic origin from tracks that are produced by upgoing muon neutrinos reaching the South Pole from the Northern Hemisphere. The Earth is used as a passive shield for the large background of cosmic-ray muons. By separating the flux of high-energy cosmic neutrinos from the lower energy flux of neutrinos of atmospheric origin, a cosmic high-energy component is isolated with a spectral index $\rm dN/dE\sim E^{-\gamma}$ with $\gamma=2.37\pm 0.09$ above an energy of $\sim 100$ TeV [16]; see Fig. 3. Also shown are the results of a second search that exclusively identifies showers initiated by electron and tau neutrinos that interact inside the detector and can be isolated from the atmospheric background to energies as low as 10 TeV [17]. Reaching us over long baselines, the cosmic neutrinos have oscillated to a $\nu_{e}:\nu_{\tau}:\nu_{\mu}$ ratio of approximately 1:1:1.

Figure 3: Current all-sky measurements of the high-energy astrophysical neutrino emission. The flux of cosmic muon neutrinos [16] inferred from the $9.5$-year upgoing-muon track analysis (solid line) with $1\sigma$ uncertainty range (shaded) is compared with the flux of showers initiated by electron and tau neutrinos [18], when assuming standard oscillations. The measurements are consistent with the expectation that each neutrino flavor contributes an identical flux to the diffuse spectrum. Both are consistent with the flux derived from the observation of a single Glashow-resonance event, shown in red.

IceCube has independently confirmed the existence of neutrinos of cosmic origin by the observation of high-energy tau neutrinos [19] and by the identification of a Glashow-resonance event where a weak intermediate $W^{-}$ boson is produced in the resonant interaction of an electron antineutrino with an atomic electron: $\bar{\nu}_{e}+e^{-}\rightarrow W^{-}\rightarrow q+\bar{q}$ [20]. Found in a dedicated search for partially contained showers of very high energy, the reconstructed energy of the Glashow event is 6.3 PeV, matching the laboratory energy to produce a $W$ boson of mass $80.37$ GeV. Given its high energy, the initial neutrino is deemed cosmic in origin and provides an independent discovery of cosmic neutrinos at the level of $5.2\sigma$.

Unlike neutrinos, gamma rays interact with photons associated with the extragalactic background light (EBL) while propagating to Earth. The gamma-rays lose energy by $e^{+}e^{-}$ pair production and the resulting electromagnetic shower subdivides the initial photon energy into multiple photons of reduced energy reaching Earth. This significantly modifies the implications of the gamma-neutrino interface, which we illustrate [4] by first parametrizing the neutrino spectrum as a power law,

$$\frac{dN}{dE_{\nu}}\propto E_{\nu}^{-\gamma_{\rm astro}},\,E_{\nu,\rm min}\leq E_{\nu}\leq E_{\nu,\rm max},\,$$ (2)

for four recent measurements taken from Refs. [21, 22, 23, 24]. Conservatively, the minimum and maximum energies are limited to the range of neutrino energies that the particular analysis is sensitive to. For each neutrino flux measurement, we derive the accompanying gamma-ray flux from the fit to the neutrino spectrum using the interface of Eq. 2. This pionic gamma-ray flux is subsequently injected into the EBL assuming that the sources follow the star-formation (SFR) history in a flat $\Lambda$CDM universe [25]. Figure 4 shows the diffuse neutrino fluxes and those of their accompanying gamma-ray counterparts. As the values of the neutrino flux spectral indices are comparable for the four measurements, the flux of the cascades is essentially determined by the value of $E_{\nu,\rm min}$ where the energy peaks. We conclude that when $E_{\nu,\rm min}\lesssim 10$ TeV, the showers initiated by the pionic photons contribute up to $30\sim 50\%$ of the diffuse extragalactic gamma-ray flux between 30 GeV and 300 GeV and nearly 100% above 500 GeV. The cascade flux exceeds the observed gamma-ray flux above $\sim 10$ GeV when $E_{\nu,\rm min}\leq 10$ TeV, assuming that the measured power-law distribution continues.

Figure 4: Gamma-ray-transparent sources? Cosmic neutrino spectra (gray curves) and gamma-ray flux resulting from the accompanying gamma rays (red curves). The data points are measurements of the diffuse cosmic neutrino flux [26, 21] and the isotropic gamma-ray background (IGRB) [27]. Neutrino spectra correspond to the best-fit single power-law models for the different IceCube analyses. Fluxes below and above the sensitivity range for the IceCube analyses are not constrained by the data and are shown as dashed curves. Pionic gamma rays from hadronic interactions are assumed to leave the sources without attenuation, propagate in the EBL, and cascade down to GeV-TeV energies.

In the end, the room left to accommodate secondary photons from TeV-PeV gamma rays and electrons is small, and it is difficult to avoid the conclusion that the extragalactic gamma-ray flux accompanying cosmic neutrinos exceeds the diffuse gamma-ray flux observed by Fermi, especially because $\sim 80\%$ of the latter is produced by blazars at the highest gamma-ray energies that are not necessarily sources of neutrinos [28]. There is no contradiction here; we infer from the calculations that the pionic gamma rays already lose energy in the target producing the neutrinos prior to propagating in the EBL. As a result, pionic gamma rays emerge below the Fermi threshold, at MeV energies or below. The observed diffuse neutrino flux originates from gamma-ray-obscured sources. The analysis presented here reinforces the conclusions of previous analyses; see Refs. [29, 30, 31, 32].

We should emphasize that the fact that powerful neutrino sources are gamma-ray-obscured should not come as a surprise. The photon and proton opacities in a neutrino-producing target are related by their cross sections (up to a kinematic factor associated with the different thresholds of the two interactions [33]),

$$\tau_{\gamma\gamma}\simeq\frac{\sigma_{\gamma\gamma}}{\kappa_{p\gamma}\sigma_{p\gamma}}\,\tau_{p\gamma}\simeq 10^{3}\,\tau_{p\gamma}\,,$$ (3)

where $\kappa_{p\gamma}\sim 0.2$ is the inelasticity in $p\gamma$ interactions. For instance, we should not expect neutrinos to be significantly produced in blazar jets that are transparent to very high-energy gamma rays. In contrast, the highly obscured dense cores close to supermassive black holes in active galaxies represent an excellent opportunity to produce neutrinos, besides providing opportunities for accelerating protons.

Recently, IceCube has identified the subdominant flux in the neutrino sky map emitted by the Milky Way. While of importance for resolving the origin of the high-energy radiation from our own galaxy, the discovery also creates new opportunities to search for the still unknown sources of Galactic cosmic rays.

2 The Discovery of Neutrinos from the Galaxy

Observing the neutrino emission from our own galaxy has been a challenge because, surprisingly, it is not a prominent feature of the neutrino sky, which is dominated by extragalactic sources. Even after combining IceCube and ANTARES data [34, 35, 36, 37], no statistically significant signal was found, achieving p-values of $\lesssim 0.02$ only. The advent of the Pass-2 data delivering a decade of superior IceCube data, and the rapid development of deep learning techniques in data science produced new opportunities to search for a Galactic component in the flux of cosmic neutrinos [38, 39, 40]. Modern neural networks allowed us to identify a larger number of neutrinos in the presence of the dominant atmospheric backgrounds and to reconstruct them with improved angular resolution. A new sample of more than 60,000 shower events achieving lower energies, improved reconstruction by a factor of two, and higher statistics by a factor of twenty made the identification of the Milky Way in the neutrino sky map possible.

Neutrino emission from the Galactic Center and the bulk of the Galactic plane originates in the southern sky, where atmospheric muons represent a challenging background for identifying neutrinos; IceCube records 100 million muons for every astrophysical neutrino observed. IceCube’s view of the Southern Hemisphere, therefore, relies on shower events starting inside the detector, which reduces the contamination of atmospheric neutrinos by an order of magnitude at TeV energies and lowers the energy threshold for observing neutrinos from more than 100 TeV to about 1 TeV. In summary, in the southern sky, the lower background, better energy resolution, and lower energy threshold of shower events compensate for their inferior angular resolution compared to secondary muon tracks. This is particularly true for searches of emission from extended objects, such as the Galactic plane.

Here, as with any search for sources in the neutrino sky map, IceCube uses the maximum likelihood method, which yields p-values that include all trials. Using this method [41], an unbinned maximum likelihood is constructed to search for spatial clustering of neutrino events. The significance of an observation of clustering is estimated by confronting the hypothesis of a source with the hypothesis that it is a fluctuation of the atmospheric background. For an event with reconstructed direction $\vec{x}_{i}=(\alpha_{i},\delta_{i})$ representing the declination and right ascension, the probability of originating from the source at $x_{s}$ is modeled as a circular two-dimensional Gaussian. The signal probability distribution function (PDF) $\mathcal{S}_{i}$ incorporates the directional uncertainty $\sigma_{i}$ and the energy $E_{i}$ for each individual event and the angular difference between the reconstructed direction of the event and the source:

$$\mathcal{S}_{i}={S}_{i}(|\vec{x}_{i}-\vec{x}_{s}|,\sigma_{i})\,\mathcal{E}_{i}(\delta_{i},E_{i},\gamma),$$ (4)

where the spatial distribution is modeled as a two-dimensional Gaussian

$$S_{i}(|\vec{x}_{i}-\vec{x}_{s}|,\sigma_{i})=\frac{1}{2\pi\sigma_{i}^{2}}\exp\Big{(}{-\frac{|\vec{x}_{i}-\vec{x}_{s}|^{2}}{2\sigma_{i}^{2}}}\Big{)}$$ (5)

and $\gamma$ is the spectral index. The energy distribution $\mathcal{E}_{i}(\delta_{i},E_{i},\gamma)$ represents the probability of observing the reconstructed muon energy $E_{i}$ given the source spectral index $\gamma$ assuming a power-law energy spectrum. A harder spectrum relative to the background can also indicate a signal.

The background PDF, $\mathcal{B}_{i}$, is assumed to be atmospheric in origin and can be calculated. Here we assume that the data in each declination/zenith range are background-dominated and the total event rate is a good approximation to the background.

The likelihood function for a point source is defined as the sum of events $N$ originating from signal and background:

$$\mathcal{L}(n_{s},x_{s},\gamma)=\prod_{i}^{events}\bigg{(}\frac{n_{s}}{N}\mathcal{S}_{i}(|x_{i}-x_{s}|,\sigma_{i},E_{i},\gamma)+\frac{N-n_{s}}{N}\mathcal{B}_{i}(\delta_{i},E_{i})\bigg{)}.$$ (6)

After determining the best fit for the number of signal events $\hat{n}_{s}$ and their spectral index $\hat{\gamma}$, the TS is calculated as

$${\rm TS}=2\log\big{[}\frac{\mathcal{L}(\hat{n}_{s},\hat{\gamma})}{\mathcal{L}(n_{s}=0)}\big{]}.$$ (7)

The significance of an observation is determined by comparing the maximized TS to the TS distribution of the background, usually obtained from data randomized in right ascension. For large sample sizes, this distribution approximately follows a chi-squared distribution, where the number of degrees of freedom corresponds to the difference in the number of free parameters between the null hypothesis and the alternate hypothesis.

The search can be modified to search for extended sources of neutrinos. For this purpose, the angular uncertainty in the spatial PDF is modified to include the extension of the source:

$$\sigma_{\rm eff}=\sqrt{\sigma_{i}^{2}+\sigma_{\rm ext}^{2}}.$$ (8)

The technique described above can be extended to any hypothesis that may represent multiple sources or any extended feature in the sky, such as the Galactic plane.

IceCube constructed a full-sky template of the Galactic plane using the GeV-energy observations by Fermi. Assuming that the observed gamma-ray flux consists of photons that are the decay products of pions produced by Galactic cosmic rays colliding with the interstellar medium, it can be transformed into a neutrino template for the Galactic plane using the multimessenger interface described above. This template is subsequently binned in equal-solid-angle bins, convolved with the detector acceptance, and smeared with a Gaussian with a width equal to each event’s angular uncertainty. Each event, therefore, has a particular spatial weight based on its position, energy, direction, estimated angular uncertainty, and the signal hypothesis under investigation. These combined spatial and energy weights are multiplied and the likelihood maximized for the analysis parameters described above.

The construction of the template requires modeling the Fermi data. Three such sky templates were used in the search using cascades. The first is a GALPROP model [42] calibrated with the Fermi observation, referred to as the $\pi^{0}$ template and illustrated in Fig. 5. The other two are KRA models that include radially dependent cosmic-ray transport properties. In particular, IceCube used the KRA models with a cosmic-ray cut-off at 5 and 50 PeV [43], denoted as the KRA${}_{\gamma}^{5}$ and KRA${}_{\gamma}^{50}$ models, respectively. In addition to these three templates, the IceCube search using tracks [44] used the CRINGE template, based on a global fit to cosmic-ray observations [45] and analytical models from [46].

Figure 5: The plane of the Milky Way galaxy in photons and neutrinos. Each panel from top to bottom shows the flux in latitude and longitude: (A) Optical flux. (B) The integrated flux in gamma rays from the Fermi-LAT 12-year survey at energies greater than 1 GeV. (C) The neutrino emission template derived from the $\pi^{0}$ template that matches the Fermi-LAT observations of the diffuse gamma-ray emission. (D) The emission template from panel (C) including the IceCube detector sensitivity to cascade-like neutrino events and the angular uncertainty of a typical signal event (7 degrees, indicated by the dotted white circle). Contours indicate the central regions that contain 20% and 50% of the predicted diffuse neutrino emission signal. (E) The pre-trial significance of the IceCube neutrino observations, calculated from the all-sky scan for point-like sources using the cascade neutrino event sample. Contours are the same as in panel (D). For more details, see Ref. [47].

The Galactic template search was performed with the same methods as previous Galactic diffuse emission searches [34, 35]. The analysis has only one unconstrained parameter, the number of signal events $n_{s}$ in the entire sky, providing the flux normalization while keeping the spectrum fixed to the template. The search with cascades rejects the background-only hypothesis with a significance of $4.71\sigma$, which is reduced to $4.5\sigma$ because of trials related to other searches performed; for details on this and the maximum likelihood method, we refer to the supplementary material in Ref. [5]. The result is shown in  Fig. 6. The search with tracks identified a Galactic contribution in the northern sky at the $2.7\sigma$ level with a consistent flux.

Figure 6: The energy spectra of neutrinos from the Galactic plane measured using three templates of the diffuse neutrino emission: the $\pi^{0}$, KRA${}_{\gamma}^{5}$ and KRA${}_{\gamma}^{50}$ models (from [5]). Dotted curves indicate the model flux while solid curves and shaded regions indicate the measured flux and $1\sigma$ uncertainties, respectively. The hatched region shows the IceCube all-sky diffuse flux [48].

The ANTARES collaboration searched for neutrino emission from the Galactic ridge region, defined as $|l|<30^{\circ}$ and $|b|<2^{\circ}$ [49]. By comparing the flux of neutrino events in the “on” region, i.e., the Galactic ridge, to that in a off-zone region and describing the neutrino spectrum as a power law, Ref. [49] identified a mild excess of neutrino emission from the Galactic ridge region at the level of $2.2\,\sigma$. The measured neutrino flux is consistent with the model of diffuse gamma rays from the $\pi^{0}$ decays from the interaction of protons with the ISM in the Galactic ridge. The best-fit spectral index of neutrinos is $dN/dE_{\nu}\propto E_{\nu}^{-2.45}$, suggesting a harder cosmic-ray spectrum in the inner Galaxy.

3 Modeling the Diffuse Neutrino Flux from the Galactic Plane

The Galactic diffuse emission (GDE) in gamma rays (GDE-$\gamma$) has been measured by Fermi-LAT from 100 MeV to 1 TeV [50]. It has been observed by ground-based gamma-ray experiments above $\sim$1 TeV [51, 52], from the parts of the Galactic plane that are accessible to the detectors, and above 100 TeV by the Tibet AS$\gamma$ [53], HAWC [52], and LHAASO-KM2A [54] observatories. Most of the photons observed by Tibet AS$\gamma$ with energy above 398 TeV do not point to known gamma-ray sources, indicating that the emission could be diffuse in nature. On the other hand, LHAASO-KM2A finds a lower GDE intensity for a similar sky region when masking known sources detected by LHAASO [54]. The comparison suggests that a significant fraction of the Tibet AS$\gamma$ flux could be contributed by unresolved sources [55].

The flux of pionic gamma rays that accompanies the flux of high-energy neutrinos from the Galactic plane is consistent with the Galactic diffuse gamma-ray emission observed by the Fermi-LAT and Tibet AS$\gamma$ experiments around 1 TeV and 0.5 PeV, respectively; see Fig. 7. This consistency suggests that the GDE-$\gamma$ could be dominantly produced by hadronic interaction above $\sim$1 TeV. In addition, the IceCube flux is comparable to the sum of the flux of individual sources not associated with pulsars and the LHAASO diffuse emission above $\sim$30 TeV (as shown in Figure 4 of [56]). This implies that the LHAASO GDE and a significant portion of the resolved and unresolved sources may dominantly originate from hadronic interactions.

Figure 7: The neutrino flux of the Galactic plane. The neutrino flux is shown in blue with a spread that includes statistical and systematic errors. The accompanying gamma-ray flux shown in red matches the lower energy Fermi data. Both are consistent with a model (dashed lines) determining the number of parent pions produced in interactions of Galactic cosmic rays interacting with the interstellar medium. The neutrino flux observed by IceCube is consistent with the one observed by the Tibet AS$\gamma$ array at the higher energies after its conversion to neutrinos [53]. The diffuse extragalactic gamma-ray (red hatched) and neutrino flux (blue hatched) are shown for comparison.

A priori, the IceCube Galactic plane flux may be explained as:

1. Unresolved emission by individual neutrino sources.

2. The truly diffuse emission produced when cosmic rays confined by the Galactic magnetic field, the “cosmic-ray sea,” interact with the interstellar medium.

3. A combination of both components.

In Scenario 1, comparing the neutrino flux converted from the gamma-ray flux of Galactic sources observed by various imaging atmospheric Cherenkov telescopes (IACTs) and air shower gamma-ray detectors with the IceCube results found that the gamma-ray sources are insufficient to accommodate the IceCube Galactic plane flux, especially considering that gamma rays may also be produced by electrons [56, 57, 58]. Therefore, sources that are currently unresolved by gamma-ray facilities need to be present to explain the IceCube observation. Particularly interesting candidate sources include young massive star clusters [59], hypernova remnants [60, 46], and giant molecular clouds [61].

Scenario 2 was long predicted by Refs. [43, 62, 46, 63, 64], for example, and more recently employed by Refs. [6, 65, 66] to explain the IceCube observation. Below we present a model illustrating this idea, realizing at the same time that given the uncertainty of the GDE-$\nu$ flux associated with the analysis templates and the potential contribution from unresolved neutrino sources, leptonic processes may still play a role, in particular between 0.1 and a few TeV.

We first model the cosmic-ray density in the Galaxy and subsequently calculate the production of photons and neutrinos in the interactions of cosmic rays with the interstellar gas.

Figure 8: Cosmic-ray proton (top curve) and helium (bottom curve) fluxes at the Earth between 1 GeV and $10^{8}$ GeV scaled by $E^{2.7}$ (adapted from Ref. [6]). The data points are from measurements by AMS-02 [67, 68], PAMELA [69], ATIC [70], CALET [71, 72], CREAM-III [73], DAMPE [74], NUCLEON [75], KASCADE [76], and IceTop [77] experiments. The errors indicate the quadratic sums of statistical and systematic uncertainties of the measurements.

In the energy range of interest, Galactic cosmic rays are dominated by proton and helium nulcei. Their spectral shape can be described by a modified power law [67] that cuts off at low and high energies and describes the break near a rigidity of 300 GV (see Fig. 8). The cosmic-ray density in the Galaxy is computed using the public simulation framework CRPropa 3.1 [78]. Stochastic differential equations are solved to describe particle propagation in the Galactic magnetic field described by the JF12 model, including both the regular and turbulent components [79]. The sources are assumed to follow the distribution of isolated radio pulsars [80, 81], which trace the spiral structure of the Galaxy. Modeling the spiral-arm source density produces a higher cosmic-ray energy density than a smooth disk when the propagated intensities for each are similarly normalized to the cosmic-ray data.

The model assumes that the spatial and spectral components of the nucleon flux are independent. This assumption could be oversimplified and does not describe the Fermi-LAT observations [82]. Variations of the cosmic-ray spectra across the Galaxy also affect the more detailed modeling of the GDE-$\nu$ flux [83]. For a detailed description of the model, we refer the reader to the supplementary material of Ref. [6].

Figure 9: Intensities of Galactic diffuse emission in gamma rays and neutrinos from two sky regions, region A (left): $25^{\circ}<l<100^{\circ}$, $|b|<5^{\circ}$, and region B (right): $50^{\circ}<l<200^{\circ}$, $|b|<5^{\circ}$ (from [6]). The data points refer to observations of various gamma-ray experiments, including Tibet AS$\gamma$ [53] (orange data points; the fainter data points at 500 TeV indicate the residual intensity after removing events near the Cygnus Cocoon), LHAASO-KM2A [84] (yellow triangle markers), Fermi-LAT [85] (pink circular markers indicating the data of the corresponding region), ARGO-YBJ [86] (green circular markers in region A), and CASA-MIA [87] (magenta square markers in region B). The error bars show $1\,\sigma$ statistical errors. The gamma-ray contribution (red) is decomposed into $\pi^{0}$ decay (grey dashed), inverse Compton (grey dash-dotted) and bremsstrahlung emission (grey dotted), with thick and thin curves indicating products from interaction with molecular (H₂) and neutral (HI) gas, respectively. The blue curves indicate the corresponding per-flavor neutrino intensities. The gamma-ray and neutrino spectra presented in this figure are from a simple model where the spatial and spectral components of the nucleon flux are assumed to be independent.

The emissivity and intensity of secondary production is subsequently computed using the HERMES code [88], with hadronic interaction cross section from Ref. 89 and the local gas density profiles from Refs. [90, 91, 92]. The calculation accounts for the cosmic-ray interactions with the neutral and molecular hydrogen assumed to have a uniform abundance ratio $n_{\rm He}/n_{\rm H}=0.1$ traced by the CO emission.

Our model, shown as the blue dashed curve in Fig. 7, agrees with the measurements using the $\pi^{0}$ and CRINGE models [93, 44] but exceed the fluxes predicted by the KRA models [47]. We note that the model presented in this work is designed to explain the IceCube $\pi^{0}$ flux as diffuse emission. When changing our source density distribution from a spiral arm model to a smooth disk model, the prediction of the secondary flux would be lowered by a factor of $\sim 2$ and comparable to the $\pi^{0}$ and KRA models at $\sim 100$ GeV.

Figure 9 further shows the pionic gamma-ray emission and the inverse Compton and synchrotron radiation of the subdominant electrons in two sky regions well observed by gamma-ray facilities. In region A, the model accommodates the IceCube and the Tibet AS$\gamma$ GDE observations but overproduces the LHAASO-KM2A measurement, which is a factor $\sim$2-3 times lower than the Tibet AS$\gamma$ measurement. The difference in multimessenger observations could be attributed to a contribution of hadronic sources that are not resolved by the Tibet AS$\gamma$ and IceCube observations. The source contribution, on the other hand, should not alter the conclusion that the hadronic emission dominates the Galactic plane emission above $\sim 30$ TeV, given the agreement between the IceCube and Tibet AS$\gamma$ fluxes.

The longitudinal and latitudinal distributions of our GDE-$\gamma$ model are roughly consistent with Fermi-LAT observations above 100 GeV. Nevertheless, this simplified model is not tailored to fit their data, known to require multiple emission components [82]. The light blue bands in Fig. 9 show the average intensity in the two Tibet regions of the $\pi^{0}$ template scaled by the full-sky normalization observed by IceCube. The model produces higher flux in region A and lower in region B, suggesting that its spatial distribution differs from the $\pi^{0}$ template.

The modeling of diffuse neutrino and gamma-ray emission is affected by several poorly known factors, including 1) cosmic-ray spectra above the rigidity $\sim$10 TV observed at the Earth, 2) the spatial dependence of the cosmic-ray density, which is determined by the source distribution in the Galaxy, timescales of the sources, and the particle diffusion in the Galactic magnetic field, and 3) the ISM gas density profile. The effects of these factors could be coupled and are not constrained by current observations [94, 95, 85, 93, 6, 93].

4 The Galaxy is a Neutrino Desert

In Fig. 7, we have also shown the measurements of the extragalactic diffuse flux of photons and neutrinos measured by Fermi and IceCube, respectively. It is striking how the neutrino flux from the Milky Way is subdominant to the aggregate flux of extragalactic sources. While the Galactic component of the neutrino flux is at the level of 10% of the extragalactic flux at $30\,\rm TeV$, it is the dominant feature in the gamma-ray flux. This contrast is expressed explicitly by the following relation:

$$\rm\frac{L_{\nu}^{EG}}{L_{\nu}^{MW}}\sim 120\,\left(\frac{\Phi_{\nu}^{EG}/\Phi_{\nu}^{MW}}{5}\right)\left(\frac{n_{0}}{0.01\,Mpc^{-3}}\right)^{-1}\left(\frac{\xi}{3}\right)^{-1}\left(\frac{F_{\epsilon}}{1}\right)\,,$$ (9)

which represents the ratio of the neutrino flux from extragalactic sources with a density $n_{0}$ and luminosity $\rm L^{EG}_{\nu}$,

$$\Phi^{\rm EG}_{\nu}=\rm n_{0}ct_{H}\,L^{EG}_{\nu}\,\,\,\xi\,,$$ (10)

and the Galactic neutrino flux with luminosity $\rm L^{MW}_{\nu}$,

$$\Phi^{\rm MW}_{\nu}=\rm\frac{3}{4\pi r_{0}^{2}}\,\,\,L_{\nu}^{MW}\,\,\,F_{\epsilon}\,.$$ (11)

In the above equations, $\rm t_{H}$ is the Hubble time and $r_{0}$ the distance to the center of the Galaxy; $\xi$ and $\rm F_{\epsilon}$ are geometric factors of order unity related to the cosmic evolution of the sources and the detailed geometry of our Galaxy, respectively. Eq.11 relates the neutrino flux of the Milky Way to the locally observed luminosity.

In Eq.9, we have assumed the source density for galaxies with a stellar mass similar to that of the Milky Way. It is clear that the extragalactic neutrino flux dominates by two orders of magnitude, even for the conservative assumption of an observed ratio of $\rm\Phi_{\nu}^{EG}/\Phi_{\nu}^{MW}\simeq 5$. We conclude that neutrino sources that do not exist in our own galaxy do exist in other galaxies. Given that the initial neutrino sources detected by IceCube are active galaxies, one may speculate that the absence of any major activity of the black hole in our own galaxy, at least for the last few million years, is responsible for the missing neutrinos.

5 Multiwavelength Search for Galactic Neutrino Sources

The softening of the cosmic-ray spectrum from $dN/dE\propto E^{-2.7}$ to $E^{-3.1}$ at several PeV, shown in Fig. 8, is known as the cosmic-ray “knee” [96]. The presence of the knee suggests that the Galactic cosmic-ray sources accelerate protons up to energies of at least a few PeV, a.k.a. “PeVatrons.” These cosmic rays produce neutrinos of 50-TeV energy and higher when interacting with the ambient matter and photons in the vicinity of their sources, forming individual, and possibly extended, Galactic neutrino sources. Confined by the Galactic magnetic field, these cosmic rays will interact with interstellar gas, forming the Galactic diffuse neutrino emission.

Identifying the Galactic sources that accelerate cosmic rays has been one of the major tasks of gamma-ray astronomy [97, 98]. Although the spectral break in the gamma-ray spectrum at sub-GeV energies may present evidence for the acceleration of cosmic rays, it may also be caused by re-acceleration of diffuse cosmic rays for a limited time period. Thus, the signature does not necessarily indicate the acceleration of freshly injected protons [99, 100].

Finding unequivocal evidence for the acceleration of protons through gamma-ray observation has been challenging. Two approaches have been pursued to identify hadronic, also referred to as pionic, gamma-ray emission.

The first approach is to directly detect the characteristic pion-decay feature at sub-GeV energies in the gamma-ray spectra. The two photons from the decay of the neutral pion ($\pi^{0}\rightarrow 2\gamma$) have an isotropic decay spectrum in the rest frame of the neutral pion, leading to a characteristic spectral feature referred to as the pion-decay bump in the vicinity of 200 MeV [101]. Using data from the Fermi-LAT, the identification of this so-called pion bump has been achieved in the spectra of the supernova remnants (SNRs) W44, IC 443 [102], W49B [103] and W51C [104]. More recently, a comprehensive search [105] for such a feature among 4FGL sources within $5^{\circ}$ from the Galactic plane identified 56 sources with a significant spectral break. This population includes 13 SNRs, 6 sources of unknown nature but overlapping with known SNRs or PWNe, 4 binaries, 2 PWNe, 1 star-formation region, the Cygnus Cocoon, along with 30 unidentified sources of unknown origin.

The second method is to use multiwavelength observations, such as the coincidence of the spatial locations of molecular clouds and gamma rays, the correlation in temporal profiles of gamma-ray and lower-energy emission of transients, and the broadband spectrum in general, to infer the origin of the gamma-ray emission. The following classes of gamma-ray sources have been the subject of intense study and have yielded promising candidate sources of cosmic-ray accelerators.

• The Galactic Center. The observation of gamma rays originating close to Sgr A* with a photon index of $\sim 2.3$ extending to beyond 10 TeV suggests the existence of a PeVatron in the Galactic Center [106, 107, 108, 109]. The gamma-ray emission may be explained by protons accelerated by the black hole Sgr A*, a nearby SNR, or the compact stellar clusters at the Galactic Center that are injected into the central molecular zone. The Galactic Center region has recently been detected by HAWC at yet higher energies [110].

• Star-formation regions and the Cygnus Cocoon. Stellar clusters, composed of massive young stars, may plausibly accelerate protons to very high energy (VHE) through the stellar winds or the shocks of supernova remnants expanding inside the bubble [111, 112, 113]. Extended gamma-ray emission with a radius of more than $2^{\circ}$ was detected from the star-forming region of Cygnus X by Fermi-LAT at 1-100 GeV, referred to as the Cygnus Cocoon [111, 114]. The gamma-ray flux cannot be explained by the Galactic diffuse emission, suggesting that either cosmic-ray electrons or protons are freshly accelerated inside the Cocoon. HAWC identified 1-100 TeV gamma rays from the Cygnus Cocoon, with a spatial morphology matching that observed at GeV energies remarkably well [115]. The gamma-ray spectrum of the Cocoon softens from $dN/dE\propto E^{-2.2}$ at 1-100 GeV to $E^{-2.64}$ at 1-100 TeV, which may indicate the termination of the accelerator or the leakage of the highest energy particles. The matching morphology of the GeV and TeV emission, the spatial profile of the gamma rays, as well as the absence of diffuse X-ray emission from the Cocoon [116, 117] all point to a hadronic source.

  The Cygnus region has also been observed by LHAASO [118]. After removing all sources from the so-called Cygnus Bubble, which they define as a region of $6^{\circ}$ radius from the center, the residual emission can be fitted with a Gaussian template with a width of $2.28^{\circ}$ and $2.17^{\circ}$ using the WCDA and KM2A data, respectively. This flux of residual emission, named LHAASO J2027$+$4119, and that from the entire Cygnus Bubble is compared to the Fermi-LAT and HAWC Cygnus Cocoon data in Fig. 10. In the extended vicinity of the Cygnus Cocoon, referred to by LHAASO as the Cygnus-X region, which is a $10^{\circ}$-radius extended region around Cygnus-X, photons with extreme energies up to 1.4 PeV have been observed [53, 118], hinting at the existence of a “superPeVatron.”

  In addition to the Cygnus Cocoon and the potential stellar cluster region at the Galactic Center, Westerlund 1 and Westerlund 2 are also detected at TeV energies [119, 120]. More candidate gamma-ray emitting stellar clusters have been identified by comparing theoretical models to compact clusters found by Gaia [121].

• Supernova remnants. While SNRs represent the dominant population of GeV gamma-ray emitters with a pion bump feature, most of them either are undetected or have a steep spectrum at VHE of 0.1-100 TeV, challenging the theory that SNRs are the accelerators of cosmic rays up to the knee [112, 113]. The scarcity of PeVatron candidates makes SNR G106.3+2.7 a particularly interesting case. It is one of the nearest middle-aged SNRs, at $\sim 800$ pc. Its emission above 10 TeV has been observed by the air shower gamma-ray detectors HAWC [122], Tibet [123], and LHAASO [124] and at multi-TeV energies by VERITAS [125] and MAGIC [126]. The direction of the highest energy photons is spatially coincident with a molecular cloud [123]. Also, the nondetection of 1-10 GeV counterparts to the VHE emission constrains the bremsstrahlung and inverse Compton emission by relativistic electrons, supporting a PeVatron scenario [127].

• Novae. An outburst of RS Ophiuchi (RS Oph), a recurrent nova with a red giant companion, was observed by the LAT, MAGIC, and H.E.S.S. telescopes [128, 129]. The similarity of the spectra and light curves at GeV and VHE energies favors a hadronic origin over the leptonic alternative. The detection of time-correlated optical and gamma-ray emission in the classical novae V906 Carinae revealed the role of radiative shocks powering cosmic-ray acceleration [130]. These observations suggest that novae as well as other types of transients powered by non-relativistic shocks could be promising neutrino emitters [131, 132, 133].

The detection of gamma-ray sources above 100 TeV, referred to as ultrahigh-energy (UHE) gamma-ray sources, provides opportunities to study both leptonic and hadronic gamma-ray emitters at the highest photon energies. To date, a total of 43 UHE gamma-ray sources have been detected at the $>4\sigma$ level [134, 135]. A good fraction of these sources are found to be associated with pulsars, challenging the traditional paradigm that $>100$-TeV photons imply the presence of hadronic PeVatrons [136].

The similarity in the spectra of pionic gamma rays and leptonic gamma rays from bremsstrahlung and inverse Compton scattering of relativistic electrons as well as the lack of multiwavelength counterparts are among the main challenges in identifying hadronic gamma-ray emitters. High-energy neutrinos, which predominantly originate from the decay of charged pions produced in hadronic interactions, provide a smoking gun for identifying the hadronic accelerators in the Milky Way. But this is also a challenging road, leaving the origin of the Galactic cosmic rays as one of the oldest puzzles in astronomy. We discuss this next.

6 Neutrino Search for HAWC and LHAASO Sources

Unlike the challenging but successful identification of the Milky Way, the search for neutrinos from the sources of the Galactic cosmic rays has come up empty so far. However, the multimessenger data suggest the presence of a component of unresolved sources in the Galactic neutrino flux observed by IceCube at the highest energies. This has led to a renewed effort to identify them.

Building on a decade-long effort to identify the Galactic cosmic-ray sources, by matching gamma-ray and neutrino spectra in order to reveal their pionic origin [137], initiated with the early IceCube data, a list of target Galactic sources has emerged in IceCube’s a priori defined source list, including 2HWC J2031$+$415, Gamma Cygni, MGRO J2019$+$37, MGRO J1908$+$06, HESS J1857$+$026, HESS J1852$-$000, HESS J1849$+$000, HESS J1843$-$033, PSR B0656$+$14, the Crab nebula, HESS J1841$-$055, and HESS J1837$-$069 [138]. Additionally, source lists of supernova remnants, pulsar wind nebulae, and unidentified objects are used for IceCube’s stacked catalogue searches [35, 138, 5]; see also Table S4 in the supplementary material of Ref. [47]. Some of the stacked Galactic source searches using the shower data sample that revealed the Galactic plane yield evidence exceeding the $>3\,\sigma$ level. However, this evidence cannot be separated from that for the diffuse emission template of the Galactic plane.

In addition to these sources, several dedicated searches for Galactic sources have been carried out starting with the early MILAGRO sources. MGRO J1908$+$06 has been one of the prominent candidate sources [139]. It has been included in IceCube’s triggered search for point-like and extended sources. A search for steady sources using 8 years of IceCube track data found this source as the second warmest source in the catalog, with a pre-trial p-value of 0.0088 [140]. It remains the most significant Galactic source in more recent searches [138].

Figure 10: Spectral energy distribution of the gamma-ray emission at the Cocoon region. The data points are from measurements by Fermi-LAT [111, 141], ARGO [142], HAWC [115], and LHAASO [118]. See Section 5 for more discussion regarding the gamma-ray observations of the Cygnus Cocoon. The solid and dashed curves denote two theory models from [115], which invoke continuous proton injection and a recent burst by a PeVatron, respectively. The magenta dotted curve indicates an upper limit converted from IceCube’s neutrino observation of the Cygnus region [143].

The IceCube and HAWC collaborations performed a joint analysis [143] to examine the correlation of neutrinos and 2HWC sources [144]. Upper limits on various Galactic sources were derived. A comparison of neutrino and gamma-ray observations of the Cygnus Cocoon is shown in Fig. 10. The joint analysis has been updated using more recent IceCube data and expanded HAWC catalogues.

A variety of stacking analysis were performed searching for the neutrino emission from classes of strong gamma-ray emitters. Using 9.5 years of all-sky data, IceCube performed a stacking search for the neutrino emission from 35 PWNe with gamma-ray emission above 1 TeV assuming a variety of weighting schemes for the sources. The analysis constrained the hadronic component of the gamma-ray emission above $\sim 10$ TeV assuming that the neutrino flux of PWNe scales as the gamma-ray flux or as the inverse of the pulsar age [145].

IceCube searched for neutrino emission from Galactic X-ray binaries in 7.5 years of data [146]. In particular, they searched for a periodic neutrino emission from 55 binaries in the northern sky as well as for a time-dependent emission using the X-ray light curves of 102 binaries in the entire sky. Finally, a search was performed for the time-integrated emission of neutrinos from four microquasars: Cyg X-3, LS 5039, LS I 61$+$303, and SS 433. No significant signal was found in the periodic and flare analysis, while in the latter, the most significant excess was found for Cyg X-3 with a post-trial p-value of 0.036, corresponding to only $1.8\,\sigma$.

With the advent of LHAASO data, IceCube performed a search for neutrinos from the 12 UHE gamma-ray sources observed by LHAASO based on 11 years of muon track events [147]. The analysis assumes an extension of the sources based on the gamma-ray observations. The neutrino spectrum is assumed to follow a power law with index between 1 and 4. Only upper limits on the hadronic component of the sources emerged. IceCube places constraints on the fraction of the gamma-ray flux originating from hadronic processes only for the Crab Nebula and LHAASO J2226+6057 at 59% and 85% of the total, when adopting a power-law spectrum as reported by LHAASO and a log-parabola spectrum by HAWC, respectively. In Ref. [147], IceCube also performed stacking searches for the LHAASO UHE sources as well as for sources lists associated with PWNe and SNRs.

Finally, using ten years of track-like events, IceCube searched for extended neutrino sources in the Galaxy [148]. They searched the neutrino sky for sources with possible extensions of 0.5, 1.0, 1.5, or 2.0 degrees. No evidence for time-integrated neutrino emission from potentially extended sources was found. In a search for the neutrino emission from 20 regions of interest coincident with TeV gamma-ray sources, the most significant region was found to be centered at 3HWC J1951$+$266, with a global significance of $2.6\,\sigma$ for an extension of $1.7^{\circ}$.

7 Extragalactic Sources

Above 10 TeV, most sources are local as a result of the absorption of gamma rays by the extragalactic background light, though a handful of extragalactic sources have been observed [149, 135]. Among them the gamma-ray burst GRB  221009A, the brightest GRB observed by Fermi-GBM to date, stands out as the first GRB observed above 10 TeV [150]. LHAASO’s WCDA observed the first main pulse of GRB 221009A at 225 s to 228 s after the GBM trigger [151]. The peak times for different energy bands agree with each other, suggesting that the peak corresponds to the onset of the afterglow. The TeV emission of GRB 221009A can be understood as synchrotron self-Compton (SSC) emission of relativistic electrons in the external shock, which is consistent with the interpretation of previous TeV afterglows [152, 153]. No gamma ray from the prompt phase was observed, likely due to the high optical depth of pair production at the early times.

Neutrino telescopes observed GRB 221009A with both real-time and follow-up searches over a wide range of energies [154, 155]. No significant deviation from background was found from MeV to PeV energies. The nondetection places tight constraints on models for neutrino emission during the prompt phase of GRBs.

In the prompt phase, neutrinos can be produced by both nonthermal and thermal protons [156, 157]. Baryons entrained in the outflow can be accelerated in internal shocks, or by magnetic reconnection in Poynting-flux-dominated outflows. The accelerated protons may interact with background photons and produce TeV–PeV neutrinos in the internal shocks [158, 159, 160, 161]. The nondetection of neutrinos from the single GRB 221009A constrains the nonthermal baryon-loading factor in the internal shock model, defined as the ratio of the fireball energy carried by protons and photons [162] $\xi_{p}=E_{p,\rm iso}/E_{\gamma,\rm iso}$, with a sensitivity that matches the constraints obtained by stacking 1172 GRBs over the history of IceCube observations [163, 164].

Neutrinos may also be produced without cosmic-ray acceleration. “Quasithermal” neutrinos may be produced by internal collisions between differentially streaming protons and neutrons as a result of the decoupling of the electrically neutral neutrons [165] (decoupling scenario) or of the transverse differences of the bulk Lorentz factor of the outflows [166, 167, 168, 169] (collision scenario). Assuming a Lorentz factor $\Gamma\sim 440$ [151] and an isotropic-equivalent energy release by the burst $E_{\gamma,\rm iso}\sim 3\times 10^{54}\,\rm erg$ [170], the nondetection of quasithermal neutrinos in the IceCube low-energy DeepCore sample (GRECO sample) [154] constrains the baryon-loading factor to $\xi_{p}\sim 0.5$ in the collision scenario.

PeV- to EeV-energy neutrino emission, which is predicted to result from ultrahigh-energy cosmic rays accelerated by the external shocks during the afterglow phase [171, 172, 162], is not well constrained. Better sensitivities at these energies provided by future telescopes such as IceCube-Gen2 [173] will enable either observing or constraining the afterglow neutrino emission.

In general, it is clear that a next generation of neutrino telescopes is required to move from discovery to astronomy. In particular, we anticipate the identification of the cosmic-ray sources in the Galaxy using the multimessenger approach described in this review.

8 Acknowledgements

The work of K.F and F.H is supported by the Office of the Vice Chancellor for Research at the University of Wisconsin–Madison with funding from the Wisconsin Alumni Research Foundation. K.F. acknowledges support from the National Science Foundation (PHY-2110821, PHY-2238916) and NASA (NMH211ZDA001N-Fermi). This work was supported by a grant from the Simons Foundation (00001470, KF). K.F acknowledges the support of the Sloan Research Fellowship. The research of F.H was also supported in part by the U.S. National Science Foundation under grants PHY-2209445 and OPP-2042807.