Probing the Sea of Cosmic Rays by Measuring Gamma-Ray Emission from Passive Giant Molecular Clouds with HAWC

Direct observations of cosmic rays measure the spectrum surrounding the vicinity of the solar system (Aguilar et al. 2015). It is generally assumed that this spectrum is representative of the cosmic-ray flux across the Galaxy, and it is usually referred to as the "sea" of cosmic rays. This assumption comes from the fact that cosmic rays, after leaving the acceleration regions, diffuse through the Galaxy owing to the deflecting effect of the interstellar magnetic field to alter the path of charged particles. If this process occurs on a timescale of millions of years, their distribution becomes homogeneous and isotropic (Strong et al. 2007).

A way to probe the paradigm of the "sea" of cosmic rays is by using gamma-ray data. Since the beginning of high-energy gamma-ray astrophysics, it has been recognized that the observation of gamma rays gives us the opportunity to measure the propagation and distribution of cosmic rays in distant parts of the Galaxy, since these can interact with the interstellar matter and radiation fields generating gamma-ray emission. This has been suggested extensively in the literature by Black & Fazio (1973), Issa & Wolfendale (1981), Aharonian (2001), and Casanova et al. (2010). The premise is to look for regions that act as targets for cosmic rays: clouds that are part of the interstellar medium. We focus on passive giant molecular clouds (GMCs), i.e., clouds with no cosmic-ray accelerators inside of them (Issa & Wolfendale 1981). These clouds are large reservoirs of gas (mainly molecular hydrogen) and dust. Masses of these complexes are of the order of ∼${10}^{5}\,{M}_{\odot }$. Their average density is 100 to 1000 times greater than the average density in the solar vicinity, which is of the order of 1–10 particles per cubic centimeter (Ferrière 2001). It is assumed that, due to the low probability of cosmic-ray sources inside the GMCs, the main gamma-ray production mechanism in these sources is the decay of neutral pions, which are produced by the collision of high-energy cosmic rays with interstellar matter in the clouds.

The flux of gamma rays is proportional to the flux of cosmic rays ${{\rm{\Phi }}}_{\mathrm{CR}}$, as well as the total mass M of the GMC, and inversely proportional to the distance square d² of the GMC (see, e.g., Casanova et al. 2010):

$$\begin{eqnarray}&&{F}_{\gamma }\propto \displaystyle \frac{M}{{d}^{2}}{{\rm{\Phi }}}_{\mathrm{CR}}.\end{eqnarray} \tag{ 1 }$$

It is then important to know the masses and distances of the GMCs. Equation (1) has been used in other analyses using data from the Fermi-LAT.

The Fermi-LAT collaboration presented their studies to understand the cosmic-ray propagation through gamma-ray data of the Orion molecular cloud (Ackermann et al. 2012b) and of the molecular clouds Chamaleon, R Coronae Australis, Cepheus, and Polaris (Ackermann et al. 2012a). In both publications the authors calculated gamma-ray emissivities and molecular mass conversion factors X_CO.

Other studies using Fermi-LAT data include the work by Neronov et al. (2012) and Neronov et al. (2017), where they look at high galactic latitude clouds that are part of the Gould Belt. They calculate the average cosmic-ray spectrum in the Galaxy.

Another example was presented in Yang et al. (2014), where gamma-ray spectra above 300 MeV were used to extract the cosmic-ray spectra from eight massive clouds. They showed that the derived spectral indices and absolute fluxes of cosmic-ray protons in the energy interval 10–100 GeV agree with the direct measurements of local cosmic rays by the PAMELA experiment (Adriani et al. 2011).

A similar study was published by Aharonian et al. (2020). With their observations, they found that the flux of cosmic rays at distances from 0.6 to 12.5 kpc also agrees with the locally measured cosmic-ray flux measured by the Alpha Magnetic Spectrometer (AMS; Aguilar et al. 2015).

A recent publication by Baghmanyan et al. (2020) also used Fermi-LAT data to find discrepancies between the measured cosmic-ray flux at Earth and the cosmic-ray flux in distant regions of the Galaxy. They report that in three GMCs (Aquila Rift, Ophiuchi, and Cepheus) they observe higher gamma-ray emission with respect to the expected emission from the "sea of cosmic rays." Most of these publications have presented results of gamma rays up to 300 GeV.

While other TeV instruments have measured the cosmic-ray flux from GMCs from regions of the Galactic plane (e.g., Abdo et al. 2008; Bartoli et al. 2012; Abramowski et al. 2016), the present work is the first attempt to measure the "sea of cosmic rays" using high-energy gamma rays above 1 TeV from high-latitude GMCs. In this paper we present measurements of gamma rays that have median energies of the order of ∼1–10 TeV using data from the High Altitude Water Cerenkov (HAWC) Observatory.

The paper outline is as follows. We present the GMCs under study in Section 2. The detector and the data set are described in Section 3.1. The analysis method is explained in Section 3.2. The results and discussion are shown in Section 4. Finally, we present our concluding remarks in Section 5.

Following a similar procedure to the one presented in Yang et al. (2014), we selected seven massive clouds identified by the CO Galactic survey of Dame et al. (1987) with the CfA 1.2 m millimeter-wave Telescope. The selected clouds are Taurus, Orion, Perseus, Ophiuchi, Monoceros, Aquila Rift, and Hercules. All of them are located less than 1 kpc away from the solar system. These clouds have high galactic latitudes ($| b| \gt 5^\circ$), and they are in the field of view of HAWC. Their masses and distances are listed in Table 1. The distances are obtained from Dame et al. (1987) and Schlafly (2014). The masses are calculated using the Planck survey (Ade et al. 2011; see Appendix A for more details).

