The LHAASO Collaboration

Introduction. — One of the most fundamental unresolved problems in astrophysics is the origin and propagation of cosmic rays (CRs). Besides the direct measurements of energy spectra and arrival directions of individual composition of CRs, the diffuse $\gamma$-ray emission plays a unique and complementary role in constraining the origin and propagation of CRs. The diffuse $\gamma$-ray emission is usually expected to be produced through interactions between CRs (nuclei and electrons/positrons) and the interstellar medium (ISM) or radiation field [1, 2, 3]. Modeling of the diffuse emission in the high-energy $\gamma$-ray bands (HE, $\lesssim$0.1 TeV) is consistent with the measurements at high and intermediate latitudes by Fermi Large Area Telescope [4], supporting the basic framework of CR propagation and interaction in the Milky Way.

Measurements of the wide-band diffuse $\gamma$-ray emission is particularly important in constraining the origin and testing the propagation models of CRs. Due either to the small field-of-view of imaging atmospheric Cherenkov telescope arrays or the limited sensitivity of extensive air shower detector arrays, the observations of diffuse $\gamma$-ray emission have been challenging for groundbased experiments for a long time. A few measurements have been reported by MILAGRO [5, 6], H.E.S.S. [7], ARGO-YBJ [8], AS$\gamma$ [9], and HAWC [10]. Particularly, the recent precise measurement by the Square Kilometer Array (KM2A) of the Large High Altitude Air Shower Observatory (LHAASO) [11] in the energy range from 10 to 1000 TeV, together with the Fermi measurement in space below TeV energy [12], gives the diffuse emission in a very wide energy range and offers new insights in understanding the origin of the diffuse $\gamma$-rays and the propagation model for CRs (e.g., [13, 14, 15, 16]). Nevertheless, the detection of diffuse emission in energy range from 1 to 10 TeV is still lacking, which hinders a good discrimination of different models. In this work, we report the measurement of diffuse $\gamma$-ray emission from the Galactic plane in the energy range from 1 to 25 TeV using the Water Cherenkov Detector Array (WCDA) of the LHAASO experiment. The energy coverage of WCDA fills the gap between Fermi and LHAASO-KM2A. Together with Fermi and KM2A, the Galactic diffuse emission covering an energy range of more than 6 orders of magnitude (from sub-GeV to PeV) will be obtained for the first time.

The LHAASO experiment. — The LHAASO [17] is a multi-component CR and $\gamma$-ray detection facility located on Haizi Mountain about 4410 meters above sea level in Daocheng, Sichuan province, China. It is designed to detect air showers induced by CR particles with energies from $\sim$TeV to $\sim$EeV and by photons from $\sim 100$ GeV to $>$PeV. The WCDA is one of the three major components of LHAASO. It is a close-packed, surface water Cherenkov detector facility with an area of 78,000 m². The WCDA consists of three ponds, two of them have an area of $150\times 150$ m² with 900 detector units each, and the third one has an area of $300\times 110$ m² with 1320 detector units. The primary goal of WCDA is to survey the northern $\gamma$-ray sky in the very-high-energy (VHE) band [17, 18]. The large field-of-view and high sensitivity of WCDA makes it a powerful detector for the study of diffuse $\gamma$-ray emission.

