© ESO 2022

1 Introduction

Cosmic rays (CRs) are believed to be well mixed in the magnetic field in our Galaxy and interact with the cold gas in the interstellar medium (ISM) to produce the diffuse γ-ray emissions. Thanks to the progress in the γ-ray astronomy, such diffuse γ-ray emissions have been well studied in recent years (Acero et al. 2016; Yang et al. 2016). Taking into account the gas content derived from molecular and atomic line spectrometry, one can obtain information on the CR distribution on the Galactic scale. By assuming a cylindrical symmetric distribution, Yang et al. (2016) and Acero et al. (2016) found that CRs are strongly inhomogeneous in our Galaxy. The radial distribution of CRs peaks at the galactocentric distance of about 4 kpc, and the derived CR spectra in the inner Galaxy are significantly harder.

The spatial distribution of CRs has direct implications on the understanding of their origin and propagation. Generally speaking, the diffusion process will make the CR distribution smoother than the distribution of CR sources, but it cannot erase all hints of the source distributions (see, e.g., Strong & Moskalenko 1998). Thus, the sharp peak of the derived CR distributions requires a significant modification to the “standard” CR propagation model in our Galaxy. As an example, Recchia et al. (2016) invoked the CR-driven wind in the Galactic plane to explain the observed CR radial profile.

However, whether such inhomogeneity is a global effect or caused by individual sources is still unclear. Aharonian et al. (2020) used the γ-ray observations of individual molecular clouds to derive the CR density therein, and unveiled a homogeneous “sea” of CRs with a constant density and spectral shape, especially for the galactocentric distances exceeding 8 kpc, as well as the Sagittarius B complex, in the region of the Galactic center (GC). However, in the region with galactocentric distances between 4 and 6 kpcs, the derived CR densities reveal significant deviation from the “sea” and each other. One possible explanation is that larger amounts of particle accelerators within this region cause the enhancement of CR density in their vicinity. The distributions of potential CR accelerators, such as supernova remnants (SNRs), pulsars, and OB stars, do reveal strong galactocentric radius dependence (Green 2015; Yusifov & Küçük 2004; Bronfman et al. 2000). Particularly, the distribution of OB stars presents a similar peak at about 4 kpc (Bronfman et al. 2000), which makes these massive star clusters an attractive source population for addressing the higher CR density and harder CR spectra in this region. Indeed, young massive star clusters (YMCs) were recently identified as a new population of γ-ray sources and are believed to be an alternative source of the Galactic CRs (Aharonian et al. 2019). In several such systems, including the Cygnus cocoon (Ackermann et al. 2011; Aharonian et al. 2019), Westerlund 1 (Abramowski et al. 2012), NGC 3603 (Yang & Aharonian 2017), Westerlund 2 (Yang et al. 2018) and 30 Dor C (H.E.S.S. Collaboration 2015), extended γ-ray emissions with hard spectra have been detected from GeV to TeV. Besides this, most of the observed γ rays reveal a good spatial correlation with the H II region, which is formed due to the photoionization of these massive stars (Yang & Aharonian 2017; Aharonian et al. 2019; Sun et al. 2020a,b). Due to strong extinction in the Galactic plane, there may be numerous unknown YMCs exist in the Galaxy. Thus, if there are additional contributions from the unknown YMCs, we expect the γ-ray emission to present similar spatial distribution as the H II gas.

To test such a hypothesis, we chose two regions centered at (l, b) = (30°, 0°) (region I) and (l, b) = (330°, 0°) (region II) to investigate the diffuse γ-ray emissions aswell as the distribution of gases in different phases. We chose these regions because they are near to the tangent point of the 4-kpc ring, and the aforementioned CR excess should have a strong impact on the total γ-ray fluxes in these directions. In the following, we give a detailed analysis of the 12-year Fermi-LAT data toward these regions and discuss the gas and CR content therein. The paper is organized as follows. In Sect. 2, we investigate the gas distribution in this region. In Sect. 3, we present details of the γ-ray data analysis. In Sect. 4, we discuss the implications of our analysis results.

2 Gas tracers