Table 1. Properties of the GMCs

GMC	Mass(±14%) $({10}^{5}\,{M}_{\odot })$	Distance (pc)	A	l (deg)	b (deg)	Decl. (deg)	R.A. (deg)
Taurus	0.11	140 ± 30	5.64	171.6	−15.8	26.2	66.5
Orion	1.62	490 ± 50	6.76	−151.0	−15.5	−3.6	87.3
Perseus	0.15	315 ± 32	1.55	159.0	−20.5	30.9	52.9
Ophiuchi	0.06	125 ± 18	3.94	−5.7	16.8	−23.5	247.5
Monoceros	0.6	830 ± 83	0.87	−145.8	−12.4	−6.8	92.3
Aquila Rift	0.77	225 ± 55	15.27	22.5	11.8	−3.72	267.6
Hercules	0.05	200 ± 30	1.16	46.0	9.0	15.7	280.6

Note. Masses are calculated using the Planck survey (Ade et al. 2011) (see Appendix A for mass calculations). Distances were obtained from Dame et al. (1987) and Schlafly (2014). The sky positions correspond to the center of the regions of interest that are analyzed.

Download table as:  ASCIITypeset image

The first four clouds that appear in Table 1 (Taurus, Orion, Perseus, and Ophiuchi) are considered as part of the Gould Belt according to Dzib et al. (2018). The Gould Belt is a region in the Galaxy composed of gas and young stars. In the sky it appears as a ring with an inclination angle of ∼20° with respect to the galactic plane, although, according to Guillout et al. (1998), it can be considered a disk-like structure. The clouds Aquila and Hercules do not line up with this structure, and the Monoceros cloud is too far to be considered part of the belt. ³¹ However, a recent study by Alves et al. (2020) has put into question the Gould Belt structure and substitutes it by the so-called Radcliff Wave. In this case, the first three clouds in Table 1 are part of the Radcliff Wave, while the fourth cloud is part of the Local arm. Assuming that these clouds have similar properties, we combine the gamma-ray emission of the clouds that belong to the Gould Belt and the Radcliff Wave in a stacked analysis to increase the sensitivity of the probe. ³² Due to the uncertainty of this classification, we also perform the stacked analysis for all the clouds. For completeness, we analyze each of them individually.

3.1. The HAWC Data

3.2. Analysis Method

The analysis is performed using the Multi-Mission Maximum Likelihood (3ML; Vianello et al. 2015) framework, together with the HAWC Accelerated Likelihood (HAL) plug-in. ³³ We perform both individual cloud analyses and a stacked analysis and assume the same spectral model for each cloud as a function of energy with the form

$$\begin{eqnarray}&&{F}_{\gamma }(E)=A\,C{\left(\displaystyle \frac{E}{{E}_{0}}\right)}^{-\alpha },\end{eqnarray} \tag{ 2 }$$

where α is the spectral index, E₀ is the pivot energy, and C is the normalization factor. The factor A, defined as ${M}_{5}/{d}_{\mathrm{kpc}}^{2}$ (${M}_{5}=M/{10}^{5}{M}_{\odot };$ ${d}_{\mathrm{kpc}}=d/1\mathrm{kpc}$), expresses the gamma-ray visibility of a cloud of mass M, located at a distance d (Aharonian 2001; Casanova et al. 2010). This factor quantifies the weighted contribution of each cloud when we perform the stacking log-likelihood analysis.

We fixed the spectral index to 2.7 assuming that gamma rays at these energies still follow the spectral shape of the "sea" of cosmic rays (Aguilar et al. 2015). We also include morphology information from each cloud based on templates generated using the Planck survey (see Appendix B on how we built the templates). This information is combined, together with HAWC's detector response, to obtain an expected number of events in each ${f}_{\mathrm{hit}}$ bin used in the likelihood calculation.

For each cloud we create a log-likelihood profile as a function of C. For individual clouds, we use this profile to find the value $\hat{C}$ that optimizes the log-likelihood. For the stacking analysis, we add the log-likelihood profiles and then find the value $\hat{C}$ that optimizes this profile. During the optimization, we constrain the normalization value to be positive.

A test statistic value TS is calculated to check for a significant excess given the source model evaluated in the entire region of interest:

$$\begin{eqnarray}&&\mathrm{TS}=2\mathrm{ln}\left[\displaystyle \frac{{ \mathcal L }(S(\hat{C})+B)}{{ \mathcal L }(B)}\right],\end{eqnarray} \tag{ 3 }$$

where S is the signal from the source model, while B is the background model, which corresponds to the background estimated from the data using direct integration (see Abeysekara et al. 2017b).

The best value $\hat{C}$ is then used as an input to a Markov Chain Monte Carlo (MCMC) to estimate a probability distribution. Since we observe a nonsignificant detection (i.e., TS < 25), we calculate the 95% credible interval upper limits from the estimated distribution of the parameter $\hat{C}$.

