Search for Neutrino Emission from the Cygnus Bubble Based on LHAASO γ-Ray Observations

The Cygnus region, which contains massive molecular and atomic clouds and young stars, is a promising Galactic neutrino source candidate. Cosmic-ray transport in the region can produce neutrinos and γ-rays. Recently, the Large High Altitude Air Shower Observatory (LHAASO) detected an ultrahigh-energy γ-ray bubble (Cygnus Bubble) in this region. Using publicly available track events detected by the IceCube Neutrino Observatory in 7 yr of full detector operation, we conduct searches for correlated neutrino signals from the Cygnus Bubble with neutrino emission templates based on LHAASO γ-ray observations. No significant signals were found for any employed templates. With the 7 TeV γ-ray flux template, we set a flux upper limit of 90% confidence level for the neutrino emission from the Cygnus Bubble to be 5.7 × 10−13 TeV−1 cm−2 s−1 at 5 TeV.


INTRODUCTION
Cosmic rays are high-energy astrophysical particles, primarily protons and atomic nuclei, while their origins have been a mystery for a century.Under the confinement of the Galactic magnetic field, the observed cosmic rays with energies up to several PeV are believed to originate from Galactic sources, called PeVatrons.Cosmic rays interact with the interstellar medium or the radiation field, generating both neutrinos (e.g., π + → µ + + ν µ ) and γ-rays (e.g., π 0 → 2γ).High energy electrons can also produce γ-rays through inverse Compton scattering.However, the cross section suffers more stringent Klein-Nishina suppression for γ-rays with energies above 100 TeV.Therefore, the coincidence between neutrinos or γ-rays (> 100 TeV) and gas clumps will provide critical evidence for the identification of hadronic PeVatrons.
The Cygnus region is an active star forming area in our Galaxy and hosts various astrophysical sources, including massive young star clusters (YMCs, e.g., Cygnus OB2), pulsar wind nebulae (PWNe, e.g., TeV J2032+4130), and supernova remnants (SNRs, e.g., γ-Cygni).Fermi-LAT detected an excess of γ-ray emission (1-100 GeV) from the direction of the Cygnus region after subtracting the interstellar background and all known sources (Ackermann et al. 2011).The hard γ-ray spectrum points to freshly accelerated cosmic rays, whether they are cosmic ray electrons or nuclei.This ∼ 2 • extended γ-ray source, known as Cygnus Cocoon, has been further observed at TeV energies by ARGO-YBJ (Bartoli et al. 2014) and HAWC (Abeysekara et al. 2021).In the latest observation of the Cygnus region, LHAASO reported the Cygnus Bubble at ultra-high energy (LHAASO Collaboration 2024), extending to more than 6 • from the core, which is much larger than the Cygnus Cocoon.The γ-ray brightness follows the distribution of the molecular gas, especially for γ-rays above 100 TeV, suggesting that these γ-rays are produced by the collision between the gas and the cosmic rays.We expect to observe high-energy neutrinos from the Cygnus Bubble.
The IceCube Neutrino Observatory has previously searched for neutrino emission from the PeVatron candidates observed by LHAASO.The hadronic components of the Crab Nebula and LHAASO J1849-0003 are constrained to be no more than ∼ 80% and ∼ 90% of the total γ-rays observed (Huang & Li 2022;Abbasi et al. 2023).As for the Cygnus region, the hadronic contribution is contained to be less than 60% (Kheirandish & Wood 2019), while the resolved sources (e.g., TeV J2032+4130) are not removed from this region.Recently, Neronov et al. (2023) claimed a 3σ excess of neutrino signals from the central region (∼ 1 • ) of the Cygnus region.However, the neutrino emission from the entire Cygnus Bubble remains unclear.
In this study, we conduct two analyses on the Cygnus Bubble using a neutrino data sample of IceCube track events from 2011 to 2018.Firstly, we search for neutrino emission from the Cygnus Bubble with a template likelihood method and set 90% C.L. upper limits on the muon neutrino flux.In addition to using the γ-ray flux maps as the neutrino emission template, we employ six other templates for comparison, testing different template radii.Secondly, we scan the region of the Cygnus Bubble and obtain neutrino hotspots, which are compared, respectively, with γ-ray hotspots, gas distribution, and sources from the TeVCat1 (Wakely & Horan 2008) and the first LHAASO catalog (Cao et al. 2023).
This paper is organized as the following structure.The LHAASO observations and the IceCube muon-track data are introduced in the section 2. The analysis methods for the template search and Cygnus Bubble scan are introduced in section 3. Results and further discussions are shown in section 4. Finally, section 5 summarizes the conclusions.