The 21 cm H I line and 2.6 mm CO line are widely used tracers for atomic hydrogen and molecular hydrogen, respectively. In this work, we used the data from the Galactic CO survey of Dame et al. (2001) with the CfA 1.2m-millimeter-wave Telescope, and the HI4PI Survey on HI gas (HI4PI Collaboration 2016). For the CO data, we used the standard assumption of a linear relationship between the velocity-integrated CO intensity, W_CO, and the column density of molecular hydrogen, N(H₂). The conversion factor X_CO may be different in different parts of the Galaxy; therefore, we used the values derived from Acero et al. (2016) for different galactocentric distances, such as $X_{\text{CO}}=0.5\times {10}^{20} {\text{cm}}^{-2}{(\text{K} \text{km}{\text{s}}^{-1})}^{-1}$ below 1.5 kpc and $X_{\text{CO}}=1.5\times {10}^{20} {\text{cm}}^{-2}{(\text{K} \text{km}{\text{s}}^{-1})}^{-1}$ above.

For the HI data, we use the following equation: $$N_{\text{HI}}(v,T_{\text{s}})=-\log \left(1-\frac{T_{\text{B}}}{T_{\text{s}}-T_{\text{bg}}}\right)T_{\text{s}}{C}_i\Delta v,$$(1)

where T_bg ≈ 2.66 K is the brightness temperature of the cosmic microwave background radiation at 21 cm, and C_i = 1.83 × 10¹⁸cm². In the case when T_B > T_s − 5 K, we truncated T_B to T_s − 5 K; T_s was chosen to be 150 K. The systematic uncertainties due to the different spin temperatures are discussed in Ackermann et al. (2012) and Acero et al. (2016). The effect, however, is quite small in most regions of the sky.

In order to investigate the galactocentric radial distribution of CRs, we need to divide the gas distribution into different rings around the GC. For this purpose we used the following relation: $$V_{\text{LSR}}=R_{\odot }\left(\frac{V(R)}{R}-\frac{V_{\odot }}{R_{\odot }}\right)\sin (l)\text{cos}(b),$$(2)

where V_LSR is the radial velocity with respect to the local standard of rest, R is the galactocentric distance, V (R) is the Galactic rotational curve, and l and b are the galactic coordinates in the line of sight (LOS). We adopted the rotational curve parameterized in Clemens (1985); V_⊙ and R_⊙ are fixed to 220 km s⁻¹ and 8.5 kpc, respectively. By applying this to both CO and HI data, we can transform the velocity information into the galactocentric distance of the gas. According to Fig. 8 of Acero et al. (2016), we split the gas into to galactocentric distance bins; for example, below and above 4.5 kpc in the analysis below. It should be noted that Eq. (2) allows emission from the forbidden velocity zones in the CO and HI data. The reason could be the non-circular motion of gas. This component contains only a small fraction of total gas; therefore, for simplicity we assigned it into the local rings.

For different reasons, the neural gas cannot always be traced by CO and H I observations (Grenier et al. 2005). In such cases (e.g. in optically thick clouds), the infrared emission from cold interstellar dust provides an alternative and independent measurement of the gas column density. In the analysis below, we also used the dust column map as gas template for cross-check. According to Eq. (4) of Planck Collaboration XIX (2011), the relation between the dust opacity and the column density can be approximated as $${\tau }_{\text{M}}(\lambda )={\left(\frac{{\tau }_{\text{D}}(\lambda )}{N_{\text{H}}}\right)}^{\text{dust}}[N_{\text{HI}}+2X_{\text{CO}}{W}_{\text{CO}}],$$(3)

where τ_M is the dust opacity as a function of the wavelength λ, ${({\tau }_{\text{D}}/N_{\text{H}})}^{\text{dust}}$ is the reference dust emissivity measured in low-N_H regions, W_CO is the integrated brightness temperature of the CO emission, and $X_{\text{CO}}=N_{{\text{H}}_2}/W_{\text{CO}}$ is the H₂∕CO conversion factor. The substitution of the latter into Eq. (3) gives $$N_{\text{H}}=N_{\text{HI}}+2N_{{\text{H}}_2}={\tau }_{\text{M}}(\lambda ){\left[{\left(\frac{{\tau }_{\text{D}}(\lambda )}{N_{\text{H}}}\right)}^{\text{dust}}\right]}^{-1}.$$(4)