We compute a limit for the ${E}^{-2.7}$ spectrum model (referred to afterward as the 2.7-model) and quasi-differential ³⁴ limits on the normalization in a similar way to previous HAWC analyses (Abeysekara et al. 2017a, 2018). The main difference between these two limits is that the quasi-differential limits are independent of any model assumption, while the limit for a broad energy range is model dependent. However, the quasi-differential limits are more conservative owing to the lowering of the statistical power by dividing the data. For the 2.7-model limit, we use the energy band going above 1 TeV. For the quasi-differential limits, we define three energy half-decade bands, between 1 and 31.6 TeV, and an overflow bin above 31.6 TeV. The range of the bands and the pivot energy, E₀, used in each range are shown in Table 2. The pivot energy is set at 10 TeV for the 2.7-model limit.

Table 2. Observed 95% Credible Upper Limits on the Gamma-Ray Emission from the GMCs (and on $\hat{C}$), Together with the Expected Median Upper Limits, as well as the 68% and 95% Containment Bands for the Expected Limits

Energy Range	1–3.16 TeV		3.16–10.0 TeV		10.0–31.6 TeV		>31.6 TeV		>1 TeV
Pivot Energy	1.78 TeV		5.62 TeV		17.8 TeV		56.2 TeV		10 TeV
	(×10−12)		(×10−13)		(×10−15)		(×10−16)		(×10−15)
	U.L.	Expected Limit	U.L.	Expected Limit	U.L.	Expected Limit	U.L.	Expected Limit	U.L.	Expected Limit
