Search for Very High Energy Gamma Rays from the Northern $\textit{Fermi}$ Bubble Region with HAWC

TL;DR

This study used HAWC data to search for very high energy gamma rays from the Northern Fermi Bubble but found no significant emission, setting upper limits that constrain certain cosmic-ray models.

Contribution

First search for very high energy gamma rays from the Northern Fermi Bubble using HAWC data, providing upper limits that challenge some hadronic emission models.

Findings

No significant gamma-ray excess detected.

Upper limits constrain proton spectra beyond 100 TeV.

Results disfavor hadronic models based on neutrino data.

Abstract

We present a search of very high energy gamma-ray emission from the Northern $\textit{Fermi}$ Bubble region using data collected with the High Altitude Water Cherenkov (HAWC) gamma-ray observatory. The size of the data set is 290 days. No significant excess is observed in the Northern $\textit{Fermi}$ Bubble region, hence upper limits above $1\,\text{TeV}$ are calculated. The upper limits are between $3\times 10^{-7}\,\text{GeV}\, \text{cm}^{-2}\, \text{s}^{-1}\,\text{sr}^{-1}$ and $4\times 10^{-8}\,\text{GeV}\,\text{cm}^{-2}\,\text{s}^{-1}\,\text{sr}^{-1}$. The upper limits disfavor a proton injection spectrum that extends beyond $100\,\text{TeV}$ without being suppressed. They also disfavor a hadronic injection spectrum derived from neutrino measurements.

Tables

Table 1. Table 1: Analysis Bins defined as the fraction of operational PMT channels triggered in an air shower event.

Analysis Bin	$f$
$f_1$	0.162 - 0.247
$f_2$	0.247 - 0.356
$f_3$	0.356 - 0.485
$f_4$	0.485 - 0.618
$f_5$	0.618 - 0.740
$f_6$	0.740 - 0.840
$f_7$	0.840 - 1.00

Table 2. Table 2: Charactersitics of the non-detection: upper limits on the differential flux in four different energy bins. Since the highest energy bin is treated as an overflow bin, only the lower boundary of that energy bin is quoted in order to be conservative.

Energy Range	Upper Limits	Sensitivity
[TeV]	[${\text{GeV}\text{cm}}^{-2}{\text{s}}^{-1}{\text{sr}}^{-1}$]	[${\text{GeV}\text{cm}}^{-2}{\text{s}}^{-1}{\text{sr}}^{-1}$]
1.2 - 3.9	3.0$\times {10}^{-7}$	3.3$\times {10}^{-7}$
3.9 - 12.4	1.0$\times {10}^{-7}$	1.1$\times {10}^{-7}$
12.4 - 39.1	0.5$\times {10}^{-7}$	0.5$\times {10}^{-7}$
$>$39.1	0.4$\times {10}^{-7}$	0.3$\times {10}^{-7}$

Table 3. Table 3: Differential flux models for the Fermi Bubbles.

Model	Description
Hadronic Model 1	$N_p\propto p^{-2.2}$
Hadronic Model 2	$N_p\propto p^{-2.1}\exp \left(-pc/14\text{TeV}\right)$
Leptonic Model 1	$N_e\propto p^{-2.17}\exp \left(-pc/1.25\text{TeV}\right)$ and IRF at 5kpc
Leptonic Model 2	$N_e\propto p^{-2.17}\exp \left(-pc/1.25\text{TeV}\right)$ and CMB
IceCube Hadronic Model	$N_p\propto p^{-2.25}\exp \left(-pc/30\text{PeV}\right)$

Table 4. Table 4: Description of the variables used in the analysis

Variable	Description
$N_i$	Number of events in sky-map pixel $i$ before gamma-hadron cuts
$N_i^{\prime }$	Number of events in sky-map pixel $i$ after gamma-hadron cuts
$⟨N_i⟩$	Number of estimated background events in sky-map pixel $i$ before gamma-hadron cuts
$⟨N_i^{\prime }⟩$	Number of estimated background events in sky-map pixel $i$ after gamma-hadron cuts
$E_i$	Excess above background in sky-map pixel $i$ before gamma-hadron cuts
$E_i^{\prime }$	Excess above background in sky-map pixel $i$ after gamma-hadron cuts

Equations

$$E_i^{\prime }=N_i^{\prime }-⟨N_i^{\prime }⟩,$$

$$N_i=G_i^T+C_i^T$$

$$G_i^T$$

$$C_i^T$$

$$N_i^{\prime }={\epsilon }_{G,i}{G}_i^T+{\epsilon }_{C,i}{C}_i^T,$$

