fermi bubbles paper

arXiv:1005.5480v3 [astro-ph.HE] 18 Oct 2010Draft version October 19, 2010 Preprint typeset using L ATEX style emulateapj v. 03/07/07 GIANT GAMMA-RAY BUBBLES FROM Fermi-LAT: AGN ACTIVITY OR BIPOLAR GALACTIC WIND? Meng Su1,3, Tracy R. Slatyer1,2, Douglas P. Finkbeiner1,2 Draft version October 19, 2010 ABSTRACT Data from the Fermi-LAT reveal two large gamma-ray bubbles, extending 50 degrees a bove and below the Galactic center, with a width of about 40 degrees in longitud e. The gamma-ray emission associated with these bubbles has a significantly harder spectrum ( dN/dE∼E−2) than the IC emission from electrons in the Galactic disk, or the gamma-rays prod uced by decay of pions from proton-ISM collisions. There is no significant spatial variation in the s pectrum or gamma-rayintensity within the bubbles, or between the north and south bubbles. The bu bbles are spatially correlated with the hard-spectrum microwave excess known as the WMAP haze; the edges of the bubbles also line up with features in the ROSATX-ray maps at 1 .5−2 keV. We argue that these Galactic gamma-ray bubbles were most likely cre ated by some large episode of energy injection in the Galactic center, such as past accretion eve nts onto the central massive black hole, or a nuclear starburst in the last ∼10 Myr. Dark matter annihilation/decay seems unlikely to generate all the features of the bubbles and the associated sig nals inWMAP andROSAT; the bubbles must be understood in order to use measurements of the d iffuse gamma-ray emission in the inner Galaxy as a probe of dark matter physics. Study of the origin a nd evolution of the bubbles also has the potential to improve our understanding of recent energe tic events in the inner Galaxy and the high-latitude cosmic ray population. Subject headings: gamma rays — ISM: jets and outflows — galaxies: active — galaxies: st arburst 1.INTRODUCTION The inner Milky Way is home to a massive black hole (MBH), surrounded by clusters of young stars and giant molecular clouds (see e.g. Morris & Serabyn 1996, for a review). The nuclear star cluster has a half-light radius of∼5pc. Although there are indications of past activity, the BH is quiescent today. Fe Kαechoes from molecular clouds around Sgr A∗ have been understood as relics of activity in the past few hundred years (Sunyaev et al. 1993; Koyama et al. 1996). On a longer timescale, one might expect relics of past activity in high energyCRs and hot gas, perhaps far off the disk. The most obvious observables would be e− CR (visible in inverse Compton gammas and microwave synchrotron) and thermal emission (X-rays). This work presents a multiwavelength study of the in- ner Galaxy and identifies several large-scale (tens of de- grees) gamma-ray features, most notably 2 large (span- ning−50◦< b <50◦) structures that we refer to as the “Fermibubbles”. We suggest that these bubbles are associated with previously discovered structures in the X-rays and microwaves, and possibly with analogous smaller-scale structures visible in the FIR. 1Institute for Theory and Computation, Harvard-Smithsonia n Center for Astrophysics, 60 Garden Street, MS-51, Cambridg e, MA 02138 USA 2Physics Department, Harvard University, Cambridge, MA 02138 USA [email protected] High-Energy Excesses Observations of gamma-ray emission in the inner Galaxyat E/lessorsimilar1GeV gobackdecadesto COS-B(Strong 1984;Strong et al.1987), and SAS-2(Fichtel et al.1975; Kniffen & Fichtel 1981) (see Bloemen 1989, for a re- view). Later data from the EGRET experiment aboard theCompton Gamma-ray Observatory extended to the high-energysideofthe π0bump(Smialkowski et al.1997; Dixon et al. 1998). However, EGRET lacked the sensi- tivity and angular resolution to reveal the detailed struc- ture of gamma-ray emission toward the inner Galaxy. TheFermi Gamma-ray Space Telescope provides greatly improved data up to ∼100GeV, with sufficient angular resolution to map out interesting structures.4 At lower energies, the ROSAT All-Sky Survey at 1.5 keV (Snowden et al. 1997) revealed a biconical X-ray structureovertheinnertensofdegreesaroundtheGalac- tic center (GC), later interpreted as a superwind bub- ble (SWB) with energetics of the order of 1054−55ergs (Sofue 2000a; Bland-Hawthorn & Cohen 2003a). On smaller scales, the Midcourse Space Experiment com- bined with IRASdata also confirms the existence of a limb-brightened bipolar structure, the so called Galactic center lobe (GCL), with origin at the GC on the de- greescale (see e.g. Law 2010, for a summary of multi- wavelength observations of GCL). The inferred energy injection of both these bipolar structures, despite their 4Seehttp://fermi.gsfc.nasa.gov/ssc/data/ 2 different scales, is ∼1054−55ergs, with an estimated age of∼106yr for the GCL and ∼107yr for the SWB. Several Galactic Center Shells (GCS), tens of pc in size, have been found with total energy of order ∼1051ergs (Sofue 2003). These shells and filaments are claimed to originate from one or more episodes of rapid energy re- lease. 1.2.Microwave Excess: the WMAP Haze Beyond direct evidence of shell structures, microwave observationsalso provide intriguing indications of energy release toward the GC. At tens of GHz, the Wilkinson Microwave Anisotropy Probe(WMAP)5providessensitivedegreeresolutionfull sky maps of diffuse microwave emission. By subtracting templates including Galactic H α, Haslam 408 MHz soft synchrotron, and thermal dust emission to remove the different known emission mechanisms in these maps, a microwaveresidualexcess(named“themicrowavehaze”) with spherical (non-disklike) morphology about ∼4 kpc in radius toward the GC (visible up to at least |b| ≈30◦) has been recognized (Finkbeiner 2004a). It has a spec- trum of about Iν∼ν−0.5, harder than typical syn- chrotron, but softer than free-free. The microwave haze waslaterinterpretedassynchrotronemissionfromahard spectrum of e−cosmic rays. Other hypotheses such as free-free, spinning dust, or thermal dust have failed to explain its morphology, spectrum, or both (Finkbeiner 2004b; Dobler & Finkbeiner 2008). However, the most recentWMAP 7-year data have not detected the haze polarization predicted by the synchrotron hypothesis (Gold et al. 2010), implying either heavily tangled mag- netic fields, or an alternative emission mechanism. With that caveat in mind, we will assume the WMAP haze is synchrotron and consider the implications. 1.3.A Hard Electron CR Spectrum A simple model, in which the electron CRs that form the haze havediffused from supernova shocks in the disk, cannot fully explain the data for standard diffusion as- sumptions. The 23 −33 GHz spectrum of the haze synchrotron is as hard as that generated from shocks, and it seems extremely unlikely that these electrons can diffuse several kpc from the disk without signifi- cant softening of the spectrum. The synchrotron cooling timescale for cosmic ray electrons emitting at frequency νisτsyn≈106B−3/2 100ν−1/2 GHzyr, where B100=B/100µG (Thompson et al. 2006). Besides the hard spectrum, it is difficult to form the distinctly non-disklike morphology of the haze with any population of sources concentrated in the disk (as is believed to be true of supernovae). The presence of a distinct component of diffuse hard e−CR far off the plane is intriguing in itself, and has motivated proposals where the haze is generated by 5http://map.gsfc.nasa.gov/pulsars, other astrophysical processes, or the annihila- tion of dark matter (Hooper et al. 2007; Cholis et al. 2009a; Zhang et al. 2009; Harding & Abazajian 2010; Kaplinghat et al. 2009; McQuinn & Zaldarriaga 2010; Malyshev et al. 2010). Other indications of excess elec- tronic activity in the Milky Way may be found in re- cent measurements of local electron and positron cos- mic rays. The ATIC, Fermiand H.E.S.S experiments have observed a hardening in the e++e−spectrum at 20−1000 GeV (Chang et al. 2008; Aharonian et al. 2009; Abdo et al. 2009), with an apparent steepening at ∼1 TeV, and the PAMELA experiment has measured a rising positron fraction above 10 GeV. Taken together, these measurementsimply a new sourceofhard electrons and positrons, which maybe related to the WMAPhaze. The coexistence of ROSAT X-ray bipolar features and theWMAP haze toward the inner Galaxy also suggests the interesting possibility of a common physical origin for these signals. 1.4.Inverse Compton Excess from Fermi-LAT Fortunately, if the WMAP haze is synchrotron radi- ation from a hard electron population located around the GC, the same CRs would also produce IC scattered gammas, allowing an independent probe of the CR pop- ulation. IC photons provide valuable complementary in- formation about the spatial distribution of the e−CR (given a model for the ISRF), which in turn can con- strain hypotheses about their origin. The unprecedented angular resolution and sensitivity of Fermi-LAT allows us toprobe the gamma-raycounterpartto the microwave haze in detail for the first time. Previous work employing the first year Fermi-LAT data isolated a spectrally hard “gamma-ray haze” with similar morphology to the WMAP microwave haze (Dobler et al. 2010). In this work, we show that the Fermi-LAT sky maps constructed from 1.6 yr data (600 days) reveal two large gamma-ray lobes, extending 50 degrees above and below the GC, with a width of about 40 degrees in longitude. These two “bubble”-like struc- tures have relatively sharp edges and are symmetric with respect to the galactic plane and the minor axis of the galactic disk. The gamma-raysignal reveals similar mor- phology to the WMAP haze, and is also suggestive of a common origin with features in the ROSATX-ray maps at 1.5 keV towards the GC. As we will discuss, the sharp edges, bilobular shape, and apparent centering on the GC of these structures suggest that they were created by some large episode of energy injection in the GC, such as a past accretion onto the central black hole, or a nuclear starburst in the last ∼10 Myr. It is well known that the GC hosts a massive black hole and massive clusters of recently formed stars (Paumard et al. 2006). Either of these could potentially provide the necessary energy injection by driving large- scale galactic winds or producing energetic jets; we will 3 outline someofthe advantagesand disadvantagesofeach scenario. 1.5.Structure of This Paper In Section 2we briefly review the Fermi-LAT data and our data analysis procedure. In Section 3, we show the 1.6yrFermidatamaps, revealthe Fermibubblefeatures and show that they are robust when different models for the expected Galactic diffuse emission are subtracted. We characterize the morphology of the bubbles in some detail and employ regression template fitting to reveal a hard, spatially uniform spectrum for the gamma-ray emission associated with the bubbles. In Section 4, we show that features in the ROSAT soft X-rays and the WMAPmicrowavehaze are spatially correlated with the Fermibubbles, and the WMAP haze and Fermibub- bles are consistent with being produced from the same electron CR population (by synchrotron and IC respec- tively). Section 5presents our understanding based on the analysis in Section 3and Section 4. Section 6dis- cusses possible scenarios to produce the Fermigamma- ray bubbles. Section 7focuses on the origin of the elec- tron CRs and the challenges in explaining the spectral and spatial profiles of the gamma-ray emission from the bubbles. We discuss the implications of the Fermibub- bles for several topics of interest in Section 8, and give our conclusions and suggest future work in Section 9. 2.FERMI DATA AND MAP MAKING The Large Area Telescope (LAT; see Gehrels & Michelson 1999; Atwood et al. 2009; as well as the Fermihomepage6) is a pair-conversion telescope consisting of 16 modules of tungsten and silicon-strip trackers, on top of a calorimeter with a thickness of 7 radiation lengths. It has a scintillating anti-coincidence detector that covers the tracker array, and a programmable trigger and data acquisition sys- tem. The LAT provides a wide field of view, and covers the energy range from about 30 MeV to 300 GeV. The spacecraft occupies a low Earth orbit with an in- clination of 25 .6◦. The field of view is so wide that the entire sky may be covered in two orbits by rocking the spacecraft north of zenith on one orbit and south of zenith on the next. This scan strategy exposes the LAT to atmospheric gammas at high zenith angles. We make use of only events designated “Class 3” (P6 V3 diffuse class) by the LAT pipeline with a zenith angle cut of 105◦. The events are binned into a full sky map using HEALPix, a convenient iso-latitude equal-area full-sky pixelizationwidelyusedintheCMBcommunity.7Spher- 6http://fermi.gsfc.nasa.gov/ 7HEALPix software and documentation can be found at http://healpix.jpl.nasa.gov , and the IDL routines used in this analysis are available as part of the IDLUTILS product a t http://sdss3data.lbl.gov/software/idlutils .ical harmonic smoothing is straightforward in this pix- elization, and we smooth each map by the appropriate kernel to obtain a Gaussian PSF of 2◦FWHM. Because the PSF of the initial map must be smaller than this, at energies below 1 GeV we use only front-converting events (which have a smaller PSF). A larger smooth- ing scale would help improve S/N, but a relatively small smoothing scale is necessary to see sharp features (such as the bubble edges). Furthermore, for the comparisons and linear combination analysis described in e.g. §3.1.3. it is necessary to smooth the maps at each energy to a common PSF. We generate maps from 1 −4◦, and find that a FWHM of 2◦works well for our purposes. See Dobler et al. (2010) for details on map construc- tion, smoothing, masking, and for instructions on how to download the maps. Our current gamma-ray maps (v2 3) constructed from theFermidata have greater signal/noise compared to the previously released v1 0 maps. They contain photon events from mission times 239557417.494176 to 291965661.204593, for about 606 days or 1.66 years of data (rather than 1 year). We refer to these as the “1.6 year maps.” As in Dobler et al. (2010), we construct maps offront-convertingand back-convertingeventssep- arately, smooth to a common PSF, and then combine them. The point source subtraction has improved: in- stead of interpolating over every source in the 3-month catalog, we use the 1-year catalog,8and subtract each point source from the maps in each energy bin, using es- timates of the PSF from the Fermiscience tools.9For the 200 brightest and 200 most variable sources, the sub- traction is noticeably imperfect, so we interpolate over the core of the PSF after subtracting the best estimate. We take care to expand the mask for very bright sources (Geminga, 3C 454.3, and LAT PSR J1836+5925). The resulting map is appropriate for diffuse work at |b|>3◦. At|b|<3◦the maps are severely compromised by the poor subtraction and interpolation over a large number of point sources. Further details of the map process- ing may be found in Appendix B of Dobler et al. (2010). The v23 maps used in this work and color versions of the map are available for download.10 3.FERMI BUBBLES 3.1.Diffuse Galactic Emission Models At low ( ∼1 GeV) energies, and close to the Galactic plane, the gamma-raysobservedby Fermiaredominated by photons from the decay of π0particles, produced by the collisions of CR protons with ambient ionized gas and dust in the ISM. Collisions of high energy CR elec- tronswith the ISM (primarily protons, but also heavier nuclei) produce bremsstrahlung radiation. The Fermi 8Available from http://fermi.gsfc.nasa.gov/ssc/data 9Seehttp://fermi.gsfc.nasa.gov/ssc/data/analysis/docume ntation/ 10Available at http://fermi.skymaps.info 4 Fermi 1 < E < 2 GeV 180 90 0 -90 -180 00 012345 012345keV cm-2 s-1 sr-1 Fermi 2 < E < 5 GeV 180 90 0 -90 -180 00 -90-4504590 00 01234 01234keV cm-2 s-1 sr-1 Fermi 5 < E < 10 GeV 180 90 0 -90 -180 00 0.00.51.01.52.02.53.0 0.00.51.01.52.02.53.0keV cm-2 s-1 sr-1 Fermi 10 < E < 20 GeV 180 90 0 -90 -180 00 -90-4504590 00 0.00.51.01.52.0 0.00.51.01.52.0keV cm-2 s-1 sr-1 Fig. 1.— All-skyFermi-LAT 1.6 year maps in 4 energy bins. Point sources have been su btracted, and large sources, including the inner disk (−2◦< b <2◦,−60◦< ℓ <60◦), have been masked. all-sky gamma-ray maps in different energy bands are shown in Figure 1. In order to uncover the Fermibubble features better, significant π0emission, bremsstrahlung, and ICemission fromthe Galactic diskmust be removed. We take three approaches for the foreground removal. One is to use the FermiDiffuse Galactic Model provided by theFermiteam11(§3.1.1). The second approach em- ploys a linear combination of templates of known emis- sion mechanisms ( §3.1.2), using existing maps from mul- tiwavelength observations and/or constructed geometric templates. The third approach is taking advantage of the lower energy band 0 .5−1.0 GeVFermimap to form a template of a diffusion emission model ( §3.1.3). 3.1.1. Fermi Diffuse Galactic Model TheFermidiffuse Galactic model12is a comprehen- sive model of Galactic gamma-ray emission from the ISM, and serves as a background estimate for point source removal. This model is based on template fits to the gamma-ray data, and includes an IC compo- nent generated by the GALPROP cosmic ray propagation code.GALPROP calculates the steady state solution to the diffusion-energy-loss equation, given the 3D gas dis- tribution, interstellar radiation field, B-field model, CR diffusion assumptions, and many other input parameters (Strong & Moskalenko 1999; Strong et al. 2009, 2007). 11Seehttp://fermi.gsfc.nasa.gov/ssc/data/access/lat/Back groundModels.html 12Available from http://fermi.gsfc.nasa.gov/ssc/data . The version of the diffuse model we use is gll.iem.v 02.The model is constrained by gamma-ray and microwave observations, locally measured CR spectra, etc. By us- ing a well motivated physical model, one can solve for the spectral and spatialdependence ofthe injection func- tion, i.e. the e−andpCR primary source spectra, as a function of position and energy. The diffuse model is the key connection between the input assumptions and the observables, and is essential for interpretation of the Fermi-LAT data. It is important to make it as complete as possible. In this model, the π0emission is modeled with maps of interstellar gas: H Ifrom the Leiden/Argentine/Bonn (LAB) Galactic Survey (Kalberla et al. 2005) and CO from the CfA composite CO survey (Dame et al. 2001). Becausethe π0emissionisafunction ofboththegasden- sity and the proton CR density, which varies with Galac- tocentric radius, it is desirable to allow the emissivity of the gas to vary. Both the H Iand CO surveys con- tain velocity information, which allows separation into six Galactocentric annuli (rings) with boundaries at 4.0, 5.5, 7.0, 10.0, 16.5, and 50 kpc. The spectrum of each is allowed to float, with the constraint that the sum of the rings along each line of sight approximates the observed signal. This freedom also allows for varying amounts of bremsstrahlung (with varying spectrum) which also scales with the ISM density. The contribution from IC is modeled with GALPROP as described above, and included in the ring fit13. 13A description of this model is available at 5 1 GeV < E < 2 GeV 180 90 0 -90 -180 00 -90-4504590 00 -1.0-0.50.00.51.0 -1.0-0.50.00.51.0keV cm-2 s-1 sr-1 2 GeV < E < 5 GeV 180 90 0 -90 -180 00 -90-4504590 00 -0.50.00.51.0 -0.50.00.51.0keV cm-2 s-1 sr-1 5 GeV < E < 10 GeV 180 90 0 -90 -180 00 -90-4504590 00 -0.4-0.20.00.20.40.60.8 -0.4-0.20.00.20.40.60.8keV cm-2 s-1 sr-1 10 GeV < E < 20 GeV 180 90 0 -90 -180 00 -90-4504590 00 -0.20.00.20.4 -0.20.00.20.4keV cm-2 s-1 sr-1 Fig. 2.— All-sky residual maps after subtracting the Fermidiffuse Galactic model from the LAT 1.6 year maps in 4 energy bi ns (see §3.1.1). Two bubble structures extending to b±50◦appear above and below the GC, symmetric about the Galactic p lane. This procedure provides a diffuse model that faithfully reproduces most of the features of the diffuse Galactic emission. One shortcoming is the existence of “dark gas” (Grenier et al. 2005), clouds with gamma-ray emission that do not appear in the H Iand CO surveys. These features are seen in dust maps (Schlegel et al. 1998) and maysimplybemolecularHcloudsunderabundantinCO. The Fermi diffuse model is primarily intended as a background for point source detection, and comes with a number of caveats. However these caveats apply mainly near the Galactic plane, and at E >50GeV. It is nev- ertheless useful for qualitatively revealing features in the diffuse emission at high latitude. In Figure 2, we show the residual maps after subtracting the Fermidiffuse Galactic model in different energy bins. A double-lobed bubble structure is clearly revealed, with similar mor- phology in the different energy bins. We note that the bubble is neither limb brightened nor centrally bright- ened, consistent with a flat projected intensity distribu- tion. 3.1.2.Simple Template-Based Diffuse Galactic Model Since the dominant foreground gamma-rays originate fromπ0gammas produced by CR protons interact- ing with the ISM, the resulting gamma-ray distribution should be morphologically correlated with other maps of spatial tracers of the ISM. A good candidate is the Schlegel, Finkbeiner, & Davis (SFD) map of Galactic fermi.gsfc.nasa.gov/ssc/data/access/lat/ring forFSSCfinal4.pdfdust, based on 100 µm far IR data (Schlegel et al. 1998). Theπ0/bremsstrahlung gamma-ray intensity is propor- tional to the ISM density ×the CR proton/electronden- sity integrated along the line of sight. As long as the CR proton/electron spectrum and density are approxi- mately spatially uniform, the ISM column density is a good tracer of π0/bremsstrahlung emission. The dust map has some advantages over gas maps: there are no problems with self absorption, no concerns about “dark gas” (Grenier et al. 2005), and the SFD dust map has sufficient spatial resolution (SFD has spatial resolution of 6’, and LAB is 36’). On the other hand, SFD con- tains no velocity information, so it is impossible to break the map into Galactocentric rings. Nevertheless, it is in- structive to employ the SFD map to build a very simple foregroundmodel. Thegoalistoremoveforegroundsina fashion that reveals the underlying structure with as few physical assumptions as possible. We will compare the resulting residuals using this simple diffuse model with those using the Fermidiffuse Galactic model. As an example, we reveal the Fermibubble structure from 1−5 GeVFermi-LAT 1.6 yr data in Figure 3. We use the SFD dust map as a template of the π0gamma foreground. The correlation between Fermiand SFD dust is striking, and the most obvious features are re- moved by this subtraction ( top row in Figure 3). This step makes the bubbles above and below the GC easily visible. The revealed bubbles are not aligned with any structures in the dust map, and cannot plausibly be an 6 Fermi 1 < E < 5 GeV 50 0 -50 -50050 012345 012345keV cm-2 s-1 sr-1 minus dust 50 0 -50 -50050 012345 012345keV cm-2 s-1 sr-1 SFD dust 50 0 -50 -50050 012345 012345keV cm-2 s-1 sr-1 no dust 50 0 -50 -50050 0.51.01.52.02.5 0.51.01.52.02.5keV cm-2 s-1 sr-1 minus disk 50 0 -50 -50050 0.51.01.52.02.5 0.51.01.52.02.5keV cm-2 s-1 sr-1 disk model 50 0 -50 -50050 0.51.01.52.02.5 0.51.01.52.02.5keV cm-2 s-1 sr-1 minus disk 50 0 -50 -50050 0.51.01.52.02.5 0.51.01.52.02.5keV cm-2 s-1 sr-1 minus bubbles 50 0 -50 -50050 0.51.01.52.02.5 0.51.01.52.02.5keV cm-2 s-1 sr-1 bubble model 50 0 -50 -50050 0.00.20.40.60.81.0 0.00.20.40.60.81.0keV cm-2 s-1 sr-1 Fig. 3.— Template decomposition of the Fermi-LAT 1.6 year 1 −5 GeV map (see §3.1.2).Top left: Point source subtracted 1 −5 GeV map, and large sources, including the inner disk ( −2◦< b <2◦,−60◦< ℓ <60◦), have been masked. Top middle: The 1−5 GeV map minus SFD dust map ( top right panel) which is used as a template of π0gammas. Middle row: Theleftpanel is the same as thetop middle panel but stretched 2 ×harder. The middlepanel subtracts a simple geometric disk template (shown in t herightpanel), representing mostly inverse Compton emission, to reveal fe atures close to the Galactic center. Two large bubbles are ap parent (spanning −50◦< b <50◦).Bottom row: Theleftpanel is the same as the middlepanel of the second row . Finally we subtract a simple bubble template ( rightpanel), with a shape derived from the edges visible in the map s, and uniform projected intensity. After subtracting the bubble template, the two bubbles features have nearly vanis hed (bottom middle panel), indicating a nearly flat intensity for the Fermi bubbles. 7 Fermi 1 < E < 5 GeV 180 90 0 -90 -180 00 -90-4504590 00 0.51.01.52.02.53.0 0.51.01.52.02.53.0keV cm-2 s-1 sr-1 Fermi 1 < E < 5 GeV 180 90 0 -90 -180 00 -90-4504590 00 0.51.01.52.02.53.0 0.51.01.52.02.53.0keV cm-2 s-1 sr-1 Fermi 5 < E < 50 GeV 180 90 0 -90 -180 00 -90-4504590 00 0.51.01.52.0 0.51.01.52.0keV cm-2 s-1 sr-1 Fermi 5 < E < 50 GeV 180 90 0 -90 -180 00 -90-4504590 00 0.51.01.52.0 0.51.01.52.0keV cm-2 s-1 sr-1 Fig. 4.— Full sky residual maps after subtracting the SFD dust and dis k templates from the Fermi-LAT 1.6 year gamma-ray maps in two energy bins. Point sources are subtracted, and large sou rces, including the inner disk ( −2◦< b <2◦,−60◦< ℓ <60◦), have been masked. Two large bubbles are seen (spanning −50◦< b <50◦) in both cases. Right panels: Apparent Fermibubble features marked in color lines, overplotted on the maps displayed in the left panels . Green dashed circles above and below the Galactic plane ind icate the approximate edges of the north and south Fermibubbles respectively. Two blue dashed arcs mark the inner (d immer) and outer (brighter) edges of the northern arc – a feature in the northern sky outside the north bubble. The r ed dotted line approximately marks the edge of Loop I. The purple dot-dashed line indicates a tentatively identi fied “donut” structure. artifact of that subtraction. Next, a simpledisk model issubtracted (Figure 3,mid- dle row). The purpose of this subtraction is to reveal the structure deeper into the plane, and allow a harder color stretch. The functional form is (csc |b|)−1 in latitude and a Gaussian ( σℓ= 30◦) in longitude. The disk model mostly removes the IC gamma-rays produced by cosmic ray electrons interacting with the ISRF including CMB, infrared, and optical photons; as discussed previously, such electrons are thought to be mostly injected in the Galactic disk by supernova shock acceleration before dif- fusing outward. Finally, we fit a simple double-lobed geometric bub- blemodel with flat gamma-ray intensity to the data, to remove the remaining large-scale residuals towards the GC (Figure 3,bottom row ). In this model, we identify the approximate edges of the two bubble-like structures towards the GC in the bottom left panel (shown with dashed green line inrightpanels of Figure 4). We then fill the identified double-lobed bubble structure with uni- form gamma ray intensity, as a template for the “ Fermi bubbles” ( bottom right panel of Figure 3). If the Fermi bubbles constitute the projection of a three dimensional two-bubble structure symmetric to the Galactic plane and the minor axis of the Galactic disk, taking the dis- tance to the GC R⊙= 8.5 kpc, the bubble centers areapproximately 10 kpc away from us and 5 kpc above and below the Galactic center, extending up to roughly 10 kpc as the most distant edge from GC has |b| ∼50◦. No structures like this appear in GALPROP models, and in factGALPROP is often run with a box-height smaller than this. Because the structures are so well centered on the GC, they are unlikely to be local. In Figure 4, we show the full sky residual maps at 1 −5 GeV and 5 −50 GeV after subtracting the SFD dust and the disk model to best reveal the Fermibubble features. Although photon Poisson noise is much greater in the 5−50 GeV map, we identify a Fermibubble structure morphologicallysimilar to the structure in the 1 −5 GeV map, present both above and below the Galactic plane. InFigure 5, weshowthefullskymapsat1 −5GeVwith the zenithal equal area (ZEA) projection with respect to both north pole and south pole. We found no interesting features appear near the poles. 3.1.3.Low Energy Fermi Map as a Diffuse Galactic Model In Figure 6, we show the 0 .5−1 GeV and 2 −50 GeV residual maps after subtracting only the SFD dust map as a template of foreground π0gammas. The residual maps should be dominated by IC emission from CR elec- trons interacting with the ISRF. We use the 0 .5−1 GeV maps as a template of IC emission from high energy elec- 8 North 1 < E < 5 GeV l = 0o l = 90o l = 180o l = 270o 012345 012345keV cm-2 s-1 sr-1 North 1 < E < 5 GeV, subtracted l = 0o l = 90o l = 180o l = 270o 0.51.01.52.02.53.0 0.51.01.52.02.53.0keV cm-2 s-1 sr-1 South 1 < E < 5 GeV l = 0ol = 90ol = 180ol = 270o 012345 012345keV cm-2 s-1 sr-1 South 1 < E < 5 GeV, subtracted l = 0ol = 90ol = 180ol = 270o 0.51.01.52.02.53.0 0.51.01.52.02.53.0keV cm-2 s-1 sr-1 Fig. 5.— The zenithal equal area (ZEA) projection with respect to bot h north pole ( upper panels ) and south pole ( lower panels ) for the 1−5 GeV energy band before ( left panels ) and after ( right panels ) subtracting the SFD dust and simple disk templates to revea l the Fermibubbles. tronsscattering starlight, andsubtractthetemplate from higher energy maps (the lower panels of Figure 6). The Fermibubble structures are clearly revealed. We thus conclude that the Fermibubbles are mostly from high energy electron CRs IC scattering on CMB photons, and IR photons at higher energies (see §5for more discus- sion). By comparing the Fermidiffuse Galactic model subtraction (Figure 2) and our simple template model subtraction(Figure 4), wefindthatthebubblestructures are robust to quite different foreground subtractions. It is difficult to see how such emission could arise – espe- cially with sharp edges Section 3.2.1– as an artifact of these subtractions. 3.2.FermiBubbles: Morphology 3.2.1.Morphological Features IntherightpanelsofFigure 4,weillustratetheedgesof theFermibubbles and someother features. We find that theFermibubbles have distinct sharp edges, rather than smoothly falling off as modeled in Dobler et al. (2010).Besides the two bubbles symmetric with respect to the Galactic plane, we find one giant northern arc that em- braces half of the north bubble, that extends from the disk up to b∼50, with ℓranging from roughly −40◦ to 0◦. It has a brighter and sharper outer edge in the 1−5 GeV map. On a even larger scale, we identify a fainter structure extended up to b∼80◦, withℓranging from roughly −80◦to 50◦. We will show in §4that this large extended structure corresponds to the North Polar Spuremission associated with Loop I(as seen for ex- ample in the Haslam 408 MHz map Haslam et al. 1982). We will discussthe possible relationof the Fermibubble, thenorthern arc , and the Loop Ifeature in §8. In the 1−5 GeV map, we also identify a “donut-like” struc- ture in the south sky with branging from roughly −35◦ to 0◦andℓfrom roughly 0◦to 40◦. The coordinates of theFermibubble edges, northern arc ,Loop I, and the “donut” feature identified from the 1 −5 GeV map are listed in Table 1. In Figure 7we compare the Fermibubble morphology 9 Fermi 0.5-1.0 GeV 180 90 0 -90 -180 00 -90-4504590 00 01234 01234keV cm-2 s-1 sr-1 Fermi 2-50 GeV 180 90 0 -90 -180 00 -90-4504590 00 0.00.51.01.5 0.00.51.01.5keV cm-2 s-1 sr-1 Fermi 2-10 GeV minus 0.5-1.0 GeV 180 90 0 -90 -180 00 -90-4504590 00 0.00.20.40.60.8 0.00.20.40.60.8keV cm-2 s-1 sr-1 Fermi 10-50 GeV minus 0.5-1.0 GeV 180 90 0 -90 -180 00 -90-4504590 00 0.00.20.40.60.8 0.00.20.40.60.8keV cm-2 s-1 sr-1 Fig. 6.— Top left: Full sky Fermi-LAT 1.6 year 0 .5−1.0 GeV map subtracts the SFD dust map as a template of π0gammas. Top right:The same as top left panel, but for energy range 2 −50 GeV (note the different gray scale for the two panels). Bottom left: The 2−10 GeV Fermigamma-ray map subtracting the top left 0.5−1.0 GeV residual map which is used as a template of ICS of starlight . Bottom right: The same as bottom left panel but for 10 −50 GeV map subtracting the top left 0.5−1.0 GeV residual map. The Fermi bubble structures are better revealed after subtracting th e lower energy 0 .5−1.0 GeV residual map with extended disk-like emission. in different energy bins. We show the difference of the 1−2 and 2−5 GeV residual maps in the upper panels ; each residual map is the result of subtracting the SFD dust map and the simple disk model to best reveal the Fermibubbles. The difference maps between the 1 −5 and 5−50 GeV maps are shown in the lower panels . The bubble features almost disappear in the difference maps, indicating that different parts of the Fermibubbles have similar spectra. To study the sharp edges of the bubbles at high lati- tude morecarefully, we examinethe (projected)intensity profiles along arcs of great circles passing through the estimated centers of the north and south bubbles, and intersecting the bubble edge (as defined in Figure 4) at |b|>28◦. Along each such ray, we define the intersection of the arc with the bubble edge to be the origin of the coordinate system; we then perform an inverse-variance- weighted average of the intensity profile along the rays (as a function of distance from the bubble edge). We subtract a constant offset from the profile along each ray, prior to averaging the rays together, to minimize aliasing of point sources onto the averaged profile, and then add the averaged offset back in at the end. The inverse variance for each data point is obtained from the Poisson errors in the original photon data, prior to any subtraction of point sources or templates (however, the smoothing of the map is taken into account). When therays are averaged together, the naive inverse variance in the result is multiplied by a factor of the annulus ra- dius (for the points being averaged together) divided by 4πσ2, whereσis the 1σvalue of the PSF, and the an- nulus width is taken to be 1◦(the spacing between the points along the rays; this is comparable to the smooth- ing scale, so there may still be unaccounted-for correla- tions between the displayed errors); this is done to take into account that the number of independent measure- ments being sampled by the rays can be far less than the number of rays, especially close to the center of the bubbles. This procedure is repeated for all the stages of the template subtraction, using the simple disk template forinverseCompton scattering(ICS) forillustration(our conclusions do not depend on this choice). The results are shown in Figure 8for the averaged (1−2)+(2−5) GeV maps, and the averaged (5 −10)+ (10−20) GeV maps. In both energy ranges the edges are clearly visible; in the south, this is true even before any templates are subtracted. The intensity profile of the north bubble is strikingly similar to profile of the south bubble. For both of the north and south bubbles, nosignificantedge-brighteningorlimb-brighteningofthe bubbles is apparent from the profiles, the flux is fairly uniform inside the bubbles. In Figure 9, we plot the intensity profile as a function of latitude from the south to the north pole. We con- 10 Fermi 1-2 GeV minus 2-5 GeV 50 0 -50 -50050 -10123 -10123keV cm-2 s-1 sr-1 Fermi 1-2 GeV minus 2-5 GeV 50 0 -50 -50050 -10123 -10123keV cm-2 s-1 sr-1 Fermi 1-5 GeV minus 5-50 GeV 50 0 -50 -50050 -10123 -10123keV cm-2 s-1 sr-1 Fermi 1-5 GeV minus 5-50 GeV 50 0 -50 -50050 -10123 -10123keV cm-2 s-1 sr-1 Fig. 7.— Comparison of the Fermibubbles between different energy bins. Top left: Subtraction of the 2 −5 GeV residual map from the 1−2 GeV residual map; each residual map is constructed from the data by regressing out the SFD dust and disk templates to best reveal the Fermibubbles. The difference map is consistent with Poisson noise away from the masked region, and the bubble features can hardly be recognized, indicating that different spatial reg ions of the Fermibubbles have the same spectrum. Top right: The same map as leftpanel, but with the Fermibubble features overplotted for comparison. The marked fea tures are the same as those plotted in Figure 4, and are listed in Table 1. Bottom row: Same as the upper panels , but subtracting the 5 −50 GeV residual map from the 1 −5 GeV residual map. 11 -20 0 20 40 Degrees from edge0123456keV cm-2 s-1 sr-1Southern map SFD subtracted SFD and simple disk subtracted SFD template Simple disk template -20 0 20 40 Degrees from edge0123456keV cm-2 s-1 sr-1Northern map SFD subtracted SFD and simple disk subtracted SFD, simple disk, Loop 1 subtracted SFD template Simple disk template Loop 1 template -20 0 20 40 Degrees from edge0123keV cm-2 s-1 sr-1Southern map SFD subtracted SFD and simple disk subtracted SFD template Simple disk template -20 0 20 40 Degrees from edge0123keV cm-2 s-1 sr-1Northern map SFD subtracted SFD and simple disk subtracted SFD, simple disk, Loop 1 subtracted SFD template Simple disk template Loop 1 template Fig. 8.— Intensity as a function of radial distance from the bubble ed ge, averaged over great circle arcs intersecting the bubble center and lying at |b|>28◦. Results are shown for ( left) the southern bubble, and ( right) the northern bubble, for ( top) the averaged 1 −2 and 2−5 GeV maps, and ( bottom) the averaged 5 −10 and 10 −20 GeV maps. Different lines show the results at different stag es of the template regression procedure and the corresponding error s are plotted (see text for an outline of the error analysis). struct great circle arcs perpendicular to the l= 0 great circle, extending 10◦in each direction (east and west), andaveragethe emissionovereachsucharc. The flatness of the bubbles with latitude (except possibly close to the Galactic plane), and the sharp edges at high latitude, are also apparent here. We note that the flat intensity of the bubbles is strik- ing. As we show in Figure 10, if we assume that the Fermibubbles are projected from spherically symmetricthree-dimensional bubbles centered above and below the GC, a non-trivial emissivity distribution in the bubble interior is required to produce a flat projected intensity distribution ( upper left panel of Figure 10). If the “bub- bles” originate from IC scattering, this suggests a rather non-uniform density distribution for the electron CRs, which – combined with a nearly uniform spectral index – presents challenges for many models for the electron injection. The expected intensity profile for a shock gen- 12 -50 0 50 b (°)110keV cm-2 s-1 sr-1Original SFD subtracted SFD, simple disk subtracted SFD, simple disk, Loop 1 subtracted 1-2+2-5 GeV -50 0 50 b (°)110keV cm-2 s-1 sr-1Original SFD subtracted SFD, simple disk subtracted SFD, simple disk, Loop 1 subtracted 5-10+10-20 GeV Fig. 9.— Intensity averaged over the central 20 degrees in longitude , as a function of latitude, for ( left) the averaged 1 −2 and 2−5 GeV maps, and ( right) the averaged 5-10 and 10-20 GeV maps. We construct great cir cle arcs perpendicular to the l= 0 great circle, extending 10◦in each direction (east and west), and average the emission o ver each such arc; the “ b” label corresponding to each arc, and the x-axis of the plot, refers to the value of batl= 0. Different lines show the results at different stages of the template subtraction process. The large oversubtraction at b∼15◦in the north, especially pronounced in the low-energy data, is associated with a bright feature in the SFD dust map. erated bubble with a compressed gas shell is shown in theupper right panel; it is noticeably limb-brightened, in contrast to observations. 3.2.2.Spectrum of Gamma-ray Emission We now attempt to estimate the spectrum of the gamma rays associated with the Fermibubbles, and the spatial variation of the spectrum. In order to reveal the hardness of the spectrum of the Fermibubbles, and quantitatively study the intensity flatness of the bub- ble interiors, we do a careful regression template fitting. First, we maximize the Poissonlikelihood of a simple dif- fuseemissionmodelinvolving5-templates(see §3.1.2). In this model, weinclude the SFD dust map (Figure 3,right panelof the top row) as a tracer of π0emission which is dominant (or nearly so) at most energies on the disk and significant even at high latitudes, the simple disk model (Figure3,right panel of thesecond row ), the bubble template (Figure 3,right panel of the bottom row), the Loop Itemplate (see e.g. Figure 11,left panel of thetop row), and a uniform background as templates to weight theFermidata properly. The significant isotropic background is due to extra- galacticemission and chargedparticle contamination, in- cluding heavy nuclei at high energies. The Fermicollab- oration has measured the extragalactic diffuse emission using additional cuts to reduce chargedparticle contami- nation Abdo et al. (2010): below ∼20 GeV the isotropic contribution in our fits is roughly a factor of 2 largerthan the extragalactic diffuse emission, but has a similar spectral slope. At energies above ∼20 GeV the isotropic contribution becomes much harder, which we attribute to charged particle contamination. Foreachsetofmodelparameters,wecomputethePois- son log likelihood, lnL=/summationdisplay ikilnµi−µi−ln(ki!), (1) whereµiis the synthetic map (i.e., linear combination of templates) at pixel i, andkis the map of observed data. The last term is a function of only the observed maps. The 1σGaussian error is calculated from the likelihood by ∆lnL= 1/2. The error bars are simply the square root of the diagonals of the covariance matrix. We refer to Appendix B of Dobler et al. (2010) for more details of the likelihood analysis. Maps of the models constructed from linear combinations of these five templates, and the residual maps between the Fermidata and the combined templates at different energy bins, are shown in Figure 11. In this fit, we mask out all pixels with Galactic lat- itude|b|<30◦(the dashed black line in the residual maps). Template-correlated spectra for the 5-template fit are shown in Figure 12. The fitting is done with regions of|b|>30◦. For a template that has units (e.g., the SFD dust map is in EB−Vmagnitudes) the correlation spectrum has obscure units (e.g. gamma-ray emission per magnitude). In such a case we multiply the correla- tion spectrum by the average SFD value in the bubble 13 Model 1 of flat intensity 50 0 -50 -50050 050100150200 050100150200keV cm-2 s-1 sr-1 Shell model 50 0 -50 -50050 50100150200 50100150200keV cm-2 s-1 sr-1 WMAP haze model 50 0 -50 -50050 10203040 10203040mK ROSAT model 50 0 -50 -50050 1020304050 102030405010-6 counts s-1 arcmin-2 Fig. 10.— Projected emission from four emissivity distributions to i llustrate the qualitative features of each. Top left: A toy model of theFermibubbles with flat projected intensity. In this model, the vol ume emissivity (assuming the ISM is optically thin to gamma r ays) is proportional to ( R2−r2)−1/2forr < R, and zero otherwise, where R=3.5 kpc is the approximate radius of the two bubbles, and ris the distance to the center of the north or south bubble. Top right: A bubble model with compressed gas shells with a thickness of 0.5 kpc; the electron CR density in the shell is a factor of 5 higher tha n in the interior of the bubbles. For this model, a limb-brigh tened edge of the bubbles is clearly visible, a feature which is not seen in theFermidata.Bottom left: An illustrative toy model for the WMAP haze. The haze synchrotron emissivity depends on the electron CR d ensity and the magnetic field; here we take B=B0e−z/z0, where z0= 2 kpc, and show the line of sight integral of this B-field throug h bubble volume. Even though the synchrotron emissivity is n ot simply the product of CR density times field strength, this panel sugges ts that the decreasing intensity of the WMAP haze at high latitudes is due to the decay of the Galactic magnetic field away from the Galac tic plane. Bottom right: A toy model for the ROSAT X-ray features. The observed soft X-rays are limb-brightened and we assume all g as is uniformly distributed within a compressed shell, with no contribution from the interior, and X-ray emission is proportional to the gas density squared. The thickness of the shell is 1 kpc. region, defined by the bottom right panel of Figure 3, masking out the |b|<30◦region. For the uniform, Loop I, and bubble templates (including inner, outer, north, and south), no renormalization is done. These templates are simply ones and zeros (smoothed to the appropriate PSF), so the outer bubble spectrum is simply the spec- trum of the bubble shell template shown in Figure 15, notthe mean of this template over the whole bubble re- gion. The normalization factors for different templates are listed in Table 2. In Figure 12, we show spectra for π0emission, bremsstrahlung and inverse Compton scattering calcu- lated using a sample GALPROP model (tuned to match locally measured protons and anti-protons as well as lo- cally measured electrons at ∼20−30 GeV), as an indi- cation of the expected spectral shapes. The spectra forthe SFD and the simple disk template reasonably match the model expectations. The dust map mostly traces theπ0emission, and the simple disk model resembles a combination of IC and bremsstrahlung emission. The spectrum for emission correlated with the Fermibubbles is clearly significantly harder than either of these com- ponents, consistent with a flatspectrum in E2dN/dE. This fact coupled with the distinct spatial morphology of theFermibubbles indicates that the IC bubbles are generated by a separate electron component. We also note that the spectrum of the bubble template falls off significantly at energy /lessorsimilar1 GeV. This feature is robust with respect to the choice of templates. The fitting co- efficients and corresponding errors of each template are listed in Table 3. We will discuss some implications of the falling spectrum in §4.3. 14 Todemonstratetherobustnessofthespectrumwehave derived for the Fermibubbles, we make use of the Fermi 0.5−1GeVresidualmap (aftersubtractingthe SFDdust map to largely remove the π0gammas) as a template of IC emission, and performa 4-templatefit ( §3.1.3). These gamma rays mostly originate from IC scattering of a rel- atively soft population ofelectrons in the disk, but might also contain gammas from IC scattering on starlight by a latitudinally extended electron population. We use the SFD dust map as a template for π0gammas as previ- ously, and include the uniform background and the bub- bletemplate asin the previous5-templatefit. The fitting is done with regions of |b|>30◦. For the SFD dust map and theFermi0.5−1 GeV IC template, the correlation coefficients are weighted by the mean of each template in the “bubble” region. The resulting model and the differ- ence maps with respect to the Fermidata, at different energy bands, are shown in Figure 13. The residuals are remarkably small. The spectrum is shown in Figure 14. TheFermi0.5−1 GeV IC template appears to contain a small fraction of π0gammas, but the spectral index is consistent with the predicted GALPROP IC component. The fitting coefficients and corresponding errors of each template are listed in Table 4. By eye, the Fermibubbles appear to possess north- south symmetry and are close to spatially uniform in intensity (with a hard cutoff at the bubble edges). To test these hypotheses more quantitatively, we split the Fermibubble template into the inner bubble andouter shelltemplates ( upper row of Figure 15), or alternatively into thenorthandsouth bubble templates ( lower row of Figure15). We then repeat the previous 5-template and 4-template fitting procedure involving either the simple disk IC template or the Fermi0.5−1 GeV IC template, but splitting the bubble template to either inner bub- ble and outer shell ornorth and south bubble templates. The goal is to identify variations in the intensity and spectral index between the bubble edge and interior, and the northern and southern bubbles. The resulting fitted spectra are shown in Figure 16. And the corresponding fitting coefficients and errors of each template are listed in Table 5, Table6, Table7, and Table 8respectively. In Figure 17, we replace the simple disk model with the Haslam 408 MHz map as the IC template, and employ the SFD map, a uniform background, the Loop Item- plate, and the double-lobed bubble ( leftpanel) or the north and south bubble (rightpanel) templates in the fit- ting. Our conclusion is that the Fermibubbles appear to be north-southsymmetric and spatially and spectrally uniform, with ahard spectrum. This statement is largely independent of our choice of template for the disk IC emission. 3.3.The “FermiBubbles” vs the “ FermiHaze” Dobler et al. (2010) employed essentially the same template regression methods, claimed that the gamma-rayemission not accounted for by the known foregrounds could be well fitted by a bivariate Gaussian with σb= 25◦,σl= 15◦. With the improved 1.6 yr data, the edges of the excess at high latitudes are seen to be quite sharp, and are not well-described by a Gaussian fall-off in in- tensity. However, the question of whether the excess is better modeled asan“egg”ora pairof“bubbles”is more subtle. The choice of the bivariate Gaussian template by Dobler et al. (2010) was intended to remove as much of the remaining gamma ray signal as possible, once the π0 and soft IC emission had been regressed out, minimizing large-scale residuals. The fits were performed with only |b|<5 masked out. In this work, on the other hand, we have masked out all emission with |b|<30◦; when attempting to subtract the disk-correlated emission, and delineating the “bubbles” template, our goal has been to isolate the sharp-edged features from the more slowly spatially varying emission, not to account for all the ob- served emission. This difference in approach can be seen in the non- negligible residuals around the inner Galaxy ( |b|/lessorsimilar20◦) in many of our maps (for example, Figure 13). An at- tempt to fit all the residual emission simultaneously with a simple template may well require a template closer to that used by Dobler et al. (2010), rather than the bubble template. However, the sharp edges now visible in the data, and their alignment with the edges of the WMAP haze and ROSAT X-ray features (as we will discuss in §4), motivate us to consider the bubbles as originating from a distinct physical mechanism. While the bubbles probably do not constitute the entire “ Fermihaze” dis- cussedbyDobler et al.(2010), they arecertainlyamajor component of it, dominating the signal at high latitudes. The question of whether the remaining residual gamma rays represent a separate physical mechanism is an inter- esting one that we defer to future work. 4.THEFERMI BUBBLES SEEN IN OTHER MAPS Inthissection, wecomparethe1 −5GeVFermibubble with the ROSAT 1.5 keV soft X-ray map, the WMAP 23 GHz microwave haze, and the Haslam 408 MHz map. Thestrikingsimilaritiesofseveralmorphologicalfeatures in these maps strongly suggest a common physical origin for theFermibubbles, WMAP haze, and X-ray edges towards the GC (Figure 18). 4.1.Comparison with ROSAT X-ray Features TheROSAT all-sky survey provides full-sky images with FWHM 12’ at energies from 0 .5−2 keV.14We com- pare the morphology of the X-ray features in ROSAT1.5 keV map with the edges of the Fermibubbles in detail in Figure 19. The limb brightened X-ray features align with the edges of both the north and south Fermibub- ble. Hints of the whole north bubble are also visible in 14http://hea-www.harvard.edu/rosat/rsdc.html 15 North bubble South bubble North arc (outer) North arc (inner ) Loop I Donut (outer) Donut (inner) (ℓ,b) [deg] ( ℓ,b) [deg] ( ℓ,b) [deg] ( ℓ,b) [deg] ( ℓ,b) [deg] ( ℓ,b) [deg] ( ℓ,b) [deg] ( 0.0, 0.0) ( 5.5, -5.0) ( 28.5, 5.0) ( 19.5, 5.0) ( 37.5, 25.0) ( 31.9, -5.0) ( 16.8, -7.6) ( -9.9, 5.0) ( 10.7, -10.0) ( 31.1, 10.0) ( 23.1, 10.0) ( 43.2, 3 0.0) ( 34.9, -9.0) ( 22.0, -9.4) ( -14.2, 10.0) ( 12.9, -15.0) ( 31.9, 15.0) ( 24.9, 15.0) ( 46.3 , 35.0) ( 37.0, -14.0) ( 24.6, -14.0) ( -14.5, 15.0) ( 15.0, -20.0) ( 34.0, 20.0) ( 26.0, 20.0) ( 47.7 , 40.0) ( 36.3, -19.0) ( 23.7, -16.5) ( -17.0, 20.0) ( 16.3, -25.0) ( 34.9, 25.0) ( 26.9, 25.0) ( 48.9 , 45.0) ( 33.5, -25.0) ( 22.4, -18.7) ( -22.3, 25.0) ( 16.0, -30.0) ( 33.7, 30.0) ( 23.7, 30.0) ( 49.8 , 50.0) ( 28.7, -29.0) ( 19.2, -21.2) ( -22.6, 30.0) ( 15.5, -35.0) ( 32.5, 35.0) ( 23.5, 35.0) ( 46.6 , 60.0) ( 11.5, -33.0) ( 13.6, -19.3) ( -21.1, 35.0) ( 15.0, -40.0) ( 30.5, 40.0) ( 20.5, 40.0) ( 40.6 , 65.0) ( 4.8, -29.0) ( 13.2, -14.3) ( -19.9, 40.0) ( 15.0, -45.0) ( 27.5, 45.0) ( 16.5, 45.0) ( 29.4 , 70.0) ( -0.4, -20.0) ( 13.6, -8.8) ( -13.6, 45.0) ( 10.6, -50.0) ( 24.7, 50.0) ( 11.7, 50.0) ( 12.5 , 75.0) ( -1.2, -15.0) ( 16.8, -7.6) ( -3.0, 50.0) ( 3.7, -52.5) ( 20.0, 52.5) ( 6.0, 52.5) ( -5.0, 78 .0) ( 0.3, -10.0) ( 1.5, 50.0) ( -6.3, -53.5) ( 14.3, 55.0) ( 0.0, 55.0) ( -19.0, 7 8.0) ( 5.5, -5.0) ( 8.7, 45.0) ( -13.8, -50.0) ( -33.0, 77.0) ( 12.3, 40.0) ( -21.8, -45.0) ( -49.1, 74.0) ( 15.4, 35.0) ( -25.3, -40.0) ( -61.5, 72.0) ( 17.0, 30.0) ( -26.7, -35.0) ( -69.2, 70.0) ( 18.3, 25.0) ( -26.3, -30.0) ( -75.2, 65.0) ( 18.5, 20.0) ( -25.6, -25.0) ( -77.6, 60.0) ( 18.4, 15.0) ( -23.0, -20.0) ( -78.3, 55.0) ( 16.0, 10.0) ( -18.8, -15.0) ( -77.6, 50.0) ( 12.0, 5.0) ( -13.8, -10.0) ( -75.5, 45.0) ( -73.5, 40.0) ( -68.3, 35.0) ( -67.0, 25.0) TABLE 1 Coordinates defining the features shown in Figure 4 and Figu re 7 Uniform SFD bubble north bubble south bubble inner bubble ou ter bubble disk loop I 0 .5−1.0 GeV - SFD 1.0 0.084 1.0 1.0 1.0 1.0 1.0 0.292 1.0 1.198 TABLE 2 Normalization factors for different templates. E range (GeV) Energy Uniform SFD dust simple disk whole bubbl e simple loop I 0.3−0.5 0.4 1.681 ±0.006 1.201 ±0.011 0.689 ±0.027 0.035 ±0.033 0.487 ±0.015 0.5−0.9 0.7 1.365 ±0.007 1.279 ±0.012 0.608 ±0.030 0.211 ±0.037 0.475 ±0.016 0.9−1.7 1.3 1.141 ±0.008 1.179 ±0.014 0.503 ±0.035 0.321 ±0.044 0.405 ±0.019 1.7−3.0 2.2 1.034 ±0.006 0.876 ±0.011 0.393 ±0.029 0.436 ±0.036 0.376 ±0.016 3.0−5.3 4.0 0.881 ±0.008 0.554 ±0.013 0.420 ±0.034 0.353 ±0.043 0.249 ±0.018 5.3−9.5 7.1 0.731 ±0.009 0.322 ±0.014 0.282 ±0.039 0.343 ±0.049 0.208 ±0.021 9.5−16.9 12.7 0.563 ±0.010 0.193 ±0.015 0.251 ±0.044 0.205 ±0.055 0.092 ±0.023 16.9−30.0 22.5 0.507 ±0.012 0.128 ±0.018 0.191 ±0.053 0.263 ±0.068 0.125 ±0.029 30.0−53.3 40.0 0.557 ±0.015 0.041 ±0.021 0.096 ±0.064 0.217 ±0.083 0.088 ±0.036 53.3−94.9 71.1 0.628 ±0.022 0.020 ±0.030 0.183 ±0.093 0.251 ±0.120 -0.043 ±0.048 94.9−168.7 126.5 0.622 ±0.030 0.080 ±0.043 0.012 ±0.125 0.319 ±0.162 -0.091 ±0.063 168.7−300.0 225.0 0.436 ±0.038 0.174 ±0.061 -0.039 ±0.158 -0.015 ±0.194 0.083 ±0.086 TABLE 3 Corresponding template fitting coefficients and errors in Figure 12 . ROSAT, as well as two sharp edges in the south that trace the south Fermibubble close to the disk. We show theROSAT 1.5 keV map overplotted with the edges of theFermibubbles, the northern arc , the“donut” and theLoop Ifeatures in the rightpanels of Figure 19. The appearance of the X-ray edges in the ROSAT 1.5 keV map, coincident with the Fermibubble edges, strongly supports the physical reality of these sharp edges. In Figure 20, we subtract the ROSAT 1.0 keV softX-ray map from the 1.5 keV map to clean up the fore- ground. We find that the extended Loop Ifeature has a softer spectrum than the X-ray features associated with the bubble edges, and is largely removed in the differ- ence map ( lower left panel of Figure 20). The residual features strikingly overlap with the edges of the Fermi bubbles ( lower right panel). No other noticeable large scale features appear in the residual X-ray map which do not appear in the gamma-rays. 16 E range (GeV) Energy Uniform SFD dust whole bubble 0 .5−1.0 GeV - SFD 0.3−0.5 0.4 1.759 ±0.006 0.883 ±0.012 -0.026 ±0.026 1.181 ±0.017 0.5−0.9 0.7 1.446 ±0.006 0.905 ±0.013 0.008 ±0.029 1.275 ±0.018 0.9−1.7 1.3 1.208 ±0.008 0.929 ±0.016 0.258 ±0.034 0.897 ±0.020 1.7−3.0 2.2 1.088 ±0.006 0.679 ±0.012 0.375 ±0.029 0.736 ±0.017 3.0−5.3 4.0 0.921 ±0.007 0.427 ±0.014 0.428 ±0.034 0.530 ±0.019 5.3−9.5 7.1 0.759 ±0.008 0.231 ±0.015 0.371 ±0.039 0.400 ±0.021 9.5−16.9 12.7 0.580 ±0.009 0.131 ±0.016 0.269 ±0.043 0.270 ±0.023 16.9−30.0 22.5 0.522 ±0.011 0.072 ±0.016 0.290 ±0.053 0.261 ±0.020 30.0−53.3 40.0 0.565 ±0.014 0.015 ±0.021 0.225 ±0.066 0.141 ±0.034 53.3−94.9 71.1 0.631 ±0.021 0.010 ±0.033 0.364 ±0.096 0.065 ±0.053 94.9−168.7 126.5 0.615 ±0.029 0.083 ±0.048 0.337 ±0.127 -0.032 ±0.076 168.7−300.0 225.0 0.441 ±0.037 0.149 ±0.068 -0.104 ±0.137 0.089 ±0.097 TABLE 4 Corresponding template fitting coefficients and errors in Figure 14 . 4.2.Comparison with WMAP Microwave Haze TheWMAPhaze is the residual remaining in WMAP microwave data after regressing out contributions from thermal dust, free-free, and “soft synchrotron” traced by the Haslam 408 MHz radio survey (Haslam et al. 1982). Therefore, it is by construction harder than the Haslam- correlated emission. We will show in this section that theWMAP synchrotron haze appears to be associated with the Fermibubbles. Figure21shows a detailed morphological comparison ofthe south Fermibubble at1 −5GeVwith the southern part of the WMAPmicrowavehaze at 23 GHz (K-band). The edge of the Fermibubbles, marked in green dashed line in the top right andlower right panels, closely traces the edge of the WMAP haze. The smaller latitudinal extension of the WMAP haze may be due to the de- cay of the magnetic field strength with latitude. These striking morphological similarities between the WMAP microwave haze and Fermigamma-ray bubble can be readily explained if the sameelectron CR population is responsible for both excesses, with the electron CRs in- teractingwith the galacticmagneticfield toproduce syn- chrotron, and interacting with the ISRF to produce IC emission. In Figure 22, we show the difference maps between the 23 GHzWMAPhaze and the 33 GHz and 41 GHz haze maps. The difference maps contain no apparent features and indicate a common spectrum of different regions in- side theWMAPhaze, suggesting a single physical origin for the bulk of the signal. We have reached the same conclusion for the Fermibubbles in Figure 7. Besides the similarity of the morphology, the relatively hardspectrumofthe Fermibubblealsomotivatesacom- mon physical origin with the WMAP haze. We now provide a simple estimate of the microwave synchrotron andgamma-rayICS signalsfromonesingle populationof hardelectronsdistributed in the inner Galaxy, to demon- stratethatthemagnitudesandspectralindicesofthetwo signals are consistent for reasonable parameter values. For a highly relativistic electron scattering on lowenergy photons, the spectrum of upscattered photons is given by Blumenthal & Gould (1970) (Cholis et al. 2009), dN dEγdǫdt=3 4σTc(mec2)2 ǫE2e/parenleftbig 2qlogq+(1+2q)(1−q) +0.5(1−q)(Γq)2/(1+Γq)/parenrightbig n(ǫ),(2) Γ = 4ǫEe/(mec2)2, q=Eγ Ee1 Γ(1−Eγ/Ee), ǫ < Eγ< EeΓ/(1+Γ). Hereǫis the initial photon energy, Eeis the electron en- ergy,Eγis the energy of the upscattered gamma ray, and n(ǫ) describes the energy distribution of the soft photons per unit volume. Where Γ ≪1, in the Thomson limit, the average energy of the upscattered photons is given by, ∝angbracketleftEγ∝angbracketright= (4/3)γ2∝angbracketleftǫ∝angbracketright, (3) whereγ=Ee/mec2is the Lorentz boost of the electron. In the Klein-Nishina (KN) limit, Γ ≫1, the spectrum insteadpeaksatthehighenergyend,andtheupscattered photon carriesawayalmost all the energyof the electron. Given a powerlaw steady-stateelectron spectrum with spectral index γ, the spectral index of the IC scattered gamma rays is ( γ+ 1)/2 in the Thomson limit, and γ+1 in the extreme KN limit (e.g. Blumenthal & Gould 1970). Photons in the Fermienergy range can be pro- duced by scattering of O(10−100) GeV electrons on starlight, which is (marginally) in the KN regime, or by Thomson scattering of much higher-energy electrons on IR or CMB photons. Consequently, the spectral index of gamma rays might be expected to vary with latitude even if the electron spectral index is uniform, becoming harderathigherlatitudeswherescatteringsintheThom- son limit dominate. Closer to the disk, where much of the ISRF energydensity is in starlight, IC scatteringsare in neither the extreme KN nor Thomson limits, and the spectrum needs to be computed carefully. Weconsiderasteady-stateelectronspectrumdescribed by a power law, dN/dE∝E−γ, with energy cutoffs at 17 Template 50 0 -50 -50050 -1012 -1012keV cm-2 s-1 sr-1 00.5 - 01.7GeV 50 0 -50 -50050 -1012 -1012keV cm-2 s-1 sr-1 Residual map 50 0 -50 -50050 -1012 -1012keV cm-2 s-1 sr-1 Template 50 0 -50 -50050 -1012 -1012keV cm-2 s-1 sr-1 01.7 - 05.3GeV 50 0 -50 -50050 -1012 -1012keV cm-2 s-1 sr-1 Residual map 50 0 -50 -50050 -1012 -1012keV cm-2 s-1 sr-1 Template 50 0 -50 -50050 -1012 -1012keV cm-2 s-1 sr-1 05.3 - 16.9GeV 50 0 -50 -50050 -1012 -1012keV cm-2 s-1 sr-1 Residual map 50 0 -50 -50050 -1012 -1012keV cm-2 s-1 sr-1 Template 50 0 -50 -50050 -1012 -1012keV cm-2 s-1 sr-1 16.9 - 53.3GeV 50 0 -50 -50050 -1012 -1012keV cm-2 s-1 sr-1 Residual map 50 0 -50 -50050 -1012 -1012keV cm-2 s-1 sr-1 Fig. 11.— The models obtained from the multi-template fits, compared w ith theFermimaps, in different energy bins. The left column shows the linear combination of the disk, Loop I, uniform, and bubble templates that provide the best fit to th eFermimaps after subtracting the best fit SFD dust template (shown in the middlecolumn). The difference maps between the combined template a nd the data are shown in the rightcolumn. The template fitting is done for the region with |b|>30◦to avoid contaminations from the Galactic disk (shown with black dashed line in the rightcolumn residual maps). The subtraction of the model largely removes the features seen in theFermimaps with |b|>30◦. We use the same gray scale for all the panels. We find that both the disk IC template and Loop Ifeatures fade off with increasing energy, but the bubble template is re quired for all the energy bands and does not fade off with incre asing energy. The oversubtraction in the residual maps, especially at low er energy bins, is due to the simple disk IC model, which is not a good template across the entire disk. However, outside the masked region, the residual maps are consistent with Poisson noise without obvious large scale features. 18 1 10 100 Photon Energy [GeV]10-810-710-610-5E2 dN/dE [GeV/cm2/s/sr]Simple disk IC template |b| > 30oUniform SFD dust Loop I Whole bubble Simple disk GALPROP π0 decay GALPROP brem GALPROP IC Fig. 12.— Correlation spectra for the 5-template fit employing a simpl e disk model for the IC (and to a lesser degree bremsstrahlung ) emission from supernova-shock-accelerated electrons (se e§3.2.2). The SFD-correlated spectrum is shown by the red shor t-dashed line which roughly traces π0emission (the gray dashed line indicates a GALPROP prediction for π0emission). The disk-correlated emission is shown by the green dashed line, which traces the soft IC (gray tripl e-dot-dashed line) and bremsstrahlung (gray dot-dashed li ne) component. The spectrum of the uniform emission, which traces the isotr opic background (including possible cosmic-ray contamina tion), is shown as a dotted brown line. The solid orange line indicates the spec trum of emission correlated with Loop I, which has a similar spectrum to the disk-correlated emission. Finally, the blue dot-dashed li ne shows the spectrum correlated with the Fermibubble template. The bubble component has a notably harder (consistent with flat) spectr um than the other template-correlated spectra, and the mode ls for the various emission mechanism generated from GALPROP, indicating that the Fermibubbles constitute a distinct component with a hard spectru m. The fitting is done over the |b|>30◦region. Note that these GALPROP “predictions” are intended only to indicate the expected sp ectral shape for these emission components, for reference. The cor relation coefficients for the SFD map and simple disk model are multiplied by the average value of these maps in the bubble region (defined b y thebottom right panel of Figure 3, with a |b|>30◦cut) to obtain the associated gamma-ray emission; see §3.2.2 for details, and Table 2 for a summary of the normalizat ion factors. 0.1 GeV and 1000 GeV. The choice of high-energy cut- off is motivated by the local measurement of the cosmic ray electron spectrum by Fermi(Abdo et al. 2009). We consider a region ∼4 kpc above the Galactic center, as an example (and since both the WMAPhaze and Fermi bubbles are reasonably well measured there), and em- ploy the model for the ISRF used in GALPROP version 50p (Porter & Strong 2005) at 4 kpc above the GC. We normalizethesynchrotrontotheapproximatevaluemea- sured by WMAP in the 23 GHz K-band (Hooper et al. 2007),∼25◦below the Galactic plane, and computethe corresponding synchrotron and IC spectra. The WMAP haze was estimated to have a spectrum Iν∝ ν−β,β= 0.39−0.67 (Dobler & Finkbeiner 2008), cor- responding approximately to an electron spectral index ofγ≈1.8−2.4; Figure 23shows our results for a mag- netic field of 10 µG and 5 µG at 4 kpc above the GC, and electron spectral indices γ= 1.8−3. We find good agreement in the case of α≈2−2.5, consistent with the spectrum of the WMAPhaze. In the default GALPROP exponential model for the Galactic magnetic field, |B|=|B0|e−z/zswith scale 19 Template 50 0 -50 -50050 -1012 -1012keV cm-2 s-1 sr-1 00.5 - 01.7GeV 50 0 -50 -50050 -1012 -1012keV cm-2 s-1 sr-1 Residual map 50 0 -50 -50050 -1012 -1012keV cm-2 s-1 sr-1 Template 50 0 -50 -50050 -1012 -1012keV cm-2 s-1 sr-1 01.7 - 05.3GeV 50 0 -50 -50050 -1012 -1012keV cm-2 s-1 sr-1 Residual map 50 0 -50 -50050 -1012 -1012keV cm-2 s-1 sr-1 Template 50 0 -50 -50050 -1012 -1012keV cm-2 s-1 sr-1 05.3 - 16.9GeV 50 0 -50 -50050 -1012 -1012keV cm-2 s-1 sr-1 Residual map 50 0 -50 -50050 -1012 -1012keV cm-2 s-1 sr-1 Template 50 0 -50 -50050 -1012 -1012keV cm-2 s-1 sr-1 16.9 - 53.3GeV 50 0 -50 -50050 -1012 -1012keV cm-2 s-1 sr-1 Residual map 50 0 -50 -50050 -1012 -1012keV cm-2 s-1 sr-1 Fig. 13.— The same as Figure 11, but using the Fermi0.5−1 GeV residual map after subtracting emission correlated wi th the SFD dust map (to remove the π0gammas) as a template for low-latitude IC emission (origina ting primarily from scatterings on starlight , and likely involving softer supernova-shock-accelerated ele ctrons). The leftcolumn shows the best fit combined template including the res idual 0.5−1 GeV map, uniform, SFD dust, and bubble templates. The middlecolumn shows the Fermi-LAT data in different energy bins. The difference maps between the combined template and the data ar e shown in the rightcolumn. The template fitting is done over the region with|b|>30◦(shown with black dashed line in the rightcolumn residual maps). We use the same grey scale for all the p anels. We find that the bubble template does not fade away with increasing e nergy. The 4-template fit works extremely well, at |b|>30◦the residual maps are consistent with Poisson noise without obvious larg e scale features, and there are no obvious sharp features in t he data closer to the disk. 20 1 10 100 Photon Energy [GeV]10-810-710-610-5E2 dN/dE [GeV/cm2/s/sr]Fermi 0.5 - 1.0 GeV IC template |b| > 30oUniform SFD dust Whole bubble 0.5 - 1.0 GeV Fermi - SFD dust GALPROP π0 decay GALPROP brem GALPROP IC Fig. 14.— Same as Figure 12, but correlation spectra for the 4-templat e fit employing the Fermi0.5−1 GeV residual map (after subtracting the SFD dust) as a template for the starlight IC. The line style is the same as Figure 12. Again, we find that t he spectrum correlated with the Fermibubble template (blue dot-dashed line) is harder (consiste nt with flat in E2dN/dE) than the spectra correlated with the other templates, and the models for the various emis sion mechanism generated from GALPROP, indicating that the Fermibubbles constitute a distinct gamma-ray component with a hard spect rum. The fitting is done for |b|>30◦. As in Figure 12, the correlation spectra have been normalized to a reference region; see §3.2.2 for details, and Table 2 for a summary of the normalizat ion factors. heightzs≈2kpc, this fieldstrengthwouldcorrespondto B0≈30−40µG or even higher. This value is consider- ably larger than commonly used (e.g. Page et al. 2007). However, models with a non-exponential halo magnetic field, as discussed by e.g. Alvarez-Mu˜ niz et al. (2002); Sun et al. (2008), can have ∼10µG fields well off the plane. We note also that the extrapolated value of B0re- quiredtoobtaingoodagreementbetweentheICandsyn- chrotron amplitudes, in the exponential model, is some- what higherthan found byDobler et al.(2010), who per- formed the comparison at 2 kpc. This apparent discrep- ancy originates from the fact that in the haze latitudi- nal profile given by Hooper et al. (2007), the emission falls off rapidly with latitude for 0 > b >−15◦, but then plateaus at b∼ −15−35◦, contrary to expectationsbased on a B-field profile exponentially falling awayfrom z= 0. This suggests either that the magnetic field inside the bubble does not fall exponentially with |z|inside the bubbles, or that the WMAP haze contains a significant free-free component at high latitude. 4.3.Evidence of a ∼700GeV Electron Excess? In Figure 24, we calculate the gamma-ray spectrum fromICscattering, usingthe standardISRF model taken fromGALPROP – as in Figure 23, but with a different en- ergy range for the electron CRs. An electron CR popula- tion with a hard low-energy cutoff at about 500 GeV can fit theFermibubble spectrum better than a single power law extending from 0 .1−1000GeV, due to the downturn in the spectrum in the lowest energy bin. Even a rather hard (dN/dE∼E−2) power law component at 300 −500 GeV produces a long tail at low energies. Interestingly, 21 Fermi bubble interior template 50 0 -50 -50050 0.00.20.40.60.81.0 0.00.20.40.60.81.0keV cm-2 s-1 sr-1 Fermi bubble shell template 50 0 -50 -50050 0.00.20.40.60.81.0 0.00.20.40.60.81.0keV cm-2 s-1 sr-1 Fermi bubble north template 50 0 -50 -50050 0.00.20.40.60.81.0 0.00.20.40.60.81.0keV cm-2 s-1 sr-1 Fermi bubble south template 50 0 -50 -50050 0.00.20.40.60.81.0 0.00.20.40.60.81.0keV cm-2 s-1 sr-1 Fig. 15.— Top row: We split our Fermibubble template (shown in the bottom right panel of Figure 3) into two components for template fitting: an interior template ( top left) and a shell template ( top right ) with uniform intensity, in order to reveal any potential sp ectrum difference with the template fitting technique. Bottom row: We split the Fermibubble template into north bubble template ( bottom left ) and south bubble template ( bottom right ). If the two bubbles have the same origin, they should not onl y have similar morphologies but also consistent spectra. this preferred 500 −700GeV energy range is rather close to the peak in local e++e−cosmic rays observed by ATIC(Chang et al. 2008)15. We note that although the estimated error bars of energy lower than 1 GeV mildly depend on the templates we use in the fitting procedure, the fall-off of the bubble intensity in the lowest energy bins is robust. Figure 23shows the same analysis for 500-900 GeV electrons with two different templates for the disk IC emission, for the bubble interior and shell separately, and for the north and south bubbles. In all cases the same cut-off in the spectrum below 1 GeV is observed. In this case, while the lowest energy bin in the gamma- rays is better fitted, the lack of low-energy CRs means 15However, note that the large peak observed by ATIC-2 and ATIC-4 appears to be in conflict with the Fermimeasurement of the spectrum, which contains no such feature.that synchrotroncan contributeto the WMAPhazeonly at the sub-percent level, unless the magnetic field in the inner galaxy is extreme (and even then the spectrum does not reproduce the observations). In this case an alternate explanation for the WMAP haze would need to be considered. As suggested initially in Finkbeiner (2004a), the WMAP haze could originate from free-free emission (thermal bremsstrahlung), and this would ex- plain the lack of a clear haze signal in WMAP polariza- tion maps. However, the spectrum of the haze is some- what softer than generally expected from free-free emis- sion; the gas temperature required is also thermally un- stable, requiringasignificantenergyinjection ( ∼1054−55 ergs) to maintain its temperature (Hooper et al. 2007; McQuinn & Zaldarriaga 2010). 4.4.Gamma-ray power and e−cosmic ray density 22 E range (GeV) Energy Uniform SFD dust simple disk inner bubbl e outer bubble simple loop I 0.3−0.5 0.4 1.681 ±0.006 1.201 ±0.011 0.683 ±0.027 0.071 ±0.042 0.004 ±0.040 0.490 ±0.015 0.5−0.9 0.7 1.365 ±0.007 1.279 ±0.012 0.607 ±0.030 0.215 ±0.048 0.207 ±0.045 0.475 ±0.016 0.9−1.7 1.3 1.142 ±0.008 1.179 ±0.014 0.498 ±0.036 0.354 ±0.057 0.293 ±0.053 0.407 ±0.019 1.7−3.0 2.2 1.034 ±0.006 0.876 ±0.011 0.403 ±0.029 0.370 ±0.047 0.491 ±0.045 0.373 ±0.016 3.0−5.3 4.0 0.880 ±0.008 0.554 ±0.013 0.426 ±0.035 0.307 ±0.056 0.393 ±0.054 0.247 ±0.018 5.3−9.5 7.1 0.731 ±0.009 0.322 ±0.014 0.284 ±0.039 0.330 ±0.064 0.354 ±0.061 0.207 ±0.021 9.5−16.9 12.7 0.562 ±0.010 0.193 ±0.015 0.265 ±0.044 0.100 ±0.070 0.293 ±0.069 0.087 ±0.023 16.9−30.0 22.5 0.506 ±0.012 0.128 ±0.018 0.201 ±0.054 0.182 ±0.088 0.335 ±0.087 0.121 ±0.029 30.0−53.3 40.0 0.556 ±0.015 0.041 ±0.021 0.098 ±0.065 0.190 ±0.106 0.243 ±0.107 0.087 ±0.036 53.3−94.9 71.1 0.627 ±0.022 0.020 ±0.030 0.187 ±0.093 0.206 ±0.151 0.294 ±0.153 -0.045 ±0.048 94.9−168.7 126.5 0.620 ±0.030 0.081 ±0.043 0.037 ±0.129 0.180 ±0.209 0.431 ±0.207 -0.098 ±0.064 168.7−300.0 225.0 0.435 ±0.038 0.179 ±0.062 -0.065 ±0.155 0.145 ±0.242 -0.178 ±0.206 0.097 ±0.086 TABLE 5 Corresponding template fitting coefficients and errors in upper left panel of Figure 16. E range (GeV) Energy Uniform SFD dust simple disk north bubbl e south bubble simple loop I 0.3−0.5 0.4 1.679 ±0.006 1.205 ±0.011 0.683 ±0.027 -0.042 ±0.050 0.069 ±0.037 0.496 ±0.015 0.5−0.9 0.7 1.363 ±0.007 1.282 ±0.012 0.603 ±0.030 0.140 ±0.056 0.241 ±0.041 0.483 ±0.017 0.9−1.7 1.3 1.138 ±0.008 1.187 ±0.015 0.493 ±0.035 0.171 ±0.065 0.385 ±0.049 0.421 ±0.020 1.7−3.0 2.2 1.030 ±0.007 0.883 ±0.012 0.384 ±0.029 0.288 ±0.054 0.503 ±0.041 0.391 ±0.016 3.0−5.3 4.0 0.881 ±0.008 0.553 ±0.013 0.421 ±0.035 0.368 ±0.064 0.346 ±0.049 0.247 ±0.019 5.3−9.5 7.1 0.731 ±0.009 0.321 ±0.014 0.283 ±0.039 0.354 ±0.072 0.337 ±0.055 0.207 ±0.022 9.5−16.9 12.7 0.562 ±0.010 0.195 ±0.016 0.248 ±0.044 0.157 ±0.078 0.230 ±0.063 0.097 ±0.024 16.9−30.0 22.5 0.509 ±0.012 0.124 ±0.018 0.196 ±0.053 0.353 ±0.100 0.216 ±0.076 0.116 ±0.030 30.0−53.3 40.0 0.556 ±0.015 0.042 ±0.021 0.094 ±0.064 0.190 ±0.117 0.232 ±0.097 0.092 ±0.038 53.3−94.9 71.1 0.628 ±0.022 0.020 ±0.030 0.184 ±0.093 0.264 ±0.164 0.243 ±0.140 -0.045 ±0.049 94.9−168.7 126.5 0.614 ±0.030 0.091 ±0.044 0.004 ±0.126 0.110 ±0.208 0.450 ±0.195 -0.068 ±0.067 168.7−300.0 225.0 0.429 ±0.039 0.186 ±0.064 -0.058 ±0.160 -0.169 ±0.254 0.061 ±0.222 0.107 ±0.092 TABLE 6 Corresponding template fitting coefficients and errors in upper right panel of Figure 16. E range (GeV) Energy Uniform SFD dust inner bubble outer bubb le 0.5−1.0 GeV - SFD 0.3−0.5 0.4 1.759 ±0.006 0.883 ±0.012 0.011 ±0.035 -0.069 ±0.037 1.180 ±0.017 0.5−0.9 0.7 1.446 ±0.006 0.905 ±0.013 0.005 ±0.039 0.012 ±0.041 1.275 ±0.018 0.9−1.7 1.3 1.208 ±0.008 0.929 ±0.016 0.284 ±0.047 0.230 ±0.048 0.896 ±0.020 1.7−3.0 2.2 1.088 ±0.006 0.680 ±0.012 0.309 ±0.039 0.444 ±0.041 0.737 ±0.017 3.0−5.3 4.0 0.921 ±0.007 0.427 ±0.014 0.399 ±0.047 0.459 ±0.049 0.530 ±0.019 5.3−9.5 7.1 0.759 ±0.008 0.231 ±0.015 0.362 ±0.053 0.382 ±0.056 0.400 ±0.021 9.5−16.9 12.7 0.580 ±0.009 0.131 ±0.016 0.195 ±0.057 0.351 ±0.063 0.271 ±0.023 16.9−30.0 22.5 0.523 ±0.011 0.069 ±0.015 0.218 ±0.071 0.368 ±0.079 0.261 ±0.018 30.0−53.3 40.0 0.565 ±0.012 0.014 ±0.015 0.202 ±0.088 0.249 ±0.096 0.144 ±0.025 53.3−94.9 71.1 0.631 ±0.021 0.010 ±0.033 0.337 ±0.128 0.394 ±0.140 0.065 ±0.053 94.9−168.7 126.5 0.614 ±0.029 0.084 ±0.048 0.232 ±0.164 0.455 ±0.190 -0.029 ±0.076 168.7−300.0 225.0 0.440 ±0.037 0.152 ±0.069 0.021 ±0.195 -0.262 ±0.175 0.088 ±0.096 TABLE 7 Corresponding template fitting coefficients and errors in lower left panel of Figure 16. In order to estimate the total gamma-ray power emit- ted by the bubbles, we must estimate the surface brightness integrated over energy, the solid angle sub- tended, and the distance (suitably averaged). From Figure12we take the intensity to be E2dN/dE = 3×10−7GeVcm−2s−1sr−1from 1−100 GeV, integrat- ingto1.4 ×10−6GeVcm−2s−1sr−1. Thebubble template used in our analysis (Figure 3) subtends 0.808 sr, yield- ing a total bubble flux of 1.13 ×10−6GeVcm−2s−1. To obtain an average distance for the emission, we approx- imate the bubbles as 2 spheres centered at b=±28◦,and directly above and below the Galactic center. For a Sun-GC distance of 8.5 kpc, this implies a distance of 9.6 kpc, and a total power (both bubbles) in the 1 −100 GeV band of 2.5 ×1040GeV/s or 4 .0×1037erg/s, which is∼5% of the total Galactic gamma-ray luminosity be- tween 0.1−100 GeV (Strong et al. 2010). Theelectroncosmic-raydensityinthebubblesrequired to generate the observed gamma rays, at any given en- ergy, depends strongly on the assumed electron spec- trum. However, typically the required values are com- parable to the locally measured electron CR density. For 23 E range (GeV) Energy Uniform SFD dust north bubble south bubb le 0.5−1.0 GeV - SFD 0.3−0.5 0.4 1.759 ±0.006 0.882 ±0.012 -0.011 ±0.047 -0.031 ±0.029 1.180 ±0.017 0.5−0.9 0.7 1.446 ±0.006 0.905 ±0.013 -0.009 ±0.052 0.014 ±0.032 1.276 ±0.018 0.9−1.7 1.3 1.207 ±0.008 0.931 ±0.016 0.201 ±0.061 0.277 ±0.038 0.899 ±0.020 1.7−3.0 2.2 1.087 ±0.006 0.681 ±0.013 0.331 ±0.050 0.391 ±0.033 0.738 ±0.017 3.0−5.3 4.0 0.922 ±0.007 0.424 ±0.014 0.500 ±0.060 0.399 ±0.039 0.527 ±0.019 5.3−9.5 7.1 0.760 ±0.008 0.228 ±0.015 0.440 ±0.067 0.342 ±0.045 0.397 ±0.021 9.5−16.9 12.7 0.580 ±0.009 0.132 ±0.016 0.231 ±0.072 0.286 ±0.051 0.272 ±0.023 16.9−30.0 22.5 0.524 ±0.011 0.068 ±0.015 0.406 ±0.091 0.236 ±0.060 0.257 ±0.008 30.0−53.3 40.0 0.565 ±0.014 0.014 ±0.015 0.227 ±0.105 0.222 ±0.077 0.143 ±0.015 53.3−94.9 71.1 0.631 ±0.021 0.011 ±0.033 0.339 ±0.151 0.376 ±0.115 0.067 ±0.054 94.9−168.7 126.5 0.608 ±0.029 0.091 ±0.049 0.067 ±0.189 0.470 ±0.162 -0.013 ±0.077 168.7−300.0 225.0 0.437 ±0.038 0.156 ±0.070 -0.216 ±0.229 -0.066 ±0.160 0.095 ±0.098 TABLE 8 Corresponding template fitting coefficients and errors in lower right panel of Figure 16. example, for the model in the first panel of Figure 24 (dN/dE∝E−2for 500 GeV ≤E≤700 GeV), the in- ferred bubble electron density is ∼10×greater than the local electron density (as measured by Fermi) at an en- ergy of 500 GeV. For a representative model from the first panel of Figure 23(dN/dE∝E−2.3for 0.1 GeV ≤E≤1000 GeV, with a 10 µGmagnetic field generat- ing the WMAP Haze via synchrotron), at 500 GeV the bubble electron density is a factor ∼2×greater than the local density. 5.INTERPRETATION AsdiscussedinDobler et al.(2010), the Fermibubbles seem most likely to originate from IC scattering, since the required electron CR population can also naturally generate the WMAP haze as a synchrotron signal. The ROSAT X-ray measurements suggest that the bubbles are hot and hence underdense regions, and thus argue againstthe gamma raysoriginatingfrom bremsstrahlung orπ0decay. Even though the material in the bubbles is likely high pressure, it is also probably very hot ( ∼107K) and has lower gas density than the ambient ISM. This would explain why the ROSAT 1.5 keV map shows a “cavity” of soft X-rays toward the center of the Fermi bubble structure, like the X-ray cavity in galaxy clus- ters (McNamara & Nulsen 2007), especially in the north Fermibubble. Furthermore, allowing the Fermibubbles to have lower density than the ambient medium means they would experience a buoyant force moving the bub- ble material awayfrom the GC (see §6.2.3for further dis- cussion), which mayhelpgeneratetheobservedmorphol- ogy. Because π0and bremsstrahlung gamma-ray emis- sion both scale as the CR density ×the gas density, an underdense region cannot be brighter unless the cosmic ray densities are greatly increased to compensate; for protons, in particular, the propagation lengths are great enough that a proton overdensity cannot reasonably ex- plain the sharp bubble edges observed in the data, if the bubbles are in a steady state. If the bubbles are expanding rapidly and highly accel- erated protons responsible for the gamma-ray emissionare trapped behind shock fronts, then sharp edges for theFermibubbles could occur naturally. However, in the presence of such a shock, electrons would also be accelerated, and would generally produce more gamma- rays than the protons via ICS (since the cooling time for electron CRs is much shorter than the cooling time for proton CRs of comparable energy). It might be thought that the presence of a bright X- ray edge could lead to a sharp edge in the gamma-ray signal, via IC scattering of electron CRs on the X-ray photons. In the Thomson limit, the energy of IC scat- teredphotonsis oforder(Γ e/2)Ee, with Γ e= 4EeEγ/m2 e (whereEγandEeare the initial photon and electron en- ergies, respectively, and meis the electron mass), and the scattering cross section is independent of the ini- tial electron and photon energies. Thus a higher-energy photon population, leading to a larger value of Γ e, al- lows IC gamma-rays at a given energy to originate from lower-energy electrons, which are much more abundant for typical electron spectra with dN/dE∼E−γ,γ/greaterorsimilar2. However, in the Klein-Nishina regime where Γ e/greaterorsimilar1 this picture changes: the energy of scattered photons is de- termined mostly by the energy of the initial electron, and the cross section scales as 1 /Γe. Scatterings of ∼50 GeV electrons already produce ∼10 GeV gamma rays; when compared to the scattering of 10 GeV electrons on X-rays in the extreme KN limit, the KN suppression of the cross section in the latter case more than counteracts the greater abundance of ∼10 GeV electrons, unless the electron spectrum is very soft (which is inconsistent with the observed signal). Thus hard UV or X-rays in the bubbles, which might be naturally expected in a high-temperature region, wouldnotmake the IC spectrum harder, and IC scat- tering on these photons is subdominant to IC scattering on the usual ISRF for the electron energies in question, unless the X-ray photon number density is much greater than the starlight photon number density. Furthermore, there is no reason to think the bubbles contains more /lessorsimilar1 eV photons, at least not in a region with a well de- fined spatial edge. Thus the sharp edges of the Fermi 24 1 10 100 Photon Energy [GeV]10-810-710-610-5E2 dN/dE [GeV/cm2/s/sr]Simple disk IC template |b| > 30oUniform SFD dust Simple disk Bubble interior Bubble shell Loop I GALPROP π0 decay GALPROP brem GALPROP IC 1 10 100 Photon Energy [GeV]10-810-710-610-5E2 dN/dE [GeV/cm2/s/sr]Simple disk IC template |b| > 30oUniform SFD dust Simple disk North bubble South bubble Loop I GALPROP π0 decay GALPROP brem GALPROP IC 1 10 100 Photon Energy [GeV]10-810-710-610-5E2 dN/dE [GeV/cm2/s/sr]Fermi 0.5 - 1.0 GeV IC template |b| > 30oUniform SFD dust Bubble interior Bubble shell 0.5 - 1.0 GeV Fermi - SFD dust GALPROP π0 decay GALPROP brem GALPROP IC 1 10 100 Photon Energy [GeV]10-810-710-610-5E2 dN/dE [GeV/cm2/s/sr]Fermi 0.5 - 1.0 GeV IC template |b| > 30oUniform SFD dust North bubble South bubble 0.5 - 1.0 GeV Fermi - SFD dust GALPROP π0 decay GALPROP brem GALPROP IC Fig. 16.— Same as Figure 12 and Figure 14, but splitting the Fermibubble template into two components for template fitting. Th e line styles are the same as Figure 12. Top row: Using the simple disk model as the IC template. In the leftpanel, we split the previous bubble template into bubble interior andbubble shell templates (see Figure 15 for the templates). The correlatio n coefficients of the 6-template fit involving the two bubble templates are shown. The purple d ash-dotted line and blue triple-dot-dashed line are for the inner bubble and the outer shell template respectively. The two templates ha ve a consistent spectrum which is significantly harder than t he other templates, indicating the bubble interior and the bubble shell have the same distinct physical origin. In the rightpanel, we split the bubble template intonorthandsouthbubbles. As we include the Loop Itemplate (which has a softer spectrum) in the north sky for re gression fitting, the north bubble has a slightly harder spectrum than the sout h bubble. Again, both of the templates have harder spectra th an any other components in the fit. Bottom row: Employing the Fermi0.5−1 GeV residual map (after subtracting the SFD dust) as a templ ate for thestarlight IC. In the leftpanel, we split the bubble template into bubble interior andbubble shell templates. In the rightpanel, we split the bubble template into north bubble andsouth bubble templates. As in Figure 12, the correlation spectra have bee n normalized to a reference region; see §3.2.2 for details, and Table 2 for a summary of the normalizat ion factors. 25 1 10 100 Photon Energy [GeV]10-810-710-610-5E2 dN/dE [GeV/cm2/s/sr]Haslam IC template |b| > 30oUniform SFD dust Haslam Loop I Whole bubble GALPROP π0 decay GALPROP brem GALPROP IC 1 10 100 Photon Energy [GeV]10-810-710-610-5E2 dN/dE [GeV/cm2/s/sr]Haslam IC template |b| > 30oUniform SFD dust Haslam North bubble South bubble Loop I GALPROP π0 decay GALPROP brem GALPROP IC Fig. 17.— Same as Figure 12, but using the Haslam 408 MHz map instead of t he simple disk model as the IC template. The Haslam map contains a bright feature associated with Loop I(see Figure 18) and is dominated by synchrotron emission fro m softer electron CRs, of energies around 1 GeV; it is not an ideal tracer of IC emissi on which depends on both electron and ISRF distribution. The resulting best-fit spectrum for the bubble template remains harder tha n the other components, but with enhanced lower energy ( /lessorsimilar2 GeV) correlation coefficients compared to Figure 12 and Figure 14. The rightpanel is the same as the leftpanel, but with the bubble template divided into north and south bubbles (see §3.2.2 for more discussion). As in Figure 12, the correlation spectra have been normalized to a reference region; see §3.2.2 for details, and Table 2 for a summary of the normalizat ion factors. bubbles, and the non-uniformity in the emissivity, most likely arise from the electron CR density rather than the photon density. The presence of similar sharp edges in theWMAPhaze (at |b| ≤30◦) supports this hypothesis, if theWMAPhazeis attributed to synchrotronradiation from the electron CRs. Similarly, the elongated shape of the Fermibubble structures perpendicular to the Galactic plane suggests that the electron CR distribution itself is extended per- pendicular to the plane. The Fermibubble morphol- ogy is a strong argument against the possibility that the WMAP haze originates from a disk-like electron distri- bution with significant longitudinal variation of the mag- netic field, as suggested by Kaplinghat et al. (2009). The limb brighteningof the X-raysin the ROSATdata (as shown in Figure 19), and the flat intensity profile of theFermibubbles, suggest the presence of a shell or shock, with increased electron CR density, coincid- ing with a hot thermal plasma. If the ambient medium several kpc above and below the GC were neutral, then bubbles of ionized gas could produce a void in the H I map (Figure 26). We see no evidence for features aligned with the bubbles in these maps, suggesting that the H I map in this part of the sky is dominated by disk emis- sion, and has nothing to do with the bubbles. If the bubbles are in a static state, the bubble edges should havelowertemperaturethanthebubbleinteriorandthus higher gas density, although shocks or MHD turbulence might lead to higher temperatures at the bubble wall.The X-rays in ROSAT may be thermal bremsstrahlung emission, so the emissivity is proportional to the ther- mal electron density ×ion density. They could also arise from charge exchange reactions occurring when the high-speed gas in the bubbles collides with the denser gas at the bubble edge (see Snowden 2009, and reference therein); this mechanism could explain the pronounced limb brightening of the X-rays. As an alternative expla- nation, the ROSAT X-ray feature might be synchrotron emission from very high energy electron CRs. Typically, though, one needs ∼50 TeV ( ∼5 TeV) electrons with ∼10µG (1 mG) magnetic field to produce ∼1 keV synchrotron photons. TheFermibubble features do notappear to be as- sociated with Loop I, a giant radio loop spanning over 100 degrees (Large et al. 1962), which is thought to be generated from the local Sco-Cen OB association. De- tections of Loop Iin high-energy gamma-rays have been claimed by Bhat et al. (1985) and also recently by Fermi (Casandjian et al. 2009); we have also discussed its pres- ence in this work (see §3.2.2). TheLoop Igamma-rays may be the IC counterpart of the synchrotron emission seen in the Haslam 408 MHz map, althoughsome ofthe emission might be π0gammas associated with H I(Figure26). We compare structures identified from the Fermi1−5 GeV maps with Loop I features in the Haslam 408 MHz map in the top row of Figure18, and see that the Fermibubbles are spatially distinct from the arcs associated with Loop I; as we have 26 Fermi 1-5 GeV 50 0 -50 -50050 0.51.01.52.02.53.0 0.51.01.52.02.53.0keV cm-2 s-1 sr-1 Haslam 408 MHz 50 0 -50 -50050 010203040 010203040K Rosat Band 6 and 7 50 0 -50 -50050 50100150200250300350400 5010015020025030035040010-6 counts s-1 arcmin-2 WMAP 23 GHz haze 50 0 -50 -50050 -0.04-0.020.000.020.040.06 -0.04-0.020.000.020.040.06mK Fig. 18.— Comparison of the Fermibubbles with features in other maps. Top left: Point-source subtracted 1 −5 GeVFermi-LAT 1.6 yr map, same as the lower left panel of Figure 3 with north and south bubble edges marked wit h green dashed line, and north arc in blue dashed line. The approximate edge of the Loop Ifeature is plotted in red dotted line, and the “donut” in purp le dot-dashed line. Top right:The Haslam 408 MHz map overplotted with the same red dotted li ne as the top left panel. The red dotted line remarkably traces the edge of the bright Loop Ifeature in the Haslam soft synchrotron map. Bottom left: theROSAT 1.5 keV X-ray map is shown together with the same color lines marking the prominent Fermibubble features. Bottom right: WMAP haze at K-band 23 GHz overplotted with Fermibubble edges. The ROSAT X-ray features and the WMAP haze trace the Fermibubbles well, suggesting a common origin for these features. shown in Figure 16, theLoop Icorrelated emission also has a softer spectrum than the Fermibubble emission. TheLoop Ifeature in the ROSAT map similarly has a softer spectrum than the limb-brightened X-ray bubble edges: as shown in Figure 20, when a low-energy map is subtracted from a higher-energy map in such a way thatLoop Ivanishes, the bubble edges remain bright. We also see additional shell structures which follow the Fermibubble edges and the northern arc in the Haslam 408 MHz map ( top rowof Figure 26). TheFermibubbles are morphologically and spectrally distinct from both the π0emission and the IC and bremsstrahlung emission from the disk electrons. As we have shown in Figure 12to Figure 17, theFermibub- bles have a distinctly hard spectrum, dNγ/dE∼E−2,with no evidence of spatial variation across the bub- bles. As shown in Figure 23, an electron population withdNe/dE∼E−2−2.5is required to produce these gamma rays by IC scattering: this is comparable to the spectrum of electrons accelerated by supernova shocks or polar cap acceleration (Biermann et al. 2010). However, diffusive propagation and cooling would be expected to soften the spectrum, making it difficult to explain the Fermibubbles by IC scattering from a steady-state pop- ulation of these electrons (a single brief injection of elec- trons with dN/dE∼E−2could generate a sufficiently hard spectrum for the bubbles ifthere was a mechanism to transport them throughout the bubble without sig- nificant cooling). The facts strongly suggest that a dis- tinct electron component with a harder spectrum than 27 Rosat Band 6 and 7 -45 -25 0 25 45 00 0306090 00 100200300400 10020030040010-6 counts s-1 arcmin-2 Rosat Band 6 and 7 -45 -25 0 25 45 00 0306090 00 100200300400 10020030040010-6 counts s-1 arcmin-2 Rosat Band 6 and 7 -45 -25 0 25 45 00 -90-60-30000 100200300400 10020030040010-6 counts s-1 arcmin-2 Rosat Band 6 and 7 -45 -25 0 25 45 00 -90-60-30000 100200300400 10020030040010-6 counts s-1 arcmin-2 Fig. 19.— Top left: TheROSAT X-ray haze features compared with the Fermibubbles’ morphology. Top row: The X-ray features of ROSAT band 6 and 7 in the north sky towards to the GC ( leftpanel) compared with Ferminorth bubble overplotted with green dashed line,northern arc feature in blue dashed line, and Loop Ifeature in red dotted line ( rightpanel).Bottom row: Same as top row but for the south bubble features. steady-state SNR generated electrons is responsible for theFermibubbles and associated signals in the WMAP andROSATdata. It has been suggested that a large population of faint millisecond pulsars (MSPs) in the Milky Way halo could contribute to the WMAP andFermihaze signals, via both pulsed gamma-ray emission and e+e−production (Malyshev et al.2010). With ahalopopulationof3 ×104 MSPs, roughly half of the spin-down power going into e+e−pairs and ∼10% going to pulsed gamma-rays, con- sistency with both the WMAP data and the first-year Fermiphoton data is possible; however, this model does not immediately explain either the rather sharp edge in the distribution of gamma-ray emission, or the features in theROSAT X-ray data (the same can of course be said for models which generate the hazes via dark mat- ter annihilation or decay). Other attempts to explain theWMAPhaze with pulsars have generally employed adisk population of pulsars, either peaking in the GC or peaking at small Galactocentric radius but going to zero in the GC (Kaplinghat et al. 2009; Harding & Abazajian 2010); such models have difficulty explaining the spheri- cal morphology of the haze. Furthermore, as pointed out by McQuinn & Zaldarriaga (2010), the regression of the 408 MHz Haslam map should remove much of the con- tribution from young pulsars in the WMAP data, since young pulsars have a similar spatial distribution to su- pernovae. How could the electron CRs possess the same hard spectrum everywhere within the bubble and extend up to 10 kpc, while experiencing a steep fall-off at the bub- ble edge? A large population of CRs might be entrained in large scale Galactic outflows from the GC and enrich the bubbles (see more discussion in §6.3). CRs could be produced along with jets, or shock accelerated CRs from magnetic reconnection inside the bubble or near its sur- 28 Rosat Band 6 and 7 -45 -25 0 25 45 00 0306090 00 100200300400 10020030040010-6 counts s-1 arcmin-2 Rosat Band 5 -45 -25 0 25 45 00 0306090 00 100200300400 10020030040010-6 counts s-1 arcmin-2 Rosat Band 6 and 7 minus Band 5 -45 -25 0 25 45 00 0306090 00 50100150200 5010015020010-6 counts s-1 arcmin-2 Rosat Band 6 and 7 minus Band 5 -45 -25 0 25 45 00 0306090 00 50100150200 5010015020010-6 counts s-1 arcmin-2 Fig. 20.— Top left: TheROSAT X-ray features in Band 6 and 7. Top right: The same region as the leftpanel, but for the ROSAT X-ray map in band 5. Bottom left: The residual X-ray features after subtracting the top right softer Band 5 map from the top leftharder Band 6 and 7 map. Bottom right: The same as the bottom left panel, overplotted with the Fermibubble features, the northern arc , andLoop I features. The residual X-ray features with harder spectrum than the diffuse Loop Ifeature align well with the Fermibubble structures. face (see more discussion in §7). However, it is challeng- ing to produce a flat intensity profile for the bubble inte- rior with a sharp edge. The ambient gas should be com- pressed to a higher density on the shell by shocks(proba- bly also enhancing the magnetic field), and brighter syn- chrotron emission on the shell would then be expected, but the haze emission observed in WMAP is not limb- brightened and shows no evidence for a shell of finite thickness (although we do see shell structure in the X- rays). A cartoon picture summarizing the morphologyof theFermibubbles and associated signals at other wave- lengths is shown in Figure 27. 6.THE ORIGIN OF THE BUBBLE STRUCTURE As we have shown in §3, the bilobular Fermibubble structures are apparently well centered on ℓ= 0◦, and symmetric with respect to reflection across the Galac- tic plane ( bottom row in Figure 3). The bubbles ex-tend down to the plane, where they appear even closer toℓ= 0◦. This alignment and north-south symmetry are unlikely unless the bubbles originate from the GC, and motivate explanations involving Massive Black Hole (MBH) activity ( §6.2) or recent starbursts towards the GC (§6.3). We note that Sco X-1 is approximately cen- tered on the north bubble at the present, but this ap- pears to be a coincidence. Given the similarity between the northern and southern bubbles and the absence of any similar feature in the south, we believe the Fermi bubbles have a GC origin. While the origin of the Fermibubbles is unknown, a rough estimate can be made for their age and the total energy required (although the latter quantity depends linearlyonthe gasdensityin thebubbles, which ispoorly constrained). From the ROSAT data, we envisage the bubblesashotlowdensity( n∼10−2cm−3)cavitiesfilled with∼2 keV gas, with (from the Fermidata) height 29 Fermi 1-5 GeV -45 -25 0 25 45 00 -90-60-30000 0.51.01.52.02.5 0.51.01.52.02.5keV cm-2 s-1 sr-1 Fermi 1-5 GeV -45 -25 0 25 45 00 -90-60-30000 0.51.01.52.02.5 0.51.01.52.02.5keV cm-2 s-1 sr-1 WMAP 23 GHz haze -45 -25 0 25 45 00 -90-60-30000 -0.04-0.020.000.020.040.06 -0.04-0.020.000.020.040.06 mK WMAP 23 GHz haze -45 -25 0 25 45 00 -90-60-30000 -0.04-0.020.000.020.040.06 -0.04-0.020.000.020.040.06mK Fig. 21.— TheFermibubbles at 1 −5 GeV (the residual map obtained by subtracting the SFD dust m ap and the disk template) compared with the WMAP K-band (23 GHz) haze (Dobler & Finkbeiner 2008). Top row: The 1−5 GeV map with ℓ= [−45◦,45◦] and b= [−90◦,0◦] (leftpanel) with the Fermisouth bubble edge overplotted in green dashed line ( rightpanel); the same as the left column third row panel of Figure 3. Bottom row: Same sky region as top row but displaying the WMAP haze at 23 GHz ( leftpanel), with the Fermisouth bubble edge overplotted in green dashed line ( rightpanel). ∼10 kpc, expanding at velocity v/lessorsimilar103km/s: thus we estimate the energy of the Fermibubbles to be about 1054−55erg, with an age of ∼107(v/1000km/s)−1yr. Energetic galactic outflows are common phenomena which have been found in both nearby and high red- shift galaxies (Veilleux et al. 2005). Galactic winds are believed to have significant impact on galaxy formation, morphology, and their environments. The main sources ofthe energyarestellarwinds, supernovae, and/oractive galactic nuclei (AGN); CRs and magnetic field pressure can also help drive galactic outflows (Everett et al. 2010; Socrates et al. 2008). In this section we present some ideas for past activity in the GC that could help to generate the shape of the Fermigamma-ray bubbles and the associated signals at other wavelengths. Section 6.1presents and discusses ev- idence of past GC activity. Section 6.2discusses the hy-pothesis that outflows from a black hole accretion event, including possible AGN jets, could form the Fermibub- bles. Section 6.3focuses on the possibility of a previous starburst phase of the Milky Way, with the Fermibub- bles being inflated by the subsequent energetic Galactic wind. 6.1.Observational Evidence of Previous GC Activity Observationsacrosstheelectromagneticspectrumhave provided constraints on dynamics and evolution of the centralgasand stellarpopulations, and currentandprior accretionactivity ofthe centralMBH. We highlightsome of the evidence of previous activities towards the GC, whichmayrelatetotheproductionofthe Fermibubbles. X-ray Reflection Nebulae in the GC —There are indica- tions of previous GC activity from X-ray echoes and time variability of reflected X-ray lines from cold iron 30 WMAP 33 minus 23 GHz haze -45 -25 0 25 45 00 -90-60-30000 -0.02-0.010.000.010.020.03 -0.02-0.010.000.010.020.03 mK WMAP 33 minus 23 GHz haze -45 -25 0 25 45 00 -90-60-30000 -0.02-0.010.000.010.020.03 -0.02-0.010.000.010.020.03mK WMAP 41 minus 23 GHz haze -45 -25 0 25 45 00 -90-60-30000 -0.02-0.010.000.010.020.03 -0.02-0.010.000.010.020.03 mK WMAP 41 minus 23 GHz haze -45 -25 0 25 45 00 -90-60-30000 -0.02-0.010.000.010.020.03 -0.02-0.010.000.010.020.03mK Fig. 22.— Difference maps of the WMAP haze.Top row: The difference map between the 23 GHz WMAP haze and 33 GHz WMAP haze. The rightpanel is the same as the left but overplotted with the south Fermibubble edge in green dashed line. Bottom row: The same as the top row, but showing the difference between the 41 GHz WMAP haze and 33 GHz WMAP haze. We see no apparent structure in the difference maps, indicating a consistent spectrum acr oss different (spatial) regions of the haze. atoms in molecular clouds around Sgr A∗including Sgr B1 and B2, Sgr C, and M0.11-0.11 (Sunyaev et al. 1993; Sunyaev & Churazov 1998). The changes in the inten- sity, spectrum, and morphology of the fluorescent iron nebulae near the GC, observed by ASCAandINTE- GRALare likely due to reflected X-rays from previ- ous activity of Sgr A∗with high luminosity ∼300 yr ago. The luminosity is ∼1.5×1039erg s−1in 2−200 keV with a power law spectrum dN/dE∝E−γwith γ=1.8±0.2 (Ponti et al. 2010; Revnivtsev et al. 2004; Nobukawa et al. 2008). The changes in the intensities and morphologies of hard X-ray nebulosities on parsec scales have been discovered (Muno et al. 2007). Outflows: Galactic Center Lobe and Expanding Molecu- lar Ring —Sofue (2000b) and Bland-Hawthorn & Cohen (2003b) have previously noted the presence of the ex- tended bipolar structure in the ROSAT data, and haveattributed it to a large-scale bipolar wind, powered by a starburst in the GC. In this picture, overpressured bubbles rise and expand adiabatically away from the injection region in the Galactic plane, driving shocks into the surrounding gas by ram pressure. Reflection of the shocks increases the density and temperature of the post-shocked gas, which has been suggested to have T≈5×106Kandn≈0.1cm−3at2 kpc abovethe plane (Bland-Hawthorn & Cohen2003b). Free-freeemissionin the resulting high-temperature plasma produces the ob- served X-ray signals. The Galactic Center Lobe (GCL) ondegreescale has estimated total kinetic energy ∼1055 erg and a dynamical time scale of ∼106yr. The size, energy,andtimescalesaresimilartothoseoftheexpand- ing molecular ring (EMR) around the GC (Kaifu et al. 1972; Scoville 1972; Totani 2006). Sofue (2000b) inter- preted the North Polar Spur (NPS) with tens of degrees scale to be an outflow from the GC with energy scale of 31 0.1 1.0 10.0 100.0 1000.0 Eγ [GeV]10-810-710-610-510-4E2 dN/dE [GeV cm-2 s-1 sr-1]γ=1.5 γ=1.8 γ=2.4 γ=3.0Inverse Compton Scattering 0.1 GeV to 1000.0 GeV 0.1 1.0 10.0 100.0 1000.0 Frequency [GHz]0.010.101.0010.00100.001000.00Iν [kJy/sr]WMAP 23 GHzSynchrotron B= 10.0 µG γ=3.0γ=2.4γ=1.8γ=1.50.1 1.0 10.0 100.0 1000.0 Eγ [GeV]10-810-710-610-510-4E2 dN/dE [GeV cm-2 s-1 sr-1]γ=1.5 γ=1.8 γ=2.4 γ=3.0Inverse Compton Scattering 0.1 GeV to 1000.0 GeV 0.1 1.0 10.0 100.0 1000.0 Frequency [GHz]0.010.101.0010.00100.001000.00Iν [kJy/sr]WMAP 23 GHzSynchrotron B= 5.0 µG γ=3.0γ=2.4γ=1.8γ=1.5 Fig. 23.— The estimated spectrum of IC gamma rays ( upper panel ) and synchrotron radiation ( lower panel ) originating from a hard electron spectrum along a line of sight 4 kpc above the Galact ic center (i.e. b≈25◦). The steady-state electron spectrum is taken to be a power law, dN/dE∝E−γ, with index γ= 1.5 (solid black ), 1.8 (blue dashed ), 2.4 (green dotted ), and 3.0 ( red dash-dotted ). In all cases the spectrum has a range of [0.1, 1000] GeV. The interstellar rad iation field model is taken from GALPROP version 50p, and the magnetic field is set to be 10 µG for the leftpanel and 5 µG for the rightpanel. The data points in the upper panels show the magnitude of the bubble emission obtained from template fitting in Figure 12 ( brown) and Figure 14 ( black) including the “whole bubble” template, as a function of energy. The lowest and highest bins contain 3 σupper limits rather than data points with 1 σerror bars, due to the large uncertainties in the haze amplitude at those energies. For reference, a rec tangular cross hatch region shows a approximate spectrum in the same place in this and subsequent figures. The data point in the lower pan el shows the magnitude of the WMAP haze averaged over b=−20◦to −30◦, for|ℓ|<10◦, in the 23 GHz K-band (the overall normalization is chosen to fit this value), and the gray area indicates the range of synchrotron spectral indices allowed for the WMAP haze by Dobler & Finkbeiner (2008). The same population of ha rd electrons can consistently generate the WMAP synchrotron haze and FermiICS bubbles. ∼1055−56erg and time scale of ∼107yr. Totani (2006) suggested that all these outflows can be attributed to the past high activity of Sgr A∗of a duration of ∼107yr, comparable to the reasonable estimation of the lifetimes of AGNs. Diffuse X-ray Emission —Muno et al. (2004) studied the diffuse X-ray emission within ∼20 pc of the GC in de- tail using Chandra observations. The hard component plasma with kT∼8 keV is spatially uniform and cor- related with a softer component with kT∼0.8 keV. Neither supernova remnants nor W-R/O stars are ob- served to produce thermal plasma hotter than kT∼3 keV. AkT∼8 keV plasma would be too hot to be confined to the GC and would form a plasma wind. A large amount of energy input ∼1040erg/s is re- quired to sustain such hot plasma. If the hot plasma is truly diffused, the required power is too large to beexplained by supernova explosions and the origin of this hot plasma, and might be explained as a result of shock heating by the wind or AGN activities. Similar diffuse hard X-ray emission has been detected from the star- burst galaxy M82 (Strickland & Heckman 2007). How- ever, Revnivtsev et al. (2009) haveresolved ∼80percent of the hard diffuse gas into faint point sources. Bipolar Hard X-ray in the GC —Chandra observations show the morphology of hot gas with a few keV seems to be bipolar with each lobe extending to ∼10 pc (see Markoff 2010, and reference therein). Three explana- tions have been suggested for the origin of the lobes: thermal wind from the central cluster of massive young stars, steady outflows from Sgr A∗, or repeated episodic outbursts from Sgr A∗(Markoff 2010). For the collective stellar wind interpretation, it is unclear why the lobes only extend up to 10 pc which has a estimated flow time 32 0.1 1.0 10.0 100.0 1000.0 Eγ [GeV]10-910-810-710-610-5E2 dN/dE [GeV cm-2 s-1 sr-1]Fermi 0.5 - 1.0 GeV IC template Simple disk IC templateICS 500.0- 700.0 GeV 0.1 1.0 10.0 100.0 1000.0 Eγ [GeV]10-910-810-710-610-5E2 dN/dE [GeV cm-2 s-1 sr-1]Fermi 0.5 - 1.0 GeV IC template Simple disk IC templateICS 1000.0-1200.0 GeV 0.1 1.0 10.0 100.0 1000.0 Eγ [GeV]10-910-810-710-610-5E2 dN/dE [GeV cm-2 s-1 sr-1]Fermi 0.5 - 1.0 GeV IC template Simple disk IC templateICS 300.0- 500.0 GeV 0.1 1.0 10.0 100.0 1000.0 Eγ [GeV]10-910-810-710-610-5E2 dN/dE [GeV cm-2 s-1 sr-1]Fermi 0.5 - 1.0 GeV IC template Simple disk IC templateICS 1000.0-1200.0 GeV Fig. 24.— The estimated spectrum of IC gamma rays originating from a ha rd electron spectrum ( dN/dE∝E−2) with a limited energy range, as in the top row of Figure 23, but with different minimum and maximum energies . The normalization of the ICS signal is fitted to the data. The three peaks are from ICS of the CMB (left peak) , FIR (middle peak), and optical/UV interstellar radiation field. The ISRF model is taken from GALPROP version 50p. A hard electron CR population with a low energy c utoff at about 500 GeV can fit the data better than a power law electron spectrum extending from 0.1 GeV to 1000 GeV (see Figure 23). ∼3×104yr, whereas the star cluster is ∼6 millon years old. Moreover, discrete blobs have been found within the lobes, and quasi-continuous winds are hard pressed to explain the origin of the blobs. The possibility of a transient jet-like feature is intriguing. The jets can be produced by accreting the debris of tidal disruption stars (see§6.2.1for more discussion). OB Stellar Disk —There are two young star disks that have been identified in the central parsec of the GC (Paumard et al.2006). In situstarformationfromdense gas accretion disks is favored over the inspiraling star cluster scenario (see e.g. Bartko et al. 2009). The gas disks could be formed as a consequence of a large inter- stellar cloud captured by the central MBH, which then cooled and fragmented to form stars. Interestingly, the two star formation events happened near-simultaneously about 6 ±2 Myr ago and the two disks are coeval to within∼1 Myr (Paumard et al. 2006) and little activity has occurred since. The two young massive star clusters ArchesandQuintuplet in the central 50 pc, with similar stellar mass, content and mass functions, were formed ∼107yr ago. It has been suggested in Paumard et al.(2006) that a global event may cause an increase in the rate of star formation, such as a passing satellite galaxy that enhanced the clouds’ collision frequency and may also provide the gas for active star formation. The unique characteristics of stellar clusters towards the GC have been used to explain the origin of the mag- netized nonthermal radio filaments, threads, and streaks (LaRosa et al. 2004; Yusef-Zadeh et al. 2004). Collec- tive winds of massive W-R and OB stars within such a dense stellar environment can produce terminal shocks that accelerate relativistic particles. The abundance and characteristics of these nonthermal radio filaments within the inner 2 degrees of the GC region can be evidence of an earlier starburst (Rosner & Bodo 1996). Yusef-Zadeh & K¨ onigl (2004) propose a jet model in which the characteristics common to both protostellar and extragalactic jets are used to explain the origin of nonthermal filaments in the Galactic center region. 6.2.Outflow from Black Hole Accretion Events The central MBH in our Milky Way, with an esti- mated mass of M BH∼4×106M⊙(e.g. Ghez et al. 2008; Gillessen et al. 2009), is currently quiescent, radiating 33 0.1 1.0 10.0 100.0 1000.0 Eγ [GeV]10-910-810-710-610-5E2 dN/dE [GeV cm-2 s-1 sr-1]Inner bubble Outer bubbleInverse Compton Scattering 500.0 GeV to 900.0 GeV Fermi 0.5 - 1.0 GeV IC template 0.1 1.0 10.0 100.0 1000.0 Eγ [GeV]10-910-810-710-610-5E2 dN/dE [GeV cm-2 s-1 sr-1]Inner bubble Outer bubbleInverse Compton Scattering 500.0 GeV to 900.0 GeV Simple disk ICS template 0.1 1.0 10.0 100.0 1000.0 Eγ [GeV]10-910-810-710-610-5E2 dN/dE [GeV cm-2 s-1 sr-1]North bubble South bubbleInverse Compton Scattering 500.0 GeV to 900.0 GeV Fermi 0.5 - 1.0 GeV IC template 0.1 1.0 10.0 100.0 1000.0 Eγ [GeV]10-910-810-710-610-5E2 dN/dE [GeV cm-2 s-1 sr-1]North bubble South bubbleInverse Compton Scattering 500.0 GeV to 900.0 GeV Simple disk ICS template Fig. 25.— The estimated spectrum of IC gamma rays originating from a ha rd electron spectrum, the same as Figure 24, but for the 500−900 GeV energy range. Top row: data points show the separately fitted spectra for the bubble interior and outer shell templates (as defined in Figure 15), with the template for the disk IC emissi on given by ( leftpanel) the 0 .5−1 GeVFermimap with dust-correlated emission subtracted, and ( rightpanel) the simple geometric disk model defined in Figure 3. Bottom row: as top row, but showing the separately fitted spectra for the north and south Fermibubbles (as defined in Figure 15). at about 8 orders of magnitude lower than the Edding- ton luminosity ( LEdd∼1044erg/s). The X-rays of Sgr A∗is weak and thermal, in contrast to the hard nonthermal power law typically observed in most low- luminosity active galactic nuclei (LLAGN; see e.g. Ho 2008). Fast X-ray flaring of Sgr A∗via nonthermal pro- cesses has been discovered by Chandra ACIS observa- tion (Baganoff et al. 2003). The observed submillime- ter/IR bump is seen to flare simultaneously with the X-rays. Due to the short timescale of the flares and the lack of evidence of any standard thin accretion disk emission (Shakura & Sunyaev 1973), a magnetic origin of the flares has been suggested, either from synchrotron or synchrotron-self Compton emission (Eckart et al. 2004; Yusef-Zadeh et al. 2006; Dodds-Eden et al. 2009). Both radiative inefficient accretion flow (RIAF) mod- els and outflow-dominated models have been invented (Falcke & Markoff 2000; Blandford & Begelman 1999; Yuan et al. 2002). Clearly the MBH has not always been so underlumi- nous: it may have experienced a long active state in the past few million or tens of million years through one or more accretionevents, driving jets out of the disk, shock-ing the ambient material, producing both gamma-rays and CRs, and appearing more similar to normal low- luminosity AGN. If the MBH were radiating at the Ed- dington luminosity, it would take only ∼103−4years to reach the estimated energy of the Fermibubbles; for a percent level accretion rate ( ∼1042erg/s), it would take ∼105−6yr, comparable to the estimated cooling time of the electron CRs. 6.2.1.Scenarios for MBH Jet Formation in the GC Currently, the X-ray flares are the only unambiguous AGN-like activity that Sgr A∗displays (Markoff 2010). However, the synchrotron radio emission is comparable to typical LLAGN with compact jets in terms of spectral characteristics. It is well known that jets associated with accretion disks surrounding black holes can efficiently generate high energy particles. BL Lacs can form rel- ativistic jets to produce TeV gamma rays, as can mi- croquasars. However, we know through multiwavelength observations that the central MBH in our Milky Way has extraordinarily low bolometric luminosity of ∼1036 ergs/s, and so is currently in its quiescent dim state, and no jets toward the GC have been physically resolved. 34 Haslam 408 MHz 50 0 -50 -50050 010203040 010203040K Haslam 408 MHz minus disk 50 0 -50 -50050 010203040 010203040K Haslam 408 MHz minus disk 50 0 -50 -50050 010203040 010203040K LAB HI 50 0 -50 -50050 0100200300400500600 0100200300400500600K LAB HI minus disk 50 0 -50 -50050 0100200300400500600 0100200300400500600K LAB HI minus disk 50 0 -50 -50050 0100200300400500600 0100200300400500600K Hα 50 0 -50 -50050 01234 01234R Hα minus disk 50 0 -50 -50050 01234 01234R Hα minus disk 50 0 -50 -50050 01234 01234R Haslam disk 50 0 -50 -50050 010203040 010203040K HI disk 50 0 -50 -50050 0100200300400500600 0100200300400500600K Hα disk 50 0 -50 -50050 01234 01234R Fig. 26.— Fermibubble features in other maps. Top row: Theleftpanel shows the half sky Haslam 408 MHz map (Haslam et al. 1982) with −90◦< ℓ <90◦, themiddlepanel subtracts a simple geometric disk template (shown in t hebottom left panel) to better reveal the structures deeper into the Galactic plane. The rightpanel is the same as the middlepanel but overplotted with the Fermibubbles, thenorthern arc , and the Loop Ifeatures identified from the 1 −5 GeVFermigamma-ray map (see Figure 4). The Loop Ifeature ( red dotted line ) align with the extended diffuse features in the Haslam 408 MH z synchrotron map (known as North Polar Spur ). The inner and outer edges of the northern arc (dashed blue lines ) overlap with two arcs in the Haslam synchrotron map. Howeve r, theFermibubbles have no apparent counterparts in this map. Second row: The same as the top row, but for the Leiden/Argentine/Bonn (LAB) Survey of Galactic H I(Kalberla et al. 2005). The middlepanel subtracts a simple disk template shown in the bottom middle panel to better reveal structures towards the GC. No apparent features have been id entified that correlate with the Fermibubbles and other features in 1 −5 GeVFermimap (there may be some faint filaments morphologically traci ng the gamma-ray features). Third row: The same as the top rowbut for the H αmap (Finkbeiner 2003). The middlepanel subtracts a simple disk template shown in the bottom right panel to reveal more structures in the inner Galaxy. No corresponding featu res have been identified morphologically similar to the stru ctures in the 1 −5 GeVFermigamma-ray maps (color line in the rightpanel). 35 Detailed examinations of the GCL have shown that the gas shell is deep into the disk, and do not support a jet origin for that structure (Law 2010). However, weak jets consistent with the spectrum of Sgr A∗can be easily hid- den by the blurring of their photosphere (Markoff et al. 2007). The jets related to the energetic accretion in the GC cangoinanydirection; theyneed notbe alignedwith the minoraxis ofthe Galaxy. This mechanismdoes not obvi- ously provide a natural explanation for the north-south symmetry of the Fermibubbles, and the relatively flat gamma-rayintensity inside the bubbles. Also, there is as yet no conclusive evidence for the presence of jets toward the GC (Muno et al. 2008). However, CR rich wide jets might have existed in the past, and the Fermibubbles might have been inflated by the jets in a relatively short time scale without significant cooling. Although the past AGN phase could be the primary heating source of the Fermibubbles, the low density and momentum associ- ated with the jets do not readily distribute the thermal energyisotropically,makingitharderto explainthemor- phologyandflatintensity ofthe Fermibubbles. Ifastar- burst phase in the GC coincided with the energetic jets, SNe in the starburst might provide a large injection of momentum and turbulence, which could help isotropize the energy distribution. Accretion of Stars —One way to form jets in the GC is for the MBH to tidally disrupt and swallow stars in the nuclear star cluster (Hills 1975; Rees 1988). The typi- cal picture is that when a star trajectory is sufficiently close to the MBH, the star is captured and disrupted by tidal forces, and after about an orbital timescale, a transient accretion disk forms from the debris of the dis- rupted star. The capture rate of the GC MBH has been estimated at ∼4.8×10−5yr−1for main-sequence stars and∼8.5×10−6yr−1for red giant stars. Sgr A∗could temporarily behave like an AGN and produce a powerful jet by accreting the debris of such stars, which may be ejected from surrounding molecular disks. If a 50 solar mass star is captured by the MBH in the GC, it gives an energy in relativistic protons as high as ∼1054−55ergs on a very short timescale ( ∼103−4yr), at a rate of about ∼1043ergs/s. Accretion of ISM —Quasi-periodic starbursts in the GC have been recently suggested as a result of the in- teractions between the stellar bar and interstellar gas (Stark et al. 2004). In this scenario, gas is driven by bar dynamics toward the inner Lindblad resonance, and ac- cumulates in a dense gas ring at 150 pc until the critical density threshold is reached. Then giant clouds can be formed on a short timescale, and move toward the cen- ter by dynamical friction. The timescale for this cycle to repeat is highly uncertain but is estimated to be of order 20 Myr (Stark et al. 2004).Accretion of IMBH —A single 104solar mass BH spiral- ing in to the GC may also trigger starbursts and change the spin of the Sgr A∗Gualandris & Merritt (2009), thereby producing precessing jets. It has been argued that one such event happens approximately every 107 years in order to create a core of old stars in the GC, of radius 0.1 pc. Chandra has detected more than 2000 hard X-ray (2 −10 keV) point sources within 23 pc of the GC (Muno et al. 2003). Some of them may harbor intermediate mass black holes. However, to our knowl- edge, no evidence of a collimated outflow has been found from the GC. There may be aligned, VLBI-scale ( ∼1 pc in the GC) radio knots that have not been discovered. Models of Past Enhanced Sgr A∗Accretion —Totani (2006) suggested a RIAF model of Sgr A∗to explain both the past high X-ray luminosity and the kinetic luminosity of the outflows inferred from observations. The required boost factor of the accretion rate ∼103−4with a time scale of∼107yr is naturally expected in the model. The induced outflow is energetic enough to support the hot (∼8 keV) plasma halo towards the GC. The sudden de- struction of the accretion flow of Sgr A∗is caused by the remnant Sgr A East passing through the MBH in the past ∼102yr. Such a model is claimed by Totani (2006) to better explain the positron production which is required to generate the 511 keV line emission to- wards the Galactic bulge observed by INTEGRAL/SPI (Weidenspointner et al. 2008). Galactic Jets in Starburst Galaxies —In other galaxies, powerful AGN radio jets interacting strongly with the hot gas have been observed (McNamara & Nulsen 2007). In relatively radio-quiet galaxies such as NGC 4636, NGC 708, and NGC 4472, (see e.g. Wang & Hu 2005, for a summary), faint X-ray “ghost” cavities appear without corresponding radio lobes. It has been suggested that capture of red giant stars and accretion onto the central MBHcan powerjets/outflowswith typical energy ∼1056 ergs, which is the required energy to form the observed cavities. Interferometric monitoring of the Seyfert galaxy NGC 3079 has found evidence of radio jet components un- dergoing compression by collision with the clumpy ISM, within a few pc of the central engine. This result sup- ports the idea that the kpc-scale superbubble originates from the spread of the momentum of jets impeded from propagating freely. The generalization of this scenario provides an explanation for why jets in Seyfert galaxies are not able to propagate to scales of kpc as do jets in radio-loud AGN (Middelberg et al. 2007). Precessing Jets —Several processes may lead to jet precession, including magnetic torques, warped discs, and gravitational torques in a binary system (see e.g. Falceta-Gon¸ calves et al. 2010). If the momenta of the MBH and the accretion disk are not perfectly aligned, 36 the Bardeen-Petterson effect could be a likely mecha- nism for precession (Bardeen & Petterson 1975), which will force the alignment ofthe disk and the MBH angular momentum. Irwin & Seaquist (1988) suggested that the “bubble” structures that have been seen in nearby starburst galax- ies couldhavebeen blownby a precessing VLBI-scalejet, a flat-spectrum radio core was recognized in NGC 3079 (Baan & Irwin 1995; Irwin & Sofue 1992). Typically, narrow, relativistic and non-precessing jets can carve out a hot gas bubble by interaction with ISM and release most of their energy far from the GC. Wide jets with large opening angles are capable of transferring momen- tumintoalargerarearesultingintheinflationoffatbub- bles (also for precessing AGN jets see Sternberg & Soker 2008). 6.2.2.Shocked Shells Driven by AGN Jet Relativistic jets dissipate their kinetic energy via in- teractions with the ISM. Jets in radio-loud AGN can in- flate a bubble composed of decelerated jet matter which is often referred to as a cocoon. Initially, the cocoon is highly overpressured against the ambient ISM and a strong shock is driven into the ambient matter. Then a thin shell isformed aroundthe cocoonbythe compressed medium. The shells are expected to be a promising site for particle acceleration. As a simple estimate of the dynamics of the expand- ing cocoon and shell, we assume the bubbles and shells are spherical, and also assume that the ambient mass density profile has a form of a power-law given by ρa(r) =ρ0(r/1kpc)−1.5, whereρ0is the mass density at r= 1 kpc. We further assume that the kinetic power of the jet,Lj, is constantin time. Under these assumptions, the dynamics can be approximately described based on the model of stellar wind bubbles. The radius of the shock is given by R(t)∼6ρ−2/7 0.01L2/7 42t6/7 7kpc, where ρ0.01=ρ0/0.01mpcm−3,L42=Lj/1042ergs s−1and t7=t/107yr. Taking the expected numbers for the Fermibubbles, we get a approximate estimate of the bubble size. The total internal energy stored in the shell can be expressed as Es∼0.1Ljt, implying that roughly 10% of the total energy released by the jet is deposited in the shell. 6.2.3.Buoyant Bubbles X-ray images have revealed shock fronts and giant cav- ities, some with bipolar structure, in the inner regions of clusters, surrounded with X-ray emitting gas. It is believed that the power in radio jets is comparable to the energy required to suppress cooling in giant ellipti- cal galaxies or even rich clusters (McNamara & Nulsen 2007). The depressions in the X-ray surface brightness of the ROSAT map may themselves indicate the presence of emptycavitiesorbubblesembeddedinhotgas,andcouldbe interpreted as a signature of previous AGN feedback with hot outflows. For adiabatic or supersonic bubbles first inflated by AGN jets, once the bubble reaches pres- sure equilibrium with the surrounding gas, it becomes buoyant and rises, because the mass density is lower in the bubble than in its surroundings. As the bubble movesawayfrom the GC, towardregionswith even lower density and pressure, it expands. The velocity at which the bubbles rise depends on the buoyancy and the drag forces. The ISM is in turn pushed by the rising bubble, which causes a upward displacement behind the bubble called “drift”. This trailing fluid can give rise to filaments of cool gas. There are indeed filamentary structures in the inner Galaxy in both H Iand Hαmaps. Their identifica- tion could support the buoyant bubble scenario. 6.3.Nuclear Starburst Another possible source of dramatic energy injection is a powerful starburst in the nucleus. Starburst induced Galactic winds are driven by the energy released by su- pernova explosions and stellar winds following an intense episode of starformation, which create an over-pressured cavity of hot gas. The galactic wind fluid is expected to have an initial temperature within the starburst re- gion in the range of 107−8K even if it has been lightly mass loaded with cold ambient gas (Chevalier & Clegg 1985). The ISM can be swept up by the mechanical en- ergy of multiple SN explosions and stellar winds. Large- scale galactic outflows have been observed in starburst galaxies both in the local Universe and at high redshifts (Veilleux et al.2005; Bland-Hawthorn et al. 2007). Star- burst episodes near the GC have been discussed in (e.g. Hartmann 1995). 6.3.1.Morphology of Outflows in Starburst Galaxies Starburst-driven galactic winds have been studied ex- tensively in both multi-waveband observations and hy- drodynamical simulations (Strickland & Stevens 2000; Veilleux et al. 2005). AGN or starburst galaxies show bipolar outflows (Gallimore et al. 2006; Sharp & Bland-Hawthorn 2010). The total energy of the superwind has been estimated as∼1055−56ergs, comparable to the estimated energy of theFermibubbles. The SpitzerSpace Telescope has found a shell-like, bipolar structure in Centaurus A , 500 pc to the north and south of the nucleus, in the mid- infrared (Quillen et al. 2006). The shell has been esti- mated to be a few million years old and its mechanical energy of 1053−55erg depends on the expanding veloc- ity. A small, few-thousand solar mass nuclear burst of star formation, or an AGN origin has been proposed to explain the formation of the shell. Recently, Westmoquette et al. (2009) showed that ion- ized gas in the starburst core of M82 is dynamically com- plexwithmanyoverlappingexpandingstructureslocated 37 at different radii, with compressed, cool, photo-ionized gas at the roots of the superwind outflow. Extra-planar warmH2knots and filaments extending more than ∼3 kpc above and below the galactic plane have also been found (e.g. Veilleux et al. 2009) NGC 253 is one of the most famous nearby starburst galaxiesandissimilartoourMilkyWayinitsoverallstar formation rate, except for a starburst region toward the center of the galaxy with spatial extent of a few hun- dred pc. A galactic wind in NGC 253 was found in Hα(McCarthy et al.1987). Strickland & Stevens(2000) discussed the spatial structure in detail in X-ray. The wind reaches out to ∼9 kpc perpendicular to the disk. Filamentary structures as part of projected conical out- flow are found in H αand near-infrared H2emission. The relatively warm gas ( ∼104K) exists close to the hot gas (∼106K). The X-ray filaments tend to be located in in- ner regions compared to H α, and are brighter where the Hαemissions are locally weak. The separation between the Hαand X-rayfilaments is ∼70pc. The spatial distri- butions of H αand X-ray indicate that the inner Galactic wind has higher temperature than the outer part. The UVemissionseemstoformashellaroundtheX-rayemis- sion (Bauer et al. 2008). VLBI and VLA observations of thenuclearregionofNGC253at22GHzshowsnodetec- tion of any compact continuum source on milliarcsecond scales, indicating no low-luminosity AGN in the central region of NGC 253. It seems that the starburst region is the most plausible explanation for the source powering the wind (Brunthaler et al. 2009). 6.3.2.The Mechanism of Galactic Winds The energy and momentum transfer of the Galactic wind could be dominated by high thermal and/or ram pressure. Materials have been swept-up and entrained as part of the wind, and the wind fluid comprises merged SN ejecta and massive star stellar wind material with ambient gas from the starburst region. Two popular ex- trinsic feedback mechanisms have been suggested: ther- mally driven winds powered by core-collapse supernovae and momentum-driven winds powered by starburst radi- ation (Cox 2005). On the other hand, the idea that CRs and magnetic fields can help to drive galacticwinds has been known for decades (Ipavich 1975; Breitschwerdt et al. 1991). For every core-collapse SN, about 10 percent of the energy release is converted into CRs ( ∼1050erg). These CRs interact with the magnetized ISM extensively and ex- change momentum through Alfvenic disturbances: the characteristic mean free path for a CR proton in the starburst phase of the GC is ∼1 pc. The effective cross section for CR protons and nuclei interacting with the ambient gas is much higher than the Thomson cross sec- tion for electrons. The luminosity at which CR collisions with gas balance gravity is about 10−6of the usual Ed- dington luminosity (Socrates et al. 2008).Momentum wind outflowing galactic supershells can also be driven by Ly αradiation pressure around star- forming regions (Dijkstra & Loeb 2009). The supershell velocity can be accelerated to 102−3km/s, and it may evenbeabletoescapefromthehostgalaxy. Theradiiare predictedtobe rsh= 0.1−10kpc, withages tsh= 1−100 Myr and energies Esh= 1053−55erg. However, the morphology of the galactic wind in nearby starburst galaxies inferred from synchrotron, H α, and H Imaps is asymmetric about both the galactic mi- nor axis and galactic plane, which may suggest the inho- mogeneous nature of the ISM. On the other hand, as we have shown, the north and south Fermibubbles are ap- proximatelysymmetric with respect to the galacticplane and the minor axis of the disk. The symmetric structure of theFermibubbles might indicate that they are not generated by subsequent interactions with ambient gas throughout the wind. Furthermore, the typical speed of galactic winds is about 200 −300 km/s. It takes about 5 ×107yr for CR electrons to reach 10 kpc, but we have not seen any evidence of cooling in the gamma ray intensity and spectrum. We probably need a faster transport mecha- nism of CRs if they were generated from the GC. How- ever, the higher the velocity of the wind which entrains the CRs, the greater the kinetic energy the wind con- tains. The estimated energy of the Fermibubbles only includes thermal energy; if they are actually kinetic en- ergydominated, then the energyrequirementto form the Fermibubbles is even larger (Gebauer & de Boer 2009; Jokipii & Morfill 1987; Lerche & Schlickeiser 1982). The estimated supernova rate in the NGC 253 star- burst region is about 0.1/yr (Engelbracht et al. 1998). Assumingeachsupernovaexplosionreleases1051erg,and 10 percent of the energy is transferred to heating the Fermibubbles, this gives a rate of energy injection com- parable to ∼1042erg/s. The star formation activity in NGC 253 has been underway for about 20 −30 Myr, and isconsideredtobe inasteadystatefortheCR transport, presumably with a smaller time scale (Engelbracht et al. 1998). It is possible that although the center of our Milky Way is currently in its quiet phase, it was recently (in the past 107yr) in a starburst phase similar to NGC 253 Heesen et al. (2009). To our knowledge, however, there is no evidence ofmassive supernova explosions( ∼104−5) in the past ∼107yr towards the GC, and no apparent Galacticwind featureshavebeen found in H α, indicating no strong recent ( ∼104yr) star formation activity. Everett et al. (2008, 2010) compared the synchrotron and soft X-ray emission from large-scale galactic wind models to ROSAT and Haslam 408 MHz maps. They show that a CR and thermally-driven wind could con- sistently fit the observations and constrain the launch- ing conditions of the wind, including the launching re- 38 gion and magnetic field strength. The comparison of the gamma-ray prediction of the wind model with Fermi- LAT diffuse emission, especially to the Fermibubbles might have important implications for the CR driven wind of our Milky Way. 6.3.3.Cosmic Rays from a Starburst A central starburst might also generate the increased population of electron CRs in the GeV-TeV energy range required to produce the gamma-raybubble signals. Star- burst galaxieshost a greatlyincreased rate ofsupernovae in the central region. The shocks from supernovae can merge and produce energetic galactic scale winds, and the enhanced population of supernova remnants is be- lieved to accelerate CRs, resulting in orders of magni- tude higher CR density than currently expected in the GC. It is likely that the GC exhibited comparable CR density in the past, with a starburst phase turned on by boosted formation of massive stars. As previously discussed, during this starburst phase the CR protons produced in the inner Galaxy scatter on the ISM with very short path lengths, producing gamma rays (via π0production and decay), and electron (and positron) CRs. Although the immediately-produced gamma-rays and high density gas and ISRF associated with the star- burstphasewouldnotbeobservabletoday,thesecondary electrons might be leftover from the past active star- burst phase, and could have been transported to ∼10 kpcbythemagneticfieldentrainedintheGalacticwinds. GeV and TeV gamma-rays have recently been detected in nearby starburst galaxies NGC 253 (Acero et al. 2009), M 82, and Large Magellanic Cloud (LMC), by Fermi-LAT (Abdo et al. 2010a), High Energy Stereo- scopic System (H.E.S.S.) (Itoh et al. 2007) and VERI- TAS (Karlsson & for the VERITAS collaboration 2009), support the starburstgalaxiesas a richsource ofhigh en- ergy gamma-rays. If a transient starburst did occur in the GC of our Milky Way, and produced a large popula- tion of CRs responsible for the observed Fermibubbles, andWMAP haze, what triggered and terminated the starburst phase is unclear. If the CRs are driven by winds, the halo magnetic field can carry CRs along the field lines from the inner disk/bulge into the halo, which could help the CR elec- trons to reach 10kpc without significant diffusive soften- ing of their spectrum. The vertical CR bulk velocity is typicallyhundredsofkm/s,whichisremarkablyconstant over the entire extent of the disk and for galaxies with different mass. In the standard picture, CRs can not stream faster than the Alfv´ en speed with respect to the static frame of the magnetic field, due to the well-known streaming instability (Kulsrud & Pearce 1969). Consid- ering nearby starburst galaxies, the vertical CR bulk speed has been measuredto be vCR= 300±30kms−1for NGC 253. For a typical magnetic field strength of B≈ 15µG and a density of the warm gas of n≈0.05cm−3leads to an Alfv´ en speed of vA=B/√4πρ≈150kms−1 (Heesen et al. 2009). The super-Alfv´ enic CR bulk speed requires that the CRs and the magnetic field which is frozen into the thermal gas of the wind are advectively transported together. The measured CR bulk speed is the superposition υCR=υwind+υA. In this picture, the spectrum of electrons injected in the starburst region might be harder than elsewhere in the Galaxy if the transient starburst phase led to a top- heavy mass function of stars. Magnetic turbulence in the GC could also be much stronger than in the Solar neighborhood: the stronger the turbulence, the faster CRs are transported, and the higher the transition en- ergy at which synchrotron/ICS losses overtake diffusive losses. It has been shown that the high star formation rate per area in the GC leads to short transport times (Becker et al. 2009). 6.3.4.Constraints on Recent Starburst Activity If theFermibubbles were generated by previous star- burst activity in the GC, we would expect to see many more supernova remnants towards the GC than have been discovered. Moreover, radioactive26Al (half-life ∼7.2×105yr) is believed to be mainly produced by massive stars, supernovae, and novae in the Galaxy (Prantzos & Diehl 1996).26Al decays into26Mg, emit- ting a gamma-ray photon of 1808.65 keV. Observations of the spectrometer (SPI) on the INTEGRAL gamma- ray observatory seem to disfavor the starburst scenario: the Galaxy-wide core-collapse supernova rate has been estimated at 1 .9±1.1 per century from the flux of26Al gamma-rays(Diehl et al. 2006). Within its half-life,26Al can only travel ∼0.1 kpc with the typical Galactic out- flow velocity of ∼100 km/s. One expects strong26Al gamma-ray emission concentrated towards the GC, with a flux comparable with the total gamma-ray flux from the disk, if the outflow was produced by a starburst (Totani 2006). However, such a strong concentration at theGalacticcenterisnotfound(Prantzos & Diehl1996), indicating that the accretion activity of Sgr A∗is more plausible as the origin of the mass outflow. Future ob- servations of the ratio of the Galactic60Fe to26Al may provide better constraints on the starburst scenario. 6.4.Other Ways to Generate the Bubbles The molecular loops in the GC could possibly be ex- plained in terms of a buoyant rise of magnetic loops due to the Parker instability. For a differentially rotating, magneticallyturbulentdisk,suchmagneticloopscaneas- ily be formed and rise out of the disk. The typical scale of such a loop is 1kpc (Machida et al. 2009). 7.THE ORIGIN OF THE COSMIC RAYS It is not necessary that the physical mechanism that creates the bubbles also injects the electron CRs respon- siblefor the Fermibubbles. It is possiblethat the bubble 39 structures were formed earlier and the electron CRs were injected by an alternative mechanism that then lights up the bubble structure with gamma-ray emission. In this section, we would like to separate the CR production from the bubble formation, and address the spatial ori- gin of the CRs. In any case, the production mechanism should gener- ate electron CRs inside the Fermibubbles, and also pre- vent them from efficiently leaving the bubbles, in order to produce the observed “sharp edge”. However, elec- trons with ∼102GeV diffuse on the order of 1 kpc before losing half their energy (McQuinn & Zaldarriaga 2010). Higher energy electrons lose energy more rapidly and so have shorterpath lengths; if the gamma-rayemis- sion fromthe bubbles is dominated by IC scattering from TeV-scale electrons injected inside the bubble, then the sharp edge of the bubbles may be natural. This would in turn imply that the gamma rays observed by Fermiare largely upscattered CMBphotons (starlight and far in- frared photons are upscattered to much higher energies), which is advantageous for generating such a latitudinally extended IC signal. If instead ∼100 GeV electrons scat- tering on starlight are primarily responsible, the electron CR density must increase markedly at high Galactic lat- itudes to compensate for the falling of the starlight den- sity to higher latitude. However, generating a very hard CR electron spectrum extending up to O(TeV) energy may be challenging for conventional CR production and acceleration mechanisms. These difficulties may be ame- liorated if the hard gamma-ray signal from the Fermi bubbles is a transient rather than a steady-state solu- tion. 7.1.CRs from the Galactic Center As discussed in §6, the electron CRs could be pro- duced in the inner Galaxy by mechanisms such as OB associations (Higdon & Lingenfelter 2005), accre- tion events, and supernova explosions, then entrained in subsequent jets or outflows and rapidly carried up to large scales, avoiding diffusive softening of the spectrum. Breitschwerdt & Schmutzler (1994) suggested that the all-sky soft X-ray emission can be explained by delayed recombination in a large-scale CR and thermal-gas pres- sure driven wind. Such a wind model has been applied to theMilkyWaywhichcouldexplaintheobservedGalactic diffuse soft X-ray emission and synchrotron. The model indicates that the Milky Way may possess a kpc-scale wind. The cooling time (denoted τ) of TeV-scale electron CRs can impose stringent constraints on such models. In the Thomson regime, the cooling time scales as (electron energy)−1, leading to estimates of τ∼105yearsfor TeV- scale electrons. However, where scattering on starlight photonsisimportant, electronenergiesof100 −1000GeV are no longer in the Thomson regime, so the scattering cross section is suppressed and the cooling time can belonger than naively expected. Figure28shows the variation in cooling time (defined as 1/(dlnE/dt)) as a function of electron energy and height above (or below) the Galactic plane, for the stan- dard ISRF model, both with and without the inclusion of synchrotron losses (for a simple exponential magnetic field profile)16. For example, at z= 2 kpc, the cool- ing time for TeV electrons is ∼3−4×105years, rising to∼7×105years at z= 5 kpc. Scatterings purely on the CMB give an upper bound on the lifetime of ∼5×106(1TeV/E) years, but even several kpc from the GC, scatterings on the infrared photons dominate at TeV-scale energies. The effect of the KN cross-section suppressionat higher electron energies can be seen in the small-zlimit where synchrotron losses are neglected, so IC scatteringofthe electrons onstarlightis an important contribution to the total energy loss rate. These relatively short lifetimes, especially close to the Galactic plane, may lead to severe difficulties in prop- agating CR electrons from the GC out to fill the bub- bles. Propagation over such large distances may also lead to significant diffusive softening ofthe electron spec- trum, whichmustbe reconciledwiththe apparentspatial uniformity of the bubbles’ (gamma-ray) spectral index. With electron injection primarily at the GC there is also no obvious natural explanation for the flat projected in- tensity profile, which seems to require sharp increases in the CR density at the bubble walls. 7.2.CRs from the Bubble Edge If the majority of the electron CRs are produced from shock acceleration within the edge of the Fermibubbles, the electron CRs in the bubble interior might be left- over CRs which undergo cooling after the shock passes through. The CRs continue to diffuse inward from the shock front while also diffusing outward; if the shock is movingfaster than the electronsdiffuse out, a sharpedge in the resulting Fermibubble gamma-rays is still ex- pected. Itisalsopossiblethat theCRs maybe secondary electrons, produced by enhanced proton-ISM interaction in shocks (within the bubble shell), where protons could be ejected from the GC and entrained in the shocks with high gas density due to shock compression. We can estimate the diffusion path length of100 −1000 GeV electrons, given the lifetimes calculated in §7.1. We use the estimated diffusion constant from GALPROP,K= 5.3×1028cm2/s at a reference rigidity of 4 GV, and take the diffusion coefficient index to be δ= 0.43 following the results of Ahn et al. (2008). Then the path length is given by, √ Kτ≈1.4/parenleftbiggE 1TeV/parenrightbigg0.43/2/parenleftbiggτ 106yr/parenrightbigg0.5 kpc.(4) 16However, note that we have treated the target photon distri- bution as isotropic. 40 Fig. 27.— A cartoon picture to summarize the observations of the Fermibubble structures. Two blue bubbles symmetric to the Galactic disk indicate the geometry of the gamma-ray bubble s observed by the Fermi-LAT. Morphologically, we see corresponding features inROSAT soft X-ray maps, shown as green arcs embracing the bubbles. T heWMAP haze shares the same edges as the Fermibubbles (the pink egg inside the blue bubbles) with smaller extensio n in latitude. These related structures may have the same phy sical origin: past AGN activities or a nuclear starburst in the GC (the yellow st ar). Thus we expect the diffusion scale to be small relative to the bubble size, although not negligible. Consequently, the electrons in the interior ofthe Fermi bubbles are unlikely to maintain a hard spectrum due to diffusion inward from the bubble walls. In this scenario, one needs to tune the electron CR distribution to get near-flatprojected intensity. Although the Fermibubble gamma-rays along any line of sight include contributions from both the bubble interior and bubble shell (due to the integration along the line of sight), the flat intensity and consistent spectrum of the inner bubble and outer shelltemplates (Figure 16)implies that theCR spectrum in the 3D bubble interior cannot be very different from the spectrum at the 3D bubble shell with electron CRs generated in situ. The edge of the Fermibubbles might contain MHD turbulent fluid with compressed gas and magnetic field. Fast magnetic reconnection rather than shocks might drive the CR acceleration. Significant magnetic reversals within the bubble shell are naturally expected, just like the heliosheath region of the solar system (Drake et al. 2006; Innes et al. 1997): the crossing of the termina- tion shock by Voyager 1 and 2 may indicate the ac- celeration within regions of fast magnetic reconnection (Lazarian & Opher2009). Itiswellknownthatmagneticfield reversal can cause magnetic reconnection. When two magnetic flux tubes of different directions become too close, rearrangement of the magnetic field lines takes place, converting the energy of the magnetic field into energy of plasma motion and heating. Such phenom- ena have been investigated extensively in e.g. the solar sphere and gamma-ray bursts. In the reconnection re- gion, the energetic particles bounce between two mag- netic tubes, undergoing first-order Fermi acceleration as the dominant acceleration process (Beck et al. 1996; Biskamp 1986; Elmegreen & Scalo 2004). On the other hand, maintaining the shape of the bub- ble, and preventing it from breaking out, is a non-trivial process. As the bubble rises through the ISM, it tends to fatten simply because the fluid moves faster on its sides than its front. The classical Rayleigh-Taylor instability appears at the leading edge of the bubble, as the inside density is lower than that of the ISM; Kelvin-Helmholtz instability occurs at the sides of the bubble due to dis- continuities in the velocity and density. Although hydro- dynamical processes may be capable of stabilizing the bubble structures against these instabilities, it has been suggestedthatalarge-scalecoherentmagneticfieldcould help to prevent disruption. We will discuss in §7.3the possibility that the magnetic field can be coherent on the 41 012345 z (kpc)104105106107108τ (years) With synchrotronICS only1 TeV100 GeV10 GeV Fig. 28.— The cooling time for electron CRs as a function of energy and h eight above the Galactic plane, for r= 0.Thin black lines:Synchrotron losses are neglected; equivalently, the B-fiel d is assumed to be negligible everywhere. Thick red lines: The cooling time calculation includes synchrotron losses in a magnetic field given by |B|= 30e−z/2kpcµG. We use the standard radiation field model from GALPROP version 50p, and define the cooling time τ= 1/(dlnE/dt). Fermibubble scale. 7.3.CR from Diffuse Production in the Bubble The hard, spatially uniform spectrum of the Fermi bubbles motivates the possibility of some diffuse injec- tion of hard CR electrons throughout the bubble volume. Such a mechanism, if present, would solve the issues of short electron cooling times, relative to the propagation time from the GC, and the hardnessof the requiredspec- trum. A uniform diffuse injection would give rise to a centrally brightened signal, so in models of this type there would most likely need to be some other mecha- nism increasing the electron intensity at the bubble wall, perhaps associated with a shock there. This requirement for multiple mechanisms may seem unwieldy, and perhaps also unnecessary, if diffusion or coolingofelectronsproducedattheshockcanexplainthe flat projected intensity profile and hard spectrum across the entire bubble. Given the large size of the bubble, however, it is unclear whether production or acceleration of electron CRs solely at the bubble walls can give rise to a sufficiently centrally bright and hard signal. Decay or annihilation of dark matter has previously been proposed as a mechanism for diffuse injection of very hard electrons and positrons; in particular, dark matterannihilationhasprovidedgoodfitstothedata(atleast in the inner ∼15−20◦) in previous studies of the WMAP haze (Hooper et al. 2007; Cholis et al. 2009a,b; Harding & Abazajian 2010; McQuinn & Zaldarriaga 2010; Lin et al. 2010; Linden et al. 2010). The IC signal from electrons produced in dark matter annihilation would naively be expected to have approximate radial symmetry about the Galactic center and fall off sharply at increasing randz, roughly tracing the distribution of dark matter density squared, but smoothed and broadened by diffusion of the electrons. This expecta- tion is in conflict with the bubble morphology, with a gamma-ray distribution elongated perpendicular to the Galactic disk, and extending to 10 kpc without much change in intensity. However, this naive picture is based on isotropic diffusion with a spatially uniform diffusion constant, in a tangled magnetic field; including the effects of an ordered magnetic field and anisotropic dif- fusion can lead to a much more bubble-like morphology (Dobler 2010). Another way to produce CRs in situinside the bub- ble is through magnetic reconnection. Electrons could be accelerated directly, or produced as secondaries from accelerated protons. The magnetic fields in the under- dense bubbles, which may be inflated by AGN outflow, may relax into an equilibrium filling the entire volume of the bubbles. The timescale depends on the magnetiza- 42 tion and helicity of the outflow and also the properties of the ISM. The magnetic field could undergo reconnec- tion on a short timescale, converting magnetic energy into heat. This mechanism can explain how the bubbles could move a large distance through the ISM without breaking up. The reconnection in the bubbles can also accelerateenergetic particles, circumventingthe problem ofsynchrotronemittershavingashorterlifetimethanthe age of the bubble they inhabit (Braithwaite 2010). The Fermibubbles might be initially highly turbulent, with a disordered magnetic field far from equilibrium (gas pres- sure and Lorentz force are not balanced); the kinetic en- ergy would then be dissipated by viscosity, especially in the low density bubble. The time scale of the relaxation to equilibrium can be estimated (Braithwaite 2010) by, τrelax≈r αvA=r√4πρ αB=r α/parenleftbigg4πr3ρ 6E/parenrightbigg1 2 (5) ≈4.4×106/parenleftBigα 0.1/parenrightBig−1/parenleftbiggr 4kpc/parenrightbigg/parenleftBigg ρ 10−5mpcm−3/parenrightBigg1 2 /parenleftbiggB 5µG/parenrightbigg−1 yr, (6) wherevA,BandEare the Alfv´ en speed, magnetic field and energy at equilibrium . We use the estimated char- acteristic values for the Fermibubbles. The resulting time scale is on the order of 106−7yr. Furthermore, the reconnection time scale is orders of magnitude shorter than the relaxation time scale: we expect to see ongoing reconnectionifequilibrium hasnot alreadybeen reached. 8.POTENTIAL CONNECTIONS TO OTHER OPEN QUESTIONS 8.1.The Cosmic Gamma-ray Background Measurement of the intensity and spectrum of cosmic gamma-ray background (Abdo et al. 2010b; Fields et al. 2010) has suggested that instead of rare but intrinsically bright active galaxies, it is numerous but individually faint normalgalaxiesthat comprisethe bulk ofthe Fermi gamma-ray background. This result infers a tighter cor- relation between cosmic star formation history and the cosmic gamma-ray background (Thompson et al. 2007; Lacki et al. 2010). Galactic outflows have been identi- fied from near infrared observations at redshift z∼2 (Alexander et al. 2010) indicating that such outflows are common features of ultraluminous infrared galaxies. Such outflowscanentrainenergeticcosmicraysto escape the galaxies, thus providing a way to contribute to the cosmic gamma-ray background at energies in the Fermi range. The Fermibubbles may provide local evidence for such a scenario. 8.2.The Origin of Hypervelocity Stars Recent surveys of hypervelocity stars (HVSs) have found 16 HVSs which are mostly B-type stars with ages>∼(1−2)×108yr (see e.g. Brown et al. 2009). These HVSs distributed in the Galactic halo are believed to originatefromtheGCinvolvinginteractionsofstarswith the MBH (Hills 1988; Yu & Tremaine 2003). Thus en- ergetic energy release from the GC which generate the Fermibubbles could be dynamically related to the ejec- tion of HVSs. Recently Lu et al. (2010) have shown that the spatial distribution of the discovered HVSs is consis- tent with being located on two thin disk planes. The ori- entationofthe planes areconsistentwith the inner/outer clockwise-rotating young stellar (CWS) disk. One pos- sibility could be that the HVSs originate from some un- knownandpreviouslyexistingdisk-likestellarstructures. HVSs might have been ejected periodically, and related accretion events produce jets which generate the Fermi bubbles. 8.3.The Future of the FermiBubbles What is the future of the Fermibubbles, are they in a “breakout” stage? An interesting possibility is that thenorthern arc and even part of the Loop Ifeature are parts of the relics of previous bubbles, and the bubble production is a periodic process. The bubbles might be fast expanding shocks which might finally expand freely into the galactic halo, thus contaminating the ISM with entrained hot gas and CRs. An intriguing possibility is that CRs and gas released from previous such bubbles to the Galactic halo may contribute to the observed dif- fuse X-ray and gamma-ray background. In any case, the study of Fermibubbles would have potential implica- tions for the understanding of the feedback mechanism from the Galaxy. 8.4.Missing Baryons and High-Velocity Clouds N-body/gas dynamical galaxy formation simulations have shown that for Milky Way like galaxies, about 70% extra baryonic mass should reside around the galaxy in form of hot gas (Sommer-Larsen 2006). Warm clouds are confined by the pressure of the ambient hot halo gas, which contains mass at least two orders of magnitude more than these warm clouds. The study of the Fermi bubbles may provide hints of hot gas feedback from the Galaxy, the search for the missing baryons (Bregman 2007) and the puzzle of high-velocity clouds. 8.5.Metallicity Although the Galactic outflow can not inject a large fraction of the ISM, a significant amount of the freshly produced metals could be channeled along the galac- tic wind. In ordinary photo-ionization it is difficult to make the [N II]/Hαratio exceed about 1.0; shock ioniza- tion/excitation is plausible once [N II]/Hαis detected. [NII]/Hαhas been estimated ∼1.5 in NGC 253, less than 1.0 in M82 where the ratio trends to increase far from the disk. 43 The presence of metals in the IGM has been inter- preted as the consequence of energetic metal-rich out- flows from galaxies. Active star formation in the in- nerGalaxymaycontaminatesthe surroundingISM, with (periodic)Galacticwindsentrainingthemetal-richgasto tens of kpc. The Fermibubbles may give some hints to understanding the feedback and metallicity of the IGM. Jets from GC in general do not imply a high metallic- ity, and detections of metal rich outflows may essentially constrain the energetic injection from jets or Galactic outflows from previous starburst towards the GC. 9.CONCLUSIONS We have identified two large gamma-ray bubbles at 1/lessorsimilarE/lessorsimilar50GeV in Fermimaps containing 1.6 years of data. They have approximately uniform surface bright- ness with sharp edges, neither limb brightened nor cen- trally brightened, and are nearly symmetric about the Galactic plane. The bubbles extend to 50◦above and below the Galactic center, with a maximum width of about 40◦in longitude. At |b|/lessorsimilar30◦, these “Fermibub- bles” have a spatial morphology similar to the WMAP microwavehaze (Finkbeiner 2004a; Dobler & Finkbeiner 2008). The ROSATsoft X-ray 1.5 keV map also reveals hard-spectrum features that align well with the edges of theFermibubbles. The similarities of the morphol- ogyand hard spectrum stronglysuggest that the WMAP haze and the Fermibubbles share a common origin. In contrast, the Fermibubble features are not aligned withLoop Ior any other feature in the Haslam 408 MHz map; while Loop Iand other shell structures appear in the gamma rays, their spectra are softer than the bubble spectrum. Furthermore, there are no convincing features spatially correlated with the bubbles in the LAB H Ior Hαmaps. To better reveal the bubble structures, we use spatial templatestoregressoutknownemissionmechanisms. To removeπ0and bremsstrahlung gammas we use the SFD dustmap, assumingthattheinterstellardusttracesgam- mas producedby cosmic rayprotons( π0decay) and elec- trons (bremsstrahlung) colliding with the ISM. To trace the gamma rays from IC scattering of disk electrons on the ISRF, we consider 3 different templates: a simple geometric disk model, the 408 MHz Haslam map, and the 0.5−1 GeVFermimap after removal of the (dust- correlated) π0gammas. We verify that our results are insensitive to this choice. As an additional cross check, we subtract the Fermidiffuse Galactic model from the data, finding that this also reveals the bubble structures. Although our template fitting technique is subject to sig- nificant uncertainties due to uncertain line of sight gas and CR distributions, these uncertainties mainly affect the intensity profile at low latitudes. For |b|>30◦and 1< E <50GeV the morphology and spectrum are com- pletely consistent for different template choices. Indeed, in the 1−5GeV maps, a significant part of the southernbubble is easily visible before anytemplate subtractions. TheFermibubbles have an energy spectrum of dN/dE∼E−2, significantly harder than other gamma- ray components. There is no apparent spatial variation in the spectrum between the bubble edge and interior, and the north and south bubbles have consistent spatial andspectralprofiles. Boththemorphologyandspectrum are consistent with the two bubbles having the same ori- gin and being the IC counterpart to the electrons which generatethe microwavehaze seen in WMAP(Finkbeiner 2004a; Dobler & Finkbeiner 2008). The spectrum of the CR electrons required to generate the Fermibubbles is harder than expected for electrons accelerated in super- novashocksinthe disk, andsuchdisk-producedelectrons would be even softer after several kpc diffusion; the mor- phology of the bubble structure is also quite different to that of the lower-energy electrons traced by the Haslam 408 MHz map. Even setting aside the WMAP haze, theFermibubbles are unlikely to originate from excess π0emission, as (by construction) they are spatially dis- tinct from the SFD dust map, their spectrum is much harderthan that of the dust-correlatedemission, and the ROSAT data suggest that the bubbles are hot and un- derdense rather than overdense. The morphology of the Fermibubbles, and the overlapofthe Fermibubbles and WMAPhaze, also strongly disfavor the hypothesis that a significant fraction of the high energy gamma rays ob- servedby Fermiinthebubbleregionarephotonsdirectly produced by dark matter annihilation. TheFermibubble structures were likely created by some large episode of energy injection in the GC, such as apastaccretioneventontothe centralMBH,oranuclear starburst in the last ∼10 Myr. We have discussed some general possibilities and considerations in this work, and found shortcomings in each scenario; it seems likely that either significant modifications to one of these ideas, or some combination of different mechanisms, will be nec- essary. Jets originating from AGN activity can potentially ac- celerate CR electrons to high energies, and transport them rapidly away from the GC; the cooling time of electrons at 100 −1000 GeV is only 105−6years, so if the CRs are injected and accelerated only in the GC, a very fast bulk transport mechanism is required to convey them throughout the bubbles before they lose a signifi- cant fraction of their energy. However, filling the bub- bles completely, with n= 10−2cm−3gas, would require a mass injection of ∼108solar masses, so in any case it is more reasonable for the bulk of the material in the bubbles to be swept up and accelerated as the bubbles expand. Energetic shocks associated with jets can have high Mach number and thus efficiently accelerate CR electrons, producing hard spectra with dN/dE∼E−2, and the total energy required to heat the bubbles is also readily achievable by accretion events onto the central 44 MBH. However,the north-southsymmetry of the bubbles has no obvious explanation in the context of an AGN jet: there is no reason for one jet to be oriented perpendicu- lar to the Galactic plane. The large width and rounded shape of the bubbles are also not typical of jets, which are generally much more collimated, although a precess- ing jet might help explain the wide opening angle of the bubbles. If the central MBH becomes active on a rela- tively short timescale, the Fermibubbles may be created by a number of past jets, which combine to give rise to the symmetric and uniform Fermibubbles. An alternate source for the large required energy injec- tionis anuclearstarburst. The wide openingangleofthe bubbles is not a problem in this case; the bubble shape is similar to that observedin NGC 3079, and the X-ray fea- tures observed by ROSAT are similar to those observed in other nearby starburst galaxies. However, no corre- sponding H αsignal of the Fermibubbles is observed, in contrast to other known starburst galaxies: this prob- lem might potentially be resolved if the H α-emitting gas has cooled in the time since the starburst phase (gas hot enough to emit the X-rays observed by ROSAT has a considerably longer cooling time). Also, generally gas filaments and clumps are observed in the X-rays in star- burst galaxies, and it would seem that a relic of a past starburst should become more clumpy with gas clouds and filaments due to cooling of the gas. However, while no such structures are obvious in the ROSAT maps, the signal-to-noise is insufficient to place strong constraints. The absence of any such filamentary structures inside theFermibubbles, on the other hand, argues against a hadronic origin for the bubble gamma ray emission. Hadronic jets might accelerate protons to high energies, and the interactions of these protons with the ISM could then produce hard π0gammas and secondary e+e−, whichwouldscatterontheISRFtoproducemoregamma rays. In this scenario, however, the gamma ray emission should trace the gas density, which we would not expect to be smooth and homogeneous. Returning to the starburst scenario, the cosmic ray ac- celeration in this case would be due to shocks at the edge ofthe bipolarwind. However,the shocksexpectedin this scenario would be relatively weak and slow-moving, and thus may not be capable of generating a sufficiently hard electron spectrum to reproduce the signal. For example, in first-order Fermi acceleration, a shock Mach number of∼3.3 is needed to obtain an electron spectral index of 2.4, as required for the synchrotron explanation of the WMAPhaze (see e.g. Jones & Ellison (1991)). TheFermibubbles have sharp edges, also suggesting the presence of a shock at the bubble walls. If the CRs producing the gamma rays have a multi-kpc diffusion length (which is not expected to be the case for 1 TeV electrons, for example), then the edges can still be sharpifthe bubbleedgeismovingoutwardfasterthantheycan diffuse. If we assume the Fermibubbles are projected structures from three dimensional symmetric blobs to- wards the GC, the flat intensity profile of the bubbles requires the emissivity to rise at the bubble walls, but remain non-negligible in the bubble interior; the lack of spatial variation in the spectral index may also constrain models where the electrons diffuse long distances from an injection point. Magnetic reconnection in the interior of the bubbles, or some other mechanism such as dark matter annihilation, may help maintain a hard spectrum throughout the bubbles by accelerating existing lower- energy electrons or injecting electrons in situ. Dark matter annihilation or decay, while an effective mechanism for injecting hard electron CRs at high lati- tudes, cannot produce the features in the ROSATX-ray maps correlated with the bubbles, and would not be ex- pected to result in sharp cutoffs in gamma-ray emission at the bubble edges. Dark matter annihilation or decay may be contributing to the bubbles, or to gamma-ray emission in the inner Galaxy that is not well subtracted by either the bubble structure template or the models for known diffuse emission mechanisms; however, under- standing the Fermibubbles will be a necessary step be- fore extracting any such dark matter signal. TheeROSITA17andPlanck18experiments will pro- vide improved measurements of the X-rays and mi- crowaves, respectively, associated with the Fermibub- bles, and so may help discriminate between these scenar- ios.eROSITA , which is expected to launch in 2012, will provide the first imaging all-sky survey of mid-energy X- rays, studying the 0 .2−12 keV energy range with ∼100 eVenergyresolutionandaPSFof20”. The Plancksatel- lite, launched in 2009, will greatly improve the measure- ments of the WMAPhaze spectrum. In addition, AMS- 0219will launch in 2011, and may significantly advance our understanding of CR acceleration and propagation, and help to refine our interpretation of the Fermibub- bles. Acknowledgments: We acknowledge helpful conversa- tions with Avery Broderick, Greg Dobler, Martin Elvis, Jim Gunn, Jill Knapp, Maxim Markivitch, David Mer- ritt, Simona Murgia, Norm Murray, Paul Nulsen, Aneta Siemiginowska, David Spergel, Anatoly Spitkovsky, and Neal Weiner. DF and TS are partially supported by NASA grant NNX10AD85G. TS is partially supported by a Sir Keith Murdoch Fellowship from the American Australian Association. This research made use of the NASA Astrophysics Data System (ADS) and the IDL Astronomy User’s Library at Goddard.20 17http://www.mpe.mpg.de/heg/www/Projects/EROSITA/main .html 18http://www.rssd.esa.int/index.php?project=Planck 19http://www.ams02.org 20Available at http://idlastro.gsfc.nasa.gov 45 REFERENCES Abdo, A. A., et al. 2009, Physical Review Letters, 102, 18110 1 —. 2010a, ApJ, 709, L152 —. 2010b, ApJ, 720, 435 Abdo, A. A., et al. 2010, Phys. Rev. Lett., 104, 101101 Acero, F., et al. 2009, Science, 326, 1080 Aharonian, F., et al. 2009, A&A, 508, 561 Ahn, H. S., et al. 2008, Astropart. Phys., 30, 133 Alexander, D. M., Swinbank, A. M., Smail, I., McDermid, R., & Nesvadba, N. P. H. 2010, MNRAS, 402, 2211 Alvarez-Mu˜ niz, J., Engel, R., & Stanev, T. 2002, ApJ, 572, 1 85 Atwood, W. B., et al. 2009, ApJ, 697, 1071 Baan, W. A., & Irwin, J. A. 1995, ApJ, 446, 602 Baganoff, F. K., et al. 2003, ApJ, 591, 891 Bardeen, J. M., & Petterson, J. A. 1975, ApJ, 195, L65+ Bartko, H., et al. 2009, ApJ, 697, 1741 Bauer, M., Pietsch, W., Trinchieri,G., Breitschwerdt, D., Ehle, M., Freyberg, M. J., & Read, A. M. 2008, A&A, 489, 1029 Beck, R., Brandenburg, A., Moss, D., Shukurov, A., & Sokoloff , D. 1996, ARA&A, 34, 155 Becker, J. K., Biermann, P. L., Dreyer, J., & Kneiske, T. M. 20 09, ArXiv e-prints Bhat, C. L., Issa, M. R., Mayer, C. J., & Wolfendale, A. W. 1985 , Nature, 314, 515 Biermann, P. L., Becker, J. K., Caceres, G., Meli, A., Seo, E. , & Stanev, T. 2010, ApJ, 710, L53 Biskamp, D. 1986, Physics of Fluids, 29, 1520 Bland-Hawthorn, J., & Cohen, M. 2003a, ApJ, 582, 246 —. 2003b, ApJ, 582, 246 Bland-Hawthorn, J., Veilleux, S., & Cecil, G. 2007, Ap&SS, 3 11, 87 Blandford, R. D., & Begelman, M. C. 1999, MNRAS, 303, L1 Bloemen, H. 1989, ARA&A, 27, 469 Blumenthal, G. R., & Gould, R. J. 1970, Reviews of Modern Physics, 42, 237 Braithwaite, J. 2010, MNRAS, 406, 705 Bregman, J. N. 2007, ARA&A, 45, 221 Breitschwerdt, D., McKenzie, J. F., & Voelk, H. J. 1991, A&A, 245, 79 Breitschwerdt, D., & Schmutzler, T. 1994, Nature, 371, 774 Brown, W. R., Geller, M. J., Kenyon, S. J., & Bromley, B. C. 200 9, ApJ, 690, L69 Brunthaler, A., Castangia, P., Tarchi, A., Henkel, C., Reid , M. J., Falcke, H., & Menten, K. M. 2009, A&A, 497, 103 Casandjian, J., Grenier, I., & for the Fermi Large Area Teles cope Collaboration. 2009, ArXiv e-prints Chang, J., et al. 2008, Nature, 456, 362 Chevalier, R. A., & Clegg, A. W. 1985, Nature, 317, 44 Cholis, I., Dobler, G., Finkbeiner, D. P., Goodenough, L., S latyer, T. R., & Weiner, N. 2009, ArXiv e-prints Cholis, I., Dobler, G., Finkbeiner, D. P., Goodenough, L., & Weiner, N. 2009a, Phys. Rev., D80, 123518 Cholis, I., Goodenough, L., & Weiner, N. 2009b, Phys. Rev., D 79, 123505 Cox, D. P. 2005, ARA&A, 43, 337 Dame, T. M., Hartmann, D., & Thaddeus, P. 2001, ApJ, 547, 792 Diehl, R., et al. 2006, Nature, 439, 45 Dijkstra, M., & Loeb, A. 2009, MNRAS, 400, 1109 Dixon, D. D., Hartmann, D. H., Kolaczyk, E. D., Samimi,J., Di ehl, R.,Kanbach, G.,Mayer-Hasselwander,H.,&Strong,A.W.199 8, New Astronomy, 3, 539 Dobler, G. 2010, Private communication Dobler, G., & Finkbeiner, D. P. 2008, ApJ, 680, 1222 Dobler, G., Finkbeiner, D. P., Cholis, I., Slatyer, T., & Wei ner, N. 2010, ApJ, 717, 825 Dodds-Eden, K., et al. 2009, ApJ, 698, 676 Drake, J. F., Swisdak, M., Che, H., & Shay, M. A. 2006, Nature, 443, 553 Eckart, A., et al. 2004, A&A, 427, 1 Elmegreen, B. G., & Scalo, J. 2004, ARA&A, 42, 211 Engelbracht, C. W., Rieke, M. J., Rieke, G. H., Kelly, D. M., & Achtermann, J. M. 1998, ApJ, 505, 639 Everett, J. E., Schiller, Q. G., & Zweibel, E. G. 2010, ApJ, 71 1, 13 Everett, J. E., Zweibel, E. G., Benjamin, R. A., McCammon, D. , Rocks, L., & Gallagher, III, J. S. 2008, ApJ, 674, 258 Falceta-Gon¸ calves, D., Caproni, A., Abraham, Z., Teixeir a, D. M., & de Gouveia Dal Pino, E. M. 2010, ApJ, 713, L74 Falcke, H., & Markoff, S. 2000, A&A, 362, 113Fichtel, C. E., Hartman, R. C., Kniffen, D. A., Thompson, D. J. , Ogelman, H., Ozel, M. E., Tumer, T., & Bignami, G. F. 1975, ApJ, 198, 163 Fields, B. D., Pavlidou, V., & Prodanovi´ c, T. 2010, ApJ, 722 , L199 Finkbeiner, D. P. 2003, ApJS, 146, 407 Finkbeiner, D. P. 2004a, Astrophys. J., 614, 186 —. 2004b Gallimore, J. F., Axon, D. J., O’Dea, C. P., Baum, S. A., & Pedl ar, A. 2006, AJ, 132, 546 Gebauer, I., & de Boer, W. 2009, ArXiv e-prints Gehrels, N., & Michelson, P. 1999, Astropart. Phys., 11, 277 Ghez, A. M., et al. 2008, ApJ, 689, 1044 Gillessen, S., Eisenhauer, F., Trippe, S., Alexander, T., G enzel, R., Martins, F., & Ott, T. 2009, ApJ, 692, 1075 Gold, B., et al. 2010 Grenier, I. A., Casandjian, J., & Terrier, R. 2005, Science, 307, 1292 Gualandris, A., & Merritt, D. 2009, ApJ, 705, 361 Harding, J. P., & Abazajian, K. N. 2010, Phys. Rev., D81, 0235 05 Hartmann, D. H. 1995, ApJ, 447, 646 Haslam, C. G. T., Salter, C. J., Stoffel, H., & Wilson, W. E. 198 2, A&AS, 47, 1 Heesen, V., Beck, R., Krause, M., & Dettmar, R. 2009, A&A, 494 , 563 Higdon, J. C., & Lingenfelter, R. E. 2005, ApJ, 628, 738 Hills, J. G. 1975, AJ, 80, 809 —. 1988, Nature, 331, 687 Ho, L. C. 2008, ARA&A, 46, 475 Hooper, D., Finkbeiner, D. P., & Dobler, G. 2007, Phys. Rev., D76, 083012 Innes, D. E., Inhester, B., Axford, W. I., & Wilhelm, K. 1997, Nature, 386, 811 Ipavich, F. M. 1975, ApJ, 196, 107 Irwin, J. A., & Seaquist, E. R. 1988, ApJ, 335, 658 Irwin, J. A., & Sofue, Y. 1992, ApJ, 396, L75 Itoh, C., et al. 2007, A&A, 462, 67 Jokipii, J. R., & Morfill, G. 1987, ApJ, 312, 170 Jones, F. C., & Ellison, D. C. 1991, Space Science Reviews, 58 , 259 Kaifu, N., Kato, T., & Iguchi, T. 1972, Nature, 238, 105 Kalberla, P. M. W., Burton, W. B., Hartmann, D., Arnal, E. M., Bajaja, E., Morras, R., & P¨ oppel, W. G. L. 2005, A&A, 440, 775 Kaplinghat, M., Phalen, D. J., & Zurek, K. M. 2009, JCAP, 0912 , 010 Karlsson, N., & for the VERITAS collaboration. 2009, ArXiv e - prints Kniffen, D. A., & Fichtel, C. E. 1981, ApJ, 250, 389 Koyama, K., Maeda, Y., Sonobe, T., Takeshima, T., Tanaka, Y. , & Yamauchi, S. 1996, PASJ, 48, 249 Kulsrud, R., & Pearce, W. P. 1969, ApJ, 156, 445 Lacki, B.C., Thompson, T. A.,Quataert, E.,Loeb, A.,& Waxma n, E. 2010, ArXiv e-prints Large, M. I., Quigley, M. J. S., & Haslam, C. G. T. 1962, MNRAS, 124, 405 LaRosa, T. N., Nord, M. E., Lazio, T. J. W., & Kassim, N. E. 2004 , ApJ, 607, 302 Law, C. J. 2010, ApJ, 708, 474 Lazarian, A., & Opher, M. 2009, ApJ, 703, 8 Lerche, I., & Schlickeiser, R. 1982, A&A, 107, 148 Lin, T., Finkbeiner, D. P., & Dobler, G. 2010, Phys. Rev. D, 82 , 023518 Linden, T., Profumo, S., & Anderson, B. 2010, Phys. Rev. D, 82 , 063529 Lu, Y., Zhang, F., & Yu, Q. 2010, ApJ, 709, 1356 Machida, M., et al. 2009, PASJ, 61, 411 Malyshev, D., Cholis, I., & Gelfand, J. D. 2010, ApJ, 722, 193 9 Markoff, S. 2010, Proceedings of the National Academy of Scie nce, 107, 7196 Markoff, S., Bower, G. C., & Falcke, H. 2007, MNRAS, 379, 1519 McCarthy, P. J., van Breugel, W., & Heckman, T. 1987, AJ, 93, 264 McNamara, B. R., & Nulsen, P. E. J. 2007, ARA&A, 45, 117 McQuinn, M., & Zaldarriaga, M. 2010, ArXiv e-prints Middelberg, E., Agudo, I., Roy, A. L., & Krichbaum, T. P. 2007 , MNRAS, 377, 731 Morris, M., & Serabyn, E. 1996, ARA&A, 34, 645 Muno, M. P., Baganoff, F. K., Brandt, W. N., Morris, M. R., & Starck, J. 2008, ApJ, 673, 251 46 Muno, M. P., Baganoff, F. K., Brandt, W. N., Park, S., & Morris, M. R. 2007, ApJ, 656, L69 Muno, M. P., et al. 2003, ApJ, 589, 225 —. 2004, ApJ, 613, 326 Nobukawa, M., et al. 2008, PASJ, 60, 191 Page, L., et al. 2007, Astrophys. J. Suppl., 170, 335 Paumard, T., et al. 2006, ApJ, 643, 1011 Ponti, G., Terrier, R., Goldwurm, A., Belanger, G., & Trap, G . 2010, ApJ, 714, 732 Porter, T. A., & Strong, A. W. 2005, 4, 77 Prantzos, N., & Diehl, R. 1996, Phys. Rep., 267, 1 Quillen, A. C., Bland-Hawthorn, J., Brookes, M. H., Werner, M. W., Smith, J. D., Stern, D., Keene, J., & Lawrence, C. R. 2006, ApJ, 641, L29 Rees, M. J. 1988, Nature, 333, 523 Revnivtsev, M., Sazonov, S., Churazov, E., Forman, W., Vikh linin, A., & Sunyaev, R. 2009, Nature, 458, 1142 Revnivtsev, M. G., et al. 2004, A&A, 425, L49 Rosner, R., & Bodo, G. 1996, ApJ, 470, L49+ Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, Astroph ys. J., 500, 525 Scoville, N. Z. 1972, ApJ, 175, L127+ Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 Sharp, R. G., & Bland-Hawthorn, J. 2010, ApJ, 711, 818 Smialkowski, A., Wolfendale, A. W., & Zhang, L. 1997, Astroparticle Physics, 7, 21 Snowden, S. L. 2009, Space Sci. Rev., 143, 253 Snowden, S. L., et al. 1997, ApJ, 485, 125 Socrates, A., Davis, S. W., & Ramirez-Ruiz, E. 2008, ApJ, 687 , 202 Sofue, Y. 2000a, ApJ, 540, 224 —. 2000b, ApJ, 540, 224 —. 2003, PASJ, 55, 445 Sommer-Larsen, J. 2006, ApJ, 644, L1 Stark, A. A., Martin, C. L., Walsh, W. M., Xiao, K., Lane, A. P. , & Walker, C. K. 2004, ApJ, 614, L41 Sternberg, A., & Soker, N. 2008, MNRAS, 389, L13 Strickland, D. K., & Heckman, T. M. 2007, ApJ, 658, 258 Strickland, D. K., & Stevens, I. R. 2000, MNRAS, 314, 511 Strong, A. W. 1984, Advances in Space Research, 3, 87Strong, A. W., Mayer-Hasselwander, H. A., Bloemen, J. B. G. M ., Hermsen, W., & Lebrun, F. 1987, A&AS, 67, 283 Strong, A. W., & Moskalenko, I. V. 1999, 4, 255 Strong, A. W., Moskalenko, I. V., Porter, T. A., J´ ohannesso n, G., Orlando, E., & Digel, S. W. 2009, ArXiv e-prints Strong, A. W., Moskalenko, I. V., & Ptuskin, V. S. 2007, Annua l Review of Nuclear and Particle Science, 57, 285 Strong, A. W., Porter, T. A., Digel, S. W., Johannesson, G., Martin, P., Moskalenko, I. V., & Murphy, E. J. 2010, ArXiv e-prints Sun, X. H., Reich, W., Waelkens, A., & Enßlin, T. A. 2008, A&A, 477, 573 Sunyaev, R., & Churazov, E. 1998, MNRAS, 297, 1279 Sunyaev, R. A., Markevitch, M., & Pavlinsky, M. 1993, ApJ, 40 7, 606 Thompson, T. A., Quataert, E., & Waxman, E. 2007, ApJ, 654, 219 Thompson, T. A., Quataert, E., Waxman, E., Murray, N., & Martin, C. L. 2006, ApJ, 645, 186 Totani, T. 2006, PASJ, 58, 965 Veilleux, S., Cecil, G., & Bland-Hawthorn, J. 2005, ARA&A, 4 3, 769 Veilleux, S., Rupke, D. S. N., & Swaters, R. 2009, ApJ, 700, L1 49 Wang, J., & Hu, C. 2005, ApJ, 630, L125 Weidenspointner, G., et al. 2008, New Astronomy Review, 52, 454 Westmoquette, M. S., Smith, L. J., Gallagher, J. S., Trancho , G., Bastian, N., & Konstantopoulos, I. S. 2009, ApJ, 696, 192 Yu, Q., & Tremaine, S. 2003, ApJ, 599, 1129 Yuan, F., Markoff, S., & Falcke, H. 2002, A&A, 383, 854 Yusef-Zadeh, F., Hewitt, J. W., & Cotton, W. 2004, ApJS, 155, 421 Yusef-Zadeh, F., & K¨ onigl, A. 2004, in Astronomical Societ y of the Pacific Conference Series,Vol.322,The Formationand Evolu tion of Massive Young Star Clusters, ed. H. J. G. L. M. Lamers, L. J. Smith, & A. Nota, 201–+ Yusef-Zadeh, F., et al. 2006, ApJ, 644, 198 Zhang, J., et al. 2009, Phys. Rev., D80, 023007