All	5.0(0.1)	5.8(${}_{-1.6}^{+2.1},{\,}_{-2.7}^{+5.0}$)	2.6(0.1)	2.7(${}_{-0.8}^{+1.1},{\,}_{-1.2}^{+2.5}$)	13.6(0.4)	14.0(${}_{-4.1}^{+5.3},{\,}_{-6.5}^{+12.0}$)	9.5(0.3)	9.6(${}_{-2.8}^{+3.7},{\,}_{-4.5}^{+8.2}$)	19.0(0.5)	23.0(${}_{-6.9}^{+8.8},{\,}_{-10.0}^{+16.0}$)
Gould Belt	3.6(0.2)	3.7(${}_{-1.0}^{+1.5},{\,}_{-1.6}^{+3.6}$)	2.3(0.1)	1.8(${}_{-0.5}^{+0.7},{\,}_{-0.8}^{+1.1}$)	13.5(0.8)	8.7(${}_{-2.4}^{+3.3},{\,}_{-4.0}^{+7.3}$)	9.4(0.7)	6.1(${}_{-1.5}^{+2.2},{\,}_{-2.8}^{+5.4}$)	18.0(1.0)	13.5(${}_{-3.9}^{+4.8},{\,}_{-6.1}^{+11.0}$)
Radcliff	2.8(0.2)	2.9(${}_{-0.8}^{+1.1},{\,}_{-1.3}^{+2.2}$)	1.8(0.1)	1.4(${}_{-0.4}^{+0.5},{\,}_{-0.7}^{+1.2}$)	10.3(0.6)	7.0(${}_{-1.9}^{+2.7},{\,}_{-3.2}^{+5.8}$)	7.4(0.5)	4.9(${}_{-1.4}^{+1.9},{\,}_{-2.3}^{+4.1}$)	13.8(1.0)	11.0(${}_{-3.5}^{+4.0},{\,}_{-5.0}^{+7.0}$)
Taurus	1.5(0.3)	1.5(${}_{-0.4}^{+0.6},{\,}_{-0.7}^{+1.4}$)	1.0(0.2)	0.8(${}_{-0.2}^{+0.3},{\,}_{-0.4}^{+0.6}$)	6.4(1.1)	4.1(${}_{-1.1}^{+1.7},{\,}_{-1.9}^{+3.4}$)	5.0(0.9)	3.1(${}_{-0.9}^{+1.3},{\,}_{-1.4}^{+2.8}$)	7.8(1.4)	5.6(${}_{-1.6}^{+2.2},{\,}_{-2.6}^{+5.3}$)
Orion	2.0(0.3)	2.6(${}_{-0.8}^{+0.9},{\,}_{-1.2}^{+2.3}$)	1.0(0.1)	1.0(${}_{-0.2}^{+0.5},{\,}_{-0.4}^{+1.0}$)	4.8(0.7)	5.2(${}_{-1.4}^{+2.0},{\,}_{-2.4}^{+4.3}$)	3.0(0.4)	3.3(${}_{-0.9}^{+1.4},{\,}_{-1.5}^{+3.1}$)	7.1(1.0)	8.2(${}_{-2.4}^{+3.1},{\,}_{-3.8}^{+7.1}$)
Perseus	1.4(0.9)	1.0(${}_{-0.3}^{+0.3},{\,}_{-0.5}^{+0.8}$)	0.74(0.5)	0.5(${}_{-0.2}^{+0.1},{\,}_{-0.3}^{+0.3}$)	4.0(2.6)	2.4(${}_{-0.7}^{+1.1},{\,}_{-1.1}^{+2.0}$)	2.9(1.9)	1.7(${}_{-0.4}^{+0.7},{\,}_{-0.7}^{+1.4}$)	5.8(3.7)	3.4(${}_{-1.0}^{+1.3},{\,}_{-1.6}^{+2.8}$)
Ophiuchi [x100]	4.7(1.2)	6.9(${}_{-1.9}^{+2.8},{\,}_{-3.1}^{+6.1}$)	1.4(0.4)	1.9(${}_{-0.5}^{+0.8},{\,}_{-0.8}^{+1.7}$)	5.2(1.3)	7.8(${}_{-2.3}^{+2.2},{\,}_{-3.9}^{+5.2}$)	2.5(0.6)	3.0(${}_{-0.7}^{+1.2},{\,}_{-1.3}^{+2.9}$	8.6(2.2)	11.0(${}_{-2.7}^{+5.0},{\,}_{-4.8}^{+9.0}$)
Monoceros	1.9(2.2)	1.5(${}_{-0.4}^{+0.6},{\,}_{-0.7}^{+1.3}$)	0.5(0.6)	0.6(${}_{-0.1}^{+0.3},{\,}_{-0.3}^{+0.6}$)	1.7(2.0)	3.1(${}_{-0.9}^{+1.1},{\,}_{-1.5}^{+2.4}$)	1.2(1.4)	1.9(${}_{-0.6}^{+0.7},{\,}_{-0.9}^{+1.7}$)	3.4(3.9)	4.6(${}_{-1.3}^{+1.7},{\,}_{-2.1}^{+3.8}$)
Aquila	3.2(0.2)	4.3(${}_{-1.2}^{+1.9},{\,}_{-2.0}^{+3.5}$)	1.3(0.1)	2.0(${}_{-0.5}^{+0.7},{\,}_{-0.9}^{+1.9}$)	6.9(0.5)	10.1(${}_{-3.2}^{+3.8},{\,}_{-4.8}^{+8.0}$)	4.2(0.3)	6.5(${}_{-1.8}^{+2.7},{\,}_{-2.9}^{+6.7}$)	9.8(0.6)	14.7(${}_{-4.3}^{+5.8},{\,}_{-6.9}^{+13.0}$)
Hercules	1.0(0.9)	1.0(${}_{-0.3}^{+0.4},{\,}_{-0.5}^{+0.9}$)	0.5(0.4)	0.5(${}_{-0.1}^{+0.2},{\,}_{-0.2}^{+0.5}$)	3.2(2.8)	3.0(${}_{-0.9}^{+1.1},{\,}_{-1.4}^{+2.5}$)	2.9(2.5)	2.3(${}_{-0.6}^{+0.9},{\,}_{-1.0}^{+1.9}$)	4.2(3.6)	4.1(${}_{-1.2}^{+1.6},{\,}_{-1.9}^{+3.5}$)

Note. The upper limits on the normalization $\hat{C}$ in parentheses are obtained by dividing the gamma-ray flux upper limits by the corresponding factor A. Flux units are in TeV⁻¹ cm⁻² s⁻¹. Gould Belt includes Taurus, Orion, Perseus, and Ophiuchi; Radcliff includes Taurus, Orion, and Perseus. Notice that the flux values of Ophiuchi are 100 times higher.

Download table as:  ASCIITypeset image

3.2.1. Significance Sky Maps

We produce significance sky maps of the GMC regions to complement the results. We use the same analysis method described in Albert et al. (2020). We calculate a likelihood ratio and then a test statistic value TS similar to Equation (3). For the signal model, we assume a ${E}^{-2.7}$ spectrum and a fixed source morphology. The morphology used is a point source. If the null hypothesis is true, the significance distribution can be approximated as a Gaussian distribution for the case of the point-source morphology. The maps and the distributions are shown in Figure 6 in Appendix C.

We did not find significant emission from the studied regions (see tables in Appendix C). Therefore, we proceed to calculate the 95% credible intervals (C.I.) as explained in Section 3.2. We also calculate the expected upper limits by performing the analysis with Poisson-fluctuated background-only maps. We obtain the 1σ and 2σ confidence bands, as well as the median of the upper limits expected from a background-only scenario.

Figure 1 shows the quasi-differential upper limits. Figure 2 shows the 2.7-model upper limits at 10 TeV. These measurements can be compared with the gamma-ray flux expected from hadronic interactions producing gamma rays from neutral pion decay. The models were produced by using the parameterization of the pp cross-section ($\tfrac{d\sigma }{{{dE}}_{\gamma }}$) process from Kafexhiu et al. (2014). We use the cosmic-ray measurements from the AMS experiment (Aguilar et al. 2015) and extrapolate the AMS fit to higher cosmic-ray energies, with the assumption that it will maintain its spectral shape. We calculate the expected gamma-ray flux by convolving the cosmic-ray flux with the pp cross section (Aharonian et al. 2020) ³⁵ :

$$\begin{eqnarray}&&{F}_{\gamma }({E}_{\gamma })=\displaystyle \frac{{\xi }_{N}}{{m}_{p}}A\int {{dE}}_{p}\displaystyle \frac{d\sigma }{{{dE}}_{\gamma }}{F}_{p}({E}_{p}),\end{eqnarray} \tag{ 4 }$$

where A is the same factor as Equation (2), ${\xi }_{N}$ is the nuclear enhancement factor and assumed to be equal to 1.8 as in Aharonian et al. (2020), and m_p is the proton mass.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As can be seen in Figure 1, none of the limits are constraining the expected gamma-ray flux in the energy range accessible to HAWC. In Figure 2 we see that if we stack all the clouds, the 2.7-model 95% C.I. upper limit excludes part of the confidence interval of the gamma-ray emission expected from a cosmic-ray flux equal to that measured at Earth. However, there is not enough evidence to reject the paradigm of an unchanging sea of cosmic rays permeating our Galaxy. The weaker upper limit constraints measured in Ophiuchi are due to the fact that this region in the sky is at high zenith angles for HAWC, where the detector has lower sensitivity. None of the limits constrain the cosmic-ray flux to a lower value than that measured locally. Table 2 contains the observed and expected upper limits on the gamma-ray flux for each of the clouds and stacked analyses.

4.1. Systematic Uncertainties

4.2. Constraints on the Cosmic-Ray Energy Density

With the assumption that the gamma-ray emission is produced by pion decay, we estimate constraints in the cosmic-ray energy density from these distant regions and compare them to those measured in the local neighborhood. Using the same expression as in Abramowski et al. (2016), and using the same energy bins as in Table 2, we can modify our spectral model, Equation (2), so that our free parameter is the cosmic-ray energy density (in $\mathrm{eV}\,{\mathrm{cm}}^{-3}$)

$$\begin{eqnarray}&&{\rho }_{\mathrm{CR}}=1.8\times {10}^{-2}\,{\left(\displaystyle \frac{{\xi }_{N}}{1.5}\right)}^{-1}\left(\displaystyle \frac{{L}_{\gamma }}{{10}^{34}\mathrm{erg}/{\rm{s}}}\right){\left(\displaystyle \frac{M}{{10}^{6}{M}_{\odot }}\right)}^{-1},\end{eqnarray} \tag{ 5 }$$

where ${\xi }_{N}$ accounts for nuclei heavier than hydrogen in both cosmic rays and interstellar matter, ${L}_{\gamma }$ is the luminosity of gamma rays, and M is the mass of the region. The luminosity is obtained as

$$\begin{eqnarray}&&{L}_{\gamma }({E}_{0}\lt {E}_{\gamma }\lt {E}_{f})=4\pi {d}^{2}{\int }_{{E}_{0}}^{{E}_{f}}{E}_{\gamma }F({E}_{\gamma }){{dE}}_{\gamma },\end{eqnarray} \tag{ 6 }$$

where d is the distance to the region. We then insert Equation (6) into Equation (5). After doing the algebra and taking care of units, we obtain

$$\begin{eqnarray}&&{\rho }_{\mathrm{CR}}=4.930\,C\,{E}_{0}^{\alpha }\,({E}_{i}^{-\alpha +2}-{E}_{f}^{-\alpha +2}),\end{eqnarray} \tag{ 7 }$$

where α is the index, E₀ is the midpoint of the energy bin, and E_i , E_f are the edges of the energy bins. Equation (2) now can be rewritten so that the free parameter is ${\rho }_{\mathrm{CR}}$. We apply the procedure in Section 3.2 to get an estimate of the average of the cosmic-ray energy density in the local Galaxy.

We compare the upper limits with our sensitivity and the cosmic-ray energy density measured by AMS. We used the measurements and fit function in Aguilar et al. (2015) to calculate the blue shaded region in Figure 3 (and Figure 5 in Appendix C), which corresponds to the ±2σ uncertainty region of the extrapolation. Table 4 contains the observed and expected upper limits on the cosmic-ray energy density for each of the clouds and stacked analyses.

Zoom In Zoom Out Reset image size

Table 4. Observed 95% Credible Upper Limits on the Cosmic-Ray Energy Density from the GMCs, Together with the Expected Median Upper Limits, as Well as the 68% and 95% Containment Bands for the Expected Limits

Energy Range	1–3.16 TeV		3.16–10.0 TeV		10.0–31.6 TeV		>31.6 TeV
Pivot Energy	1.78 TeV		5.62 TeV		17.8 TeV		56.2 TeV
	U.L.	Expected Limit	U.L.	Expected Limit	U.L.	Expected Limit	U.L.	Expected Limit
All	1.9	1.1(${}_{-0.3}^{+0.5},{\,}_{-0.6}^{+1.0}$)	0.8	0.5(${}_{-0.1}^{+0.2},{\,}_{-0.2}^{+0.4}$)	0.4	0.2(${}_{-0.07}^{+0.1},{\,}_{-0.1}^{+0.2}$)	0.2	0.1(${}_{-0.03}^{+0.05},{\,}_{-0.07}^{+0.1}$)
Gould Belt	2.7	2.7(${}_{-0.8}^{+1.0},{\,}_{-1.3}^{+2.0}$)	1.8	1.4(${}_{-0.4}^{+0.6},{\,}_{-0.6}^{+1.2}$)	1.0	0.8(${}_{-0.3}^{+0.3},{\,}_{-0.4}^{+1.0}$)	0.7	0.5(${}_{-0.1}^{+0.2},{\,}_{-0.2}^{+0.6}$)
Radcliff	2.7	2.7(${}_{-0.7}^{+1.2},{\,}_{-1.2}^{+2.0}$)	1.5	1.3(${}_{-0.4}^{+0.4},{\,}_{-0.6}^{+1.0}$)	0.9	0.7(${}_{-0.2}^{+0.3},{\,}_{-0.3}^{+0.7}$)	0.9	0.5(${}_{-0.2}^{+0.2},{\,}_{-0.3}^{+0.5}$)
Taurus	3.6	3.5(${}_{-1.0}^{+1.3},{\,}_{-1.7}^{+3.2}$)	2.2	1.8(${}_{-0.5}^{+0.6},{\,}_{-0.8}^{+1.6}$)	1.4	1.0(${}_{-0.3}^{+0.4},{\,}_{-0.5}^{+0.7}$)	1.2	0.7(${}_{-0.2}^{+0.3},{\,}_{-0.3}^{+0.6}$)
Orion	3.7	5.0(${}_{-1.5}^{+1.8},{\,}_{-2.3}^{+3.9}$)	2.0	2.1(${}_{-0.6}^{+0.8},{\,}_{-0.9}^{+1.7}$)	0.9	1.0(${}_{-0.3}^{+0.4},{\,}_{-0.5}^{+0.8}$)	0.6	0.7(${}_{-0.2}^{+0.3},{\,}_{-0.3}^{+0.5}$)
Perseus	11.3	8.3(${}_{-2.4}^{+2.9},{\,}_{-3.9}^{+6.6}$)	6.0	4.0(${}_{-1.0}^{+1.5},{\,}_{-1.8}^{+3.3}$)	3.4	2.2(${}_{-0.6}^{+0.7},{\,}_{-1.0}^{+1.8}$)	2.4	1.6(${}_{-0.4}^{+0.5},{\,}_{-0.7}^{+1.2}$)
Ophiuchi [×1000]	1.6	2.5(${}_{-0.8}^{+0.8},{\,}_{-1.3}^{+2.5}$)	0.5	0.6(${}_{-0.2}^{+0.2},{\,}_{-0.3}^{+0.5}$)	0.2	0.2(${}_{-0.06}^{+0.08},{\,}_{-0.1}^{+0.2}$)	0.08	0.1(${}_{-0.02}^{+0.03},{\,}_{-0.03}^{+0.08}$)
Monoceros	28.6	24.2(${}_{-6.5}^{+8.3},{\,}_{-11.6}^{+26.8}$)	7.9	8.4(${}_{-2.2}^{+3.2},{\,}_{-3.6}^{+6.0}$)	2.4	3.9(${}_{-1.1}^{+1.8},{\,}_{-1.7}^{+3.5}$)	2.0	2.6(${}_{-0.7}^{+1.1},{\,}_{-1.2}^{+2.8}$)
Aquila	2.8	3.7(${}_{-1.1}^{+1.4},{\,}_{-1.7}^{+3.1}$)	1.2	1.8(${}_{-0.5}^{+0.6},{\,}_{-0.8}^{+1.3}$)	0.6	0.9(${}_{-0.2}^{+0.3},{\,}_{-0.4}^{+0.7}$)	0.4	0.6(${}_{-0.2}^{+0.2},{\,}_{-0.3}^{+0.5}$)
Hercules	11.8	12.0(${}_{-3.2}^{+4.5},{\,}_{-5.3}^{+12.8}$)	5.6	6.2(${}_{-1.8}^{+2.4},{\,}_{-2.9}^{+5.0}$)	3.7	3.2(${}_{-0.9}^{+1.4},{\,}_{-1.4}^{+3.0}$)	3.2	2.5(${}_{-0.7}^{+1.0},{\,}_{-1.1}^{+2.2}$)

Note. Units are in ×10⁻³ eV cm⁻³. Gould Belt includes Taurus, Orion, Perseus, and Ophiuchi; Radcliff Wave includes Taurus, Orion, and Perseus. Notice that the energy density values of Ophiuchi are 1000 times higher.

Download table as:  ASCIITypeset image

From these limits we do not have evidence that there is any deviation from the paradigm of the sea of cosmic rays permeating our Galaxy. It is also important to remark again that the gamma-ray expectation is based on an extrapolation of cosmic-ray data, so these estimated limits give a constraint on the cosmic-ray energy density at energies of ${E}_{\mathrm{CR}}\sim 10\,{E}_{\gamma }$ (see, e.g., Kelner et al. 2006), assuming purely pion decay photons.

With the purpose to test the cosmic-ray sea paradigm and probe whether its density is independent of the location, we performed measurements of gamma-ray flux from high-latitude passive GMCs using data from the HAWC Observatory. Since no significant excess was observed, we calculated upper limits at the 95% credible interval, for individual and stacked gamma-ray emission of clouds that are part of the Gould Belt as described in Section 2. The gamma-ray flux expected from pure hadronic interactions of the cosmic-ray flux with passive molecular clouds is below 10⁻¹¹ TeV⁻¹ cm⁻² s⁻¹ above 10 TeV. The quasi-differential limits are less than a factor of 10 higher than the prediction. The most stringent constraint that we observe is in the 2.7-model limit when stacking all clouds. However, with the current data, we do not find enough evidence to reject the paradigm of the sea of cosmic rays.

Using the same analysis method for the differential gamma-ray limits, we also estimated the cosmic-ray energy density of the clouds. Using the data from AMS and extrapolating the results to HAWC energies, we see that the HAWC limits are higher by a factor of less than 10 for the case of the stacked clouds, as well as Taurus and Aquila.

With current detector settings, as described in Section 3.1, doubling the exposure time would lower the constraints by a factor of 0.7, which will allow us to make a definitive statement about the cosmic-ray flux in distant regions of the Galaxy, at least in the case of the stacked analysis. The sensitivity of this analysis will also improve with the help of new installed outrigger detectors at the HAWC site, as well as the current development of new algorithms for event reconstruction, gamma–hadron separation, and background estimation.

We acknowledge the support from the US National Science Foundation (NSF); the US Department of Energy Office of High-Energy Physics; the Laboratory Directed Research and Development (LDRD) program of Los Alamos National Laboratory; Consejo Nacional de Ciencia y Tecnología (CONACyT), México, grants 271051, 232656, 260378, 179588, 254964, 258865, 243290, 132197, A1-S-46288, A1-S-22784, cátedras 873, 1563, 341, 323, Red HAWC, México; DGAPA-UNAM grants IG101320, IN106521, IN111716-3, IN111419, IA102019, IN110621; VIEP-BUAP; PIFI 2012, 2013, PROFOCIE 2014, 2015; the University of Wisconsin Alumni Research Foundation; the Institute of Geophysics, Planetary Physics, and Signatures at Los Alamos National Laboratory; Polish Science Centre grant, DEC-2017/27/B/ST9/02272; Coordinación de la Investigación Científica de la Universidad Michoacana; Royal Society—Newton Advanced Fellowship 180385; Generalitat Valenciana, grant CIDEGENT/2018/034; and Chulalongkorn University's CUniverse (CUAASC) grant. Thanks to Scott Delay, Luciano Díaz, and Eduardo Murrieta for technical support.

Although some of the mass values of the clouds can be found in the literature, we calculate them ourselves owing to the fact that we are building the regions of interest for the analysis, focusing on the high-density areas of the clouds where we assume most of the gamma rays originate. As mentioned in the main text, we use data from the Planck survey. We start by calculating the column density:

$$\begin{eqnarray}&&{N}_{{\rm{H}}}={\tau }_{D}/{\left(\frac{{\tau }_{D}}{{N}_{{\rm{H}}}}\right)}^{\mathrm{ref}},\end{eqnarray} \tag{ A1 }$$

where the reference value used is ${({\tau }_{D}/{N}_{{\rm{H}}})}^{\mathrm{ref}}\,=(1.18\pm 0.17)\times {10}^{-26}\,{\mathrm{cm}}^{2}$ for 353 GHz Ade et al. (2011). ${\tau }_{D}$ is the measured dust opacity. The mass of the cloud is then calculated as

$$\begin{eqnarray}&&{M}_{\mathrm{dust}}={N}_{{\rm{H}}}{\rm{\Omega }}{d}^{2}{m}_{{\rm{H}}},\end{eqnarray} \tag{ A2 }$$

where Ω is the angular area of the cloud and d is the distance to the cloud.

We apply a cut on the dust opacity value to select the high-density regions of 5 × 10⁻⁵ (for Hercules, in which the density is lower, we use a cut of $2.5\times {10}^{-5}$) as done in Yang et al. (2014).

We generated our molecular cloud templates using the Planck survey. We apply the same methodology as in Yang et al. (2014). First, we calculate the column density from the dust optical depth map at 353 GHz from Equation (A1). We apply the corresponding cut on the opacity value as mentioned in Appendix A, and then we normalize the map to the integral of the column density divided by the size of the spatial bin. This procedure ensures that the correct units are derived when the 3ML spatial model is used. Figure 4 shows the results after applying this procedure.

Zoom In Zoom Out Reset image size

Table 5 shows the TS results, together with the corresponding flux measurements of Table 2. As can be seen, none of the results are significant (i.e., TS > 25). Similarly, Table 6 shows the TS results and cosmic-ray energy density measurements of Table 4. Figure 5 shows the results of the cosmic-ray energy density of the rest of the clouds that were not presented in the main text. Finally, Figure 6 shows the HAWC significance maps of the GMC regions, as well as the distribution of the significance of the regions. We note that, except for two clouds, the distributions of the significance behave as expected: a Gaussian distribution with mean of 0 and width of 1. The Perseus region has a shifted mean of 0.27, while the Ophiuchi region has a width of 0.71. This could be related owing to systematic effects in the background estimation at this particular decl. or statistical fluctuations. To check this, we rotated the two affected regions of interest in R.A. and confirmed that the significance distributions behaved as expected. The systematic uncertainties due to the deviation in the width and mean from the Gaussian expectation are not affecting the expected limits reported in Section 4, for which we use Poisson-fluctuated background-only maps.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 5. Gamma-Ray Flux Measurements and 95% C.I. Upper Limits (in Parentheses)

Energy Range	1–3.16 TeV		3.16–10.0 TeV		10.0–31.6 TeV		>31.6 TeV		>1 TeV
Pivot Energy	1.78 TeV		5.62 TeV		17.8 TeV		56.2 TeV		10 TeV
	TS	Flux $[\times {10}^{-12}]$	TS	Flux [$\times {10}^{-13}]$	TS	Flux [$\times {10}^{-15}]$	TS	Flux [$\times {10}^{-16}]$	TS	Flux [×10−15]
Gould Belt	0.0	${1.3}_{-0.9}^{+1.4}$(3.6)	0.8	${1.0}_{-0.6}^{+0.8}$(2.3)	1.4	${6.5}_{-3.8}^{+4.2}$(13.5)	1.5	${4.2}_{-2.6}^{+2.9}$(9.4)	0.79	${7.5}_{-4.9}^{+5.9}$(18.0)
Radcliff	0.0	${1.0}_{-0.7}^{+1.1}$(2.8)	0.8	${0.8}_{-0.5}^{+0.6}$(1.8)	1.4	${4.6}_{-2.8}^{+3.3}$(10.3)	1.5	${3.4}_{-1.9}^{+2.2}$(7.4)	0.79	${5.9}_{-3.8}^{+4.5}$(13.8)
Taurus	0.0	${0.6}_{-0.4}^{+0.6}$(1.5)	0.5	${0.4}_{-0.3}^{+0.3}$(1.0)	1.5	${2.9}_{-1.7}^{+2.0}$(6.4)	2.0	${2.5}_{-1.4}^{+1.5}$(5.0)	0.8	${3.4}_{-2.2}^{+2.6}$(7.8)
Orion	0.0	${0.6}_{-0.5}^{+0.8}$(2.0)	0.0	${0.1}_{-0.3}^{+0.4}$(1.0)	0.0	${1.5}_{-1.0}^{+1.9}$(4.8)	0.0	${1.0}_{-7.4}^{+1.2}$(3.0)	0.0	${2.3}_{-1.6}^{+2.6}$(7.1)
Perseus	1.3	${0.6}_{-0.4}^{+0.5}$(1.4)	1.9	${0.4}_{-0.2}^{+0.2}$(0.74)	2.3	${2.0}_{-1.0}^{+1.2}$(4.0)	2.3	${1.5}_{-0.8}^{+0.9}$(2.9)	2.8	${3.0}_{-1.6}^{+1.7}$(5.8)
Ophiuchi [x100]	0.0	${1.4}_{-1.0}^{+1.8}$(4.7)	0.0	${0.4}_{-0.3}^{+0.6}$(1.4)	0.0	${1.6}_{-1.1}^{+2.0}$(5.3)	0.0	${0.8}_{-0.5}^{+1.0}$(2.5)	0.0	${2.7}_{-2.0}^{+3.5}$(8.6)
Monoceros	0.5	${0.8}_{-0.5}^{+0.7}$(1.9)	0.0	${0.2}_{-0.1}^{+0.2}$(0.52)	0.0	${0.4}_{-0.3}^{+0.7}$(1.7)	0.0	${0.8}_{-0.7}^{+3.0}$(1.2)	0.0	${1.0}_{-0.7}^{+1.3}$(3.4)
Aquila	0.0	${1.0}_{-0.7}^{+1.1}$(3.2)	0.0	${0.4}_{-0.3}^{+0.5}$(1.3)	0.0	${1.9}_{-1.4}^{+2.7}$(6.9)	0.0	${1.2}_{-0.8}^{+1.6}$(4.2)	0.0	${2.6}_{-1.9}^{+3.7}$(9.8)
Hercules	0.0	${0.4}_{-0.2}^{+0.4}$(1.0)	0.0	${0.2}_{-0.1}^{+0.2}$(0.5)	0.0	${1.1}_{-0.8}^{+1.1}$(3.2)	0.44	${1.2}_{-0.8}^{+0.9}$(2.9)	0.0	${1.5}_{-1.0}^{+1.6}$(4.2)

Note. Flux units are in TeV⁻¹ cm⁻² s⁻¹. Gould Belt includes Taurus, Orion, Perseus, and Ophiuchi; Radcliff Wave includes Taurus, Orion, and Perseus.

Download table as:  ASCIITypeset image

Table 6. Cosmic-Ray Energy Density Quasi-differential 95% C.I. Upper Limits $[\times {10}^{-3}$ eV cm⁻³] in Each Energy Bin

Energy Range	1–3.16 TeV		3.16–10.0 TeV		10.0–31.6 TeV		>31.6 TeV
Pivot Energy	1.78 TeV		5.62 TeV		17.8 TeV		56.2 TeV
	TS		TS		TS		TS
All	0.0	${0.7}_{-0.5}^{+0.6}$(1.9)	0.0	0.4 ± 0.3(0.8)	0.0	0.2 ± 0.2(0.4)	0.0	0.1 ± 0.1(0.2)
Gould Belt	0.0	${1.1}_{-0.8}^{+0.9}$(2.7)	0.9	0.8 ± 0.5(1.8)	1.4	0.5 ± 0.3(1.0)	1.4	0.3 ± 0.2(0.7)
Radcliff	0.0	1.1 ± 0.8(2.7)	0.9	0.8 ± 0.5(1.5)	1.5	0.5 ± 0.3(0.9)	1.4	0.3 ± 0.2(0.9)
Taurus	0.0	1.5 ± 1.1(3.6)	0.5	0.9 ± 0.7(2.2)	1.5	${0.7}_{-0.3}^{+0.4}$(1.4)	2.0	${0.6}_{-0.3}^{+0.4}$(1.2)
Orion	0.0	1.4 ± 1.1(3.7)	0.0	0.8 ± 0.6(2.0)	0.0	0.4 ± 0.7(0.9)	0.0	0.2 ± 0.2(0.6)
Perseus	1.3	${5.6}_{-3.5}^{+3.4}$(11.3)	1.9	${3.1}_{-1.8}^{+1.7}$(6.0)	2.3	1.7 ± 0.9(3.4)	2.3	1.2 ± 0.7(2.4)
Ophiuchi [×1000]	0.0	0.6 ± 0.5(1.6)	0.0	0.2 ± 0.1(0.5)	0.0	${0.07}_{-0.05}^{+0.06}$(0.2)	0.0	0.03 ± 0.02(0.08)
Monoceros	0.5	13.0 ± 9.0(28.6)	0.0	3.1 ± 2.4(7.9)	0.0	0.9 ± 0.7(2.4)	0.0	0.7 ± 0.6(2.0)
Aquila	0.0	1.1 ± 0.8(2.8)	0.0	${0.4}_{0.3}^{+0.4}$(1.2)	0.0	0.2 ± 0.2(0.6)	0.0	0.1 ± 0.1(0.4)
Hercules	0.0	5.0 ± 4.0(11.8)	0.0	${2.2}_{-1.6}^{+1.7}$(5.6)	0.1	1.6 ± 1.1(3.7)	0.4	0.1 ± 0.1(3.2)

Download table as:  ASCIITypeset image