FERMI BUBBLE γ-RAYS AS A RESULT OF DIFFUSIVE INJECTION OF GALACTIC COSMIC RAYS

Recent detailed analysis of Fermi-LAT data has discovered two giant γ-ray emission regions extending up to ∼50° (∼10 kpc) in Galactic latitude above and below the Galactic center (GC) with a width of ∼40° in longitude (Su et al. 2010). The γ-ray emission regions, now popularly known as the "Fermi bubbles" (FBs), coincide well with the Wilkinson Microwave Anisotropy Probe (WMAP) haze at low latitudes (Finkbeiner 2004; Dobler & Finkbeiner 2008), and share their edges with the ROSAT X-ray map (Snowden et al. 1997). More recently, the Planck satellite experiments have shown that the morphology of the FBs is highly consistent with that of the microwave haze (Planck Collaboration 2013), and radio measurements have also found that the bubbles are coincident with two giant radio lobes that appear to originate from the GC (Carretti et al. 2013). This correlation between the multi-wavelength observations seem to suggest that the FB γ-rays (measured in the range of ∼1–100 GeV) are produced by high-energy electrons via an inverse Compton scattering process, as the same electrons can also simultaneously produce synchrotron radiations in the presence of magnetic fields (Dobler et al. 2010). Moreover, the fact that the FBs are symmetric across the Galactic plane and also centered on the GC encourages the assumption that the bubbles, and the associated γ-ray emissions, have their origin at the GC.

However, because of severe radiative losses, it would be extremely difficult to transport high-energy electrons from the GC to the far edges of the bubbles. For electrons relevant for producing the FB γ-rays, convection by Galactic winds would require a wind speed as high as ∼10⁴ km s⁻¹ which is more than an order of magnitude larger than the typical Galactic wind speed of ∼200–300 km s⁻¹, and diffusive transport would need a diffusion coefficient which is ∼2–3 orders of magnitude larger than the standard Galactic value. This problem can be overcome using models based on the jet activity of the central active galactic nucleus where jet speeds of over 10⁴ km s⁻¹ are readily achievable (Guo & Mathews 2012; Yang et al. 2012). An alternative solution is to consider the production of electrons inside the bubbles itself, for instance, by multiple shock waves generated by periodic star capture by the central supermassive black hole (Cheng et al. 2011) or by second-order Fermi acceleration by the plasma wave turbulence present inside the bubbles (Mertsch & Sarkar 2011).

On the other hand, hadronic models suffer less constraint at least from the particle injection point of view. Crocker & Aharonian (2011) showed that the FB γ-rays can be explained if cosmic-ray (CR) nuclei from the GC are injected into the bubbles by fast winds, and if the particles remain trapped inside the bubbles for over ∼10¹⁰ yr. In almost all the models proposed so far (see Su et al. 2010 for more scenarios), an additional process of high-energy particle production, either inside or outside the bubbles, has been considered. In this letter, we present a simple model which does not invoke any additional sources or particle production processes other than those responsible for the production of the bulk of the Galactic CRs.

In our model, it is assumed that CRs (mainly protons), after leaving their sources, undergo diffusive propagation through the Galaxy. If the FBs are absent of sources, the diffusive streaming of CRs in the direction of the density gradient can result in a net flux of CRs injected into the bubbles. In addition, if the bubbles are expanding, there can also be an additional CR injection that scales linearly with the expansion velocity of the bubbles. For uniform injection throughout the bubble surface, the total CR power injected into each bubble (assuming spherical shape) is given by L = 4πR²uε, where R represents the bubble radius, u the streaming or injection velocity of CRs, and ε the CR energy density. For CRs streaming at Alfvén speed, typically ∼100 km s⁻¹, we obtain L ∼ 4 × 10⁴⁰ erg s⁻¹ for R = 4.5 kpc and ε = 1 eV cm⁻³, the locally measured energy density. This amount of injected power is ∼2 orders of magnitude larger than the power required in CR protons to produce the measured γ-ray luminosity of ∼2 × 10³⁷ erg s⁻¹ from each bubble. Even if the CR density in the halo is lower by one order of magnitude with respect to the local value, this rough estimate shows that the injection of some fraction of Galactic CRs can easily account for the measured γ-ray luminosity.