Data analysis. — The data used in this work is collected by WCDA between March 5, 2021 and July 31, 2023. After quality cuts, the total live time is 763 days. The event selection criteria are the same as in Ref. [19]. The parameter PINCness, defined as ${\mathcal{P}}=\frac{1}{N}\sum_{i=1}^{N}\frac{(\zeta_{i}-\langle\zeta_{i}\rangle)^{2}}{\sigma_{\zeta_{i}}^{2}}$, where $\zeta_{i}$ is the logarithm of the $i$-th PMT charge, $\langle\zeta_{i}\rangle$ and $\sigma_{\zeta_{i}}$ are the average value and dispersion of a sample of gamma-like events, is applied to exclude hadronic showers [20]. In this work we set ${\mathcal{P}}<1.1$. The region of interesting (ROI) is same as that used in the KM2A analysis [11], i.e., the inner Galaxy plane ($15^{\circ}<l<125^{\circ}$, $|b|<5^{\circ}$) and the outer Galactic plane ($125^{\circ}<l<235^{\circ}$, $|b|<5^{\circ}$. The analysis in different longitude regions are also performed. The events are binned by number of hits in six bins, i.e., 60-100, 100-200, 200-300, 300-500, 500-800, $\geq$800, which represent increasing energies from $\sim 1$ TeV to $\sim 25$ TeV. For the outer region, the measurement is performed for the last four bins to keep a relatively high significance. We also update the KM2A results for $25<E<1000$ TeV, following the analysis of [11], but using the full array data from July 20, 2021 to July 31, 2023 and the updated source mask compared to the one used in this work.

The background, mainly from residual CR events, is estimated using the “direct integral method” [21]. This method assumes that the detection efficiency is slowly varying and the time-average within a relatively short time window can properly reflect the background content. The time step to calculate the background is chosen to be 1 hour, and the time window is chosen to be $\pm 5$ hours of each step. Some spurious large-scale structures of the background may exist, which have been corrected in the analysis (see Fig. S1 of the Supplemental Material) [11]. Quite a number of pointlike and extended sources have been detected [22, 23, 19]. To reduce the impact of $\gamma$-rays sources on the background estimate, we mask the Galactic plane with latitudes $|b|\leq 10^{\circ}$ for declinations $\delta\leq 50^{\circ}$ and $|b|\leq 5^{\circ}$ for $\delta>50^{\circ}$. The sky regions around sources out of the Galactic plane, based on the TeVCat as of December 2022 [24] and WCDA/KM2A source list [19], are also masked. The mask radius is chosen to be 5 times of $\sigma\equiv\sqrt{\sigma_{\rm ext}^{2}+\sigma_{\rm psf}^{2}}$, where $\sigma_{\rm ext}$ is the fitted Gaussian extension of the source and $\sigma_{\rm psf}$ is the Gaussian width of the point spread function of WCDA. In this work we adopt $\sigma_{\rm psf}=0.5^{\circ}$ for all energy bins to keep the same sky regions.

To measure the diffuse emission, we also mask known sources detected by LHAASO [19] and those listed in TeVCat [24], with different mask radii, to balance the remaining sky region and the impacts from sources. We adopt a mask radius of 2.5 times of $\sigma$ as defined above for each source. For pulsar halos Geminga and Monogem, whose morphologies are not well described by simple Gaussian functions, the mask radius is set to be $8^{\circ}$ which can largely remove the emission from these two sources. The regions of interest (ROI) can be seen in Fig. 1, where blank regions are masked out. Since the source list and source extensions of WCDA are different from those of KM2A, the ROIs also differ slightly from the previous KM2A analysis [11]. Specifically, the solid angle of the ROI is 0.172 (0.255) for the inner (outer) region in this work, compared with 0.206 (0.268) of the previous KM2A ROI [11]. We also estimate the residual contamination from resolved sources, through a fitting of the morphological distribution of the data for different rings surrounding resolved sources (represented by the multiplicative factor of source extension $\sigma$), using simulated templates of resolved sources and the diffuse emission. See Fig. S2 of the Supplemental Material for illustration. The contamination fraction, given in Table S1 of the Supplemental Material, will be subtracted when evaluating the diffuse fluxes.

A likelihood fitting method is employed to derive the flux of the diffuse emission, in the whole energy band or individual energy bin. To properly account for the spectral shape of the emission in a wide energy band, from TeV to PeV, a smoothly broken power-law (SBPL) function, $\phi(E)=\phi_{0}\left(E/10\,{\rm TeV}\right)^{-\alpha}\left[1+(E/E_{\rm br})^{s}\right]^{(\alpha-\beta)/s}$, is assumed, where $E_{\rm br}$ is the break energy, $\alpha$ and $\beta$ are the photon indices before and after the break, and $s$ characterizes the smoothness of the break which is fixed to be $5$ in this work. For the fitting of the WCDA data alone in a limited energy range, we assume a power-law spectrum $\phi(E)=\phi_{0}\left(E/10\,{\rm TeV}\right)^{-\alpha}$ instead. The test statistic (TS) is defined as two times of the logarithmic likelihood ratio of the signal plus background hypothesis ($H_{1}$) and the background only hypothesis ($H_{0}$), i.e., ${\rm TS}=2\ln({\mathcal{L}}_{s+b}/{\mathcal{L}}_{b})$. The TS value follows approximately a $\chi^{2}$ distribution with two (four) degrees of freedom for PL (SBPL) model. The forward-folding method is used to determine the model parameters, i.e., for a given set of model parameters, we convolve the instrumental response functions to get the expected number of events, and then construct the Poisson likelihood according to the observations.

Results on diffuse emission. — Fig. 1 shows the significance maps for $N_{\rm hit}\geq 200$ in Galactic coordinate for the inner (top panel) and outer (bottom panel) regions, respectively. The one-dimensional significance distributions of these two regions and the control regions shifted by $20^{\circ}$ in Galactic latitudes are given in Fig. S3 of the Supplemental Material. The detection significance of diffuse emission is about $24.6\sigma$ ($9.1\sigma$) in the inner (outer) region. For the first time, the multi-TeV diffuse emission in the outer Galactic plane is detected by LHAASO-WCDA.

The spectral energy distributions (SED) of the diffuse emission from $\sim$TeV to PeV in the two regions are given in Fig. 2. The shorter error bars represent the statistical errors, and the longer ones are the quadratic sum of the statistical errors and the systematical uncertainties (see below for details). The data of the fluxes together with the statistical and systematic uncertainties can be found in Tables S2 and S3 of the Supplemental Material. The diffuse $\gamma$-ray fluxes are lower by a factor of several than those of CR electrons and positrons [25, 26, 27, 28], as shown in Fig. S4 of the Supplemental Material. However, in our data-driven background estimate method, the electron and positron background is largely suppressed given their nearly isotropic spatial distribution. For a better comparison with the neutrino emission [29], we also present the total fluxes without source masks in Fig. S5 of the Supplemental Material.

The wide-band spectra are fitted with the SBPL function, with parameters being presented in Table 1. The fitting results indicate that breaks of the spectra, around $30$ TeV, exist. Quantitatively, a single power-law fitting gives a TS value of 1099.10 (272.29) for the inner (outer) region, while the SBPL fitting gives 1131.25 (273.40). Therefore, the significance of the spectral break in the inner region is about $5.3\sigma$, given two more free parameters. For the outer region the break is not statistically significant. Within the statistical uncertainties, the spectral parameters in the inner and outer Galactic plane are consistent with each other. We further divide the inner region into three sub-regions, with Galactic longitudes of $l\in[15^{\circ},50^{\circ}]$, $l\in[50^{\circ},90^{\circ}]$ and $l\in[90^{\circ},125^{\circ}]$, and re-fit their spectra. The results are shown in Fig. 3. The spectral index of the third sub-region is found to be slightly harder ($\sim 2\sigma$ significance) than the other two sub-regions, as can be seen in Table 2. Such differences remain for the background estimate with different mask region ($|b|<10^{\circ}$) and time window (24 hours). These results indicate that spectral variations may exist across the Galactic plane.

We also derive the spatial distributions of the WCDA diffuse emission along Galactic latitudes and longitudes, as shown in Fig. 4. The results are compared with the gas distributions obtained from the sum of atomic hydrogen from the HI4PI survey [30], molecular hydrogen from the CO survey [31] with a CO-to-H₂ conversion factor of $X_{\rm CO}=1.6\times 10^{20}$ cm² $({\rm K\,km\,s^{-1}})^{-1}$ [4], and an ionized hydrogen component from the model of [32]. An alternative way is to use the dust opacity to trace the gas column density [33]. Differences from these two methods vary for different regions of the Galaxy. On average, we find that the inferred gas column densities in the ROIs defined in this work are smaller than 10%. It is shown that the latitude distributions are roughly consistent with the gas distributions in both the inner and outer regions. The fitting of the gas template to the data gives $p$-values of $0.003-0.013$ of the two regions. However, the longitude distribution deviates clearly from the gas distribution. The gas template fitting shows a $p$-value of $6.0\times 10^{-21}$. Similar results have also been obtained for the KM2A measurements at higher energies [11] (see Fig. S6 of the Supplemental Material for the comparison between WCDA and KM2A).

Systematic uncertainties. — One of the major sources of systematic uncertainties is the background estimate. We vary the time window from $\pm 5$ hours to $\pm 12$ hours, and also change the mask regions for background estimate to $|b|\leq 10^{\circ}$ for all declinations to obtain the variations of the results. The other effect related with the background is the spurious large-scale efficiency correction (see Sec. A of the Supplemental Material), for which we adopt different smoothing angles of $15^{\circ}$, $20^{\circ}$, $25^{\circ}$, and $30^{\circ}$. We find that the impacts of those changes on parameters $\phi_{0}$ are about $6\%~{}(16\%)$ and on $\alpha$ are about 0.05 (0.13), for the inner (outer) region. The variation of the atmosphere condition during the operation affects the detection efficiency, which is not properly modelled in the simulation. For point-like sources this effect is estimated to be about 8% for the flux ($\phi_{0}$) and 0.02 for the spectral index ($\alpha$) [19]. Similar impacts on the diffuse emission is expected. Combining these results, we obtain the total systematic uncertainties as 10% (18%) for $\phi_{0}$ and 0.05 (0.13) for $\alpha$, for the inner (outer) region, respectively. The discussion of systematic uncertainties of the KM2A results can be referred to Ref. [11]. The systematic uncertainty of the flux in each energy bin has been shown in Fig. 2 and also in Tables S2 and S3 of the Supplemental Material.

Discussion. — We compare the data with the predictions of the hadronic interactions between CRs and the ISM. We adopt the parameterizations of the local proton and helium spectra as given in Ref. [11], which include the uncertainties of the current measurements, especially for those measured by indirect detection experiments above $\sim 100$ TeV. We assume that the CR intensity is uniform in the Galaxy¹¹1 Note that this assumption is over-simplified since variations of the CR intensities across the Galaxy were found to be present [34, 35]. Assuming a GALPROP model, which approximates the propagation halo as a cylinder, employs a numerical method to solve the transport equations, and includes experimental/observational inputs about the nuclear cross sections and gas distributions [2], we find that the expected diffuse fluxes are higher (lower) by 10% (40%) than the results obtained here for the inner (outer) region after the same source mask., and calculate the $\gamma$-ray emission using the local CR spectra and the line-of-sight integrated gas content. The $\gamma$-ray production cross section used is from the AAfrag package [36]. The expected diffuse $\gamma$-ray fluxes are shown by grey shaded bands in Fig. 2. We can see that in the inner region the measurements are higher by a factor of $\sim 2.7$ than the prediction, and in the outer region the measurements are higher by $\sim 1.5$ times. These results are consistent with those found in the LHAASO-KM2A analysis [11]. The HAWC measurements in a different sky region of the inner Galaxy also showed a factor of 2 higher fluxes compared with a DRAGON model prediction [10]. In the outer region, the predicted spectrum is consistent with the data within uncertainties. We test the consistency of the spectral shape between the prediction and the measurement in the inner region through a fit with free normalization of the model curve, and get $\chi^{2}=86.14$ and 19.16 for 13 degrees of freedom, for the “high” and “low” predictions (corresponding to the upper and lower edges of the band shown in Fig. 2). The $p$-values of the fittings are $7.61\times 10^{-13}$ and 0.12, respectively. Given the current large uncertainties of the spectral measurements of individual CR species, we cannot draw a firm conclusion on whether the diffuse $\gamma$-ray data is consistent with the knee of the light CR components or not based on the spectral shape only.

If we compare the data (including those measured by Fermi-LAT between 1 and 500 GeV [12], updated with the new mask of this work) with the GALPROP model prediction²²2 The propagation and source injection parameters were adjusted [37] to fit the up-to-date measurements of CR primary and secondary nuclei [38, 39, 40]. One should be aware that, a homogeneous, isotropic diffusion is assumed when calculating the CR propagation, which may be too simple [41]. in a wider energy range, from GeV to PeV, we find that the excess mainly exists from several GeV to tens of TeV (see Fig. 5). Such an excess may come from a population of unresolved, low surface brightness sources such as pulsar wind nebulae or pulsar halos [1, 42, 43, 13, 16]. The spatial variations of the spectra of the diffuse emission may further supports the hypothesis that unresolved sources may exist, although other possibilities such as spatially-dependent spectral slopes of CRs can not be excluded. Following Ref. [12], we add an empirical component with spectrum $\propto E^{-2.45}\exp\left(-E/20\,{\rm TeV}\right)$ to the CR propagation model prediction, and find a good agreement with the wide-band measurements, as shown in Fig. 5. Other sources like young massive star clusters [44], the confinement and interaction of CRs around acceleration sources [45, 46, 47, 48], as well as the change of the canonical propagation scenario [49, 50, 51, 52] were also discussed in literature. The multi-messenger analysis of the diffuse $\gamma$ rays and neutrinos [29] would be helpful in understanding the nature of the diffuse emission.

Summary. — In this work we measure the diffuse $\gamma$-ray emission in an energy range of about $1-25$ TeV from the Galactic plane using the LHAASO-WCDA data. The detection significance of the diffuse emission, after masking sky regions around resolved sources, is $24.6\sigma$ for the inner Galactic plane ($15^{\circ}<l<125^{\circ}$, $|b|<5^{\circ}$) and $9.1\sigma$ for the outer Galactic plane ($125^{\circ}<l<235^{\circ}$, $|b|<5^{\circ}$) regions, respectively. The diffuse emission from the outer Galactic plane is for the first time detected in the multi-TeV energy range. The WCDA measurements fill the gap between direct detection by space telescopes and ultra-high-energy measurements by groundbased experiments, and offer a full energy coverage from GeV to PeV for the Galactic diffuse $\gamma$-ray emission. The SEDs of the diffuse emission in the two regions have been measured with high precision. The WCDA spectra are consistent with power-law models, with indices of $-2.67\pm 0.05_{\rm stat}$ and $-2.83\pm 0.19_{\rm stat}$. The joint WCDA-KM2A spectrum shows a break around 30 TeV in the inner region, with indices changing by about $0.5$. The break in the outer region is not significant. We also show that there are slight spectral variations across the Galactic plane in the inner region. The measured fluxes are higher by a factor of $1.5-2.7$ than the simple prediction of hadronic interactions between CRs (with locally measured fluxes) and the ISM, indicating that there are unmodelled source components or spatial variations of CR distribution.

Acknowledgements. — We would like to thank all staff members who work at the LHAASO site above 4400 meter above the sea level year-round to maintain the detector and keep the water recycling system, electricity power supply and other components of the experiment operating smoothly. We are grateful to Chengdu Management Committee of Tianfu New Area for the constant financial support for research with LHAASO data. We appreciate the computing and data service support provided by the National High Energy Physics Data Center for the data analysis in this paper. This work is supported by the following grants: the Department of Science and Technology of Sichuan Province, China (No. 2024ZYD0111, 24NSFSC2319), the National Natural Science Foundation of China (Nos. 12220101003, 12173039, 12393851, 12393852, 12393853, 12393854, 12205314, 12105301, 12305120, 12261160362, 12105294, U1931201, 12375107, 12342502), the Project for Young Scientists in Basic Research of Chinese Academy of Sciences (No. YSBR-061), the National Science and Technology Development Agency (NSTDA) of Thailand, and the National Research Council of Thailand (NRCT) under the High-Potential Research Team Grant Program (N42A650868).