$$⟨N_i⟩=G_i^I+C_i^I$$

$$⟨N_i^{\prime }⟩={\epsilon }_{G,i}{G}_i^I+{\epsilon }_{C,i}{C}_i^I.$$

$${\epsilon }_{G,i}=\frac{\int h^{\prime }(e)de}{\int h(e)de}.$$

$${\epsilon }_{C,i}=⟨N_i^{\prime }⟩/⟨N_i⟩.$$

$$E_i^{\prime }={\epsilon }_{C,i}{C}_i+{\epsilon }_{G,i}{G}_i.$$

$$E_i=C_i+G_i.$$

$$G_i=\frac{E_i^{\prime }-{\epsilon }_{C,i}{E}_i}{{\epsilon }_{G,i}-{\epsilon }_{C,i}}.$$

$$G_i^{\prime }={\epsilon }_{G,i}{G}_i$$

$$w_l^k=\frac{M{E}_l^k}{⟨N^{\prime }{⟩}_l},$$

$$R_k=\frac{\sum _{l=1}^f{w}_l^k{G}_l^{\prime }}{\sum _{l=0}^f{w}_l^{k}M{E}_l^k},$$

$$\delta R_k=\frac{\sqrt{\sum _{l=1}^f(w_l^k\delta G_l^{\prime }{)}^2}}{\sum _{l=1}^f{w}_l^{k}M{E}_l^k}.$$

$$F_k=(R\pm \delta R{)}_{k}F(E_k),$$

$$S_k=E_k/\sqrt{{⟨N^{\prime }⟩}_k},$$

$$G_i^{\prime }={\epsilon }_{G,i}{G}_i,$$

$$G_i=\frac{E_i^{\prime }-{\epsilon }_{C,i}{E}_i}{{\epsilon }_{G,i}-{\epsilon }_{C,i}}.$$

$${\epsilon }_{C,i}=\sum _j⟨N_j^{\prime }⟩/\sum _j⟨N_j⟩,$$

$${\left(\frac{\delta G}{|G|}\right)}^2=\frac{1}{(E_i^{\prime }-{\epsilon }_{C,i}{E}_i{)}^2}\left[\delta E_i^{{}^{\prime }2}+({\epsilon }_{C,i}{E}_i{)}^2[(\frac{\delta {\epsilon }_{C,i}}{{\epsilon }_{C,i}}{)}^2+(\frac{\delta E_i}{E_i}{)}^2]\right]+\frac{\delta {\epsilon }_{C,i}^2}{({\epsilon }_{G,i}-{\epsilon }_{C,i}{)}^2}.$$

$$\frac{\delta {\epsilon }_{C,i}}{{\epsilon }_{C,i}^t}=\frac{|{\epsilon }_{C,i}-{\epsilon }_{C,i}^t|}{{\epsilon }_{C,i}^t}=\frac{{\sum }_j({\epsilon }_{G,i}-{\epsilon }_{C,i}^t)G_j^I}{{\epsilon }_{C,i}^t{\sum }_j(C_j^I+G_j^I)}\approx \frac{{\sum }_j({\epsilon }_{G,i})G_j^I}{{\epsilon }_{C,i}^t{\sum }_j(C_j^I)}<\frac{{\sum }_j({\epsilon }_{G,i})G_j}{{\epsilon }_{C,i}^t{\sum }_j(C_j^I)},$$

$$\frac{N^{\prime }-⟨N^{\prime }⟩}{⟨N^{\prime }⟩}=\frac{{\epsilon }_{G}G+{\epsilon }_C{C}^I}{{\epsilon }_C{C}^I}.$$

$$\frac{N^{\prime }-⟨N^{\prime }⟩}{⟨N^{\prime }⟩}=\frac{{\epsilon }_{G}G}{{\epsilon }_C{C}_I}=0.6\times 10^{-3}.$$

$$(\delta G_i{)}_{sys}=\frac{\partial G_i}{\partial {\epsilon }_{C,i}}\delta {\epsilon }_{C,i}$$

Full text

Astroparticle physics — gamma rays — Fermi Bubbles — Gamma-Ray Astronomy

\AuthorCollaborationLimit

=500

1 Introduction

The search for a counterpart of the microwave haze  (Dobler et al., 2010) in gamma-ray data, using the Fermi Large Area Telescope (LAT), revealed the existence of two large structures extending up to $55\degree$ above and below the Galactic Plane (Dobler et al., 2010; Su et al., 2010). Due to their bubble-like shape they received the name of Fermi Bubbles.