Here for the dust emissivity at 353 GHz, we used ${({\tau }_{\text{D}}/N_{\text{H}})}_{353 \text{GHz}}^{\text{dust}}=(1.18\pm 0.17)\times {10}^{-26}$ cm² taken from Table 3 of Planck Collaboration XIX (2011). The derived total gas column density distributions in region I and region II are shown in the top panel of Fig. 1.

We also used the dust data to derive the residual templates for the missing “dark” gas component. We fit the dust opacity maps as a linear combination of HI and CO maps, and, in this way, find the residual map to be the dark gas template. Then, we iterated this fit by includingthe dark gas template until a convergence is achieved (Yang et al. 2016). This method is similar to the derivation of the $E{(B-V)}_{\text{res}}$ templates used by the Fermi-LAT collaboration (Acero et al. 2016), where instead of the extinction maps we used the dust opacity maps. Indeed, E(B − V) has nearly perfect linear correlation with the dust opacity, especially in higher column regions (Planck Collaboration XI 2014). Thus, our method and the method used in Acero et al. (2016) should give similar results.

We used the free-free emission map obtained from the joint analysis of Planck, WMAP, and 408 MHz observations (Planck Collaboration X 2016) to derive the map of H II column density. First, we converted the emission measure into free-free intensity (I_ν) by using the conversion factor at 353-GHz in Table 1 of Finkbeiner (2003). Then, we applied Eq. (5) from Sodroski et al. (1997), $$\text{nnbegin}aligned{N}_{\text{ion}Hii}=\text{ \&}1.2\times {10}^{15}{\text{cm}}^{-2}{\left(\frac{T_{\text{e}}}{1\text{K}}\right)}^{0.35}{\left(\frac{\nu }{1\text{GHz}}\right)}^{0.1}{\left(\frac{n_{\text{e}}}{1{\text{cm}}^{\text{-3}}}\right)}^{-1}\text{nn\&}\times \frac{I_{\nu }}{1\text{Jy}{\text{sr}}^{\text{-1}}},\text{nnend}aligned$$(5)

to convert the free-free intensity into column density in each pixel, with a frequency at ν = 353 GHz, and an electron temperature of T_e = 8000 K. This equation also shows that the H II column density is inverse proportional to the effective density of electrons n_e. Here, we chose n_e = 2 cm⁻³ (Sodroski et al. 1997) to estimate the upper limit of the H II column density. We note that n_e varies in different positions of the Galaxy; however, in this study, most of the H II column densities are from the inner Galaxy. As discussed in Sodroski et al. (1997), the assumption of using n_e = 2 cm⁻³ as the effective electron density in the H II regions in the inner Galaxy produce the same gas-to-dust mass ratio for the H II component as the mean values of the gas-to-dust mass ratio for the H2 and HI components within the same ranges of radial distance. Thus, in the following study we fixed n_e = 2 cm⁻³ to derive the column density of the H II component. The derived H II column density maps of both regions are shown in the bottom panel of Fig. 1.

Fig. 1 Dust and H II templates used in Fermi-LAT data analysis. Top panel: derived total gas column density maps (in units of cm−2) from Planck dust opacity for region I (left) and region II (right). Bottom panel: derived H II gas density maps (in units of cm−2) from the joint analysis of Planck, WMAP, and 408 MHz observations for region I (left) and region II (right).

3 Fermi-LAT data analysis

We collected the Fermi-LAT Pass 8 data from August 4, 2008 (MET 239557417) until October 28, 2020 (MET 625621368), and used the Fermitools from the Conda distribution¹ together with the latest version of the instrument response functions (IRFs) P8R3_SOURCE_V3 for the γ-ray emission analysis. Two 20° × 20° square regions centered at the position of (l = 30°, b = 0°) and (l = 330°, b = 0°) were chosen as two regions of interest (ROI). Here, we selected the “source” class events, then applied the recommended data cut expression (DATA_QUAL > 0)&&(LAT_CONFIG == 1) to exclude time periods when some spacecraft events affected the data quality. Moreover, to reduce the background contamination from the Earth’s albedo, only the events with zenith angles under 90° were included for the analysis.