For CRs with an equilibrium number density N_g in the Galaxy and propagating with diffusion coefficient D_g, their diffusive injection flux into the bubbles is given by F_dif = D_g∇N_g∝D_gN_g, calculated at the edges of the bubbles. Since N_g is related to the CR source spectrum Q as N_g∝Q/D_g, we get F_dif∝Q. Thus, CRs streaming diffusively into the bubbles will follow the same spectral shape as the source CR spectrum in the Galaxy.

Once injected, CRs undergo diffusive propagation inside the bubbles, while at the same time also be convected radially outward by the expanding plasma. If the plasma inside the bubbles is extremely turbulent as suggested by X-ray observations of the Galactic bulge (Yao & Wang 2007) or if the magnetic field lines are highly tangled (McQuinn & Zaldarriaga 2011), then the diffusion is expected to be much slower, and the CR transport will be dominated by the convection process. This will result in a CR distribution that peaks toward the bubble edges, in agreement with the observations (Su et al. 2010). Moreover, for an energy-independent diffusion (which will be justified later), the CR spectrum inside the bubbles will follow the injection spectrum. This is important because for the hadronic origin γ-rays, their spectrum should mimic that of the parent CRs. The observed γ-ray spectral index of ∼2.0–2.1 above ∼1 GeV is quite close to the CR source index of Γ ∼ 2.1–2.3 required to produce the locally observed CR spectrum.

Inside the bubbles, secondary electrons (and positrons) can be continuously produced from inelastic collisions of CR protons with the bubble plasma. These electrons can produce synchrotron radiations, and can account for the observed radio and microwave emissions. The production spectrum of such electrons follows the proton spectrum, and has the same index as the source CR spectrum in the Galaxy. If the electrons suffer continuous synchrotron losses, their spectrum at age t_age of the bubbles will follow E^−Γ for E < E_c, and E^{−Γ − 1} for E > E_c, where E denotes the electron energy, and E_c∝1/t_age is the energy at which the energy loss time equals t_age. Detailed analysis of the combined Planck-WMAP data in the frequency range of 23–61 GHz have inferred that the electron spectral index in the region is ∼2.1 (Planck Collaboration 2013). In the present model, this hard spectrum can be explained if the electrons responsible for the microwave emissions have energies less than E_c.

2.1. Cosmic-ray Spectrum

Considering that the bubbles expand outward with constant velocity, the CR proton distribution inside a bubble can be described by a one-dimensional diffusion-loss equation written in the comoving frame of expansion:

$$\begin{equation} \frac{\partial }{\partial x}\left(D_\mathrm{b}\frac{\partial N_\mathrm{b}}{\partial x}\right)-\frac{N_\mathrm{b}}{\tau }=\frac{\partial N_\mathrm{b}}{\partial t}, \end{equation} \tag{ 1 }$$

where N_b(x, E, t) represents the differential number density of particles with kinetic energy E at time t and position x measured from the bubble boundary with x > 0 (< 0) representing the region inside (outside) the bubble, D_b represents the diffusion coefficient inside the bubbles, and τ = 1/(n_bvσ) represents the CR inelastic collision time with the bubble plasma of density n_b, v the CR velocity, and σ the collision cross-section. Energetic outflows such as unstable large-scale Galactic winds in the inner region of the Galaxy (e.g., Everett et al. 2008) can generate a turbulence wave spectrum inside the bubbles that follows k⁻² in the short wavelength regime k ≫ 1/L, where k denotes the wave number and L is the characteristic length of turbulence injection (e.g., Bykov & Toptygin 1987). Then, for turbulence presumably injected at scales of several parsecs, which is much larger than the gyro-radii of CRs relevant for the present study, the CR diffusion coefficient is expected to be energy independent. We take D_b = K × 10²⁸ cm² s⁻¹, where K is a constant. Moreover, since we are mainly interested in CRs with energies above ∼1 GeV, we neglect the very weak energy dependence of σ (Kelner et al. 2006), and take a constant value of σ = 32 mb.

For the Q(E) spectrum of particles injected at x = 0 at time t', the solution of Equation (1) is obtained as,

$$\begin{eqnarray} N_\mathrm{b}(x,E,t)=&\frac{Q(E)}{\sqrt{\pi D_\mathrm{b}(t-t^\prime)}}\exp \left[\frac{-x^2}{4D_\mathrm{b}(t-t^\prime)}-\frac{(t-t^\prime)}{\tau }\right].\nonumber\\ \end{eqnarray} \tag{ 2 }$$