The gamma-ray emission of the Fermi Bubbles presents a hard spectrum —$dN/dE\sim E^{-2}$— in the energy range from approximately 1 GeV to 100 GeV. The surface brightness is roughly uniform in both bubbles —except for a structure inside the South Bubble called the cocoon— and the total luminosity of the bubbles for galactic longitude $|b|>10\degree$ and between 100 MeV and 500 GeV was found to be $4.4^{+2.4}_{-0.9}\times 10^{37}$ erg  s*-1*  (Ackermann et al., 2014).

The origin of the Fermi Bubbles is still uncertain. Different models have been proposed to explain their formation. Most of the models revolve around the idea of outflows from the galactic center which then interact with the interstellar medium. The outflow can be generated by activity of the nucleus in our galaxy producing a jet (Guo & Mathews, 2012; Guo et al., 2012), wind from long time-scale star formation (Crocker & Aharonian, 2011), periodic star capture processes by the supermassive black hole in the Galactic Center (Cheng et al., 2011), or by winds produced by the hot accretion flow in Sgr A*∗*  (Mou et al., 2015).

The production of gamma rays is also under dispute. Hadronic and leptonic models are the main mechanisms to explain the gamma-ray production. Photons of hadronic origin are due to the decay of neutral pions that are produced in the interaction of protons with nuclei in the interstellar medium (ISM). These protons are injected in the bubble regions by the outflow processes mentioned before or they can be accelerated inside the bubble as proposed by  Fujita et al. (2013, 2014). Some of these models,  (Crocker & Aharonian, 2011; Fujita et al., 2013), predict the possibility of high-energy gamma rays. In the leptonic model, high-energy photons are produced by inverse Compton scattering from the interaction of energetic electrons with photons from the interstellar radiation fields (IRF) or cosmic microwave background (CMB). The division between hadronic and leptonic models should not be strict, but rather, a combination of both models can be possible as shown in  Cheng et al. (2011); Ackermann et al. (2014). Observations at other wavelengths, specifically at lower energies, have helped to constrain some models. For instance, the microwave haze —produced by synchrotron radiation— can help to constrain the electron population, which can also radiate in gamma rays  (Dobler et al., 2010; Su et al., 2010; Crocker & Aharonian, 2011; Mou et al., 2015; Guo et al., 2012).

The same principle can apply at very high energies (VHE; $>$100 GeV), where observations can constrain the population of the highest-energy cosmic rays. Considering that the Northern Fermi Bubble region is in the field of view of the High Altitude Water Cherenkov (HAWC), a search for VHE gamma rays (above 1 TeV) is presented. The paper is divided as follows: The HAWC observatory and the data set used in the analysis are defined in Section 2, the analysis procedure is described in Section 3, and the results are discussed in Section 4 and summarized in Section 5.

2 The HAWC Observatory & the Data Set

HAWC is a ground-based gamma-ray observatory. It is located between Volcán Sierra Negra and Pico de Orizaba near Puebla, Mexico, at an altitude of 4100 m a.s.l. and at $(18\degree 59^{\prime}41"$N, $97\degree 18^{\prime}30"$W). The observatory has a duty cycle of $>$95% and a large field of view of $\sim$2 sr, which allows it to cover $8.4$ sr in a day (Abeysekara et al., 2013). The instrument consists of an array of 300 water Cherenkov detectors (WCDs). The WCDs are steel tanks of 7.3 m in diameter and 5 m in height, filled with water up to 4.5 m. Each WCD is filled with $\sim$ 200,000  L of purified water. The array provides an effective area of $\sim$22,000 m2. Inside the WCDs, four photomultiplier tubes (PMTs) facing upward are attached to the bottom.

A simple multiplicity trigger is applied to find extensive air showers in the data. For the present analysis, the trigger requires 28 PMTs detecting Cherenkov light within a 150 ns time window to be activated in order to identify a shower event. After the processing and calibration of the events, the air shower cores, footprint brightness in the array, and gamma- and cosmic-ray directions are reconstructed. More information on the detector, calibration, and reconstruction is presented in  Abeysekara et al. (2017b).

The HAWC observatory began science operations in August 2013, when it was still under construction. The analysis described in this paper uses data between 2014 November 27th to 2016 February 11th.

The data set is divided into seven event-size bins represented by the fraction f of functioning PMT channels triggered in an air shower event. The energy of the observed gamma rays is related to the shower event size that is measured in the HAWC array. The range of f for this analysis goes from 0.162 to 1.00. Table 1 shows the ranges for each analysis bin.