For each region, the background model includes the sources in the Fermi-LAT ten-year Source Catalog (4FGL-DR2, Ballet et al. 2020) within the ROI enlarged by 10°, and the normalizations and spectral parameters were left free for sources with significance > 5σ within the ROI.

3.1 Diffuse γ-ray emission

As for the diffuse background components, we firstly used the Fermi-LAT Galactic diffuse background model (gll_iem_v07.fits), hereafter referred to as the “Fermi” model. Since we are interested in the diffuse γ-ray emissions, we also built our own Galactic diffuse background models for cross-checking. The diffuse γ-ray emission mainly come from the pion-decay process induced by the inelastic collision between CR protons and ambient gas, as well as the inverse Compton (IC) scattering of CR electrons in the interstellar radiation fields (ISRFs). We calculated the IC component using GALPROP² (Vladimirov et al. 2011), which uses information regarding CR electrons and ISRFs. For the GALPROP webrun we used the parameter set ^SY^Z4^R30^T150^C2 in Ackermann et al. (2012) as fiducial model. We also used the 128 parameter sets in Ackermann et al. (2012) to generate the corresponding IC emission templates; and for each region, we chose the best template by comparing the maximum likelihood achieved by each likelihood fitting process.

For the pion-decay component, as a zero-order approximation, the CRs are homogeneously distributed in the CR sea in our Galaxy. Thus, we can use the gas distribution map to model the pion-decay component of the diffuse γ-ray emission. Asmentioned above, to avoid the uncertainties of dark gas and H I spin temperatures, we used the gas column density maps derived from the Planck dust opacity. Instead of the template gll_iem_v07.fits in the Fermi models, models using total gas distribution maps derived from Planck dust opacity (as shown in the top panels of Fig. 1) and IC emission templates are referred to as “dust” models.

However, since the metallicity varies with galactocentric distance, the dust-to-gas ratio may also vary in the Galaxy. Furthermore, as derived in Acero et al. (2016), the spectra and density of CRs vary in galactocentric distance. According to the results in Acero et al. (2016), both the spectra and density of CRs and the X_CO factor vary significantly below and above 4.5 kpc. Thus, we split the gas into to rings with galactocentric radius below and above 4.5 kpc and include both of them in the pion-decay component, and the gas column density distributions within different galactocentric distances are shown in Fig. 2. In this model, as mentioned above, a dark gas component must be invoked to account for the gas that cannnot be traced by CO and 21 cm observations. Thus, we also added the dark gas component as described in the last section. Models that replace the total gas distribution map with split gas distribution maps and the dark gas component maps (as shown in Fig. 3) are referred to as “gas” models. Finally, as we investigated the possible contribution of an H II component in the diffuse γ-ray emissions, we also added the H II column density map (as shown in the bottom panel of Fig. 1) as a diffuse template in the following analysis. In the following analysis, we used a log-parabola function ($\text{d}N/\text{d}E=N_0{(E/E_{\text{b}})}^{-\alpha -\beta \log (E/E_{\text{b}})}$) as the energy spectrum of all the gas template and left both the normalization and spectral parameters free in the likelihood fitting below.

3.2 Results

Using photons of energies within 0.2–200 GeV, we performed the standard likelihood analysis of Fermi-LAT data with the 4FGL-DR2 sources and diffuse components mentioned above. To check the possible contribution from the H II component, we performed the analysis with and without a template derived from H II column density, respectively. To compare the goodness of fit in these two models, we also calculated the Akaike information criterion (AIC) value for each model. AIC was defined as AIC = −2log(likelihood) + 2k, where k is the number of free parameters in the model. The log(likelihood) and corresponding AIC value for each fit are listed in Table 1. In addition, we compared the overall maximum likelihood of the model with H II template (alternative hypothesis, L) and that of the model without (null hypothesis, L₀) to obtain the test statistics (TS) of the H II component which is defined as − 2log(L₀∕L) following Lande et al. (2012). Then, the corresponding significance (σ_H II) of the H II component can be estimated as $~\sqrt{{\text{TS}}_{H\text{textsc}ii}}$. The TS value and the significance of the additional H II component for each model are also listed in Table 1.