The above solution does not include the adiabatic energy loss resulting from the spherical expansion of the bubble. This is taken into account as follows. The adiabatic energy loss rate is given by,

$$\begin{equation} \frac{dE}{dt}=-\frac{2}{3}\frac{E}{t}. \end{equation} \tag{ 3 }$$

Therefore, particles with energy E at the present time t had energy E' = E(t/t')^2/3 at an earlier time t', and distributed within an energy interval dE' = (t/t')^2/3dE. Then, the particle distribution in the presence of adiabatic losses can be obtained by replacing Q(E) → Q(E')(t/t')^2/3 in Equation (2). And for continuous injection of particles with flux F(E), the distribution at time t is obtained by writing Q(E') → F(E')dt', and integrating Equation (2) over t' as,

$$\begin{eqnarray} N_\mathrm{b}(x,E,t)=&\frac{1}{\sqrt{\pi D_\mathrm{b}}}\int _0^t dt^\prime \;\frac{F(E^\prime)}{\sqrt{(t-t^\prime)}}\left(\frac{t}{t^\prime }\right)^{2/3}\nonumber \\ &\times \exp \left[\frac{-x^2}{4D_\mathrm{b}(t-t^\prime)}-\frac{(t-t^\prime)}{\tau }\right]. \end{eqnarray} \tag{ 4 }$$

Note that in the present model, the injection flux is

$$\begin{equation} F=\left[D_\mathrm{g}\frac{dN_\mathrm{g}}{dx}+ UN_\mathrm{g}\right]_{x=0}, \end{equation} \tag{ 5 }$$

where the first term represents the flux of Galactic CRs injected due to their diffusive motion in the Galaxy F_dif, and the second term represents the injection flux due to the expansion of the bubble in the interstellar medium F_exp. Finally, the solutions in terms of the radial coordinate r measured from the bubble center can be found by simply replacing x with R − r in Equation (4).

To estimate the injection fluxes, we write the Galactic CR density in the pure diffusion model as function of perpendicular distance z to the Galactic plane as (Thoudam 2008),

$$\begin{equation} N_\mathrm{g}(z,E)\propto \frac{Q(E)}{D_\mathrm{g}(E)} f(z) \end{equation} \tag{ 6 }$$

where Q(E) represents the source spectrum and D_g is the diffusion coefficient in the Galaxy. The function

$$\begin{equation} f(z)=\int ^\infty _0\frac{\sinh [K(H-z)]\times \mathrm{{\textit {J}}_1}(K\Re)dK}{\sinh (KH)\left[K\coth (KH)+\eta v\sigma /(2D_\mathrm{g})\right]} \end{equation} \tag{ 7 }$$

has a very weak energy dependence, where J₁ is the Bessel function of order 1, H represents the halo boundary taken to be large enough to contain the FBs, ℜ = 20 kpc the radial size of the source distribution, and η the averaged surface density of interstellar gas in the Galactic plane. For H = 10 kpc, we obtain D_g = D₀(E/3GeV)^0.6 based on the measured boron-to-carbon ratio, where D₀ = 6 × 10²⁸ cm² s⁻¹.

Then, the CR density gradient along z follows,

$$\begin{equation} \frac{dN_\mathrm{g}}{dz}\propto \frac{Q(E)}{D_\mathrm{g}}\frac{df}{dz}. \end{equation} \tag{ 8 }$$

From Equation (8), it can be noticed that

$$\begin{equation} D_\mathrm{g}\frac{dN_\mathrm{g}}{dz}\propto Q(E) \end{equation} \tag{ 9 }$$

which shows that the diffusive injection flux F_dif follows the source CR spectrum as mentioned before. On the other hand, the injection flux due to the expansion F_exp follows the ambient CR spectrum, and hence is steeper than F_dif.

To calculated F_exp, we estimate the expansion velocity of the bubble based on the combined Planck-WMAP measurements. Assuming that the electrons radiate at critical frequency, the highest measured frequency of 61 GHz corresponds to an electron energy of 65.6 GeV for an estimated magnetic field strength of B = 1.3 μG inside the bubbles. The latter is calculated at the position of the bubble center (taken at z = 5 kpc above the GC) using the relation B(z) = 7e^−z/3 kpc μG used in the GALPROP CR propagation code (Strong et al. 2010). The fact that Planck-WMAP data has not found any break or steepening in the spectrum up to 61 GHz, implies that the age of the bubbles must be less than the synchrotron loss time of 65.6 GeV electrons. This gives an upper limit for the bubble age of t_age = 1.1 × 10⁸ yr. For the present radius of R = 4.5 kpc for the bubbles, the corresponding lower limit of the expansion velocity is obtained as U = 39.6 km s⁻¹. Choosing a different scale height, say 5 kpc, will give a slightly larger magnetic field of ∼2.5 μG inside the bubbles, lowering the upper limit of t_age and increasing the lower limit of U by factor of around two. These values of magnetic field are less than the estimated equipartition values of ∼6 μG inside the bubbles (Carretti et al. 2013), and ∼15 μG of the bubble walls (Jones et al. 2012).

The CR injection flux might vary for different positions, and also for different directions in the Galaxy, but, in the present study, we neglect such variations and assume uniform injection over the entire bubble surface. We assume that the injection flux corresponds to that at z = 5 kpc as C[F_dif + F_exp]_{z = 5 kpc}, where we have introduced a constant C, hereafter referred to as the injection fraction, in order to take care of the unknown actual fraction of CRs injected. Its value will be determined based on the measured γ-ray data. If the diffusion coefficient in the Galaxy has a spatial dependence that scales inversely with the magnetic field, we can write D₀(z) = 6 × 10²⁸e^z/3 kpc cm² s⁻¹ which gives a value of D₀ = 31.7 × 10²⁸ cm² s⁻¹ at z = 5 kpc. For this value of D₀ and U = 39.6 km s⁻¹ obtained above, the injection flux is dominated by F_dif for CR energies above 1 GeV. For larger values of U, F_exp will become significant especially at lower energies. This effect will be discussed later in Section 3.

The CR source spectrum is chosen to be a broken power-law with Γ = 2.2 and 2.09 for energies below and above 300 GeV. For the assumed form of D_g, this choice of spectral index reproduces well the proton spectrum recently measured by the ATIC (Panov et al. 2007), CREAM (Yoon et al. 2011), and PAMELA (Adriani et al. 2011) experiments which exhibit a slight hardening above ∼300 GeV.

2.2. Gamma-ray Emission

For the proton distribution N_b(r, E, t), the γ-ray emissivity q_γ(r, E_γ, t) of energy E_γ is calculated using the inelastic interaction cross-section given by Kelner et al. (2006). Then, the γ-ray intensity in a given direction characterized by the Galactic longitude l and the latitude b is calculated as

$$\begin{equation} I_\gamma (l,b,E_\gamma,t)=\frac{1}{4\pi }\int ^{y_2}_{y_1} q_\gamma (y,E_\gamma,t)dy, \end{equation} \tag{ 10 }$$

where the integration is performed along the line of sight distance y, and the integration limits are determined from the points of intersection of the line of sight with the bubble surface. In Equation (10), the emissivity previously written as a function of r has been carefully expressed in terms of y.

At any given time, the CR distribution inside the bubbles is governed mainly by the competition between diffusion and convection. Diffusion tends to make the distribution uniform while convection (and also, the other subdominant effects: inelastic collision and adiabatic losses) does the opposite. Thus, faster diffusion will produce a more uniform distribution due to the increasing CR diffusion distance in the given time. This is shown in Figure 1 where we have plotted the normalized distributions of 10 GeV CRs inside a bubble for different values of D_b taking K = (0.05–0.6). For the calculation (and also, in the following), we take t = t_age = 1.1 × 10⁸ yr, U = 39.6 km s⁻¹, and n_b = 3 × 10⁻³ cm⁻³. The latter is the averaged value in the region taken from Everett et al. (2008) that explains the diffuse soft X-ray emission. It can be seen that the distribution becomes flatter toward the bubble center with increasing K. The distribution is expected to be similar at all energies because of the energy-independent nature of D_b.

Zoom In Zoom Out Reset image size

Once the values of t, U and n_b are fixed, the CR distribution inside the bubbles is determined by the choice of D_b. Here, its value is chosen such that the resulting projected γ-ray intensity distribution matches the measured profile. It is found that choosing K = 0.26 produces a good fit to the measured data as shown in Figure 2, where the data corresponds to the southern bubble. The model predictions for both the averaged 1–5 GeV (solid line) and 5–20 GeV (dashed line) are added with backgrounds obtained by fitting horizontal lines to the respective data between 5° and 40°. It can be mentioned that neither the hadronic model presented in Crocker & Aharonian (2011) nor the leptonic model based on diffusive shock acceleration (Cheng et al. 2011) can satisfactorily explain the measured sharp edges shown in Figure 2. Both these models predicted a constant volume emissivity throughout the bubbles which will produce softer edges on the projected profile. Also shown in Figure 2 are the predictions for 100 GeV (dotted line) and 1000 GeV (dot-dashed line) γ-rays which can be tested in future. As expected, their profiles look very similar to those at low energies. Our result at high energies, say at 1000 GeV, is clearly different from that expected from the leptonic stochastic acceleration model presented in Mertsch & Sarkar (2011) which predicted a significant edge brightening at high energies.

Zoom In Zoom Out Reset image size

The γ-ray spectrum averaged over a whole bubble is shown in Figure 3 (top). The calculation assumes an injection fraction of C = 0.8 which is also the same injection value used in Figure 2. The model prediction is found to be in good agreement with the data in the energy range of ∼1–100 GeV where the measurement uncertainties are small. It is interesting to see that the same source index required to explain the measured CR spectrum in the pure diffusion model also reproduces the FB γ-ray spectrum. This would be difficult if one uses the re-acceleration model of CR propagation which requires a steeper source index of ∼2.4.

Zoom In Zoom Out Reset image size

Spectra for the inner and outer regions of a bubble separated as in Su et al. (2010) are shown in Figure 3 (bottom). The outer region is taken as a shell with thickness 1 kpc. Our model predictions are in good agreement with the measured data between ∼1–100 GeV, and reproduce the measured spectral uniformity in the two regions. Even the slightly higher intensity measured in the outer region seems to be explained.

The results presented above consider the lower limit of expansion velocity inferred from the combined Planck-WMAP data. Choosing higher expansion velocity will increase the contribution of F_exp, thereby making the total injection spectrum steeper. This will result into a steeper γ-ray spectrum. This is shown in Figure 4 (top) for different velocities in the range of U = 39.6–700 km s⁻¹, where the spectra are normalized to that of U = 39.6 km s⁻¹ at 10 GeV. It can be noticed that the model spectrum becomes steeper as U increases. This is shown explicitly in Figure 4 (bottom), where the spectral indices for the different results are plotted as function of U. The indices (represented by points) are obtained by fitting the model spectra between 1–50 GeV. The solid line together with the shaded region represents the measured index of −2.09 ± 0.06 obtained by fitting the measured spectrum in the same energy range. It can be noticed that already at U = 180 km s⁻¹, the model prediction becomes inconsistent with the measured index. This sets an upper limit on the expansion velocity at U < 180 km s⁻¹, which translates into a lower limit of the bubble age at t_age > 2.44 × 10⁷ yr. Consequently, this implies a synchrotron break in the electron spectrum at energies below 298 GeV, and also a corresponding break in the synchrotron emission spectrum at frequencies below 1258 GHz, which can be checked by future measurements.

Zoom In Zoom Out Reset image size

For a typical sound speed of ∼100 km s⁻¹ (for a gas temperature of ∼10⁶ K) in the Galactic halo, the limit of U < 180 km s⁻¹ implies a Mach number M < 1.8, even smaller at M < 1.2 if the CR pressure is also taken into account. Such weak shocks are known to be rather inefficient for particle acceleration, justifying our neglect of CR acceleration by the expanding bubbles. Moreover, for a Galactic gravitational acceleration value of ∼10⁻⁸ cm s⁻² in the halo (Ferrière 1998), this upper limit of U gives an extent of the plasma flow of ∼6–7 kpc, roughly the size of the FBs. In addition, the present values of n_b and B for the bubbles give an Alfvén speed of 52 km s⁻¹, which corresponds to a plasma beta value of β ∼ 4, indicating a low magnetic pressure bubbles.

We have shown that the γ-rays from the FBs can be a result of diffusive injection of Galactic CR protons during their propagation through the Galaxy. Some important observed properties of the bubbles are explained. Unlike other existing models, our proposed model does not consider any additional particle production processes or sources other than those responsible for the production of Galactic CRs.

I wish to thank Jörg Rachen for very detailed and insightful discussions on the model. I also thank David Jones for various helpful conversations on the topic.