Standard selection cuts are applied to the data that pass the trigger condition. The signals in each PMT are required to have $>$1 photoelectrons (PEs) and are required to be between 150 ns before and 400 ns after the trigger. In addition, it is required that more than 90% of the PMT channels are functioning during the observation. Finally, cuts are applied to distinguish between gamma rays and hadronic cosmic rays, the latter being the main background of measurements with the HAWC observatory. All the cuts are optimized by studying the Crab Nebula in the HAWC data  (Abeysekara et al., 2017b).

3 Analysis

The analysis is focused on measuring the flux of gamma rays in the Northern Fermi Bubble Region, since this region is inside the field of view of HAWC. The main challenge of the analysis is to estimate the background of the data set. The procedure to deal with the background is divided in three steps:

Distinguish the air shower signatures between cosmic rays and gamma rays. The gamma-hadron cuts select the gamma-like showers. 2. 2.

Find the isotropic flux of cosmic rays and gamma rays. This is found by using direct integration and it is explained in Section 3.1. 3. 3.

Remove effects of the large-scale anisotropy as seen in (Abeysekara et al., 2014). The procedure is shown in Section 3.2.

3.1 Direct Integration Isotropic Background Estimation

The positions of the events are binned in equatorial coordinates using the HEALPix scheme (Gorski et al., 2005). These are referred to as sky maps. For the analysis we set the pixel size to be $\sim 0.11\degree$.

The isotropic background is estimated using the direct integration (DI) technique described in Atkins et al. (2003). The background is integrated over 24 hours and therefore only data were used when the detector performance was stable for 24 hours, since this is a requirement for the integration technique. This results in a lifetime for the analysis of 290 days.

As explained in  Abeysekara et al. (2014), an analysis based on a background integration period of $\delta t$ is sensitive to potential signal excesses of an RA size smaller than $\delta t\,\cdot 15\degree$ hour*-1*. Using a 24 hour integration period ensures that the analysis is sensitive to the Fermi Bubbles which extend to $\sim 50\degree$ in RA.

Since the estimation of the isotropic background can be biased by strong known sources in the data, a region of interest (ROI) masking is used, as shown in Figure 1. The ROI masking covers the galactic plane [$\pm 6\degree$], as well as circular regions of radius $3\degree$, $1.3\degree$,$1\degree$ and $1\degree$, respectively, for Geminga, the Crab Nebula, Mrk 421 and 501. Region A and B from the small-scale cosmic-ray anisotropy are also masked. Their shapes are obtained from the results in Abeysekara et al. (2014), by requiring that the significances in the sky map without gamma-hadron cuts are greater than 4$\sigma$. Finally, the ROI for the Northern Fermi Bubble was obtained from the Fermi Diffuse Model pass 7 version 6 111See http://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html.

The shape of the Northern Fermi Bubble above 1 TeV is unknown. We perform a gamma-ray flux excess search within the boundaries of the Northern Bubble as detected by Fermi below TeV energies.

3.2 Gamma-Ray Excess Calculation

For each analysis bin, sky maps are created after applying the gamma-hadron cuts. The isotropic background is then estimated and an excess sky map is obtained. The excess in each pixel $i$ is given by the following equation,

[TABLE]

where $N^{\prime}_{i}$ is the observed data after gamma-hadron cuts, and $\langle N^{\prime}_{i}\rangle$ is the isotropic background estimated after gamma-hadron cuts. However, the excess sky map in the lower analysis bins reveals the large-scale anisotropy as seen in  (Abeysekara et al., 2014). This is because there is enough statistics to calculate the background with an accuracy of one part per mille in these bins. The Northern Fermi Bubble is located at a deficit region, therefore this systematic effect needs to be removed. The subtraction of this cosmic-ray feature is achieved by using the data without gamma-hadron cuts.

The data without gamma-hadron separation is composed of a total number of gamma rays and cosmic rays,

[TABLE]

where $N_{i}$ is the data without gamma-hadron cuts in the pixel $i$, $G^{T}_{i}$ is the number of gamma rays in the pixel $i$, and $C^{T}_{i}$ is the number of cosmic rays in the pixel $i$. $G^{T}_{i}$ and $C^{T}_{i}$ can also be decomposed in terms of an isotropic component and an excess(or deficit). This is expressed as

[TABLE]

where $G^{I}_{i}$, $C^{I}_{i}$ are the isotropic components; and $G_{i}$, $C_{i}$ are the gamma-ray and cosmic-ray excesses or deficits.