We found that in both regions the inclusion of the H II templates for dust and gas models improves the fit dramatically, and the results are even improved for the Fermi models. However, the overall results reveal that the gas models cannot fit the data as well as the Fermi or dust models. Furthermore, by subtracting the fit model map calculated using gtmodel from the observed γ-ray counts map (as shown in Fig. 4), we generated the residual maps of these different models. Then, we divided the residual counts by the square root of the observed γ-ray counts in each pixel, which should be equal to the significance (in σ) of each pixel approximately. We used these residual significance (signal-to-noise, S/N) maps to compare the fitting quality of these models. The improvement can be also seen from the comparisons of the results with and without the H II templates, which are shown in Figs. A.1 and A.2. In addition, we used a Gaussian function to fit the distribution of the S/N ratios of each model to check whether the model is a statistically good description of the data. As shown in the bottom panel of Figs. A.1 and A.2, each distribution can be well described by the Gaussian function with a mean of ~ 0.0 and standard deviation of ~1.0, which indicates an unbiased significance distribution of the background.

Fig. 2 Gas column density maps (in units of cm−2) below 4.5 kpc (top panel) and above 4.5 kpc (bottom panel) for region I (left) and region II (right), which are derived from Planck dust opacity using method described in Sect. 2.

3.3 Spectral analysis

The TS values of the H II templates in the Fermi model are much smaller than those in the other two models. This is not unexpected, since in the standard Fermi diffuse background model the artificial “patches” are added to absorb the residual in the Galactic plane³, and the standard diffuse model is not suggested to study the diffuse emission regarding the interstellar medium. And the IC contribution is also fixed. Thus, in the spectral analysis we used only the gas and dust templates described above.

We then divided the energy range 0.2–200 GeV into nine logarithmically spaced energy bins and extracted the spectral energy distribution (SED) of diffuse γ rays in these regions via the maximum likelihood analysis in each energy bin. The uncertainties include 68% statistical errors for the energy flux and the systematic errors due to the uncertainties in LAT effective areas. The results for the dust model and dust+H II models, as well as gas model and gas+H II models, are shown in Figs. 5–8, respectively.

In Figs. 6–8, we have divided the SEDs with their corresponding gas column densities derived from Sect. 2 to obtain the γ-ray emissivities per H atom, which is proportional to the CR density. The derived emissivities are compared with the one predicted by using the local interstellar spectrum (LIS) of CRs (black line, Casandjian 2015). We found that, in the dust model without including H II templates, the derived CR spectrum from the dust model of each region is significantly harder than the LIS, which is consistent with the results for the inner Galaxy in Yang et al. (2016) and Acero et al. (2016). However, for the dust+H II models, the derived γ-ray emissivities associated with the dust templates are in good agreement with the local emissivities, which imply similar CR spectra in these regions and the solar neighbourhood. Meanwhile, the derived emissivities associated with the H II component are much harder, with an index of about − 2.1.

For the gas model without the H II templates, the derived γ-ray emissivities for the inner rings (below 4.5 kpc) are significantly higher than the LIS value and the outer rings (above 4.5 kpc), while the spectra in the inner rings is significantly harder. These results are also in good agreement with the results of Yang et al. (2016) and Acero et al. (2016). When including H II templates, for region I the γ-ray emissivities for both inner and outer gas rings are consistent with the LIS value in both the spectra and normalizations. For region II, the γ-ray emissivities’ spectral indices in both rings are consistent with the LIS value, while the normalization in the inner ring is higher. The γ-ray emissivities from H II templates are also significantly harder in both regions, and are consistent with the results derived from dust+H IImodels in both spectra and normalizations.

It should be noted that the normalization of the derived emissivities associated with H II component is uncertain. This is because the ionized gas only contributes a small portion to the total gas column in our Galaxy. Moreover, it is possible that the CRs associated with ionized gas also illuminate other gas components. In Figs. 6 and 8, we assumed only H II gas are responsible for the γ rays associated with the H II templates. Thus, the emissivities derived here can be regarded as the upper limits.