IceCube Neutrino Sample
The IceCube Neutrino Observatory is a cubic kilometer detector located at the South Pole (Abbasi et al. 2009).Installed between 1.45 and 2.45 km below the surface of ice, IceCube consists of 86 strings equipped with digital optical modules (DOMs), which can detect Cherenkov light emitted by the secondary charged particles (Aartsen et al. 2017a).Muon neutrinos propagating inside the Earth can produce ultrarelativistic muons via charge-current (CC) interactions.These muons, when traversing the detector, will leave a track-like signature.With high statistics and a typical angular resolution of ≲ 1 • at ∼ TeV (Aartsen et al. 2020), track-like events are adequately used for neutrino source searches.
In the analyses, we use 7 years of all-sky muon track data collected by the completed 86-string detector, namely IC86-2011 (IC86-I) and IC86-2012-18 (IC86-II) (IceCube Collaboration 2021).The data consists of three components: (i) Experimental data events, including the reconstructed direction with the R.A. (α) and decl.(δ), angular uncertainty (σ), and reconstructed muon energy (E rec ) for each event.(ii) Instrument response functions, including the effective area A eff (E ν , δ ν ) and the smearing function M (E rec |E ν , δ ν ).The smearing function gives the fraction count of simulated signal events within the reconstructed muon energy (E ν , δ ν , E rec ) bin relative to all events in the (E ν , δ ν ) bin.With the smearing function, the signal energy probability density function (PDF) of the likelihood can be derived under a source spectrum assumption.(iii) The detector uptime, which records the periods of data taking.

LHAASO γ-Ray Data
LHAASO comprises composite detection arrays that aim to study cosmic rays and γ-rays (Cao 2010).Located at ∼ 29 • North in Sichuan Province, China, LHAASO covers a large sky region spanning from −21 • to 79 • in declination.Two arrays of LHAASO, the Kilometer Square Array (KM2A) and the Water Cherenkov Detector Array (WCDA), are used for γ-ray detection.The ∼ 1.3 km 2 KM2A is able to detect photons with energies from 10 TeV to several PeV, while the 0.078 km 2 WCDA probes lower energy photons ranging from 100 GeV to 20 TeV (He 2018).The angular resolution of KM2A is 0.4 • at 30 TeV and can reach 0.2 • at 1 PeV (Addazi et al. 2022).For the WCDA, the angular resolution is better than 0.2 • at 10 TeV.
The Cygnus Bubble, recently reported by LHAASO, was measured by both KM2A and WCDA (LHAASO Collaboration 2024).The residual structure extends to ∼ 10 • after the removal of all resolved γ-ray sources and applying a circular mask with a radius of 2.5 • around LHAASO J2018+3651.The γ-ray excess within a radius of 6 • is still clear after accounting for the diffuse γ-ray background.In the 6 • radius region, the energy spectrum of Cygnus Bubble is fitted by a log-parabola function with energies from 2 TeV to 2 PeV.The fitted photon index of Cygnus Bubble is Γ = (2.71± 0.02) + (0.11 ± 0.02) × log 10 (E/10 TeV).Eight photons with energies above 1 PeV are detected in this region, indicating the existence of super PeVtron(s).The significance maps of different energy bands show a brightening in the center associated with massive molecular clouds.The γ-rays from the Cygnus Bubble are characterized by four components: two diffuse components with the γ-ray emission proportional to the column density of atomic (H I, HI4PI Collaboration et al. 2016) and molecular clouds (MCs, Dame et al. 2001), and two extended sources, LHAASO J2031+4057 and LHAASO J2027+4119.LHAASO J2031+4057 is only observed by WCDA at energy ranges below 20 TeV, while the other three components are observed both by WCDA and KM2A.The γ-ray flux map measured by LHAASO can be obtained with the spatial and the spectral information of these components.