The data after gamma-hadron separation also contains gamma rays and cosmic rays but the composition is different due to the rejection efficiency of the gamma-hadron separation cuts,

[TABLE]

where $\varepsilon_{G,i}$ and $\varepsilon_{C,i}$ are the gamma and hadron efficiencies after applying the gamma-hadron cuts.

For completeness the isotropic background for the data before and after gamma-hadron cuts are written as follows,

[TABLE]

The gamma passing rate efficiency $\varepsilon_{G,i}$ is obtained using simulations. The detector response is simulated in each of the seven analysis bins and for 5$\degree$ declination bands between $-37.5\degree$ and $77.5\degree$. Each bin contains an energy histogram that is expected for the simulated signal. We compute the number of events in the energy histograms and the ratio for the events with gamma-hadron cuts $h^{\prime}(e)$ over the events with no gamma-hadron cuts $h(e)$, where $e$ is the energy. Therefore, the efficiency can be written as

[TABLE]

The hadron passing rate efficiency $\varepsilon_{C,i}$ is estimated from the data since the total number of cosmic rays outnumbers the total number of gamma rays. In order to avoid bright sources, we use equations 5 and 6 to estimate $\varepsilon_{C,i}$,

[TABLE]

The isotropic gamma-ray emission can be safely neglected in the sums for equations 5 and 6. To account for declination dependence, the data in the same declination (or HEALPix ring) as in pixel $i$ is used so that $\varepsilon_{C,i}=\sum_{j}\langle N^{\prime}_{j}\rangle/\sum_{j}\langle N_{j}\rangle$, where $j$ is the pixel in the ring containing $i$.

By combining Equations 4 and 6, Equation 1 can be re-written as

[TABLE]

Using Equations 2 and 5 give an equation for the excess in pixel $i$ for the data without gamma-hadron cuts,

[TABLE]

Finally, the number of gamma rays is obtained from Equations 9 and 10,

[TABLE]

The efficiency $\varepsilon_{G,i}$ is applied to the number of gamma rays $G_{i}$ to obtain the number of excess events measured by the detector.

[TABLE]

The previous equation is used to calculate the number of gamma rays G in each pixel inside the Northern Bubble region as defined in Figure 1 and then summed to get a total excess $G^{\prime}=\sum_{i}G^{\prime}_{i}$ in each analysis bin.

As mentioned in Section 3.1, the shape of the Fermi Bubbles at high energies is unknown, though some authors suggest that the size of the bubbles increases with energy (Fujita et al., 2013; Yang et al., 2014; Mou et al., 2015). In this case, calculating the flux in the smaller region of the MeV-GeV excess is the more conservative approach.

The description of the variables is presented in appendix A.

The uncertainty calculation for $G_{i}$ is shown in appendix B.

3.3 Testing the Analysis Method

The analysis method is tested on simulated sky maps containing a dipole distribution as shown in Figure 2 assuming no sources are present.

A rate map in the local coordinates of HAWC containing a snapshot of 24 sec of data is generated. Since HAWC observations cover a local sky of zenith angles $0\degree<\theta<45\degree$, the rate map is generated for this zenith angle range. Using the dipole distribution given in Figure 2, the total sky event rate from HAWC data, and information from the detector response, a rate in each pixel is obtained. After the 24 sec period the rate map is reset and the procedure is started again. In this way a simulated data set is generated that is of the same size as the real data set analyzed in this paper for both cases of without and with gamma-hadron cuts. An example of a resulting simulated sky map $f_{1}$ is shown in Figure 3. The upper panel presents the resulting map after simply subtracting the estimated background from the data, the lower panel shows the excess map after applying the procedure described in Section 3.2.

Figure 4 shows the resulting simulated excesses in each $f$ bin. A comparison is made between the event excesses derived from simply subtracting the estimated DI background from the simulated data (blue points) and the event excess obtained by the method from Section 3.2 (red points). The effect of the simulated cosmic-ray anisotropy results in systematically lower excesses for the lowest two $f$ bins if the method from Section 3.2 is not applied.

The method is also tested by adding a strong Fermi Bubble-like gamma-ray emission. The spectrum is assumed to be a power-law with spectral index $\gamma=2.0$ and normalization of $5.03\times 10^{-7}$ GeV*-1* cm*-2* s*-1* sr*-1*, both values obtained by fitting the Fermi data points in the range of 1 GeV to 150 GeV.

Using this assumption and extending the spectrum to TeV energies, the analysis procedure was tested. If the Northern Fermi Bubbles had this spectrum, the HAWC observatory would have detected it.