Fig. 3 Dark gas templates for region I (left) and region II (right), which are derived from dust data using method described in Sect. 2.

Fig. 4 Observed counts map of γ rays in the energy range of 0.2–200 GeV for region I (left) and region II (right).

Fig. 5 SEDs of diffuse γ-ray emissions associated with different components for region I (left panel) and region II (right panel). The black and red points show the γ rays associated with total gas (the dust template) for both dust models and dust+H II models, respectively. The blue points represent the γ rays related to H II gas (the H II template) obtained for dust+H II models.

Fig. 6 Same as Fig. 5, but the SEDs have been normalized to the γ-ray emissivity per H atom, which is proportional to the CR density. Also plotted as black curves are the predicted γ-ray emissivities for the local CR spectra (LIS).

Fig. 7 γ-ray emissivities per H atom of diffuse γ-ray emission associated with to the two gas rings and for gas models for region I (left panel) and region II (right panel), respectively.

4 Discussion

In this work, we found a hard γ-ray component associated H II in two regions in the Galactic plane. The derived intensity and spectrum of such a γ-ray component is stable when adopting different gas tracers and IC models. Including such a component will also alter the γ-ray emissivities associated with the gas in the Galactic plane. This can be seen in both the dust and gas models in our analysis. In the dust model, the inclusion of the H II component will deduce a similar γ-ray emissivities in the whole line of sight as the local value, which implies that the average CR density and spectra in the line of sight is similar to the local CR density and spectrum. Such a conclusion is different from those in Acero et al. (2016) and Yang et al. (2016). Moreover, in the gas model, the gases are split into two rings with galactocentric radius below and above 4.5 kpc. Without the H II component, we reproduced the results in Acero et al. (2016) and Yang et al. (2016); thus, the γ-ray emissivities per H atom are higher in the inner ring and reveal a significantly harder spectrum. Including the H II template, we found the γ-ray emissivities in both rings become consistent with the local value for region I. For region II, the spectra are similar to the local value, while for the inner rings the normalizations are still significantly higher. However, we note that the normalization of γ-ray emissivities per H atom is inversely proportional to the X_co conversion factor, which is fixed in our analysis. Indeed, X_co depends on the metallicity and can be different in different part of the Galaxy (Acero et al. 2016).

To conclude, we found a robust hard γ-ray emission component associated with the H II gas, and the γ-ray emissions associated with other gas reveal a significant softer spectrum with similar spectral shape. In the following discussion, we label the γ-ray emission and corresponding CR population related to the dust and gas templates and H II templates as “soft” and “hard” components, respectively. As for the hard component, the absolute normalization can not be determined decisively, as we discussed above. By using the upper limit of the emissivities in the most extreme case, the hard component are only produced in the H II region, the derived CR density can be 100 times higher than the LIS value (see Figs. 6 and 8). Such enhancement of CR densities is expected near CR sources such as SNRs and YMCs. The dimension of the regions with enhanced CR density can be of more than 100 pc (Aharonian et al. 2019). Thus in this case, the hard γ-ray component associated with H II gas can be interpreted as the enhanced γ-ray emission in the vicinity of sources. To illustrate such a scenario, we plotted in Fig. 9 the γ-ray spectrum assuming a power-law parent proton spectrum with an index of − 2.2, which is similar to the CR spectrum derived in the vicinity of YMCs (Aharonian et al. 2019). On the other hand, the hard γ-ray emission can also have a leptonic origin, especially considering the enhanced radiation fields in these H II regions. The interstellar radiation fields (ISRFs), which provide the low-energy photon targets for the IC process, can be divided into three components in our Galaxy (Popescu et al. 2017), including the optical and UV fields from the star light, the infrared fields from the dust emission, and the cosmic microwave backgrounds. In the H II regions, the ionizing massive stars will increase the optical and UV fields significantly and thus produce additional IC emissions. We calculated such a component assuming the electron spectrum is the same as the electron spectrum measured in the solar neighborhood and from radio observations (Strong et al. 2011), and the optical and UV components are described by a grey-body spectrum with a temperature of 5000 K. The derived IC γ-ray fluxes are also shown in Fig. 9. In the calculated γ-ray flux from both pion-decay and the IC component, the total normalization is rather arbitrary, since both the average gas density and the energy densities of the enhanced optical and UV radiation fields are unknown. But the derived spectral shapes already reveal significant differences. The IC components, due to the Klein-Ninisha effects, are softened above dozens of GeVs. The current data cannot rule out the leptonic origin, but further observations above 100 GeV may distinguish two such possibilities. If there is another hard electron spectrum in these H II regions, the IC component can be harder than the curve calculated here, but the Klein-Nishina effects for the optical and UV components will inevitably play a role above several hundred GeV in the γ-ray spectrum (Popescu et al. 2017). Thus, the observations above several hundred GeV may distinguish such two scenarios.