Template Search
An unbinned maximum likelihood method is widely used in neutrino point source searches (Braun et al. 2008(Braun et al. , 2010)).A detailed description of the point source likelihood can be found in Appendix A. Considering the large extension (∼ 6 • ) of the Cygnus Bubble, point source likelihood is not suitable for this analysis.Here, we search for neutrino emission associated with the Cygnus Bubble following the ps-template method (Aartsen et al. 2017b).There are two modifications compared to the point source likelihood.Firstly, a spatial template, rather than a two-dimensional (2D) Gaussian function, is used to describe the spatial distribution of signal events.The ps-template method accounts for the extension of the source by mapping the changing detector acceptance and convolving the template with the angular uncertainty of the events.Secondly, unlike in the point source likelihood where the background is estimated using scrambled data with negligible point source signal contribution, for a large extended source, the signal events in the data should be subtracted.Therefore, a signal-subtracted likelihood is constructed, and the background is estimated using scrambled data with the signal contamination subtracted (Pinat & Sánchez 2018).
The event-wise template likelihood (Aartsen et al. 2017b) is defined as where n s is the number of signal events under the source spectra assumption with a spectral index of γ, and N is the total number of events.S i is the signal PDF, which is related to the location x i , angular uncertainty σ i , and muon energy proxy E i of the i-th event.D i is the scrambled background PDF, which is obtained from the data and therefore contains the scrambled signal component, while S i is the scrambled signal PDF.Each PDF consists of a spatial term and an energy term.In the template likelihood, the neutrino spectrum is fixed and we only fit the number of signal events ns by maximizing the likelihood.The neutrino spectrum is calculated with the parameterized energy distribution of secondary particles produced in p-p interactions, assuming that all the γ-rays are originate from hadronic processes (details can be found in Appendix B).
The construction of the signal spatial PDF in the template likelihood can be described as We start with the Cygnus Bubble template T spat , which is treated as the neutrino spatial template.By convolving it with IceCube acceptance M acc , we can obtain the true neutrino direction after accounting for the detector efficiency.Then this map is smoothed with a 2D Gaussian of width σ i to account for the angular uncertainty of the events.Finally, this map is normalized to unity.The background PDF D i is constructed as the same method in point source likelihood, while the scrambled signal PDF S i is constructed following (Pinat 2017).
We use the γ-ray emission and the gas column density to weigh the neutrino emission within a 6 • radius from the bubble center.In the γ-ray flux template, neutrino emission is assumed to follow the γ-ray flux map at 7 TeV (WCDA) and 50 TeV (KM2A).In the H I and MC templates, neutrino emission follows the column densities of H I and MC, respectively, derived from the H I and CO emission.In the hydrogen (MC+H I) template, neutrino emission follows the total hydrogen column density (2N H2 + N HI ).In the Gaussian templates, neutrino emission follows the 2D Gaussian distribution with σ = 0.33 • for LHAASO J2031+4057 and σ = 2.28 • for LHAASO J2027+4119.Finally, we employ a uniform template for comparison with the other templates, assuming a uniform spatial distribution of the neutrino flux.The above templates cover the region of LHAASO 2018+3651.In addition to the 6 • radius, we explore templates with radii of 0.7 • and 1.2 • based on recent findings (Neronov et al. 2023), as well as a 10 • radius according to LHAASO's measurements.The center of each template is LHAASO J2032+4102 (R.A. = 308.05• , decl.= 41.05 • ).

Cygnus Bubble Scan
We scan a 28 • × 22 • region, extending to ∼ 10 • to include the entire bubble, to investigate the relation between the neutrino and γ-ray hotspots and sources, as well as the distribution of MC and H I gas.The region is divided into a grid of points, each occupying an area of 0.1 • × 0.1 • .Because the γ-ray significance maps are smoothed with a Gaussian kernel of σ = 0.3 • , we scan this region with sources having a matching σ s = 0.3 • extension for consistency.In the likelihood, the signal spatial PDF for extended source is modified from a 2D Gaussian used for the point source likelihood, as shown in Equation (A6).The modified signal spatial PDF is defined as where σ s = 0.3 • .Other PDFs remain the same as in the point source likelihood.For each grid we maximize the likelihood (see Appendix A) by fitting two parameters: the number of signal events ns and the spectral index γ assuming a power-law energy spectrum.