4 Results and Discussion

4.1 Gamma-Ray Excess Results

Figure 6 shows skymaps of the result of the first analysis bin, f1. The figure shows a sky map without gamma-hadron cuts, and sky maps with gamma-hadron cuts before and after applying our procedure.

Figure 7 shows the results of the summed excess inside the bubble region after applying our procedure. No significant excess inside the bubble region is observed, therefore upper limits on the differential flux are calculated. The upper limits are also compared to the differential sensitivity of the HAWC observatory. The upper limits give the maximum flux intensity that is plausible given the observed counts in the HAWC data. The differential sensitivity quantifies the power of the detection procedure and is based on finding an $\alpha$-level threshold (related to background fluctuations claimed as detections) and the probability $\beta$ to detect a source222The definitions of upper limit and sensitivity are the same as upper bound and upper limit in Kashyap et al. (2010)..

4.1.1 Calculating the Upper Limits

The differential flux is calculated from the measured excess by comparing the signal observed in the data to an expected signal obtained for each of the fractional $f$ analysis bins using simulations. Since the energy response histograms for each analysis bin overlap (see Figure 8), the excesses measured in the analysis bins are combined in a weighted sum.

The procedure is as follows: a differential flux is assumed in an energy bin of width $\Delta\log(E/1\,\text{TeV})$. The width of differential energy bins is defined such that the results are independent of spectral assumptions. Using the HAWC detector response, an expected signal for the Northern Bubble region is obtained for each fractional $f$ bin. Taking into account the previous values, the weight in the energy bin $k$ for the fractional bin $f_{l}$ is calculated as

[TABLE]

where $w^{k}_{l}$ is the weight in the energy bin $k$ for the fractional $f_{l}$; $ME^{k}_{l}$ is the expected signal in the energy bin $k$ for the fractional $f_{l}$, and $\langle N^{\prime}\rangle_{l}$ is the background estimated in fractional bin $f_{l}$. This procedure results in a matrix that allows “projection” of the fractional $f$ analysis bin space onto the energy space.

Using the weights, the ratio of the observed signal and the expected signal is calculated

[TABLE]

and the uncertainty in the ratio as

[TABLE]

The ratio is used to obtain an estimation of the flux in the energy bin $k$

[TABLE]

where $F(E_{k})$ is the flux assumption at the energy bin $k$ used to obtain the expected excess $ME^{k}_{l}$.

The upper limit calculation is then performed in the energy bins. The prescription of  Helene (1983) is used to calculate an upper limit on the differential flux derived from equation 16. A 95% confidence level is chosen.

4.1.2 Calculating the Sensitivity

The sensitivity is calculated based on Kashyap et al. (2010)333Named upper limit in the reference instead of sensitivity. The procedure consists of setting a small probability for false positives (Type I error) and setting a probability of detection when there is a source (related to Type II error). The probability for false positives is set to $\alpha=0.05$ and the probability of detection is set to $\beta=0.5$. This is to compare the detection power of the observatory to the calculated upper limit at 95% confidence level.

The calculation is performed by using the measured background and doing simulations for a simulated Fermi Bubble of varying flux. For the simulation we assume a power-law with an index of -2.75 in the differential energy bin. For each analysis bin, the total background counts and the total expected number of events from the simulated source are calculated inside the bubble region. Following the same procedure as in Section 4.1.1, the analysis bins are combined to get the total number of events for each energy bin. In each energy bin, a null hypothesis histogram and an alternative hypothesis histogram are created for the quantity

[TABLE]

where $E_{k}$ is obtained by Poisson-fluctuating $\left\langle N^{\prime}\right\rangle_{k}$ for the null hypothesis, and then subtracting $\left\langle N^{\prime}\right\rangle_{k}$ from this value; or by Poisson-fluctuating $\left\langle N^{\prime}\right\rangle_{k}+ME_{k}$ for the alternative hypothesis, and then subtracting $\left\langle N^{\prime}\right\rangle_{k}$ from this value. The Poisson fluctuations are performed 10000 times to fill the histograms. The null hypothesis histogram is used to find the $\alpha$-level detection threshold and the alternative hypothesis histogram is used to find the flux normalization that is required to obtain a probability of detection of 0.5.

4.2 Differential Flux of the Fermi Bubbles

The first energy bin is centered at 2.2 TeV, which is the median energy of fraction f1 assuming a power law spectrum of index $\gamma=2.75$ (see Figure 8). The energy bin width is set to $\Delta\log(E/1\,\text{TeV})=0.5$ which is comparable to the width of the energy histograms. The energy range covers up to the highest energy at which HAWC is sensitive ($\sim$100 TeV).