Interestingly, a two-component CR model was used in Yang & Aharonian (2019) to explain the hardening of CR spectrum above 200 GV, as well as the increasing positron/electron and antiproton/proton ratios. In such a scenario, CRs are injected from two population of sources with different spectral properties (more specifically, a 1st source population with soft spectra and a 2nd source population with hard spectra). For the second source population, CRs accumulate the grammage of about 1 g cm⁻² in the vicinity of the sources. Accordingly, they should be effective γ-ray emitters. Yang & Aharonian (2019) calculated diffuse γ-ray emissions associated with the second (hard) source population and found it can dominate over the first (soft) source population above about 100 GeV, which is in agreement with our measurements shown in Fig. 5. As mentioned in Yang & Aharonian (2019), one possible candidate for the second source population can be YMC. Indeed, the OB star distributions do have a strong peak at about 4 kpc away from the GC (Bronfman et al. 2000). Therefore, we expect a strong enhancement of the hard γ-ray emissions associated with these sources in this region, which can naturally explain the high energy density and hard spectrum of CRs derived from the 4–6 kpc rings in our Galaxy (Acero et al. 2016; Yang et al. 2016).

On the other hand, the similar density and spectral shape of the soft CR component with the LIS reveal the possible CR sea all over the Galaxy. Such a possibility is also unveiled in Aharonian et al. (2020), in which the γ-ray observations on individual molecular clouds are used to measure the CR density in different positions of our Galaxy. Such a uniform CR distribution is not trivial, since if the CR propagation is dominated by diffusion, the spatial distribution of the sources should always leave a trace in the spatial distribution of CRs (Strong & Moskalenko 1998). However, all the promising candidates of CR accelerators, such as SNRs (Green 2015), pulsars (Yusifov & Küçük 2004), and OB stars (Bronfman et al. 2000), reveal a significant inhomogeneous distribution in our Galaxy. Thus, a larger CR halo or strong re-acceleration may be needed to explain the uniform distribution of the CRs, which will lead to a significant modification to the current paradigm of CR propagation in our Galaxy.

Finally, we note that although the fitting of γ-ray data is significantly improved by introducing the H II gas component, there are still strong residuals (see Figs. A.1 and A.2). This is not unexpected if we assume the hard γ-ray component is from the vicinity of CR sources. Since the CR distribution should be significantly inhomogeneous near the source, we expect no perfect spatial correlation between the gas and the γ-ray emissions. In this regard, a detailed analysis in each compact H II region may be helpful. These regions are often extremely crowded and different sources can be confused and overlapping in the line of sight (Liu et al. 2019). However, as shown in Aharonian et al. (2019) these sources can also be PeVatrons. If this is the case, the corresponding γ-ray emission should be higher than 10 TeV. In this energy band, the background due to the diffuse γ-ray emission associated with soft CR components and individual γ-ray sources suchas pulsars, SNRs, and pulsar wind nebulae can be strongly suppressed. The extensive shower arrays such as LHAASO (Sciascio 2016) and HAWC (Abeysekara et al. 2013), as well as the Cerenkov telescope arrays such as H.E.S.S., MAGIC, VERITAS, and CTA would be ideal instruments to explore these regions.