Template Search Results
The results of template searches using the γ-ray flux maps of the Cygnus Bubble are summarized in Table 1.Although some excess from the Cygnus Bubble are observed, the results are not statistically significant.The γ-ray flux template at 7 TeV yields the lower pretrial p-value of 0.176 (0.9σ).We set upper limits on the muon neutrino flux, which are found to be ∼ 3 times higher than the theoretically expected neutrino flux and shown in Figure 1.In the previous search in Cygnus region, no neutrino excess (n s = 0) was found with a pretrial p-value of 0.80 (Kheirandish & Wood 2019), which was probably due to the use of different data samples and templates.The results for the other six templates are summarized in Table 2.Among them, the Gaussian template for LHAASO J2031+4057 (σ = 0.33 • ) yields the lowest pretrial p-value of 0.007 (2.4σ), while the H I template yields the largest p-value of 0.291 (0.6σ).The lower significant result for the γ-ray flux map at 50 TeV is probably due to the lack of LHAASO J2031+4057 component.
The template search results with different template radii of 0.7 • , 1.2 • , 6.0 • , and 10.0 • are shown in Figure 2. The MC template with a radius of 1.2 • gives the most significant result, with a pretrial p-value of 2.6 × 10 −3 (2.8σ).At larger radii of 6 • and 10 • , the neutrino excess of γ-ray flux template at 7 TeV is more significant.At the radius of 0.7 • , the significance of neutrino excess is not sensitive to templates.The best fit number of signal event (n s = 46.6)for the MC template (1.2 • ) is much higher than the expected signal event number (n exp = 4.2) within the central 1.2 • region.It seems challenging to attribute all the observed neutrino excess solely to the neutrinos accompanying the γ-rays from the Cygnus Bubble observed.
While the excess of the neutrino signal is not substantial, our results are still consistent with the hadronic origin of the γ-ray emission from the Cygnus Bubble, as the upper limits of the flux exceed the theoretically predicted neutrino flux.To obtain more significant results, additional through-going track events are required.Furthermore, cascade events might be more suitable for measuring neutrinos originating from the extensive region of the Cygnus Bubble.This is because the high angular resolution of track events does not provide a significant advantage in reducing background, and cascade events have lower atmospheric background compared to track events.] flux derived from LHAASO (6°, mask J2018) flux derived from LHAASO (6°) 90% C.L. upper limit ( -ray flux map at 7 TeV) Figure 1.Upper limits (90% C.L.) on the muon neutrino flux (red) resulting from the template searches using the γ-ray flux map at 7 TeV for the Cygnus Bubble.The bold red solid line shows the center energy range contributing to 90% significance.The expected muon neutrino flux (black and gray), derived from LHAASO γ-ray observations assuming hadronuclear interactions, is also shown.The 2.5 • region centered on LHAASO J2018+3651 is masked for the gray line.

Cygnus Bubble Scan Results
The results of the Cygnus Bubble scan are shown in Figure 3, with the upper panels (A-C) illustrating the neutrino significance map and the lower panels (D-F) illustrating the neutrino excess map.In the entire scan region, the most significant point, indicated by the white cross, is found at R.A. = 303.35• and decl.= 43.75• with a pretrial p-value of 2.2 × 10 −3 (2.9σ).This point is located 4.4 • away from the template center, with the best fit parameters being ns = 22.2 and γ = 2.3.In the central 2 • region, the most significant point is found at R.A. = 308.25 • and decl.= 40.45• with a pretrial p-value of 6.3 × 10 −3 (2.5σ).This point is located ∼ 0.6 • away from the template center, with the best fit parameters being ns = 31.7 and γ = 4.0.We further conduct Monte Carlo simulations by scrambling the R.A.s of IceCube events to compute the post-trial probability of the neutrino hotspot being more significant than the observed one.The post-trial p-value is 0.96 in the 28 • × 22 • region and 0.84 (0.18) within the 10 • (2 • ) region.
The γ-ray significance map observed by LHAASO partly correlates with the H I distribution and clearly associates with the dense MC clumps.If these γ-rays are produced by cosmic rays interacting with the surrounding gas, the accompanying neutrinos are expected to follow the γ-ray and gas distributions.The neutrino hotspot in the bubble center is spatially associated with the γ-ray hotspot below 20 TeV.However, the neutrino significance map in the larger region (see the upper panels A-C in Figure 3) doesn't exhibit an obvious association with the γ-ray significance map.Similarly, the neutrino excess map (see the lower panels D-F in Figure 3) also lacks a clear correlation with the gas distribution or the γ-ray distribution.LHAASO has observed eight PeV photons from the Cygnus Bubble.Two of these PeV photons are associated with the neutrino hotspot in the bubble center.One PeV photon is located close the neutrino excess around (R.A. = 305.05• , decl.= 45.95 • ), and another PeV photon is located close to the neutrino excess around (R.A. = 300.55• , decl.= 43.55• ), which is close to the blazar MAGIC J2001+435 (R.A. = 300.32• , decl.= 43.88• ).IceCube had previously reported the upper limit on neutrino flux from this source, MG4 J200112+4352, in source-list searches using ten years of track data, yielding a pretrial p-value of 0.21 (0.8σ) (Aartsen et al. 2020).