Table 2 shows the values of the upper limits and sensitivities for each energy bin. The upper limits obtained from data are consistent with the detection power of HAWC. Figure 9 shows the upper limits together with the flux measurement of the Fermi Bubbles made by the Fermi Collaboration (Ackermann et al., 2014). Different leptonic and hadronic models are also present in Figure 9.

The two leptonic models are obtained from Ackermann et al. (2014). In these models, the emission is due to inverse Compton scattering. Two radiation fields are used: the IRF at 5 kpc above the Galactic Plane and photons from the CMB. The electron spectrum interacting with the radiation fields is modeled as a power-law with an exponential cutoff. The spectral index has a value of $2.17\pm 0.5^{+0.33}_{-0.89}$ and the cutoff energy is $1.25\pm 0.13^{+1.73}_{-0.68}$ TeV.

The two cyan hadronic models, also obtained from Ackermann et al. (2014), assume a power-law and a power-law with cutoff for the injection spectrum of the hadrons. These protons interact with the ISM producing neutral pions that decay into gamma rays. The spectrum was obtained using the library cparamlib444https://github.com/niklask/cparamlib, which implements the cross sections from Kamae et al. (2006), for the production of gamma rays through hadronic interactions. The spectral index for the power-law is 2.2; the spectral index for the power-law with cutoff is $2.13\pm 0.01^{+0.15}_{-0.52}$ with a cutoff energy of $14\pm 7^{+6}_{-13}$ TeV. Using the fit results obtained in  Ackermann et al. (2014), we extrapolate the results for the hadronic models above 100 TeV. The upper limits derived from HAWC data exclude the hadronic injection without a cutoff, that best fits the GeV gamma-ray data, above 3.9 TeV.

The hadronic model represented by the red line is obtained from  Lunardini et al. (2015). This model is the counterpart of a neutrino flux model that best fits the IceCube data. The IceCube data corresponds to five events that are spatially correlated with the Fermi Bubbles. The differential flux model was obtained by taking into account the flux from both bubbles. Above 10 TeV, the HAWC upper limits exclude the parent proton spectrum predicted from IceCube data.

Table 3 gives a summary of the different models.

Early reports such as  Crocker & Aharonian (2011); Fujita et al. (2013), presented the possibility of observing TeV gamma rays. The intensity was predicted to be $\leq E^{2}\,F(\text{TeV})\,\sim\,10^{-6}\,\,\text{GeV}\text{cm}^{-2}\text{s}^{-1}\text{sr}^{-1}$. The result presented here sets a stricter upper limit.

The result is not constraining the main contribution to the spectrum of the Fermi Bubbles. Nevertheless our result may imply, for a hadronic model, that there is a cutoff in the proton spectrum.  Ackermann et al. (2014) showed that GeV gamma-ray spectrum cuts off around 100 GeV. The cutoff for the parent proton spectrum in this case could be around 1 TeV  (Cheng et al., 2015).

As mentioned in Section 3.2, Fujita et al. (2013); Yang et al. (2014); Mou et al. (2015) propose that the size of the bubbles increases with energy. While defining the search region to be the same as the excess detected at GeV energies is a more conservative approach, it may be interesting to increase the size of the latter in a follow-up analysis.

Increasing the sensitivity at energies $<$1 TeV is another objective for future analysis. Compared to recent(Abeysekara et al., 2017b, a) or future (HAWC Collaboration, 2017, in preparation) publications of the analysis of HAWC data, this analysis uses only the seven highest event-size bins. At energies $\lesssim$1 TeV, the large-scale anisotropy signal (or any significant, spatially-extended feature) causes signal contamination in the estimation of the background because the structure takes up a large portion of the field-of-view of HAWC, significantly altering the all-sky rate. An iterative procedure for the DI method will be followed as explained in Ahlers et al. (2016) and has been shown to remove this artifact.

5 Conclusions