Fig. 8 γ-ray emissivities per H atom of diffuse γ-ray emission associated with to the two gas rings obtained for dust +H II models for region I (left panel) and region II (right panel), respectively.

Fig. 9 γ-ray flux associated with H II gas obtained for dust+H II models (blue) and gas+H II models, for region I (left panel) and region II (right panel), respectively. Also plotted are the possible contributions from IC and pion decay processes. The details can be found in the text.

Acknowledgements

B.L. is supported by the NSFC under grant 12103049. R.Y. is supported by the NSFC under grants 11421303, 12041305 and the national youth thousand talents program in China.

Appendix A Residual significance (signal-to-noise, S/N) maps for regions I and II

Fig. A.1 Comparison of the residual significance (signal-to-noise, S/N) maps in the energy range of 0.2–200 GeV for models without an H II component (top) and with an H II component (middle) of region I. The left panels show the results of Fermi and Fermi+H II models, the middle panels show the results of dust and dust+H II models, and the right panels illustrate the results of gas and gas+H II models. The color scales are from − 5σ to 5σ. The bottom panel show the distribution of significance (bin size is 0.05 σ of each model without or without the HII component, in which each solid line illustrates the Gaussian function fitted to the corresponding distribution and μ and σ represent the mean and the standard deviation of each fit.

Fig. A.2 Same as Fig. A.1, but for region II.

References

All Tables

All Figures

	Fig. 1 Dust and H II templates used in Fermi-LAT data analysis. Top panel: derived total gas column density maps (in units of cm−2) from Planck dust opacity for region I (left) and region II (right). Bottom panel: derived H II gas density maps (in units of cm−2) from the joint analysis of Planck, WMAP, and 408 MHz observations for region I (left) and region II (right).
In the text

	Fig. 2 Gas column density maps (in units of cm−2) below 4.5 kpc (top panel) and above 4.5 kpc (bottom panel) for region I (left) and region II (right), which are derived from Planck dust opacity using method described in Sect. 2.
In the text

	Fig. 3 Dark gas templates for region I (left) and region II (right), which are derived from dust data using method described in Sect. 2.
In the text

	Fig. 4 Observed counts map of γ rays in the energy range of 0.2–200 GeV for region I (left) and region II (right).
In the text

	Fig. 5 SEDs of diffuse γ-ray emissions associated with different components for region I (left panel) and region II (right panel). The black and red points show the γ rays associated with total gas (the dust template) for both dust models and dust+H II models, respectively. The blue points represent the γ rays related to H II gas (the H II template) obtained for dust+H II models.
In the text

	Fig. 6 Same as Fig. 5, but the SEDs have been normalized to the γ-ray emissivity per H atom, which is proportional to the CR density. Also plotted as black curves are the predicted γ-ray emissivities for the local CR spectra (LIS).
In the text

	Fig. 7 γ-ray emissivities per H atom of diffuse γ-ray emission associated with to the two gas rings and for gas models for region I (left panel) and region II (right panel), respectively.
In the text

	Fig. 8 γ-ray emissivities per H atom of diffuse γ-ray emission associated with to the two gas rings obtained for dust +H II models for region I (left panel) and region II (right panel), respectively.
In the text

	Fig. 9 γ-ray flux associated with H II gas obtained for dust+H II models (blue) and gas+H II models, for region I (left panel) and region II (right panel), respectively. Also plotted are the possible contributions from IC and pion decay processes. The details can be found in the text.
In the text

Fig. A.1 Comparison of the residual significance (signal-to-noise, S/N) maps in the energy range of 0.2–200 GeV for models without an H II component (top) and with an H II component (middle) of region I. The left panels show the results of Fermi and Fermi+H II models, the middle panels show the results of dust and dust+H II models, and the right panels illustrate the results of gas and gas+H II models. The color scales are from − 5σ to 5σ. The bottom panel show the distribution of significance (bin size is 0.05 σ of each model without or without the HII component, in which each solid line illustrates the Gaussian function fitted to the corresponding distribution and μ and σ represent the mean and the standard deviation of each fit.

In the text

	Fig. A.2 Same as Fig. A.1, but for region II.
In the text