CONCLUSION
In this study, we conducted template searches and a scan of the Cygnus Bubble observed by LHAASO to investigate a potential correlation between neutrinos and γ-rays using 7 years (2011-2018) of IceCube muon neutrino data observed by the full detector.Using various spatial templates of neutrino emission based on different assumptions, we found no significant neutrino signals in the Cygnus Bubble.The most significant result of the template searches is obtained with the MC template in a 1.2 • radius, centered at LHAASO J2032+4102 (R.A. = 308.05• , decl.= 41.05 • ), yielding a pretrial p-value of 2.6 × 10 −3 (2.8σ).As for the signals from larger regions with radii of 6 • and 10 • , the γ-ray flux template at 7 TeV yields more significant results than other templates.We obtained 90% C.L. upper limits on the neutrino flux for each template.By comparing the resulting upper limits with the theoretically predicted neutrino flux based on γ-ray observations assuming hadronuclear interactions, we conclude that the neutrino result is consistent with the pure hadronic origin of the γ-ray emission from the Cygnus Bubble.The neutrino signals exhibit a stronger tendency to follow the MC distribution compared to the γ-ray flux distribution in the central region (∼ 1 • ) of the Cygnus Bubble.The p-value is the probability of background TS being greater than the observed one.The distribution of background TS can be derived from the simulations by scrambling the IceCube events in R.A.s.The distribution consists of two parts: the over-fluctuating part with T S > 0 and the under-fluctuating part with T S = 0.If the shape of neutrino spectrum is fixed, the over-fluctuating part is expected to follow a χ 2 1 distribution, and the fraction of the under-fluctuating part is 0.5.Thus, the pretrial p-value can be expressed as where σ pp is the cross section of p-p interactions, F ν (F γ ) is the energy distribution probability of neutrinos (γrays) generated from a cosmic-ray proton with the energy of E p (Kelner et al. 2006), V is the volume of emission region, n H is the number density of hydrogen, d is the distance to observer, and J p is the energy spectrum of cosmicray density.The neutrino spectrum can be obtained from the best fit γ-ray spectrum ϕ γ with a photon index of Γ = 2.71 + 0.11 × log 10 (E γ /10 TeV), as the integral over the emission region is the same for both neutrino and γ-ray emissions.The observed neutrinos are assumed to be equal flavor ratio due to neutrino oscillation.

Figure 2 .
Figure 2. The results (significance) of template searches are shown as a function of the template radius.Results from various templates with radii of 0.7 • , 1.2 • , 6.0 • , and 10.0 • are presented.The results using the LHAASO γ-ray flux template at 7 TeV are represented as red squares.

Figure 3 .
Figure 3. Neutrino significance map with the pretrial p-value −log10p (panels A-C) and neutrino excess map with the best fit number of signal events ns (panels D-F) for the Cygnus Bubble scan.The γ-ray significance map with energies above 100 TeV (panel A), 25 − 100 TeV (panel B), 2 − 20 TeV (panel C) are shown by contours starting from 3σ and increasing in steps of 3σ.The spatial distribution of the MC template (panel D), the H I template (panel E), and the γ-ray flux template at 7 TeV (panel F) are indicated by contours smoothed with a Gaussian kernel of σ = 0.5 • for comparison with the neutrino excess map.Eight photons with energies beyond 1 PeV are shown as gold triangles.Sources from TeVCat and the first LHAASO catalog located in this region are indicated by green plus signs.The most significant point in the entire scan region (R.A. = 303.35• , decl.= 43.75• ) and in the central 2 • region (R.A. = 308.25 • , decl.= 40.45• ) are indicated by white crosses.

Table 1 .
Results of template searches using the γ-ray flux templates of the Cygnus Bubble for a 6 • region.The spatial template, the best fit number of signal events ns, the 90% C.L. muon neutrino upper limit flux at 5 TeV with units of TeV −1 cm −2 s −1 , and the pretrial p-value of each search are listed.

Table 2 .
Results of template searches.Same as Table1but for the MC, H I, hydrogen, uniform, and two 2D Gaussian templates within a 6 • radius region.
THE NEUTRINO AND GAMMA-RAY CONNECTIONIn p-p interactions, the energy spectra of neutrinos and γ-rays can be expressed as ϕ ν,γ (E) = B.