A search of high-energy gamma rays in the Northern Fermi Bubble region has been presented by using 290 days of data from the HAWC observatory. No significant excess is found above 1.2 TeV in the search area and the 95% C.L. flux upper limits are calculated and compared to the differential sensitivity with $\alpha=0.05$ and $\beta=0.5$. The upper limits are between 3$\times 10^{-7}$ $\text{GeV}\,\text{cm}^{-2}\,\text{s}^{-1}\,\text{sr}^{-1}$ and 4$\times 10^{-8}$ $\text{GeV}\,\text{cm}^{-2}\,\text{s}^{-1}\,\text{sr}^{-1}$ . The upper limits, for gamma-ray energies between 3.9 TeV and 120 TeV, disfavor the emission of hadronic models that try to explain the GeV gamma-ray emission detected by the Fermi LAT. This makes a continuation of the proton injection above 100 TeV highly unlikely (solid cyan line in Figure 9). The HAWC upper limits also disfavor a hadronic injection spectrum derived from IceCube measurements. The present result does not allow unequivocal conclusions about the hadronic or leptonic origin of the Fermi bubbles though. A future analysis of HAWC data will include a better sensitivity, especially at lower energies and possibly larger search regions according to the predictions of some theoretical models.

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, 239762, 254964, 271737, 258865, 243290, 132197), Laboratorio Nacional HAWC de rayos gamma; L’OREAL Fellowship for Women in Science 2014; Red HAWC, México; DGAPA-UNAM (grants RG100414, IN111315, IN111716-3, IA102715, 109916, IA102917); 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-2014/13/B/ST9/945; Coordinación de la Investigación Científica de la Universidad Michoacana. Thanks to Luciano Díaz and Eduardo Murrieta for technical support. We thank Cecilia Lunardini for providing us with the gamma-ray flux model from IceCube data.

Appendix A Table of variables

Appendix B Uncertainty Calculation of G’

The number of gamma rays, as presented in Section 3.2, is given by

[TABLE]

where $G_{i}$ is given by

[TABLE]

The value of $\varepsilon_{C,i}$ is obtained by the equation

[TABLE]

where $j$ are the pixels in the same HEALPix ring as pixel $i$.

The uncertainty $\delta G_{i}$ is calculated as

[TABLE]

The uncertainties of the different terms are $\delta E_{i}^{{}^{\prime}}=\sqrt{\langle N^{\prime}_{i}\rangle}$; $\delta E_{i}=\sqrt{\langle N_{i}\rangle}$; and $\delta\varepsilon_{C,i}=|\varepsilon_{C,i}|\sqrt{\alpha_{i}\left(\frac{1}{\sum_{j}\langle N^{\prime}_{j}\rangle}+\frac{1}{\sum_{j}\langle N_{j}\rangle}\right)}$, where $\alpha_{i}$ is the relative exposure of the observed sky map to the direct integration background. It is calculated as $\alpha_{i}=\Delta\Omega/(\Delta\theta\,\Delta t\,15^{\degree}\,hr^{-1}\,\cos\delta$), where $\Delta\Omega$ is the pixel area, $\Delta\theta$ is the pixel size, $\Delta t$ is the integration time and $\delta$ is the declination.

The systematic uncertainty on $G=\sum_{i}G_{i}$ due to the gamma-ray content in the variable $\varepsilon_{C,i}$ is estimated. First we calculate the relative error of the measured $\varepsilon_{C,i}$ to the true value $\varepsilon^{t}_{C,i}$ where the superscript is for “true”.

[TABLE]

where the numerator is close to the gamma-ray signal after gamma-hadron cuts and the denominator is close to the isotropic background after gamma-hadron cuts. By using the information from Abeysekara et al. (2017b), the relative error is estimated.

The total trigger rate for the HAWC observatory is 24kHz for the 2sr field of view. Assuming that most of the Crab events come from a $1^{\degree}$ radius, we can obtain an estimation of the background rate events from the Crab. This background rate is defined as $BR=24000(2\pi(1-\cos(1^{\degree}))/2sr$. This background rate is proportional to the isotropic component after gamma-hadron cuts $\varepsilon_{C}C_{I}$ The total number of observed events from the Crab is $166.85\text{ events/transit}$ or $0.00772\text{Hz}$ for 6 hour/transit. This is proportional to the excess gamma rays after gamma-hadron cuts $\varepsilon_{G}G^{T}$. We calculate the following ratio

[TABLE]

The ratio $C/C_{I}$ is $O(10^{-4})$(Abeysekara et al., 2014), so the ratio $\frac{N^{\prime}-\langle N^{\prime}\rangle}{\langle N^{\prime}\rangle}$ can be approximated as

[TABLE]

The systematic error in $G_{i}$ can then be written as

[TABLE]

Assuming gaussian regime, $(\delta G)_{sys}=\sum_{i}(\delta G_{i})^{2}_{sys}$, where $i$ is pixel number. The ratio $(\delta G)_{sys.}/(\delta G)_{stat.}$ is of order $O(10^{-4})$, so the systematic uncertainty is ignored.